1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 /*
  22  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  23  */
  24 /*
  25  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  26  * Use is subject to license terms.
  27  */
  28 
  29 /*
  30  * __vsinf: single precision vector sin
  31  *
  32  * Algorithm:
  33  *
  34  * For |x| < pi/4, approximate sin(x) by a polynomial x+x*z*(S0+
  35  * z*(S1+z*S2)) and cos(x) by a polynomial 1+z*(-1/2+z*(C0+z*(C1+
  36  * z*C2))), where z = x*x, all evaluated in double precision.
  37  *
  38  * Accuracy:
  39  *
  40  * The largest error is less than 0.6 ulps.
  41  */
  42 
  43 #include <sys/isa_defs.h>
  44 
  45 #ifdef _LITTLE_ENDIAN
  46 #define HI(x)   *(1+(int *)&x)
  47 #define LO(x)   *(unsigned *)&x
  48 #else
  49 #define HI(x)   *(int *)&x
  50 #define LO(x)   *(1+(unsigned *)&x)
  51 #endif
  52 
  53 #ifdef __RESTRICT
  54 #define restrict _Restrict
  55 #else
  56 #define restrict
  57 #endif
  58 
  59 extern int __vlibm_rem_pio2m(double *, double *, int, int, int);
  60 
  61 static const double C[] = {
  62         -1.66666552424430847168e-01,    /* 2^ -3 * -1.5555460000000 */
  63         8.33219196647405624390e-03,     /* 2^ -7 *  1.11077E0000000 */
  64         -1.95187909412197768688e-04,    /* 2^-13 * -1.9956B60000000 */
  65         1.0,
  66         -0.5,
  67         4.16666455566883087158e-02,     /* 2^ -5 *  1.55554A0000000 */
  68         -1.38873036485165357590e-03,    /* 2^-10 * -1.6C0C1E0000000 */
  69         2.44309903791872784495e-05,     /* 2^-16 *  1.99E24E0000000 */
  70         0.636619772367581343075535,     /* 2^ -1  * 1.45F306DC9C883 */
  71         6755399441055744.0,             /* 2^ 52  * 1.8000000000000 */
  72         1.570796326734125614166,        /* 2^  0  * 1.921FB54400000 */
  73         6.077100506506192601475e-11,    /* 2^-34  * 1.0B4611A626331 */
  74 };
  75 
  76 #define S0      C[0]
  77 #define S1      C[1]
  78 #define S2      C[2]
  79 #define one     C[3]
  80 #define mhalf   C[4]
  81 #define C0      C[5]
  82 #define C1      C[6]
  83 #define C2      C[7]
  84 #define invpio2 C[8]
  85 #define c3two51 C[9]
  86 #define pio2_1  C[10]
  87 #define pio2_t  C[11]
  88 
  89 #define PREPROCESS(N, index, label)                                     \
  90         hx = *(int *)x;                                                 \
  91         ix = hx & 0x7fffffff;                                               \
  92         t = *x;                                                         \
  93         x += stridex;                                                   \
  94         if (ix <= 0x3f490fdb) { /* |x| < pi/4 */                  \
  95                 if (ix == 0) {                                          \
  96                         y[index] = t;                                   \
  97                         goto label;                                     \
  98                 }                                                       \
  99                 y##N = (double)t;                                       \
 100                 n##N = 0;                                               \
 101         } else if (ix <= 0x49c90fdb) { /* |x| < 2^19*pi */                \
 102                 y##N = (double)t;                                       \
 103                 medium = 1;                                             \
 104         } else {                                                        \
 105                 if (ix >= 0x7f800000) { /* inf or nan */             \
 106                         y[index] = t / t;                               \
 107                         goto label;                                     \
 108                 }                                                       \
 109                 z##N = y##N = (double)t;                                \
 110                 hx = HI(y##N);                                          \
 111                 n##N = ((hx >> 20) & 0x7ff) - 1046;                   \
 112                 HI(z##N) = (hx & 0xfffff) | 0x41600000;                     \
 113                 n##N = __vlibm_rem_pio2m(&z##N, &y##N, n##N, 1, 0);     \
 114                 if (hx < 0) {                                                \
 115                         y##N = -y##N;                                   \
 116                         n##N = -n##N;                                   \
 117                 }                                                       \
 118                 z##N = y##N * y##N;                                     \
 119                 if (n##N & 1) { /* compute cos y */                 \
 120                         f##N = (float)(one + z##N * (mhalf + z##N *     \
 121                             (C0 + z##N * (C1 + z##N * C2))));           \
 122                 } else { /* compute sin y */                            \
 123                         f##N = (float)(y##N + y##N * z##N * (S0 +       \
 124                             z##N * (S1 + z##N * S2)));                  \
 125                 }                                                       \
 126                 y[index] = (n##N & 2)? -f##N : f##N;                        \
 127                 goto label;                                             \
 128         }
 129 
 130 #define PROCESS(N)                                                      \
 131         if (medium) {                                                   \
 132                 z##N = y##N * invpio2 + c3two51;                        \
 133                 n##N = LO(z##N);                                        \
 134                 z##N -= c3two51;                                        \
 135                 y##N = (y##N - z##N * pio2_1) - z##N * pio2_t;          \
 136         }                                                               \
 137         z##N = y##N * y##N;                                             \
 138         if (n##N & 1) { /* compute cos y */                         \
 139                 f##N = (float)(one + z##N * (mhalf + z##N * (C0 +       \
 140                     z##N * (C1 + z##N * C2))));                         \
 141         } else { /* compute sin y */                                    \
 142                 f##N = (float)(y##N + y##N * z##N * (S0 + z##N * (S1 +  \
 143                     z##N * S2)));                                       \
 144         }                                                               \
 145         *y = (n##N & 2)? -f##N : f##N;                                      \
 146         y += stridey
 147 
 148 void
 149 __vsinf(int n, float *restrict x, int stridex, float *restrict y,
 150     int stridey)
 151 {
 152         double          y0, y1, y2, y3;
 153         double          z0, z1, z2, z3;
 154         float           f0, f1, f2, f3, t;
 155         int             n0, n1, n2, n3, hx, ix, medium;
 156 
 157         y -= stridey;
 158 
 159         for (;;) {
 160 begin:
 161                 y += stridey;
 162 
 163                 if (--n < 0)
 164                         break;
 165 
 166                 medium = 0;
 167                 PREPROCESS(0, 0, begin);
 168 
 169                 if (--n < 0)
 170                         goto process1;
 171 
 172                 PREPROCESS(1, stridey, process1);
 173 
 174                 if (--n < 0)
 175                         goto process2;
 176 
 177                 PREPROCESS(2, (stridey << 1), process2);
 178 
 179                 if (--n < 0)
 180                         goto process3;
 181 
 182                 PREPROCESS(3, (stridey << 1) + stridey, process3);
 183 
 184                 if (medium) {
 185                         z0 = y0 * invpio2 + c3two51;
 186                         z1 = y1 * invpio2 + c3two51;
 187                         z2 = y2 * invpio2 + c3two51;
 188                         z3 = y3 * invpio2 + c3two51;
 189 
 190                         n0 = LO(z0);
 191                         n1 = LO(z1);
 192                         n2 = LO(z2);
 193                         n3 = LO(z3);
 194 
 195                         z0 -= c3two51;
 196                         z1 -= c3two51;
 197                         z2 -= c3two51;
 198                         z3 -= c3two51;
 199 
 200                         y0 = (y0 - z0 * pio2_1) - z0 * pio2_t;
 201                         y1 = (y1 - z1 * pio2_1) - z1 * pio2_t;
 202                         y2 = (y2 - z2 * pio2_1) - z2 * pio2_t;
 203                         y3 = (y3 - z3 * pio2_1) - z3 * pio2_t;
 204                 }
 205 
 206                 z0 = y0 * y0;
 207                 z1 = y1 * y1;
 208                 z2 = y2 * y2;
 209                 z3 = y3 * y3;
 210 
 211                 hx = (n0 & 1) | ((n1 & 1) << 1) | ((n2 & 1) << 2) |
 212                     ((n3 & 1) << 3);
 213                 switch (hx) {
 214                 case 0:
 215                         f0 = (float)(y0 + y0 * z0 * (S0 + z0 * (S1 + z0 * S2)));
 216                         f1 = (float)(y1 + y1 * z1 * (S0 + z1 * (S1 + z1 * S2)));
 217                         f2 = (float)(y2 + y2 * z2 * (S0 + z2 * (S1 + z2 * S2)));
 218                         f3 = (float)(y3 + y3 * z3 * (S0 + z3 * (S1 + z3 * S2)));
 219                         break;
 220 
 221                 case 1:
 222                         f0 = (float)(one + z0 * (mhalf + z0 * (C0 +
 223                             z0 * (C1 + z0 * C2))));
 224                         f1 = (float)(y1 + y1 * z1 * (S0 + z1 * (S1 + z1 * S2)));
 225                         f2 = (float)(y2 + y2 * z2 * (S0 + z2 * (S1 + z2 * S2)));
 226                         f3 = (float)(y3 + y3 * z3 * (S0 + z3 * (S1 + z3 * S2)));
 227                         break;
 228 
 229                 case 2:
 230                         f0 = (float)(y0 + y0 * z0 * (S0 + z0 * (S1 + z0 * S2)));
 231                         f1 = (float)(one + z1 * (mhalf + z1 * (C0 +
 232                             z1 * (C1 + z1 * C2))));
 233                         f2 = (float)(y2 + y2 * z2 * (S0 + z2 * (S1 + z2 * S2)));
 234                         f3 = (float)(y3 + y3 * z3 * (S0 + z3 * (S1 + z3 * S2)));
 235                         break;
 236 
 237                 case 3:
 238                         f0 = (float)(one + z0 * (mhalf + z0 * (C0 +
 239                             z0 * (C1 + z0 * C2))));
 240                         f1 = (float)(one + z1 * (mhalf + z1 * (C0 +
 241                             z1 * (C1 + z1 * C2))));
 242                         f2 = (float)(y2 + y2 * z2 * (S0 + z2 * (S1 + z2 * S2)));
 243                         f3 = (float)(y3 + y3 * z3 * (S0 + z3 * (S1 + z3 * S2)));
 244                         break;
 245 
 246                 case 4:
 247                         f0 = (float)(y0 + y0 * z0 * (S0 + z0 * (S1 + z0 * S2)));
 248                         f1 = (float)(y1 + y1 * z1 * (S0 + z1 * (S1 + z1 * S2)));
 249                         f2 = (float)(one + z2 * (mhalf + z2 * (C0 +
 250                             z2 * (C1 + z2 * C2))));
 251                         f3 = (float)(y3 + y3 * z3 * (S0 + z3 * (S1 + z3 * S2)));
 252                         break;
 253 
 254                 case 5:
 255                         f0 = (float)(one + z0 * (mhalf + z0 * (C0 +
 256                             z0 * (C1 + z0 * C2))));
 257                         f1 = (float)(y1 + y1 * z1 * (S0 + z1 * (S1 + z1 * S2)));
 258                         f2 = (float)(one + z2 * (mhalf + z2 * (C0 +
 259                             z2 * (C1 + z2 * C2))));
 260                         f3 = (float)(y3 + y3 * z3 * (S0 + z3 * (S1 + z3 * S2)));
 261                         break;
 262 
 263                 case 6:
 264                         f0 = (float)(y0 + y0 * z0 * (S0 + z0 * (S1 + z0 * S2)));
 265                         f1 = (float)(one + z1 * (mhalf + z1 * (C0 +
 266                             z1 * (C1 + z1 * C2))));
 267                         f2 = (float)(one + z2 * (mhalf + z2 * (C0 +
 268                             z2 * (C1 + z2 * C2))));
 269                         f3 = (float)(y3 + y3 * z3 * (S0 + z3 * (S1 + z3 * S2)));
 270                         break;
 271 
 272                 case 7:
 273                         f0 = (float)(one + z0 * (mhalf + z0 * (C0 +
 274                             z0 * (C1 + z0 * C2))));
 275                         f1 = (float)(one + z1 * (mhalf + z1 * (C0 +
 276                             z1 * (C1 + z1 * C2))));
 277                         f2 = (float)(one + z2 * (mhalf + z2 * (C0 +
 278                             z2 * (C1 + z2 * C2))));
 279                         f3 = (float)(y3 + y3 * z3 * (S0 + z3 * (S1 + z3 * S2)));
 280                         break;
 281 
 282                 case 8:
 283                         f0 = (float)(y0 + y0 * z0 * (S0 + z0 * (S1 + z0 * S2)));
 284                         f1 = (float)(y1 + y1 * z1 * (S0 + z1 * (S1 + z1 * S2)));
 285                         f2 = (float)(y2 + y2 * z2 * (S0 + z2 * (S1 + z2 * S2)));
 286                         f3 = (float)(one + z3 * (mhalf + z3 * (C0 +
 287                             z3 * (C1 + z3 * C2))));
 288                         break;
 289 
 290                 case 9:
 291                         f0 = (float)(one + z0 * (mhalf + z0 * (C0 +
 292                             z0 * (C1 + z0 * C2))));
 293                         f1 = (float)(y1 + y1 * z1 * (S0 + z1 * (S1 + z1 * S2)));
 294                         f2 = (float)(y2 + y2 * z2 * (S0 + z2 * (S1 + z2 * S2)));
 295                         f3 = (float)(one + z3 * (mhalf + z3 * (C0 +
 296                             z3 * (C1 + z3 * C2))));
 297                         break;
 298 
 299                 case 10:
 300                         f0 = (float)(y0 + y0 * z0 * (S0 + z0 * (S1 + z0 * S2)));
 301                         f1 = (float)(one + z1 * (mhalf + z1 * (C0 +
 302                             z1 * (C1 + z1 * C2))));
 303                         f2 = (float)(y2 + y2 * z2 * (S0 + z2 * (S1 + z2 * S2)));
 304                         f3 = (float)(one + z3 * (mhalf + z3 * (C0 +
 305                             z3 * (C1 + z3 * C2))));
 306                         break;
 307 
 308                 case 11:
 309                         f0 = (float)(one + z0 * (mhalf + z0 * (C0 +
 310                             z0 * (C1 + z0 * C2))));
 311                         f1 = (float)(one + z1 * (mhalf + z1 * (C0 +
 312                             z1 * (C1 + z1 * C2))));
 313                         f2 = (float)(y2 + y2 * z2 * (S0 + z2 * (S1 + z2 * S2)));
 314                         f3 = (float)(one + z3 * (mhalf + z3 * (C0 +
 315                             z3 * (C1 + z3 * C2))));
 316                         break;
 317 
 318                 case 12:
 319                         f0 = (float)(y0 + y0 * z0 * (S0 + z0 * (S1 + z0 * S2)));
 320                         f1 = (float)(y1 + y1 * z1 * (S0 + z1 * (S1 + z1 * S2)));
 321                         f2 = (float)(one + z2 * (mhalf + z2 * (C0 +
 322                             z2 * (C1 + z2 * C2))));
 323                         f3 = (float)(one + z3 * (mhalf + z3 * (C0 +
 324                             z3 * (C1 + z3 * C2))));
 325                         break;
 326 
 327                 case 13:
 328                         f0 = (float)(one + z0 * (mhalf + z0 * (C0 +
 329                             z0 * (C1 + z0 * C2))));
 330                         f1 = (float)(y1 + y1 * z1 * (S0 + z1 * (S1 + z1 * S2)));
 331                         f2 = (float)(one + z2 * (mhalf + z2 * (C0 +
 332                             z2 * (C1 + z2 * C2))));
 333                         f3 = (float)(one + z3 * (mhalf + z3 * (C0 +
 334                             z3 * (C1 + z3 * C2))));
 335                         break;
 336 
 337                 case 14:
 338                         f0 = (float)(y0 + y0 * z0 * (S0 + z0 * (S1 + z0 * S2)));
 339                         f1 = (float)(one + z1 * (mhalf + z1 * (C0 +
 340                             z1 * (C1 + z1 * C2))));
 341                         f2 = (float)(one + z2 * (mhalf + z2 * (C0 +
 342                             z2 * (C1 + z2 * C2))));
 343                         f3 = (float)(one + z3 * (mhalf + z3 * (C0 +
 344                             z3 * (C1 + z3 * C2))));
 345                         break;
 346 
 347                 default:
 348                         f0 = (float)(one + z0 * (mhalf + z0 * (C0 +
 349                             z0 * (C1 + z0 * C2))));
 350                         f1 = (float)(one + z1 * (mhalf + z1 * (C0 +
 351                             z1 * (C1 + z1 * C2))));
 352                         f2 = (float)(one + z2 * (mhalf + z2 * (C0 +
 353                             z2 * (C1 + z2 * C2))));
 354                         f3 = (float)(one + z3 * (mhalf + z3 * (C0 +
 355                             z3 * (C1 + z3 * C2))));
 356                 }
 357 
 358                 *y = (n0 & 2)? -f0 : f0;
 359                 y += stridey;
 360                 *y = (n1 & 2)? -f1 : f1;
 361                 y += stridey;
 362                 *y = (n2 & 2)? -f2 : f2;
 363                 y += stridey;
 364                 *y = (n3 & 2)? -f3 : f3;
 365                 continue;
 366 
 367 process1:
 368                 PROCESS(0);
 369                 continue;
 370 
 371 process2:
 372                 PROCESS(0);
 373                 PROCESS(1);
 374                 continue;
 375 
 376 process3:
 377                 PROCESS(0);
 378                 PROCESS(1);
 379                 PROCESS(2);
 380         }
 381 }