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 /* 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24 */ 25 /* 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27 * Use is subject to license terms. 28 */ 29 30 #include <sys/isa_defs.h> 31 #include <sys/ccompile.h> 32 33 #ifdef _LITTLE_ENDIAN 34 #define HI(x) *(1+(int*)x) 35 #define LO(x) *(unsigned*)x 36 #else 37 #define HI(x) *(int*)x 38 #define LO(x) *(1+(unsigned*)x) 39 #endif 40 41 #ifdef __RESTRICT 42 #define restrict _Restrict 43 #else 44 #define restrict 45 #endif 46 47 /* 48 * vcos.1.c 49 * 50 * Vector cosine function. Just slight modifications to vsin.8.c, mainly 51 * in the primary range part. 52 * 53 * Modification to primary range processing. If an argument that does not 54 * fall in the primary range is encountered, then processing is continued 55 * in the medium range. 56 * 57 */ 58 59 extern const double __vlibm_TBL_sincos_hi[], __vlibm_TBL_sincos_lo[]; 60 61 static const double 62 half[2] = { 0.5, -0.5 }, 63 one = 1.0, 64 invpio2 = 0.636619772367581343075535, /* 53 bits of pi/2 */ 65 pio2_1 = 1.570796326734125614166, /* first 33 bits of pi/2 */ 66 pio2_2 = 6.077100506303965976596e-11, /* second 33 bits of pi/2 */ 67 pio2_3 = 2.022266248711166455796e-21, /* third 33 bits of pi/2 */ 68 pio2_3t = 8.478427660368899643959e-32, /* pi/2 - pio2_3 */ 69 pp1 = -1.666666666605760465276263943134982554676e-0001, 70 pp2 = 8.333261209690963126718376566146180944442e-0003, 71 qq1 = -4.999999999977710986407023955908711557870e-0001, 72 qq2 = 4.166654863857219350645055881018842089580e-0002, 73 poly1[2]= { -1.666666666666629669805215138920301589656e-0001, 74 -4.999999999999931701464060878888294524481e-0001 }, 75 poly2[2]= { 8.333333332390951295683993455280336376663e-0003, 76 4.166666666394861917535640593963708222319e-0002 }, 77 poly3[2]= { -1.984126237997976692791551778230098403960e-0004, 78 -1.388888552656142867832756687736851681462e-0003 }, 79 poly4[2]= { 2.753403624854277237649987622848330351110e-0006, 80 2.478519423681460796618128289454530524759e-0005 }; 81 82 static const unsigned thresh[2] = { 0x3fc90000, 0x3fc40000 }; 83 84 /* Don't __ the following; acomp will handle it */ 85 extern double fabs(double); 86 extern void __vlibm_vcos_big(int, double *, int, double *, int, int); 87 88 /* 89 * y[i*stridey] := cos( x[i*stridex] ), for i = 0..n. 90 * 91 * Calls __vlibm_vcos_big to handle all elts which have abs >~ 1.647e+06. 92 * Argument reduction is done here for elts pi/4 < arg < 1.647e+06. 93 * 94 * elts < 2^-27 use the approximation 1.0 ~ cos(x). 95 */ 96 void 97 __vcos(int n, double * restrict x, int stridex, double * restrict y, 98 int stridey) 99 { 100 double x0_or_one[4], x1_or_one[4], x2_or_one[4]; 101 double y0_or_zero[4], y1_or_zero[4], y2_or_zero[4]; 102 double x0, x1, x2, *py0 = 0, *py1 = 0, *py2, *xsave, *ysave; 103 unsigned hx0, hx1, hx2, xsb0, xsb1 = 0, xsb2; 104 int i, biguns, nsave, sxsave, sysave; 105 volatile int v __unused; 106 nsave = n; 107 xsave = x; 108 sxsave = stridex; 109 ysave = y; 110 sysave = stridey; 111 biguns = 0; 112 113 x0 = *x; /* 'x0' may be used uninitialized */ 114 do /* MAIN LOOP */ 115 { 116 /* Gotos here so _break_ exits MAIN LOOP. */ 117 LOOP0: /* Find first arg in right range. */ 118 xsb0 = HI(x); /* get most significant word */ 119 hx0 = xsb0 & ~0x80000000; /* mask off sign bit */ 120 if (hx0 > 0x3fe921fb) { 121 /* Too big: arg reduction needed, so leave for second part */ 122 biguns = 1; 123 goto MEDIUM; 124 } 125 if (hx0 < 0x3e400000) { 126 /* Too small. cos x ~ 1. */ 127 v = *x; 128 *y = 1.0; 129 x += stridex; 130 y += stridey; 131 i = 0; 132 if (--n <= 0) 133 break; 134 goto LOOP0; 135 } 136 x0 = *x; 137 py0 = y; 138 x += stridex; 139 y += stridey; 140 i = 1; 141 if (--n <= 0) 142 break; 143 144 LOOP1: /* Get second arg, same as above. */ 145 xsb1 = HI(x); 146 hx1 = xsb1 & ~0x80000000; 147 if (hx1 > 0x3fe921fb) 148 { 149 biguns = 2; 150 goto MEDIUM; 151 } 152 if (hx1 < 0x3e400000) 153 { 154 v = *x; 155 *y = 1.0; 156 x += stridex; 157 y += stridey; 158 i = 1; 159 if (--n <= 0) 160 break; 161 goto LOOP1; 162 } 163 x1 = *x; 164 py1 = y; 165 x += stridex; 166 y += stridey; 167 i = 2; 168 if (--n <= 0) 169 break; 170 171 LOOP2: /* Get third arg, same as above. */ 172 xsb2 = HI(x); 173 hx2 = xsb2 & ~0x80000000; 174 if (hx2 > 0x3fe921fb) 175 { 176 biguns = 3; 177 goto MEDIUM; 178 } 179 if (hx2 < 0x3e400000) 180 { 181 v = *x; 182 *y = 1.0; 183 x += stridex; 184 y += stridey; 185 i = 2; 186 if (--n <= 0) 187 break; 188 goto LOOP2; 189 } 190 x2 = *x; 191 py2 = y; 192 193 /* 194 * 0x3fc40000 = 5/32 ~ 0.15625 195 * Get msb after subtraction. Will be 1 only if 196 * hx0 - 5/32 is negative. 197 */ 198 i = (hx0 - 0x3fc40000) >> 31; 199 i |= ((hx1 - 0x3fc40000) >> 30) & 2; 200 i |= ((hx2 - 0x3fc40000) >> 29) & 4; 201 switch (i) 202 { 203 double a0, a1, a2, w0, w1, w2; 204 double t0, t1, t2, z0, z1, z2; 205 unsigned j0, j1, j2; 206 207 case 0: /* All are > 5/32 */ 208 j0 = (xsb0 + 0x4000) & 0xffff8000; 209 j1 = (xsb1 + 0x4000) & 0xffff8000; 210 j2 = (xsb2 + 0x4000) & 0xffff8000; 211 HI(&t0) = j0; 212 HI(&t1) = j1; 213 HI(&t2) = j2; 214 LO(&t0) = 0; 215 LO(&t1) = 0; 216 LO(&t2) = 0; 217 x0 -= t0; 218 x1 -= t1; 219 x2 -= t2; 220 z0 = x0 * x0; 221 z1 = x1 * x1; 222 z2 = x2 * x2; 223 t0 = z0 * (qq1 + z0 * qq2); 224 t1 = z1 * (qq1 + z1 * qq2); 225 t2 = z2 * (qq1 + z2 * qq2); 226 w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 227 w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 228 w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 229 j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 230 j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 231 j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 232 xsb0 = (xsb0 >> 30) & 2; 233 xsb1 = (xsb1 >> 30) & 2; 234 xsb2 = (xsb2 >> 30) & 2; 235 a0 = __vlibm_TBL_sincos_hi[j0+1]; /* cos_hi(t) */ 236 a1 = __vlibm_TBL_sincos_hi[j1+1]; 237 a2 = __vlibm_TBL_sincos_hi[j2+1]; 238 /* cos_lo(t) sin_hi(t) */ 239 t0 = __vlibm_TBL_sincos_lo[j0+1] - (__vlibm_TBL_sincos_hi[j0+xsb0]*w0 - a0*t0); 240 t1 = __vlibm_TBL_sincos_lo[j1+1] - (__vlibm_TBL_sincos_hi[j1+xsb1]*w1 - a1*t1); 241 t2 = __vlibm_TBL_sincos_lo[j2+1] - (__vlibm_TBL_sincos_hi[j2+xsb2]*w2 - a2*t2); 242 243 *py0 = a0 + t0; 244 *py1 = a1 + t1; 245 *py2 = a2 + t2; 246 break; 247 248 case 1: 249 j1 = (xsb1 + 0x4000) & 0xffff8000; 250 j2 = (xsb2 + 0x4000) & 0xffff8000; 251 HI(&t1) = j1; 252 HI(&t2) = j2; 253 LO(&t1) = 0; 254 LO(&t2) = 0; 255 x1 -= t1; 256 x2 -= t2; 257 z0 = x0 * x0; 258 z1 = x1 * x1; 259 z2 = x2 * x2; 260 t0 = z0 * (poly3[1] + z0 * poly4[1]); 261 t1 = z1 * (qq1 + z1 * qq2); 262 t2 = z2 * (qq1 + z2 * qq2); 263 t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0)); 264 w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 265 w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 266 j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 267 j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 268 xsb1 = (xsb1 >> 30) & 2; 269 xsb2 = (xsb2 >> 30) & 2; 270 a1 = __vlibm_TBL_sincos_hi[j1+1]; 271 a2 = __vlibm_TBL_sincos_hi[j2+1]; 272 t1 = __vlibm_TBL_sincos_lo[j1+1] - (__vlibm_TBL_sincos_hi[j1+xsb1]*w1 - a1*t1); 273 t2 = __vlibm_TBL_sincos_lo[j2+1] - (__vlibm_TBL_sincos_hi[j2+xsb2]*w2 - a2*t2); 274 *py0 = one + t0; 275 *py1 = a1 + t1; 276 *py2 = a2 + t2; 277 break; 278 279 case 2: 280 j0 = (xsb0 + 0x4000) & 0xffff8000; 281 j2 = (xsb2 + 0x4000) & 0xffff8000; 282 HI(&t0) = j0; 283 HI(&t2) = j2; 284 LO(&t0) = 0; 285 LO(&t2) = 0; 286 x0 -= t0; 287 x2 -= t2; 288 z0 = x0 * x0; 289 z1 = x1 * x1; 290 z2 = x2 * x2; 291 t0 = z0 * (qq1 + z0 * qq2); 292 t1 = z1 * (poly3[1] + z1 * poly4[1]); 293 t2 = z2 * (qq1 + z2 * qq2); 294 w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 295 t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1)); 296 w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 297 j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 298 j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 299 xsb0 = (xsb0 >> 30) & 2; 300 xsb2 = (xsb2 >> 30) & 2; 301 a0 = __vlibm_TBL_sincos_hi[j0+1]; 302 a2 = __vlibm_TBL_sincos_hi[j2+1]; 303 t0 = __vlibm_TBL_sincos_lo[j0+1] - (__vlibm_TBL_sincos_hi[j0+xsb0]*w0 - a0*t0); 304 t2 = __vlibm_TBL_sincos_lo[j2+1] - (__vlibm_TBL_sincos_hi[j2+xsb2]*w2 - a2*t2); 305 *py0 = a0 + t0; 306 *py1 = one + t1; 307 *py2 = a2 + t2; 308 break; 309 310 case 3: 311 j2 = (xsb2 + 0x4000) & 0xffff8000; 312 HI(&t2) = j2; 313 LO(&t2) = 0; 314 x2 -= t2; 315 z0 = x0 * x0; 316 z1 = x1 * x1; 317 z2 = x2 * x2; 318 t0 = z0 * (poly3[1] + z0 * poly4[1]); 319 t1 = z1 * (poly3[1] + z1 * poly4[1]); 320 t2 = z2 * (qq1 + z2 * qq2); 321 t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0)); 322 t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1)); 323 w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 324 j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 325 xsb2 = (xsb2 >> 30) & 2; 326 a2 = __vlibm_TBL_sincos_hi[j2+1]; 327 t2 = __vlibm_TBL_sincos_lo[j2+1] - (__vlibm_TBL_sincos_hi[j2+xsb2]*w2 - a2*t2); 328 *py0 = one + t0; 329 *py1 = one + t1; 330 *py2 = a2 + t2; 331 break; 332 333 case 4: 334 j0 = (xsb0 + 0x4000) & 0xffff8000; 335 j1 = (xsb1 + 0x4000) & 0xffff8000; 336 HI(&t0) = j0; 337 HI(&t1) = j1; 338 LO(&t0) = 0; 339 LO(&t1) = 0; 340 x0 -= t0; 341 x1 -= t1; 342 z0 = x0 * x0; 343 z1 = x1 * x1; 344 z2 = x2 * x2; 345 t0 = z0 * (qq1 + z0 * qq2); 346 t1 = z1 * (qq1 + z1 * qq2); 347 t2 = z2 * (poly3[1] + z2 * poly4[1]); 348 w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 349 w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 350 t2 = z2 * (poly1[1] + z2 * (poly2[1] + t2)); 351 j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 352 j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 353 xsb0 = (xsb0 >> 30) & 2; 354 xsb1 = (xsb1 >> 30) & 2; 355 a0 = __vlibm_TBL_sincos_hi[j0+1]; 356 a1 = __vlibm_TBL_sincos_hi[j1+1]; 357 t0 = __vlibm_TBL_sincos_lo[j0+1] - (__vlibm_TBL_sincos_hi[j0+xsb0]*w0 - a0*t0); 358 t1 = __vlibm_TBL_sincos_lo[j1+1] - (__vlibm_TBL_sincos_hi[j1+xsb1]*w1 - a1*t1); 359 *py0 = a0 + t0; 360 *py1 = a1 + t1; 361 *py2 = one + t2; 362 break; 363 364 case 5: 365 j1 = (xsb1 + 0x4000) & 0xffff8000; 366 HI(&t1) = j1; 367 LO(&t1) = 0; 368 x1 -= t1; 369 z0 = x0 * x0; 370 z1 = x1 * x1; 371 z2 = x2 * x2; 372 t0 = z0 * (poly3[1] + z0 * poly4[1]); 373 t1 = z1 * (qq1 + z1 * qq2); 374 t2 = z2 * (poly3[1] + z2 * poly4[1]); 375 t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0)); 376 w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 377 t2 = z2 * (poly1[1] + z2 * (poly2[1] + t2)); 378 j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 379 xsb1 = (xsb1 >> 30) & 2; 380 a1 = __vlibm_TBL_sincos_hi[j1+1]; 381 t1 = __vlibm_TBL_sincos_lo[j1+1] - (__vlibm_TBL_sincos_hi[j1+xsb1]*w1 - a1*t1); 382 *py0 = one + t0; 383 *py1 = a1 + t1; 384 *py2 = one + t2; 385 break; 386 387 case 6: 388 j0 = (xsb0 + 0x4000) & 0xffff8000; 389 HI(&t0) = j0; 390 LO(&t0) = 0; 391 x0 -= t0; 392 z0 = x0 * x0; 393 z1 = x1 * x1; 394 z2 = x2 * x2; 395 t0 = z0 * (qq1 + z0 * qq2); 396 t1 = z1 * (poly3[1] + z1 * poly4[1]); 397 t2 = z2 * (poly3[1] + z2 * poly4[1]); 398 w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 399 t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1)); 400 t2 = z2 * (poly1[1] + z2 * (poly2[1] + t2)); 401 j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 402 xsb0 = (xsb0 >> 30) & 2; 403 a0 = __vlibm_TBL_sincos_hi[j0+1]; 404 t0 = __vlibm_TBL_sincos_lo[j0+1] - (__vlibm_TBL_sincos_hi[j0+xsb0]*w0 - a0*t0); 405 *py0 = a0 + t0; 406 *py1 = one + t1; 407 *py2 = one + t2; 408 break; 409 410 case 7: /* All are < 5/32 */ 411 z0 = x0 * x0; 412 z1 = x1 * x1; 413 z2 = x2 * x2; 414 t0 = z0 * (poly3[1] + z0 * poly4[1]); 415 t1 = z1 * (poly3[1] + z1 * poly4[1]); 416 t2 = z2 * (poly3[1] + z2 * poly4[1]); 417 t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0)); 418 t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1)); 419 t2 = z2 * (poly1[1] + z2 * (poly2[1] + t2)); 420 *py0 = one + t0; 421 *py1 = one + t1; 422 *py2 = one + t2; 423 break; 424 } 425 426 x += stridex; 427 y += stridey; 428 i = 0; 429 } while (--n > 0); /* END MAIN LOOP */ 430 431 /* 432 * CLEAN UP last 0, 1, or 2 elts. 433 */ 434 if (i > 0) /* Clean up elts at tail. i < 3. */ 435 { 436 double a0, a1, w0, w1; 437 double t0, t1, z0, z1; 438 unsigned j0, j1; 439 440 if (i > 1) 441 { 442 if (hx1 < 0x3fc40000) 443 { 444 z1 = x1 * x1; 445 t1 = z1 * (poly3[1] + z1 * poly4[1]); 446 t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1)); 447 t1 = one + t1; 448 *py1 = t1; 449 } 450 else 451 { 452 j1 = (xsb1 + 0x4000) & 0xffff8000; 453 HI(&t1) = j1; 454 LO(&t1) = 0; 455 x1 -= t1; 456 z1 = x1 * x1; 457 t1 = z1 * (qq1 + z1 * qq2); 458 w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 459 j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 460 xsb1 = (xsb1 >> 30) & 2; 461 a1 = __vlibm_TBL_sincos_hi[j1+1]; 462 t1 = __vlibm_TBL_sincos_lo[j1+1] 463 - (__vlibm_TBL_sincos_hi[j1+xsb1]*w1 - a1*t1); 464 *py1 = a1 + t1; 465 } 466 } 467 if (hx0 < 0x3fc40000) 468 { 469 z0 = x0 * x0; 470 t0 = z0 * (poly3[1] + z0 * poly4[1]); 471 t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0)); 472 t0 = one + t0; 473 *py0 = t0; 474 } 475 else 476 { 477 j0 = (xsb0 + 0x4000) & 0xffff8000; 478 HI(&t0) = j0; 479 LO(&t0) = 0; 480 x0 -= t0; 481 z0 = x0 * x0; 482 t0 = z0 * (qq1 + z0 * qq2); 483 w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 484 j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 485 xsb0 = (xsb0 >> 30) & 2; 486 a0 = __vlibm_TBL_sincos_hi[j0+1]; 487 t0 = __vlibm_TBL_sincos_lo[j0+1] - (__vlibm_TBL_sincos_hi[j0+xsb0]*w0 - a0*t0); 488 *py0 = a0 + t0; 489 } 490 } /* END CLEAN UP */ 491 492 return; 493 494 /* 495 * Take care of BIGUNS. 496 * 497 * We have jumped here in the middle of processing after having 498 * encountered a medium range argument. Therefore things are in a 499 * bit of a tizzy. 500 */ 501 502 MEDIUM: 503 504 x0_or_one[1] = 1.0; 505 x1_or_one[1] = 1.0; 506 x2_or_one[1] = 1.0; 507 x0_or_one[3] = -1.0; 508 x1_or_one[3] = -1.0; 509 x2_or_one[3] = -1.0; 510 y0_or_zero[1] = 0.0; 511 y1_or_zero[1] = 0.0; 512 y2_or_zero[1] = 0.0; 513 y0_or_zero[3] = 0.0; 514 y1_or_zero[3] = 0.0; 515 y2_or_zero[3] = 0.0; 516 517 if (biguns == 3) 518 { 519 biguns = 0; 520 xsb0 = xsb0 >> 31; 521 xsb1 = xsb1 >> 31; 522 goto loop2; 523 } 524 else if (biguns == 2) 525 { 526 xsb0 = xsb0 >> 31; 527 biguns = 0; 528 goto loop1; 529 } 530 biguns = 0; 531 532 do 533 { 534 double fn0, fn1, fn2, a0, a1, a2, w0, w1, w2, y0, y1, y2; 535 unsigned hx; 536 int n0, n1, n2; 537 538 /* 539 * Find 3 more to work on: Not already done, not too big. 540 */ 541 542 loop0: 543 hx = HI(x); 544 xsb0 = hx >> 31; 545 hx &= ~0x80000000; 546 if (hx > 0x413921fb) /* (1.6471e+06) Too big: leave it. */ 547 { 548 if (hx >= 0x7ff00000) /* Inf or NaN */ 549 { 550 x0 = *x; 551 *y = x0 - x0; 552 } 553 else 554 biguns = 1; 555 x += stridex; 556 y += stridey; 557 i = 0; 558 if (--n <= 0) 559 break; 560 goto loop0; 561 } 562 x0 = *x; 563 py0 = y; 564 x += stridex; 565 y += stridey; 566 i = 1; 567 if (--n <= 0) 568 break; 569 570 loop1: 571 hx = HI(x); 572 xsb1 = hx >> 31; 573 hx &= ~0x80000000; 574 if (hx > 0x413921fb) 575 { 576 if (hx >= 0x7ff00000) 577 { 578 x1 = *x; 579 *y = x1 - x1; 580 } 581 else 582 biguns = 1; 583 x += stridex; 584 y += stridey; 585 i = 1; 586 if (--n <= 0) 587 break; 588 goto loop1; 589 } 590 x1 = *x; 591 py1 = y; 592 x += stridex; 593 y += stridey; 594 i = 2; 595 if (--n <= 0) 596 break; 597 598 loop2: 599 hx = HI(x); 600 xsb2 = hx >> 31; 601 hx &= ~0x80000000; 602 if (hx > 0x413921fb) 603 { 604 if (hx >= 0x7ff00000) 605 { 606 x2 = *x; 607 *y = x2 - x2; 608 } 609 else 610 biguns = 1; 611 x += stridex; 612 y += stridey; 613 i = 2; 614 if (--n <= 0) 615 break; 616 goto loop2; 617 } 618 x2 = *x; 619 py2 = y; 620 621 n0 = (int) (x0 * invpio2 + half[xsb0]); 622 n1 = (int) (x1 * invpio2 + half[xsb1]); 623 n2 = (int) (x2 * invpio2 + half[xsb2]); 624 fn0 = (double) n0; 625 fn1 = (double) n1; 626 fn2 = (double) n2; 627 n0 = (n0 + 1) & 3; /* Add 1 (before the mod) to make sin into cos */ 628 n1 = (n1 + 1) & 3; 629 n2 = (n2 + 1) & 3; 630 a0 = x0 - fn0 * pio2_1; 631 a1 = x1 - fn1 * pio2_1; 632 a2 = x2 - fn2 * pio2_1; 633 w0 = fn0 * pio2_2; 634 w1 = fn1 * pio2_2; 635 w2 = fn2 * pio2_2; 636 x0 = a0 - w0; 637 x1 = a1 - w1; 638 x2 = a2 - w2; 639 y0 = (a0 - x0) - w0; 640 y1 = (a1 - x1) - w1; 641 y2 = (a2 - x2) - w2; 642 a0 = x0; 643 a1 = x1; 644 a2 = x2; 645 w0 = fn0 * pio2_3 - y0; 646 w1 = fn1 * pio2_3 - y1; 647 w2 = fn2 * pio2_3 - y2; 648 x0 = a0 - w0; 649 x1 = a1 - w1; 650 x2 = a2 - w2; 651 y0 = (a0 - x0) - w0; 652 y1 = (a1 - x1) - w1; 653 y2 = (a2 - x2) - w2; 654 a0 = x0; 655 a1 = x1; 656 a2 = x2; 657 w0 = fn0 * pio2_3t - y0; 658 w1 = fn1 * pio2_3t - y1; 659 w2 = fn2 * pio2_3t - y2; 660 x0 = a0 - w0; 661 x1 = a1 - w1; 662 x2 = a2 - w2; 663 y0 = (a0 - x0) - w0; 664 y1 = (a1 - x1) - w1; 665 y2 = (a2 - x2) - w2; 666 xsb0 = HI(&x0); 667 i = ((xsb0 & ~0x80000000) - thresh[n0&1]) >> 31; 668 xsb1 = HI(&x1); 669 i |= (((xsb1 & ~0x80000000) - thresh[n1&1]) >> 30) & 2; 670 xsb2 = HI(&x2); 671 i |= (((xsb2 & ~0x80000000) - thresh[n2&1]) >> 29) & 4; 672 switch (i) 673 { 674 double t0, t1, t2, z0, z1, z2; 675 unsigned j0, j1, j2; 676 677 case 0: 678 j0 = (xsb0 + 0x4000) & 0xffff8000; 679 j1 = (xsb1 + 0x4000) & 0xffff8000; 680 j2 = (xsb2 + 0x4000) & 0xffff8000; 681 HI(&t0) = j0; 682 HI(&t1) = j1; 683 HI(&t2) = j2; 684 LO(&t0) = 0; 685 LO(&t1) = 0; 686 LO(&t2) = 0; 687 x0 = (x0 - t0) + y0; 688 x1 = (x1 - t1) + y1; 689 x2 = (x2 - t2) + y2; 690 z0 = x0 * x0; 691 z1 = x1 * x1; 692 z2 = x2 * x2; 693 t0 = z0 * (qq1 + z0 * qq2); 694 t1 = z1 * (qq1 + z1 * qq2); 695 t2 = z2 * (qq1 + z2 * qq2); 696 w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 697 w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 698 w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 699 j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 700 j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 701 j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 702 xsb0 = (xsb0 >> 30) & 2; 703 xsb1 = (xsb1 >> 30) & 2; 704 xsb2 = (xsb2 >> 30) & 2; 705 n0 ^= (xsb0 & ~(n0 << 1)); 706 n1 ^= (xsb1 & ~(n1 << 1)); 707 n2 ^= (xsb2 & ~(n2 << 1)); 708 xsb0 |= 1; 709 xsb1 |= 1; 710 xsb2 |= 1; 711 a0 = __vlibm_TBL_sincos_hi[j0+n0]; 712 a1 = __vlibm_TBL_sincos_hi[j1+n1]; 713 a2 = __vlibm_TBL_sincos_hi[j2+n2]; 714 t0 = (__vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+n0]; 715 t1 = (__vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+n1]; 716 t2 = (__vlibm_TBL_sincos_hi[j2+((n2+xsb2)&3)] * w2 + a2 * t2) + __vlibm_TBL_sincos_lo[j2+n2]; 717 *py0 = ( a0 + t0 ); 718 *py1 = ( a1 + t1 ); 719 *py2 = ( a2 + t2 ); 720 break; 721 722 case 1: 723 j0 = n0 & 1; 724 j1 = (xsb1 + 0x4000) & 0xffff8000; 725 j2 = (xsb2 + 0x4000) & 0xffff8000; 726 HI(&t1) = j1; 727 HI(&t2) = j2; 728 LO(&t1) = 0; 729 LO(&t2) = 0; 730 x0_or_one[0] = x0; 731 x0_or_one[2] = -x0; 732 y0_or_zero[0] = y0; 733 y0_or_zero[2] = -y0; 734 x1 = (x1 - t1) + y1; 735 x2 = (x2 - t2) + y2; 736 z0 = x0 * x0; 737 z1 = x1 * x1; 738 z2 = x2 * x2; 739 t0 = z0 * (poly3[j0] + z0 * poly4[j0]); 740 t1 = z1 * (qq1 + z1 * qq2); 741 t2 = z2 * (qq1 + z2 * qq2); 742 t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0)); 743 w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 744 w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 745 j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 746 j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 747 xsb1 = (xsb1 >> 30) & 2; 748 xsb2 = (xsb2 >> 30) & 2; 749 n1 ^= (xsb1 & ~(n1 << 1)); 750 n2 ^= (xsb2 & ~(n2 << 1)); 751 xsb1 |= 1; 752 xsb2 |= 1; 753 a1 = __vlibm_TBL_sincos_hi[j1+n1]; 754 a2 = __vlibm_TBL_sincos_hi[j2+n2]; 755 t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0); 756 t1 = (__vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+n1]; 757 t2 = (__vlibm_TBL_sincos_hi[j2+((n2+xsb2)&3)] * w2 + a2 * t2) + __vlibm_TBL_sincos_lo[j2+n2]; 758 *py0 = t0; 759 *py1 = ( a1 + t1 ); 760 *py2 = ( a2 + t2 ); 761 break; 762 763 case 2: 764 j0 = (xsb0 + 0x4000) & 0xffff8000; 765 j1 = n1 & 1; 766 j2 = (xsb2 + 0x4000) & 0xffff8000; 767 HI(&t0) = j0; 768 HI(&t2) = j2; 769 LO(&t0) = 0; 770 LO(&t2) = 0; 771 x1_or_one[0] = x1; 772 x1_or_one[2] = -x1; 773 x0 = (x0 - t0) + y0; 774 y1_or_zero[0] = y1; 775 y1_or_zero[2] = -y1; 776 x2 = (x2 - t2) + y2; 777 z0 = x0 * x0; 778 z1 = x1 * x1; 779 z2 = x2 * x2; 780 t0 = z0 * (qq1 + z0 * qq2); 781 t1 = z1 * (poly3[j1] + z1 * poly4[j1]); 782 t2 = z2 * (qq1 + z2 * qq2); 783 w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 784 t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1)); 785 w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 786 j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 787 j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 788 xsb0 = (xsb0 >> 30) & 2; 789 xsb2 = (xsb2 >> 30) & 2; 790 n0 ^= (xsb0 & ~(n0 << 1)); 791 n2 ^= (xsb2 & ~(n2 << 1)); 792 xsb0 |= 1; 793 xsb2 |= 1; 794 a0 = __vlibm_TBL_sincos_hi[j0+n0]; 795 a2 = __vlibm_TBL_sincos_hi[j2+n2]; 796 t0 = (__vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+n0]; 797 t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1); 798 t2 = (__vlibm_TBL_sincos_hi[j2+((n2+xsb2)&3)] * w2 + a2 * t2) + __vlibm_TBL_sincos_lo[j2+n2]; 799 *py0 = ( a0 + t0 ); 800 *py1 = t1; 801 *py2 = ( a2 + t2 ); 802 break; 803 804 case 3: 805 j0 = n0 & 1; 806 j1 = n1 & 1; 807 j2 = (xsb2 + 0x4000) & 0xffff8000; 808 HI(&t2) = j2; 809 LO(&t2) = 0; 810 x0_or_one[0] = x0; 811 x0_or_one[2] = -x0; 812 x1_or_one[0] = x1; 813 x1_or_one[2] = -x1; 814 y0_or_zero[0] = y0; 815 y0_or_zero[2] = -y0; 816 y1_or_zero[0] = y1; 817 y1_or_zero[2] = -y1; 818 x2 = (x2 - t2) + y2; 819 z0 = x0 * x0; 820 z1 = x1 * x1; 821 z2 = x2 * x2; 822 t0 = z0 * (poly3[j0] + z0 * poly4[j0]); 823 t1 = z1 * (poly3[j1] + z1 * poly4[j1]); 824 t2 = z2 * (qq1 + z2 * qq2); 825 t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0)); 826 t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1)); 827 w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 828 j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 829 xsb2 = (xsb2 >> 30) & 2; 830 n2 ^= (xsb2 & ~(n2 << 1)); 831 xsb2 |= 1; 832 a2 = __vlibm_TBL_sincos_hi[j2+n2]; 833 t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0); 834 t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1); 835 t2 = (__vlibm_TBL_sincos_hi[j2+((n2+xsb2)&3)] * w2 + a2 * t2) + __vlibm_TBL_sincos_lo[j2+n2]; 836 *py0 = t0; 837 *py1 = t1; 838 *py2 = ( a2 + t2 ); 839 break; 840 841 case 4: 842 j0 = (xsb0 + 0x4000) & 0xffff8000; 843 j1 = (xsb1 + 0x4000) & 0xffff8000; 844 j2 = n2 & 1; 845 HI(&t0) = j0; 846 HI(&t1) = j1; 847 LO(&t0) = 0; 848 LO(&t1) = 0; 849 x2_or_one[0] = x2; 850 x2_or_one[2] = -x2; 851 x0 = (x0 - t0) + y0; 852 x1 = (x1 - t1) + y1; 853 y2_or_zero[0] = y2; 854 y2_or_zero[2] = -y2; 855 z0 = x0 * x0; 856 z1 = x1 * x1; 857 z2 = x2 * x2; 858 t0 = z0 * (qq1 + z0 * qq2); 859 t1 = z1 * (qq1 + z1 * qq2); 860 t2 = z2 * (poly3[j2] + z2 * poly4[j2]); 861 w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 862 w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 863 t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2)); 864 j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 865 j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 866 xsb0 = (xsb0 >> 30) & 2; 867 xsb1 = (xsb1 >> 30) & 2; 868 n0 ^= (xsb0 & ~(n0 << 1)); 869 n1 ^= (xsb1 & ~(n1 << 1)); 870 xsb0 |= 1; 871 xsb1 |= 1; 872 a0 = __vlibm_TBL_sincos_hi[j0+n0]; 873 a1 = __vlibm_TBL_sincos_hi[j1+n1]; 874 t0 = (__vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+n0]; 875 t1 = (__vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+n1]; 876 t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2); 877 *py0 = ( a0 + t0 ); 878 *py1 = ( a1 + t1 ); 879 *py2 = t2; 880 break; 881 882 case 5: 883 j0 = n0 & 1; 884 j1 = (xsb1 + 0x4000) & 0xffff8000; 885 j2 = n2 & 1; 886 HI(&t1) = j1; 887 LO(&t1) = 0; 888 x0_or_one[0] = x0; 889 x0_or_one[2] = -x0; 890 x2_or_one[0] = x2; 891 x2_or_one[2] = -x2; 892 y0_or_zero[0] = y0; 893 y0_or_zero[2] = -y0; 894 x1 = (x1 - t1) + y1; 895 y2_or_zero[0] = y2; 896 y2_or_zero[2] = -y2; 897 z0 = x0 * x0; 898 z1 = x1 * x1; 899 z2 = x2 * x2; 900 t0 = z0 * (poly3[j0] + z0 * poly4[j0]); 901 t1 = z1 * (qq1 + z1 * qq2); 902 t2 = z2 * (poly3[j2] + z2 * poly4[j2]); 903 t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0)); 904 w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 905 t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2)); 906 j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 907 xsb1 = (xsb1 >> 30) & 2; 908 n1 ^= (xsb1 & ~(n1 << 1)); 909 xsb1 |= 1; 910 a1 = __vlibm_TBL_sincos_hi[j1+n1]; 911 t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0); 912 t1 = (__vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+n1]; 913 t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2); 914 *py0 = t0; 915 *py1 = ( a1 + t1 ); 916 *py2 = t2; 917 break; 918 919 case 6: 920 j0 = (xsb0 + 0x4000) & 0xffff8000; 921 j1 = n1 & 1; 922 j2 = n2 & 1; 923 HI(&t0) = j0; 924 LO(&t0) = 0; 925 x1_or_one[0] = x1; 926 x1_or_one[2] = -x1; 927 x2_or_one[0] = x2; 928 x2_or_one[2] = -x2; 929 x0 = (x0 - t0) + y0; 930 y1_or_zero[0] = y1; 931 y1_or_zero[2] = -y1; 932 y2_or_zero[0] = y2; 933 y2_or_zero[2] = -y2; 934 z0 = x0 * x0; 935 z1 = x1 * x1; 936 z2 = x2 * x2; 937 t0 = z0 * (qq1 + z0 * qq2); 938 t1 = z1 * (poly3[j1] + z1 * poly4[j1]); 939 t2 = z2 * (poly3[j2] + z2 * poly4[j2]); 940 w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 941 t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1)); 942 t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2)); 943 j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 944 xsb0 = (xsb0 >> 30) & 2; 945 n0 ^= (xsb0 & ~(n0 << 1)); 946 xsb0 |= 1; 947 a0 = __vlibm_TBL_sincos_hi[j0+n0]; 948 t0 = (__vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+n0]; 949 t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1); 950 t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2); 951 *py0 = ( a0 + t0 ); 952 *py1 = t1; 953 *py2 = t2; 954 break; 955 956 case 7: 957 j0 = n0 & 1; 958 j1 = n1 & 1; 959 j2 = n2 & 1; 960 x0_or_one[0] = x0; 961 x0_or_one[2] = -x0; 962 x1_or_one[0] = x1; 963 x1_or_one[2] = -x1; 964 x2_or_one[0] = x2; 965 x2_or_one[2] = -x2; 966 y0_or_zero[0] = y0; 967 y0_or_zero[2] = -y0; 968 y1_or_zero[0] = y1; 969 y1_or_zero[2] = -y1; 970 y2_or_zero[0] = y2; 971 y2_or_zero[2] = -y2; 972 z0 = x0 * x0; 973 z1 = x1 * x1; 974 z2 = x2 * x2; 975 t0 = z0 * (poly3[j0] + z0 * poly4[j0]); 976 t1 = z1 * (poly3[j1] + z1 * poly4[j1]); 977 t2 = z2 * (poly3[j2] + z2 * poly4[j2]); 978 t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0)); 979 t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1)); 980 t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2)); 981 t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0); 982 t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1); 983 t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2); 984 *py0 = t0; 985 *py1 = t1; 986 *py2 = t2; 987 break; 988 } 989 990 x += stridex; 991 y += stridey; 992 i = 0; 993 } while (--n > 0); 994 995 if (i > 0) 996 { 997 double fn0, fn1, a0, a1, w0, w1, y0, y1; 998 double t0, t1, z0, z1; 999 unsigned j0, j1; 1000 int n0, n1; 1001 1002 if (i > 1) 1003 { 1004 n1 = (int) (x1 * invpio2 + half[xsb1]); 1005 fn1 = (double) n1; 1006 n1 = (n1 + 1) & 3; /* Add 1 (before the mod) to make sin into cos */ 1007 a1 = x1 - fn1 * pio2_1; 1008 w1 = fn1 * pio2_2; 1009 x1 = a1 - w1; 1010 y1 = (a1 - x1) - w1; 1011 a1 = x1; 1012 w1 = fn1 * pio2_3 - y1; 1013 x1 = a1 - w1; 1014 y1 = (a1 - x1) - w1; 1015 a1 = x1; 1016 w1 = fn1 * pio2_3t - y1; 1017 x1 = a1 - w1; 1018 y1 = (a1 - x1) - w1; 1019 xsb1 = HI(&x1); 1020 if ((xsb1 & ~0x80000000) < thresh[n1&1]) 1021 { 1022 j1 = n1 & 1; 1023 x1_or_one[0] = x1; 1024 x1_or_one[2] = -x1; 1025 y1_or_zero[0] = y1; 1026 y1_or_zero[2] = -y1; 1027 z1 = x1 * x1; 1028 t1 = z1 * (poly3[j1] + z1 * poly4[j1]); 1029 t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1)); 1030 t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1); 1031 *py1 = t1; 1032 } 1033 else 1034 { 1035 j1 = (xsb1 + 0x4000) & 0xffff8000; 1036 HI(&t1) = j1; 1037 LO(&t1) = 0; 1038 x1 = (x1 - t1) + y1; 1039 z1 = x1 * x1; 1040 t1 = z1 * (qq1 + z1 * qq2); 1041 w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 1042 j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 1043 xsb1 = (xsb1 >> 30) & 2; 1044 n1 ^= (xsb1 & ~(n1 << 1)); 1045 xsb1 |= 1; 1046 a1 = __vlibm_TBL_sincos_hi[j1+n1]; 1047 t1 = (__vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+n1]; 1048 *py1 = ( a1 + t1 ); 1049 } 1050 } 1051 n0 = (int) (x0 * invpio2 + half[xsb0]); 1052 fn0 = (double) n0; 1053 n0 = (n0 + 1) & 3; /* Add 1 (before the mod) to make sin into cos */ 1054 a0 = x0 - fn0 * pio2_1; 1055 w0 = fn0 * pio2_2; 1056 x0 = a0 - w0; 1057 y0 = (a0 - x0) - w0; 1058 a0 = x0; 1059 w0 = fn0 * pio2_3 - y0; 1060 x0 = a0 - w0; 1061 y0 = (a0 - x0) - w0; 1062 a0 = x0; 1063 w0 = fn0 * pio2_3t - y0; 1064 x0 = a0 - w0; 1065 y0 = (a0 - x0) - w0; 1066 xsb0 = HI(&x0); 1067 if ((xsb0 & ~0x80000000) < thresh[n0&1]) 1068 { 1069 j0 = n0 & 1; 1070 x0_or_one[0] = x0; 1071 x0_or_one[2] = -x0; 1072 y0_or_zero[0] = y0; 1073 y0_or_zero[2] = -y0; 1074 z0 = x0 * x0; 1075 t0 = z0 * (poly3[j0] + z0 * poly4[j0]); 1076 t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0)); 1077 t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0); 1078 *py0 = t0; 1079 } 1080 else 1081 { 1082 j0 = (xsb0 + 0x4000) & 0xffff8000; 1083 HI(&t0) = j0; 1084 LO(&t0) = 0; 1085 x0 = (x0 - t0) + y0; 1086 z0 = x0 * x0; 1087 t0 = z0 * (qq1 + z0 * qq2); 1088 w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 1089 j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 1090 xsb0 = (xsb0 >> 30) & 2; 1091 n0 ^= (xsb0 & ~(n0 << 1)); 1092 xsb0 |= 1; 1093 a0 = __vlibm_TBL_sincos_hi[j0+n0]; 1094 t0 = (__vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+n0]; 1095 *py0 = ( a0 + t0 ); 1096 } 1097 } 1098 1099 if (biguns) 1100 __vlibm_vcos_big(nsave, xsave, sxsave, ysave, sysave, 0x413921fb); 1101 } 1102