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