125c28e83SPiotr Jasiukajtis /* 225c28e83SPiotr Jasiukajtis * CDDL HEADER START 325c28e83SPiotr Jasiukajtis * 425c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the 525c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License"). 625c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License. 725c28e83SPiotr Jasiukajtis * 825c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 925c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing. 1025c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions 1125c28e83SPiotr Jasiukajtis * and limitations under the License. 1225c28e83SPiotr Jasiukajtis * 1325c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each 1425c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 1525c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the 1625c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying 1725c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner] 1825c28e83SPiotr Jasiukajtis * 1925c28e83SPiotr Jasiukajtis * CDDL HEADER END 2025c28e83SPiotr Jasiukajtis */ 2125c28e83SPiotr Jasiukajtis 2225c28e83SPiotr Jasiukajtis /* 2325c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 2425c28e83SPiotr Jasiukajtis */ 2525c28e83SPiotr Jasiukajtis /* 2625c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 2725c28e83SPiotr Jasiukajtis * Use is subject to license terms. 2825c28e83SPiotr Jasiukajtis */ 2925c28e83SPiotr Jasiukajtis 3025c28e83SPiotr Jasiukajtis #include <sys/isa_defs.h> 3125c28e83SPiotr Jasiukajtis #include <sys/ccompile.h> 3225c28e83SPiotr Jasiukajtis 3325c28e83SPiotr Jasiukajtis #ifdef _LITTLE_ENDIAN 3425c28e83SPiotr Jasiukajtis #define HI(x) *(1+(int*)x) 3525c28e83SPiotr Jasiukajtis #define LO(x) *(unsigned*)x 3625c28e83SPiotr Jasiukajtis #else 3725c28e83SPiotr Jasiukajtis #define HI(x) *(int*)x 3825c28e83SPiotr Jasiukajtis #define LO(x) *(1+(unsigned*)x) 3925c28e83SPiotr Jasiukajtis #endif 4025c28e83SPiotr Jasiukajtis 4125c28e83SPiotr Jasiukajtis #ifdef __RESTRICT 4225c28e83SPiotr Jasiukajtis #define restrict _Restrict 4325c28e83SPiotr Jasiukajtis #else 4425c28e83SPiotr Jasiukajtis #define restrict 4525c28e83SPiotr Jasiukajtis #endif 4625c28e83SPiotr Jasiukajtis 4725c28e83SPiotr Jasiukajtis extern const double __vlibm_TBL_sincos_hi[], __vlibm_TBL_sincos_lo[]; 4825c28e83SPiotr Jasiukajtis 4925c28e83SPiotr Jasiukajtis static const double 5025c28e83SPiotr Jasiukajtis half[2] = { 0.5, -0.5 }, 5125c28e83SPiotr Jasiukajtis one = 1.0, 5225c28e83SPiotr Jasiukajtis invpio2 = 0.636619772367581343075535, 5325c28e83SPiotr Jasiukajtis pio2_1 = 1.570796326734125614166, 5425c28e83SPiotr Jasiukajtis pio2_2 = 6.077100506303965976596e-11, 5525c28e83SPiotr Jasiukajtis pio2_3 = 2.022266248711166455796e-21, 5625c28e83SPiotr Jasiukajtis pio2_3t = 8.478427660368899643959e-32, 5725c28e83SPiotr Jasiukajtis pp1 = -1.666666666605760465276263943134982554676e-0001, 5825c28e83SPiotr Jasiukajtis pp2 = 8.333261209690963126718376566146180944442e-0003, 5925c28e83SPiotr Jasiukajtis qq1 = -4.999999999977710986407023955908711557870e-0001, 6025c28e83SPiotr Jasiukajtis qq2 = 4.166654863857219350645055881018842089580e-0002, 6125c28e83SPiotr Jasiukajtis poly1[2]= { -1.666666666666629669805215138920301589656e-0001, 6225c28e83SPiotr Jasiukajtis -4.999999999999931701464060878888294524481e-0001 }, 6325c28e83SPiotr Jasiukajtis poly2[2]= { 8.333333332390951295683993455280336376663e-0003, 6425c28e83SPiotr Jasiukajtis 4.166666666394861917535640593963708222319e-0002 }, 6525c28e83SPiotr Jasiukajtis poly3[2]= { -1.984126237997976692791551778230098403960e-0004, 6625c28e83SPiotr Jasiukajtis -1.388888552656142867832756687736851681462e-0003 }, 6725c28e83SPiotr Jasiukajtis poly4[2]= { 2.753403624854277237649987622848330351110e-0006, 6825c28e83SPiotr Jasiukajtis 2.478519423681460796618128289454530524759e-0005 }; 6925c28e83SPiotr Jasiukajtis 7025c28e83SPiotr Jasiukajtis static const unsigned thresh[2] = { 0x3fc90000, 0x3fc40000 }; 7125c28e83SPiotr Jasiukajtis 7225c28e83SPiotr Jasiukajtis /* Don't __ the following; acomp will handle it */ 7325c28e83SPiotr Jasiukajtis extern double fabs(double); 7425c28e83SPiotr Jasiukajtis extern void __vlibm_vsin_big(int, double *, int, double *, int, int); 7525c28e83SPiotr Jasiukajtis 7625c28e83SPiotr Jasiukajtis void 7725c28e83SPiotr Jasiukajtis __vsin(int n, double * restrict x, int stridex, double * restrict y, 7825c28e83SPiotr Jasiukajtis int stridey) 7925c28e83SPiotr Jasiukajtis { 8025c28e83SPiotr Jasiukajtis double x0_or_one[4], x1_or_one[4], x2_or_one[4]; 8125c28e83SPiotr Jasiukajtis double y0_or_zero[4], y1_or_zero[4], y2_or_zero[4]; 8225c28e83SPiotr Jasiukajtis double x0, x1, x2, *py0 = 0, *py1 = 0, *py2, *xsave, *ysave; 8325c28e83SPiotr Jasiukajtis unsigned hx0, hx1, hx2, xsb0, xsb1 = 0, xsb2; 8425c28e83SPiotr Jasiukajtis int i, biguns, nsave, sxsave, sysave; 85c81a25e9SToomas Soome volatile int v __unused; 8625c28e83SPiotr Jasiukajtis nsave = n; 8725c28e83SPiotr Jasiukajtis xsave = x; 8825c28e83SPiotr Jasiukajtis sxsave = stridex; 8925c28e83SPiotr Jasiukajtis ysave = y; 9025c28e83SPiotr Jasiukajtis sysave = stridey; 9125c28e83SPiotr Jasiukajtis biguns = 0; 9225c28e83SPiotr Jasiukajtis 93*8b8bfb2aSToomas Soome x0 = *x; /* error: 'x0' may be used uninitialized */ 9425c28e83SPiotr Jasiukajtis do 9525c28e83SPiotr Jasiukajtis { 9625c28e83SPiotr Jasiukajtis LOOP0: 9725c28e83SPiotr Jasiukajtis xsb0 = HI(x); 9825c28e83SPiotr Jasiukajtis hx0 = xsb0 & ~0x80000000; 9925c28e83SPiotr Jasiukajtis if (hx0 > 0x3fe921fb) 10025c28e83SPiotr Jasiukajtis { 10125c28e83SPiotr Jasiukajtis biguns = 1; 10225c28e83SPiotr Jasiukajtis goto MEDIUM; 10325c28e83SPiotr Jasiukajtis } 10425c28e83SPiotr Jasiukajtis if (hx0 < 0x3e400000) 10525c28e83SPiotr Jasiukajtis { 10625c28e83SPiotr Jasiukajtis v = *x; 10725c28e83SPiotr Jasiukajtis *y = *x; 10825c28e83SPiotr Jasiukajtis x += stridex; 10925c28e83SPiotr Jasiukajtis y += stridey; 11025c28e83SPiotr Jasiukajtis i = 0; 11125c28e83SPiotr Jasiukajtis if (--n <= 0) 11225c28e83SPiotr Jasiukajtis break; 11325c28e83SPiotr Jasiukajtis goto LOOP0; 11425c28e83SPiotr Jasiukajtis } 11525c28e83SPiotr Jasiukajtis x0 = *x; 11625c28e83SPiotr Jasiukajtis py0 = y; 11725c28e83SPiotr Jasiukajtis x += stridex; 11825c28e83SPiotr Jasiukajtis y += stridey; 11925c28e83SPiotr Jasiukajtis i = 1; 12025c28e83SPiotr Jasiukajtis if (--n <= 0) 12125c28e83SPiotr Jasiukajtis break; 12225c28e83SPiotr Jasiukajtis 12325c28e83SPiotr Jasiukajtis LOOP1: 12425c28e83SPiotr Jasiukajtis xsb1 = HI(x); 12525c28e83SPiotr Jasiukajtis hx1 = xsb1 & ~0x80000000; 12625c28e83SPiotr Jasiukajtis if (hx1 > 0x3fe921fb) 12725c28e83SPiotr Jasiukajtis { 12825c28e83SPiotr Jasiukajtis biguns = 2; 12925c28e83SPiotr Jasiukajtis goto MEDIUM; 13025c28e83SPiotr Jasiukajtis } 13125c28e83SPiotr Jasiukajtis if (hx1 < 0x3e400000) 13225c28e83SPiotr Jasiukajtis { 13325c28e83SPiotr Jasiukajtis v = *x; 13425c28e83SPiotr Jasiukajtis *y = *x; 13525c28e83SPiotr Jasiukajtis x += stridex; 13625c28e83SPiotr Jasiukajtis y += stridey; 13725c28e83SPiotr Jasiukajtis i = 1; 13825c28e83SPiotr Jasiukajtis if (--n <= 0) 13925c28e83SPiotr Jasiukajtis break; 14025c28e83SPiotr Jasiukajtis goto LOOP1; 14125c28e83SPiotr Jasiukajtis } 14225c28e83SPiotr Jasiukajtis x1 = *x; 14325c28e83SPiotr Jasiukajtis py1 = y; 14425c28e83SPiotr Jasiukajtis x += stridex; 14525c28e83SPiotr Jasiukajtis y += stridey; 14625c28e83SPiotr Jasiukajtis i = 2; 14725c28e83SPiotr Jasiukajtis if (--n <= 0) 14825c28e83SPiotr Jasiukajtis break; 14925c28e83SPiotr Jasiukajtis 15025c28e83SPiotr Jasiukajtis LOOP2: 15125c28e83SPiotr Jasiukajtis xsb2 = HI(x); 15225c28e83SPiotr Jasiukajtis hx2 = xsb2 & ~0x80000000; 15325c28e83SPiotr Jasiukajtis if (hx2 > 0x3fe921fb) 15425c28e83SPiotr Jasiukajtis { 15525c28e83SPiotr Jasiukajtis biguns = 3; 15625c28e83SPiotr Jasiukajtis goto MEDIUM; 15725c28e83SPiotr Jasiukajtis } 15825c28e83SPiotr Jasiukajtis if (hx2 < 0x3e400000) 15925c28e83SPiotr Jasiukajtis { 16025c28e83SPiotr Jasiukajtis v = *x; 16125c28e83SPiotr Jasiukajtis *y = *x; 16225c28e83SPiotr Jasiukajtis x += stridex; 16325c28e83SPiotr Jasiukajtis y += stridey; 16425c28e83SPiotr Jasiukajtis i = 2; 16525c28e83SPiotr Jasiukajtis if (--n <= 0) 16625c28e83SPiotr Jasiukajtis break; 16725c28e83SPiotr Jasiukajtis goto LOOP2; 16825c28e83SPiotr Jasiukajtis } 16925c28e83SPiotr Jasiukajtis x2 = *x; 17025c28e83SPiotr Jasiukajtis py2 = y; 17125c28e83SPiotr Jasiukajtis 17225c28e83SPiotr Jasiukajtis i = (hx0 - 0x3fc90000) >> 31; 17325c28e83SPiotr Jasiukajtis i |= ((hx1 - 0x3fc90000) >> 30) & 2; 17425c28e83SPiotr Jasiukajtis i |= ((hx2 - 0x3fc90000) >> 29) & 4; 17525c28e83SPiotr Jasiukajtis switch (i) 17625c28e83SPiotr Jasiukajtis { 17725c28e83SPiotr Jasiukajtis double a0, a1, a2, w0, w1, w2; 17825c28e83SPiotr Jasiukajtis double t0, t1, t2, z0, z1, z2; 17925c28e83SPiotr Jasiukajtis unsigned j0, j1, j2; 18025c28e83SPiotr Jasiukajtis 18125c28e83SPiotr Jasiukajtis case 0: 18225c28e83SPiotr Jasiukajtis j0 = (xsb0 + 0x4000) & 0xffff8000; 18325c28e83SPiotr Jasiukajtis j1 = (xsb1 + 0x4000) & 0xffff8000; 18425c28e83SPiotr Jasiukajtis j2 = (xsb2 + 0x4000) & 0xffff8000; 18525c28e83SPiotr Jasiukajtis HI(&t0) = j0; 18625c28e83SPiotr Jasiukajtis HI(&t1) = j1; 18725c28e83SPiotr Jasiukajtis HI(&t2) = j2; 18825c28e83SPiotr Jasiukajtis LO(&t0) = 0; 18925c28e83SPiotr Jasiukajtis LO(&t1) = 0; 19025c28e83SPiotr Jasiukajtis LO(&t2) = 0; 19125c28e83SPiotr Jasiukajtis x0 -= t0; 19225c28e83SPiotr Jasiukajtis x1 -= t1; 19325c28e83SPiotr Jasiukajtis x2 -= t2; 19425c28e83SPiotr Jasiukajtis z0 = x0 * x0; 19525c28e83SPiotr Jasiukajtis z1 = x1 * x1; 19625c28e83SPiotr Jasiukajtis z2 = x2 * x2; 19725c28e83SPiotr Jasiukajtis t0 = z0 * (qq1 + z0 * qq2); 19825c28e83SPiotr Jasiukajtis t1 = z1 * (qq1 + z1 * qq2); 19925c28e83SPiotr Jasiukajtis t2 = z2 * (qq1 + z2 * qq2); 20025c28e83SPiotr Jasiukajtis w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 20125c28e83SPiotr Jasiukajtis w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 20225c28e83SPiotr Jasiukajtis w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 20325c28e83SPiotr Jasiukajtis j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 20425c28e83SPiotr Jasiukajtis j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 20525c28e83SPiotr Jasiukajtis j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 20625c28e83SPiotr Jasiukajtis xsb0 = (xsb0 >> 30) & 2; 20725c28e83SPiotr Jasiukajtis xsb1 = (xsb1 >> 30) & 2; 20825c28e83SPiotr Jasiukajtis xsb2 = (xsb2 >> 30) & 2; 20925c28e83SPiotr Jasiukajtis a0 = __vlibm_TBL_sincos_hi[j0+xsb0]; 21025c28e83SPiotr Jasiukajtis a1 = __vlibm_TBL_sincos_hi[j1+xsb1]; 21125c28e83SPiotr Jasiukajtis a2 = __vlibm_TBL_sincos_hi[j2+xsb2]; 21225c28e83SPiotr Jasiukajtis t0 = (__vlibm_TBL_sincos_hi[j0+1] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+xsb0]; 21325c28e83SPiotr Jasiukajtis t1 = (__vlibm_TBL_sincos_hi[j1+1] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+xsb1]; 21425c28e83SPiotr Jasiukajtis t2 = (__vlibm_TBL_sincos_hi[j2+1] * w2 + a2 * t2) + __vlibm_TBL_sincos_lo[j2+xsb2]; 21525c28e83SPiotr Jasiukajtis *py0 = a0 + t0; 21625c28e83SPiotr Jasiukajtis *py1 = a1 + t1; 21725c28e83SPiotr Jasiukajtis *py2 = a2 + t2; 21825c28e83SPiotr Jasiukajtis break; 21925c28e83SPiotr Jasiukajtis 22025c28e83SPiotr Jasiukajtis case 1: 22125c28e83SPiotr Jasiukajtis j1 = (xsb1 + 0x4000) & 0xffff8000; 22225c28e83SPiotr Jasiukajtis j2 = (xsb2 + 0x4000) & 0xffff8000; 22325c28e83SPiotr Jasiukajtis HI(&t1) = j1; 22425c28e83SPiotr Jasiukajtis HI(&t2) = j2; 22525c28e83SPiotr Jasiukajtis LO(&t1) = 0; 22625c28e83SPiotr Jasiukajtis LO(&t2) = 0; 22725c28e83SPiotr Jasiukajtis x1 -= t1; 22825c28e83SPiotr Jasiukajtis x2 -= t2; 22925c28e83SPiotr Jasiukajtis z0 = x0 * x0; 23025c28e83SPiotr Jasiukajtis z1 = x1 * x1; 23125c28e83SPiotr Jasiukajtis z2 = x2 * x2; 23225c28e83SPiotr Jasiukajtis t0 = z0 * (poly3[0] + z0 * poly4[0]); 23325c28e83SPiotr Jasiukajtis t1 = z1 * (qq1 + z1 * qq2); 23425c28e83SPiotr Jasiukajtis t2 = z2 * (qq1 + z2 * qq2); 23525c28e83SPiotr Jasiukajtis t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0)); 23625c28e83SPiotr Jasiukajtis w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 23725c28e83SPiotr Jasiukajtis w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 23825c28e83SPiotr Jasiukajtis j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 23925c28e83SPiotr Jasiukajtis j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 24025c28e83SPiotr Jasiukajtis xsb1 = (xsb1 >> 30) & 2; 24125c28e83SPiotr Jasiukajtis xsb2 = (xsb2 >> 30) & 2; 24225c28e83SPiotr Jasiukajtis a1 = __vlibm_TBL_sincos_hi[j1+xsb1]; 24325c28e83SPiotr Jasiukajtis a2 = __vlibm_TBL_sincos_hi[j2+xsb2]; 24425c28e83SPiotr Jasiukajtis t0 = x0 + x0 * t0; 24525c28e83SPiotr Jasiukajtis t1 = (__vlibm_TBL_sincos_hi[j1+1] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+xsb1]; 24625c28e83SPiotr Jasiukajtis t2 = (__vlibm_TBL_sincos_hi[j2+1] * w2 + a2 * t2) + __vlibm_TBL_sincos_lo[j2+xsb2]; 24725c28e83SPiotr Jasiukajtis *py0 = t0; 24825c28e83SPiotr Jasiukajtis *py1 = a1 + t1; 24925c28e83SPiotr Jasiukajtis *py2 = a2 + t2; 25025c28e83SPiotr Jasiukajtis break; 25125c28e83SPiotr Jasiukajtis 25225c28e83SPiotr Jasiukajtis case 2: 25325c28e83SPiotr Jasiukajtis j0 = (xsb0 + 0x4000) & 0xffff8000; 25425c28e83SPiotr Jasiukajtis j2 = (xsb2 + 0x4000) & 0xffff8000; 25525c28e83SPiotr Jasiukajtis HI(&t0) = j0; 25625c28e83SPiotr Jasiukajtis HI(&t2) = j2; 25725c28e83SPiotr Jasiukajtis LO(&t0) = 0; 25825c28e83SPiotr Jasiukajtis LO(&t2) = 0; 25925c28e83SPiotr Jasiukajtis x0 -= t0; 26025c28e83SPiotr Jasiukajtis x2 -= t2; 26125c28e83SPiotr Jasiukajtis z0 = x0 * x0; 26225c28e83SPiotr Jasiukajtis z1 = x1 * x1; 26325c28e83SPiotr Jasiukajtis z2 = x2 * x2; 26425c28e83SPiotr Jasiukajtis t0 = z0 * (qq1 + z0 * qq2); 26525c28e83SPiotr Jasiukajtis t1 = z1 * (poly3[0] + z1 * poly4[0]); 26625c28e83SPiotr Jasiukajtis t2 = z2 * (qq1 + z2 * qq2); 26725c28e83SPiotr Jasiukajtis w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 26825c28e83SPiotr Jasiukajtis t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1)); 26925c28e83SPiotr Jasiukajtis w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 27025c28e83SPiotr Jasiukajtis j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 27125c28e83SPiotr Jasiukajtis j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 27225c28e83SPiotr Jasiukajtis xsb0 = (xsb0 >> 30) & 2; 27325c28e83SPiotr Jasiukajtis xsb2 = (xsb2 >> 30) & 2; 27425c28e83SPiotr Jasiukajtis a0 = __vlibm_TBL_sincos_hi[j0+xsb0]; 27525c28e83SPiotr Jasiukajtis a2 = __vlibm_TBL_sincos_hi[j2+xsb2]; 27625c28e83SPiotr Jasiukajtis t0 = (__vlibm_TBL_sincos_hi[j0+1] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+xsb0]; 27725c28e83SPiotr Jasiukajtis t1 = x1 + x1 * t1; 27825c28e83SPiotr Jasiukajtis t2 = (__vlibm_TBL_sincos_hi[j2+1] * w2 + a2 * t2) + __vlibm_TBL_sincos_lo[j2+xsb2]; 27925c28e83SPiotr Jasiukajtis *py0 = a0 + t0; 28025c28e83SPiotr Jasiukajtis *py1 = t1; 28125c28e83SPiotr Jasiukajtis *py2 = a2 + t2; 28225c28e83SPiotr Jasiukajtis break; 28325c28e83SPiotr Jasiukajtis 28425c28e83SPiotr Jasiukajtis case 3: 28525c28e83SPiotr Jasiukajtis j2 = (xsb2 + 0x4000) & 0xffff8000; 28625c28e83SPiotr Jasiukajtis HI(&t2) = j2; 28725c28e83SPiotr Jasiukajtis LO(&t2) = 0; 28825c28e83SPiotr Jasiukajtis x2 -= t2; 28925c28e83SPiotr Jasiukajtis z0 = x0 * x0; 29025c28e83SPiotr Jasiukajtis z1 = x1 * x1; 29125c28e83SPiotr Jasiukajtis z2 = x2 * x2; 29225c28e83SPiotr Jasiukajtis t0 = z0 * (poly3[0] + z0 * poly4[0]); 29325c28e83SPiotr Jasiukajtis t1 = z1 * (poly3[0] + z1 * poly4[0]); 29425c28e83SPiotr Jasiukajtis t2 = z2 * (qq1 + z2 * qq2); 29525c28e83SPiotr Jasiukajtis t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0)); 29625c28e83SPiotr Jasiukajtis t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1)); 29725c28e83SPiotr Jasiukajtis w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 29825c28e83SPiotr Jasiukajtis j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 29925c28e83SPiotr Jasiukajtis xsb2 = (xsb2 >> 30) & 2; 30025c28e83SPiotr Jasiukajtis a2 = __vlibm_TBL_sincos_hi[j2+xsb2]; 30125c28e83SPiotr Jasiukajtis t0 = x0 + x0 * t0; 30225c28e83SPiotr Jasiukajtis t1 = x1 + x1 * t1; 30325c28e83SPiotr Jasiukajtis t2 = (__vlibm_TBL_sincos_hi[j2+1] * w2 + a2 * t2) + __vlibm_TBL_sincos_lo[j2+xsb2]; 30425c28e83SPiotr Jasiukajtis *py0 = t0; 30525c28e83SPiotr Jasiukajtis *py1 = t1; 30625c28e83SPiotr Jasiukajtis *py2 = a2 + t2; 30725c28e83SPiotr Jasiukajtis break; 30825c28e83SPiotr Jasiukajtis 30925c28e83SPiotr Jasiukajtis case 4: 31025c28e83SPiotr Jasiukajtis j0 = (xsb0 + 0x4000) & 0xffff8000; 31125c28e83SPiotr Jasiukajtis j1 = (xsb1 + 0x4000) & 0xffff8000; 31225c28e83SPiotr Jasiukajtis HI(&t0) = j0; 31325c28e83SPiotr Jasiukajtis HI(&t1) = j1; 31425c28e83SPiotr Jasiukajtis LO(&t0) = 0; 31525c28e83SPiotr Jasiukajtis LO(&t1) = 0; 31625c28e83SPiotr Jasiukajtis x0 -= t0; 31725c28e83SPiotr Jasiukajtis x1 -= t1; 31825c28e83SPiotr Jasiukajtis z0 = x0 * x0; 31925c28e83SPiotr Jasiukajtis z1 = x1 * x1; 32025c28e83SPiotr Jasiukajtis z2 = x2 * x2; 32125c28e83SPiotr Jasiukajtis t0 = z0 * (qq1 + z0 * qq2); 32225c28e83SPiotr Jasiukajtis t1 = z1 * (qq1 + z1 * qq2); 32325c28e83SPiotr Jasiukajtis t2 = z2 * (poly3[0] + z2 * poly4[0]); 32425c28e83SPiotr Jasiukajtis w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 32525c28e83SPiotr Jasiukajtis w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 32625c28e83SPiotr Jasiukajtis t2 = z2 * (poly1[0] + z2 * (poly2[0] + t2)); 32725c28e83SPiotr Jasiukajtis j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 32825c28e83SPiotr Jasiukajtis j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 32925c28e83SPiotr Jasiukajtis xsb0 = (xsb0 >> 30) & 2; 33025c28e83SPiotr Jasiukajtis xsb1 = (xsb1 >> 30) & 2; 33125c28e83SPiotr Jasiukajtis a0 = __vlibm_TBL_sincos_hi[j0+xsb0]; 33225c28e83SPiotr Jasiukajtis a1 = __vlibm_TBL_sincos_hi[j1+xsb1]; 33325c28e83SPiotr Jasiukajtis t0 = (__vlibm_TBL_sincos_hi[j0+1] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+xsb0]; 33425c28e83SPiotr Jasiukajtis t1 = (__vlibm_TBL_sincos_hi[j1+1] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+xsb1]; 33525c28e83SPiotr Jasiukajtis t2 = x2 + x2 * t2; 33625c28e83SPiotr Jasiukajtis *py0 = a0 + t0; 33725c28e83SPiotr Jasiukajtis *py1 = a1 + t1; 33825c28e83SPiotr Jasiukajtis *py2 = t2; 33925c28e83SPiotr Jasiukajtis break; 34025c28e83SPiotr Jasiukajtis 34125c28e83SPiotr Jasiukajtis case 5: 34225c28e83SPiotr Jasiukajtis j1 = (xsb1 + 0x4000) & 0xffff8000; 34325c28e83SPiotr Jasiukajtis HI(&t1) = j1; 34425c28e83SPiotr Jasiukajtis LO(&t1) = 0; 34525c28e83SPiotr Jasiukajtis x1 -= t1; 34625c28e83SPiotr Jasiukajtis z0 = x0 * x0; 34725c28e83SPiotr Jasiukajtis z1 = x1 * x1; 34825c28e83SPiotr Jasiukajtis z2 = x2 * x2; 34925c28e83SPiotr Jasiukajtis t0 = z0 * (poly3[0] + z0 * poly4[0]); 35025c28e83SPiotr Jasiukajtis t1 = z1 * (qq1 + z1 * qq2); 35125c28e83SPiotr Jasiukajtis t2 = z2 * (poly3[0] + z2 * poly4[0]); 35225c28e83SPiotr Jasiukajtis t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0)); 35325c28e83SPiotr Jasiukajtis w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 35425c28e83SPiotr Jasiukajtis t2 = z2 * (poly1[0] + z2 * (poly2[0] + t2)); 35525c28e83SPiotr Jasiukajtis j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 35625c28e83SPiotr Jasiukajtis xsb1 = (xsb1 >> 30) & 2; 35725c28e83SPiotr Jasiukajtis a1 = __vlibm_TBL_sincos_hi[j1+xsb1]; 35825c28e83SPiotr Jasiukajtis t0 = x0 + x0 * t0; 35925c28e83SPiotr Jasiukajtis t1 = (__vlibm_TBL_sincos_hi[j1+1] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+xsb1]; 36025c28e83SPiotr Jasiukajtis t2 = x2 + x2 * t2; 36125c28e83SPiotr Jasiukajtis *py0 = t0; 36225c28e83SPiotr Jasiukajtis *py1 = a1 + t1; 36325c28e83SPiotr Jasiukajtis *py2 = t2; 36425c28e83SPiotr Jasiukajtis break; 36525c28e83SPiotr Jasiukajtis 36625c28e83SPiotr Jasiukajtis case 6: 36725c28e83SPiotr Jasiukajtis j0 = (xsb0 + 0x4000) & 0xffff8000; 36825c28e83SPiotr Jasiukajtis HI(&t0) = j0; 36925c28e83SPiotr Jasiukajtis LO(&t0) = 0; 37025c28e83SPiotr Jasiukajtis x0 -= t0; 37125c28e83SPiotr Jasiukajtis z0 = x0 * x0; 37225c28e83SPiotr Jasiukajtis z1 = x1 * x1; 37325c28e83SPiotr Jasiukajtis z2 = x2 * x2; 37425c28e83SPiotr Jasiukajtis t0 = z0 * (qq1 + z0 * qq2); 37525c28e83SPiotr Jasiukajtis t1 = z1 * (poly3[0] + z1 * poly4[0]); 37625c28e83SPiotr Jasiukajtis t2 = z2 * (poly3[0] + z2 * poly4[0]); 37725c28e83SPiotr Jasiukajtis w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 37825c28e83SPiotr Jasiukajtis t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1)); 37925c28e83SPiotr Jasiukajtis t2 = z2 * (poly1[0] + z2 * (poly2[0] + t2)); 38025c28e83SPiotr Jasiukajtis j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 38125c28e83SPiotr Jasiukajtis xsb0 = (xsb0 >> 30) & 2; 38225c28e83SPiotr Jasiukajtis a0 = __vlibm_TBL_sincos_hi[j0+xsb0]; 38325c28e83SPiotr Jasiukajtis t0 = (__vlibm_TBL_sincos_hi[j0+1] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+xsb0]; 38425c28e83SPiotr Jasiukajtis t1 = x1 + x1 * t1; 38525c28e83SPiotr Jasiukajtis t2 = x2 + x2 * t2; 38625c28e83SPiotr Jasiukajtis *py0 = a0 + t0; 38725c28e83SPiotr Jasiukajtis *py1 = t1; 38825c28e83SPiotr Jasiukajtis *py2 = t2; 38925c28e83SPiotr Jasiukajtis break; 39025c28e83SPiotr Jasiukajtis 39125c28e83SPiotr Jasiukajtis case 7: 39225c28e83SPiotr Jasiukajtis z0 = x0 * x0; 39325c28e83SPiotr Jasiukajtis z1 = x1 * x1; 39425c28e83SPiotr Jasiukajtis z2 = x2 * x2; 39525c28e83SPiotr Jasiukajtis t0 = z0 * (poly3[0] + z0 * poly4[0]); 39625c28e83SPiotr Jasiukajtis t1 = z1 * (poly3[0] + z1 * poly4[0]); 39725c28e83SPiotr Jasiukajtis t2 = z2 * (poly3[0] + z2 * poly4[0]); 39825c28e83SPiotr Jasiukajtis t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0)); 39925c28e83SPiotr Jasiukajtis t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1)); 40025c28e83SPiotr Jasiukajtis t2 = z2 * (poly1[0] + z2 * (poly2[0] + t2)); 40125c28e83SPiotr Jasiukajtis t0 = x0 + x0 * t0; 40225c28e83SPiotr Jasiukajtis t1 = x1 + x1 * t1; 40325c28e83SPiotr Jasiukajtis t2 = x2 + x2 * t2; 40425c28e83SPiotr Jasiukajtis *py0 = t0; 40525c28e83SPiotr Jasiukajtis *py1 = t1; 40625c28e83SPiotr Jasiukajtis *py2 = t2; 40725c28e83SPiotr Jasiukajtis break; 40825c28e83SPiotr Jasiukajtis } 40925c28e83SPiotr Jasiukajtis 41025c28e83SPiotr Jasiukajtis x += stridex; 41125c28e83SPiotr Jasiukajtis y += stridey; 41225c28e83SPiotr Jasiukajtis i = 0; 41325c28e83SPiotr Jasiukajtis } while (--n > 0); 41425c28e83SPiotr Jasiukajtis 41525c28e83SPiotr Jasiukajtis if (i > 0) 41625c28e83SPiotr Jasiukajtis { 41725c28e83SPiotr Jasiukajtis double a0, a1, w0, w1; 41825c28e83SPiotr Jasiukajtis double t0, t1, z0, z1; 41925c28e83SPiotr Jasiukajtis unsigned j0, j1; 42025c28e83SPiotr Jasiukajtis 42125c28e83SPiotr Jasiukajtis if (i > 1) 42225c28e83SPiotr Jasiukajtis { 42325c28e83SPiotr Jasiukajtis if (hx1 < 0x3fc90000) 42425c28e83SPiotr Jasiukajtis { 42525c28e83SPiotr Jasiukajtis z1 = x1 * x1; 42625c28e83SPiotr Jasiukajtis t1 = z1 * (poly3[0] + z1 * poly4[0]); 42725c28e83SPiotr Jasiukajtis t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1)); 42825c28e83SPiotr Jasiukajtis t1 = x1 + x1 * t1; 42925c28e83SPiotr Jasiukajtis *py1 = t1; 43025c28e83SPiotr Jasiukajtis } 43125c28e83SPiotr Jasiukajtis else 43225c28e83SPiotr Jasiukajtis { 43325c28e83SPiotr Jasiukajtis j1 = (xsb1 + 0x4000) & 0xffff8000; 43425c28e83SPiotr Jasiukajtis HI(&t1) = j1; 43525c28e83SPiotr Jasiukajtis LO(&t1) = 0; 43625c28e83SPiotr Jasiukajtis x1 -= t1; 43725c28e83SPiotr Jasiukajtis z1 = x1 * x1; 43825c28e83SPiotr Jasiukajtis t1 = z1 * (qq1 + z1 * qq2); 43925c28e83SPiotr Jasiukajtis w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 44025c28e83SPiotr Jasiukajtis j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 44125c28e83SPiotr Jasiukajtis xsb1 = (xsb1 >> 30) & 2; 44225c28e83SPiotr Jasiukajtis a1 = __vlibm_TBL_sincos_hi[j1+xsb1]; 44325c28e83SPiotr Jasiukajtis t1 = (__vlibm_TBL_sincos_hi[j1+1] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+xsb1]; 44425c28e83SPiotr Jasiukajtis *py1 = a1 + t1; 44525c28e83SPiotr Jasiukajtis } 44625c28e83SPiotr Jasiukajtis } 44725c28e83SPiotr Jasiukajtis if (hx0 < 0x3fc90000) 44825c28e83SPiotr Jasiukajtis { 44925c28e83SPiotr Jasiukajtis z0 = x0 * x0; 45025c28e83SPiotr Jasiukajtis t0 = z0 * (poly3[0] + z0 * poly4[0]); 45125c28e83SPiotr Jasiukajtis t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0)); 45225c28e83SPiotr Jasiukajtis t0 = x0 + x0 * t0; 45325c28e83SPiotr Jasiukajtis *py0 = t0; 45425c28e83SPiotr Jasiukajtis } 45525c28e83SPiotr Jasiukajtis else 45625c28e83SPiotr Jasiukajtis { 45725c28e83SPiotr Jasiukajtis j0 = (xsb0 + 0x4000) & 0xffff8000; 45825c28e83SPiotr Jasiukajtis HI(&t0) = j0; 45925c28e83SPiotr Jasiukajtis LO(&t0) = 0; 46025c28e83SPiotr Jasiukajtis x0 -= t0; 46125c28e83SPiotr Jasiukajtis z0 = x0 * x0; 46225c28e83SPiotr Jasiukajtis t0 = z0 * (qq1 + z0 * qq2); 46325c28e83SPiotr Jasiukajtis w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 46425c28e83SPiotr Jasiukajtis j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 46525c28e83SPiotr Jasiukajtis xsb0 = (xsb0 >> 30) & 2; 46625c28e83SPiotr Jasiukajtis a0 = __vlibm_TBL_sincos_hi[j0+xsb0]; 46725c28e83SPiotr Jasiukajtis t0 = (__vlibm_TBL_sincos_hi[j0+1] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+xsb0]; 46825c28e83SPiotr Jasiukajtis *py0 = a0 + t0; 46925c28e83SPiotr Jasiukajtis } 47025c28e83SPiotr Jasiukajtis } 47125c28e83SPiotr Jasiukajtis 47225c28e83SPiotr Jasiukajtis return; 47325c28e83SPiotr Jasiukajtis 47425c28e83SPiotr Jasiukajtis /* 47525c28e83SPiotr Jasiukajtis * MEDIUM RANGE PROCESSING 47625c28e83SPiotr Jasiukajtis * Jump here at first sign of medium range argument. We are a bit 47725c28e83SPiotr Jasiukajtis * confused due to the jump.. fix up several variables and jump into 47825c28e83SPiotr Jasiukajtis * the nth loop, same as was being processed above. 47925c28e83SPiotr Jasiukajtis */ 48025c28e83SPiotr Jasiukajtis 48125c28e83SPiotr Jasiukajtis MEDIUM: 48225c28e83SPiotr Jasiukajtis 48325c28e83SPiotr Jasiukajtis x0_or_one[1] = 1.0; 48425c28e83SPiotr Jasiukajtis x1_or_one[1] = 1.0; 48525c28e83SPiotr Jasiukajtis x2_or_one[1] = 1.0; 48625c28e83SPiotr Jasiukajtis x0_or_one[3] = -1.0; 48725c28e83SPiotr Jasiukajtis x1_or_one[3] = -1.0; 48825c28e83SPiotr Jasiukajtis x2_or_one[3] = -1.0; 48925c28e83SPiotr Jasiukajtis y0_or_zero[1] = 0.0; 49025c28e83SPiotr Jasiukajtis y1_or_zero[1] = 0.0; 49125c28e83SPiotr Jasiukajtis y2_or_zero[1] = 0.0; 49225c28e83SPiotr Jasiukajtis y0_or_zero[3] = 0.0; 49325c28e83SPiotr Jasiukajtis y1_or_zero[3] = 0.0; 49425c28e83SPiotr Jasiukajtis y2_or_zero[3] = 0.0; 49525c28e83SPiotr Jasiukajtis 49625c28e83SPiotr Jasiukajtis if (biguns == 3) 49725c28e83SPiotr Jasiukajtis { 49825c28e83SPiotr Jasiukajtis biguns = 0; 49925c28e83SPiotr Jasiukajtis xsb0 = xsb0 >> 31; 50025c28e83SPiotr Jasiukajtis xsb1 = xsb1 >> 31; 50125c28e83SPiotr Jasiukajtis goto loop2; 50225c28e83SPiotr Jasiukajtis } 50325c28e83SPiotr Jasiukajtis else if (biguns == 2) 50425c28e83SPiotr Jasiukajtis { 50525c28e83SPiotr Jasiukajtis xsb0 = xsb0 >> 31; 50625c28e83SPiotr Jasiukajtis biguns = 0; 50725c28e83SPiotr Jasiukajtis goto loop1; 50825c28e83SPiotr Jasiukajtis } 50925c28e83SPiotr Jasiukajtis biguns = 0; 51025c28e83SPiotr Jasiukajtis 51125c28e83SPiotr Jasiukajtis do 51225c28e83SPiotr Jasiukajtis { 51325c28e83SPiotr Jasiukajtis double fn0, fn1, fn2, a0, a1, a2, w0, w1, w2, y0, y1, y2; 51425c28e83SPiotr Jasiukajtis unsigned hx; 51525c28e83SPiotr Jasiukajtis int n0, n1, n2; 51625c28e83SPiotr Jasiukajtis 51725c28e83SPiotr Jasiukajtis loop0: 51825c28e83SPiotr Jasiukajtis hx = HI(x); 51925c28e83SPiotr Jasiukajtis xsb0 = hx >> 31; 52025c28e83SPiotr Jasiukajtis hx &= ~0x80000000; 52125c28e83SPiotr Jasiukajtis if (hx < 0x3e400000) 52225c28e83SPiotr Jasiukajtis { 52325c28e83SPiotr Jasiukajtis v = *x; 52425c28e83SPiotr Jasiukajtis *y = *x; 52525c28e83SPiotr Jasiukajtis x += stridex; 52625c28e83SPiotr Jasiukajtis y += stridey; 52725c28e83SPiotr Jasiukajtis i = 0; 52825c28e83SPiotr Jasiukajtis if (--n <= 0) 52925c28e83SPiotr Jasiukajtis break; 53025c28e83SPiotr Jasiukajtis goto loop0; 53125c28e83SPiotr Jasiukajtis } 53225c28e83SPiotr Jasiukajtis if (hx > 0x413921fb) 53325c28e83SPiotr Jasiukajtis { 53425c28e83SPiotr Jasiukajtis if (hx >= 0x7ff00000) 53525c28e83SPiotr Jasiukajtis { 53625c28e83SPiotr Jasiukajtis x0 = *x; 53725c28e83SPiotr Jasiukajtis *y = x0 - x0; 53825c28e83SPiotr Jasiukajtis } 53925c28e83SPiotr Jasiukajtis else 54025c28e83SPiotr Jasiukajtis biguns = 1; 54125c28e83SPiotr Jasiukajtis x += stridex; 54225c28e83SPiotr Jasiukajtis y += stridey; 54325c28e83SPiotr Jasiukajtis i = 0; 54425c28e83SPiotr Jasiukajtis if (--n <= 0) 54525c28e83SPiotr Jasiukajtis break; 54625c28e83SPiotr Jasiukajtis goto loop0; 54725c28e83SPiotr Jasiukajtis } 54825c28e83SPiotr Jasiukajtis x0 = *x; 54925c28e83SPiotr Jasiukajtis py0 = y; 55025c28e83SPiotr Jasiukajtis x += stridex; 55125c28e83SPiotr Jasiukajtis y += stridey; 55225c28e83SPiotr Jasiukajtis i = 1; 55325c28e83SPiotr Jasiukajtis if (--n <= 0) 55425c28e83SPiotr Jasiukajtis break; 55525c28e83SPiotr Jasiukajtis 55625c28e83SPiotr Jasiukajtis loop1: 55725c28e83SPiotr Jasiukajtis hx = HI(x); 55825c28e83SPiotr Jasiukajtis xsb1 = hx >> 31; 55925c28e83SPiotr Jasiukajtis hx &= ~0x80000000; 56025c28e83SPiotr Jasiukajtis if (hx < 0x3e400000) 56125c28e83SPiotr Jasiukajtis { 56225c28e83SPiotr Jasiukajtis v = *x; 56325c28e83SPiotr Jasiukajtis *y = *x; 56425c28e83SPiotr Jasiukajtis x += stridex; 56525c28e83SPiotr Jasiukajtis y += stridey; 56625c28e83SPiotr Jasiukajtis i = 1; 56725c28e83SPiotr Jasiukajtis if (--n <= 0) 56825c28e83SPiotr Jasiukajtis break; 56925c28e83SPiotr Jasiukajtis goto loop1; 57025c28e83SPiotr Jasiukajtis } 57125c28e83SPiotr Jasiukajtis if (hx > 0x413921fb) 57225c28e83SPiotr Jasiukajtis { 57325c28e83SPiotr Jasiukajtis if (hx >= 0x7ff00000) 57425c28e83SPiotr Jasiukajtis { 57525c28e83SPiotr Jasiukajtis x1 = *x; 57625c28e83SPiotr Jasiukajtis *y = x1 - x1; 57725c28e83SPiotr Jasiukajtis } 57825c28e83SPiotr Jasiukajtis else 57925c28e83SPiotr Jasiukajtis biguns = 1; 58025c28e83SPiotr Jasiukajtis x += stridex; 58125c28e83SPiotr Jasiukajtis y += stridey; 58225c28e83SPiotr Jasiukajtis i = 1; 58325c28e83SPiotr Jasiukajtis if (--n <= 0) 58425c28e83SPiotr Jasiukajtis break; 58525c28e83SPiotr Jasiukajtis goto loop1; 58625c28e83SPiotr Jasiukajtis } 58725c28e83SPiotr Jasiukajtis x1 = *x; 58825c28e83SPiotr Jasiukajtis py1 = y; 58925c28e83SPiotr Jasiukajtis x += stridex; 59025c28e83SPiotr Jasiukajtis y += stridey; 59125c28e83SPiotr Jasiukajtis i = 2; 59225c28e83SPiotr Jasiukajtis if (--n <= 0) 59325c28e83SPiotr Jasiukajtis break; 59425c28e83SPiotr Jasiukajtis 59525c28e83SPiotr Jasiukajtis loop2: 59625c28e83SPiotr Jasiukajtis hx = HI(x); 59725c28e83SPiotr Jasiukajtis xsb2 = hx >> 31; 59825c28e83SPiotr Jasiukajtis hx &= ~0x80000000; 59925c28e83SPiotr Jasiukajtis if (hx < 0x3e400000) 60025c28e83SPiotr Jasiukajtis { 60125c28e83SPiotr Jasiukajtis v = *x; 60225c28e83SPiotr Jasiukajtis *y = *x; 60325c28e83SPiotr Jasiukajtis x += stridex; 60425c28e83SPiotr Jasiukajtis y += stridey; 60525c28e83SPiotr Jasiukajtis i = 2; 60625c28e83SPiotr Jasiukajtis if (--n <= 0) 60725c28e83SPiotr Jasiukajtis break; 60825c28e83SPiotr Jasiukajtis goto loop2; 60925c28e83SPiotr Jasiukajtis } 61025c28e83SPiotr Jasiukajtis if (hx > 0x413921fb) 61125c28e83SPiotr Jasiukajtis { 61225c28e83SPiotr Jasiukajtis if (hx >= 0x7ff00000) 61325c28e83SPiotr Jasiukajtis { 61425c28e83SPiotr Jasiukajtis x2 = *x; 61525c28e83SPiotr Jasiukajtis *y = x2 - x2; 61625c28e83SPiotr Jasiukajtis } 61725c28e83SPiotr Jasiukajtis else 61825c28e83SPiotr Jasiukajtis biguns = 1; 61925c28e83SPiotr Jasiukajtis x += stridex; 62025c28e83SPiotr Jasiukajtis y += stridey; 62125c28e83SPiotr Jasiukajtis i = 2; 62225c28e83SPiotr Jasiukajtis if (--n <= 0) 62325c28e83SPiotr Jasiukajtis break; 62425c28e83SPiotr Jasiukajtis goto loop2; 62525c28e83SPiotr Jasiukajtis } 62625c28e83SPiotr Jasiukajtis x2 = *x; 62725c28e83SPiotr Jasiukajtis py2 = y; 62825c28e83SPiotr Jasiukajtis 62925c28e83SPiotr Jasiukajtis n0 = (int) (x0 * invpio2 + half[xsb0]); 63025c28e83SPiotr Jasiukajtis n1 = (int) (x1 * invpio2 + half[xsb1]); 63125c28e83SPiotr Jasiukajtis n2 = (int) (x2 * invpio2 + half[xsb2]); 63225c28e83SPiotr Jasiukajtis fn0 = (double) n0; 63325c28e83SPiotr Jasiukajtis fn1 = (double) n1; 63425c28e83SPiotr Jasiukajtis fn2 = (double) n2; 63525c28e83SPiotr Jasiukajtis n0 &= 3; 63625c28e83SPiotr Jasiukajtis n1 &= 3; 63725c28e83SPiotr Jasiukajtis n2 &= 3; 63825c28e83SPiotr Jasiukajtis a0 = x0 - fn0 * pio2_1; 63925c28e83SPiotr Jasiukajtis a1 = x1 - fn1 * pio2_1; 64025c28e83SPiotr Jasiukajtis a2 = x2 - fn2 * pio2_1; 64125c28e83SPiotr Jasiukajtis w0 = fn0 * pio2_2; 64225c28e83SPiotr Jasiukajtis w1 = fn1 * pio2_2; 64325c28e83SPiotr Jasiukajtis w2 = fn2 * pio2_2; 64425c28e83SPiotr Jasiukajtis x0 = a0 - w0; 64525c28e83SPiotr Jasiukajtis x1 = a1 - w1; 64625c28e83SPiotr Jasiukajtis x2 = a2 - w2; 64725c28e83SPiotr Jasiukajtis y0 = (a0 - x0) - w0; 64825c28e83SPiotr Jasiukajtis y1 = (a1 - x1) - w1; 64925c28e83SPiotr Jasiukajtis y2 = (a2 - x2) - w2; 65025c28e83SPiotr Jasiukajtis a0 = x0; 65125c28e83SPiotr Jasiukajtis a1 = x1; 65225c28e83SPiotr Jasiukajtis a2 = x2; 65325c28e83SPiotr Jasiukajtis w0 = fn0 * pio2_3 - y0; 65425c28e83SPiotr Jasiukajtis w1 = fn1 * pio2_3 - y1; 65525c28e83SPiotr Jasiukajtis w2 = fn2 * pio2_3 - y2; 65625c28e83SPiotr Jasiukajtis x0 = a0 - w0; 65725c28e83SPiotr Jasiukajtis x1 = a1 - w1; 65825c28e83SPiotr Jasiukajtis x2 = a2 - w2; 65925c28e83SPiotr Jasiukajtis y0 = (a0 - x0) - w0; 66025c28e83SPiotr Jasiukajtis y1 = (a1 - x1) - w1; 66125c28e83SPiotr Jasiukajtis y2 = (a2 - x2) - w2; 66225c28e83SPiotr Jasiukajtis a0 = x0; 66325c28e83SPiotr Jasiukajtis a1 = x1; 66425c28e83SPiotr Jasiukajtis a2 = x2; 66525c28e83SPiotr Jasiukajtis w0 = fn0 * pio2_3t - y0; 66625c28e83SPiotr Jasiukajtis w1 = fn1 * pio2_3t - y1; 66725c28e83SPiotr Jasiukajtis w2 = fn2 * pio2_3t - y2; 66825c28e83SPiotr Jasiukajtis x0 = a0 - w0; 66925c28e83SPiotr Jasiukajtis x1 = a1 - w1; 67025c28e83SPiotr Jasiukajtis x2 = a2 - w2; 67125c28e83SPiotr Jasiukajtis y0 = (a0 - x0) - w0; 67225c28e83SPiotr Jasiukajtis y1 = (a1 - x1) - w1; 67325c28e83SPiotr Jasiukajtis y2 = (a2 - x2) - w2; 67425c28e83SPiotr Jasiukajtis xsb0 = HI(&x0); 67525c28e83SPiotr Jasiukajtis i = ((xsb0 & ~0x80000000) - thresh[n0&1]) >> 31; 67625c28e83SPiotr Jasiukajtis xsb1 = HI(&x1); 67725c28e83SPiotr Jasiukajtis i |= (((xsb1 & ~0x80000000) - thresh[n1&1]) >> 30) & 2; 67825c28e83SPiotr Jasiukajtis xsb2 = HI(&x2); 67925c28e83SPiotr Jasiukajtis i |= (((xsb2 & ~0x80000000) - thresh[n2&1]) >> 29) & 4; 68025c28e83SPiotr Jasiukajtis switch (i) 68125c28e83SPiotr Jasiukajtis { 68225c28e83SPiotr Jasiukajtis double t0, t1, t2, z0, z1, z2; 68325c28e83SPiotr Jasiukajtis unsigned j0, j1, j2; 68425c28e83SPiotr Jasiukajtis 68525c28e83SPiotr Jasiukajtis case 0: 68625c28e83SPiotr Jasiukajtis j0 = (xsb0 + 0x4000) & 0xffff8000; 68725c28e83SPiotr Jasiukajtis j1 = (xsb1 + 0x4000) & 0xffff8000; 68825c28e83SPiotr Jasiukajtis j2 = (xsb2 + 0x4000) & 0xffff8000; 68925c28e83SPiotr Jasiukajtis HI(&t0) = j0; 69025c28e83SPiotr Jasiukajtis HI(&t1) = j1; 69125c28e83SPiotr Jasiukajtis HI(&t2) = j2; 69225c28e83SPiotr Jasiukajtis LO(&t0) = 0; 69325c28e83SPiotr Jasiukajtis LO(&t1) = 0; 69425c28e83SPiotr Jasiukajtis LO(&t2) = 0; 69525c28e83SPiotr Jasiukajtis x0 = (x0 - t0) + y0; 69625c28e83SPiotr Jasiukajtis x1 = (x1 - t1) + y1; 69725c28e83SPiotr Jasiukajtis x2 = (x2 - t2) + y2; 69825c28e83SPiotr Jasiukajtis z0 = x0 * x0; 69925c28e83SPiotr Jasiukajtis z1 = x1 * x1; 70025c28e83SPiotr Jasiukajtis z2 = x2 * x2; 70125c28e83SPiotr Jasiukajtis t0 = z0 * (qq1 + z0 * qq2); 70225c28e83SPiotr Jasiukajtis t1 = z1 * (qq1 + z1 * qq2); 70325c28e83SPiotr Jasiukajtis t2 = z2 * (qq1 + z2 * qq2); 70425c28e83SPiotr Jasiukajtis w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 70525c28e83SPiotr Jasiukajtis w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 70625c28e83SPiotr Jasiukajtis w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 70725c28e83SPiotr Jasiukajtis j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 70825c28e83SPiotr Jasiukajtis j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 70925c28e83SPiotr Jasiukajtis j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 71025c28e83SPiotr Jasiukajtis xsb0 = (xsb0 >> 30) & 2; 71125c28e83SPiotr Jasiukajtis xsb1 = (xsb1 >> 30) & 2; 71225c28e83SPiotr Jasiukajtis xsb2 = (xsb2 >> 30) & 2; 71325c28e83SPiotr Jasiukajtis n0 ^= (xsb0 & ~(n0 << 1)); 71425c28e83SPiotr Jasiukajtis n1 ^= (xsb1 & ~(n1 << 1)); 71525c28e83SPiotr Jasiukajtis n2 ^= (xsb2 & ~(n2 << 1)); 71625c28e83SPiotr Jasiukajtis xsb0 |= 1; 71725c28e83SPiotr Jasiukajtis xsb1 |= 1; 71825c28e83SPiotr Jasiukajtis xsb2 |= 1; 71925c28e83SPiotr Jasiukajtis a0 = __vlibm_TBL_sincos_hi[j0+n0]; 72025c28e83SPiotr Jasiukajtis a1 = __vlibm_TBL_sincos_hi[j1+n1]; 72125c28e83SPiotr Jasiukajtis a2 = __vlibm_TBL_sincos_hi[j2+n2]; 72225c28e83SPiotr Jasiukajtis t0 = (__vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+n0]; 72325c28e83SPiotr Jasiukajtis t1 = (__vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+n1]; 72425c28e83SPiotr Jasiukajtis t2 = (__vlibm_TBL_sincos_hi[j2+((n2+xsb2)&3)] * w2 + a2 * t2) + __vlibm_TBL_sincos_lo[j2+n2]; 72525c28e83SPiotr Jasiukajtis *py0 = ( a0 + t0 ); 72625c28e83SPiotr Jasiukajtis *py1 = ( a1 + t1 ); 72725c28e83SPiotr Jasiukajtis *py2 = ( a2 + t2 ); 72825c28e83SPiotr Jasiukajtis break; 72925c28e83SPiotr Jasiukajtis 73025c28e83SPiotr Jasiukajtis case 1: 73125c28e83SPiotr Jasiukajtis j0 = n0 & 1; 73225c28e83SPiotr Jasiukajtis j1 = (xsb1 + 0x4000) & 0xffff8000; 73325c28e83SPiotr Jasiukajtis j2 = (xsb2 + 0x4000) & 0xffff8000; 73425c28e83SPiotr Jasiukajtis HI(&t1) = j1; 73525c28e83SPiotr Jasiukajtis HI(&t2) = j2; 73625c28e83SPiotr Jasiukajtis LO(&t1) = 0; 73725c28e83SPiotr Jasiukajtis LO(&t2) = 0; 73825c28e83SPiotr Jasiukajtis x0_or_one[0] = x0; 73925c28e83SPiotr Jasiukajtis x0_or_one[2] = -x0; 74025c28e83SPiotr Jasiukajtis y0_or_zero[0] = y0; 74125c28e83SPiotr Jasiukajtis y0_or_zero[2] = -y0; 74225c28e83SPiotr Jasiukajtis x1 = (x1 - t1) + y1; 74325c28e83SPiotr Jasiukajtis x2 = (x2 - t2) + y2; 74425c28e83SPiotr Jasiukajtis z0 = x0 * x0; 74525c28e83SPiotr Jasiukajtis z1 = x1 * x1; 74625c28e83SPiotr Jasiukajtis z2 = x2 * x2; 74725c28e83SPiotr Jasiukajtis t0 = z0 * (poly3[j0] + z0 * poly4[j0]); 74825c28e83SPiotr Jasiukajtis t1 = z1 * (qq1 + z1 * qq2); 74925c28e83SPiotr Jasiukajtis t2 = z2 * (qq1 + z2 * qq2); 75025c28e83SPiotr Jasiukajtis t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0)); 75125c28e83SPiotr Jasiukajtis w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 75225c28e83SPiotr Jasiukajtis w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 75325c28e83SPiotr Jasiukajtis j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 75425c28e83SPiotr Jasiukajtis j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 75525c28e83SPiotr Jasiukajtis xsb1 = (xsb1 >> 30) & 2; 75625c28e83SPiotr Jasiukajtis xsb2 = (xsb2 >> 30) & 2; 75725c28e83SPiotr Jasiukajtis n1 ^= (xsb1 & ~(n1 << 1)); 75825c28e83SPiotr Jasiukajtis n2 ^= (xsb2 & ~(n2 << 1)); 75925c28e83SPiotr Jasiukajtis xsb1 |= 1; 76025c28e83SPiotr Jasiukajtis xsb2 |= 1; 76125c28e83SPiotr Jasiukajtis a1 = __vlibm_TBL_sincos_hi[j1+n1]; 76225c28e83SPiotr Jasiukajtis a2 = __vlibm_TBL_sincos_hi[j2+n2]; 76325c28e83SPiotr Jasiukajtis t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0); 76425c28e83SPiotr Jasiukajtis t1 = (__vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+n1]; 76525c28e83SPiotr Jasiukajtis t2 = (__vlibm_TBL_sincos_hi[j2+((n2+xsb2)&3)] * w2 + a2 * t2) + __vlibm_TBL_sincos_lo[j2+n2]; 76625c28e83SPiotr Jasiukajtis *py0 = t0; 76725c28e83SPiotr Jasiukajtis *py1 = ( a1 + t1 ); 76825c28e83SPiotr Jasiukajtis *py2 = ( a2 + t2 ); 76925c28e83SPiotr Jasiukajtis break; 77025c28e83SPiotr Jasiukajtis 77125c28e83SPiotr Jasiukajtis case 2: 77225c28e83SPiotr Jasiukajtis j0 = (xsb0 + 0x4000) & 0xffff8000; 77325c28e83SPiotr Jasiukajtis j1 = n1 & 1; 77425c28e83SPiotr Jasiukajtis j2 = (xsb2 + 0x4000) & 0xffff8000; 77525c28e83SPiotr Jasiukajtis HI(&t0) = j0; 77625c28e83SPiotr Jasiukajtis HI(&t2) = j2; 77725c28e83SPiotr Jasiukajtis LO(&t0) = 0; 77825c28e83SPiotr Jasiukajtis LO(&t2) = 0; 77925c28e83SPiotr Jasiukajtis x1_or_one[0] = x1; 78025c28e83SPiotr Jasiukajtis x1_or_one[2] = -x1; 78125c28e83SPiotr Jasiukajtis x0 = (x0 - t0) + y0; 78225c28e83SPiotr Jasiukajtis y1_or_zero[0] = y1; 78325c28e83SPiotr Jasiukajtis y1_or_zero[2] = -y1; 78425c28e83SPiotr Jasiukajtis x2 = (x2 - t2) + y2; 78525c28e83SPiotr Jasiukajtis z0 = x0 * x0; 78625c28e83SPiotr Jasiukajtis z1 = x1 * x1; 78725c28e83SPiotr Jasiukajtis z2 = x2 * x2; 78825c28e83SPiotr Jasiukajtis t0 = z0 * (qq1 + z0 * qq2); 78925c28e83SPiotr Jasiukajtis t1 = z1 * (poly3[j1] + z1 * poly4[j1]); 79025c28e83SPiotr Jasiukajtis t2 = z2 * (qq1 + z2 * qq2); 79125c28e83SPiotr Jasiukajtis w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 79225c28e83SPiotr Jasiukajtis t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1)); 79325c28e83SPiotr Jasiukajtis w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 79425c28e83SPiotr Jasiukajtis j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 79525c28e83SPiotr Jasiukajtis j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 79625c28e83SPiotr Jasiukajtis xsb0 = (xsb0 >> 30) & 2; 79725c28e83SPiotr Jasiukajtis xsb2 = (xsb2 >> 30) & 2; 79825c28e83SPiotr Jasiukajtis n0 ^= (xsb0 & ~(n0 << 1)); 79925c28e83SPiotr Jasiukajtis n2 ^= (xsb2 & ~(n2 << 1)); 80025c28e83SPiotr Jasiukajtis xsb0 |= 1; 80125c28e83SPiotr Jasiukajtis xsb2 |= 1; 80225c28e83SPiotr Jasiukajtis a0 = __vlibm_TBL_sincos_hi[j0+n0]; 80325c28e83SPiotr Jasiukajtis a2 = __vlibm_TBL_sincos_hi[j2+n2]; 80425c28e83SPiotr Jasiukajtis t0 = (__vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+n0]; 80525c28e83SPiotr Jasiukajtis t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1); 80625c28e83SPiotr Jasiukajtis t2 = (__vlibm_TBL_sincos_hi[j2+((n2+xsb2)&3)] * w2 + a2 * t2) + __vlibm_TBL_sincos_lo[j2+n2]; 80725c28e83SPiotr Jasiukajtis *py0 = ( a0 + t0 ); 80825c28e83SPiotr Jasiukajtis *py1 = t1; 80925c28e83SPiotr Jasiukajtis *py2 = ( a2 + t2 ); 81025c28e83SPiotr Jasiukajtis break; 81125c28e83SPiotr Jasiukajtis 81225c28e83SPiotr Jasiukajtis case 3: 81325c28e83SPiotr Jasiukajtis j0 = n0 & 1; 81425c28e83SPiotr Jasiukajtis j1 = n1 & 1; 81525c28e83SPiotr Jasiukajtis j2 = (xsb2 + 0x4000) & 0xffff8000; 81625c28e83SPiotr Jasiukajtis HI(&t2) = j2; 81725c28e83SPiotr Jasiukajtis LO(&t2) = 0; 81825c28e83SPiotr Jasiukajtis x0_or_one[0] = x0; 81925c28e83SPiotr Jasiukajtis x0_or_one[2] = -x0; 82025c28e83SPiotr Jasiukajtis x1_or_one[0] = x1; 82125c28e83SPiotr Jasiukajtis x1_or_one[2] = -x1; 82225c28e83SPiotr Jasiukajtis y0_or_zero[0] = y0; 82325c28e83SPiotr Jasiukajtis y0_or_zero[2] = -y0; 82425c28e83SPiotr Jasiukajtis y1_or_zero[0] = y1; 82525c28e83SPiotr Jasiukajtis y1_or_zero[2] = -y1; 82625c28e83SPiotr Jasiukajtis x2 = (x2 - t2) + y2; 82725c28e83SPiotr Jasiukajtis z0 = x0 * x0; 82825c28e83SPiotr Jasiukajtis z1 = x1 * x1; 82925c28e83SPiotr Jasiukajtis z2 = x2 * x2; 83025c28e83SPiotr Jasiukajtis t0 = z0 * (poly3[j0] + z0 * poly4[j0]); 83125c28e83SPiotr Jasiukajtis t1 = z1 * (poly3[j1] + z1 * poly4[j1]); 83225c28e83SPiotr Jasiukajtis t2 = z2 * (qq1 + z2 * qq2); 83325c28e83SPiotr Jasiukajtis t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0)); 83425c28e83SPiotr Jasiukajtis t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1)); 83525c28e83SPiotr Jasiukajtis w2 = x2 * (one + z2 * (pp1 + z2 * pp2)); 83625c28e83SPiotr Jasiukajtis j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 83725c28e83SPiotr Jasiukajtis xsb2 = (xsb2 >> 30) & 2; 83825c28e83SPiotr Jasiukajtis n2 ^= (xsb2 & ~(n2 << 1)); 83925c28e83SPiotr Jasiukajtis xsb2 |= 1; 84025c28e83SPiotr Jasiukajtis a2 = __vlibm_TBL_sincos_hi[j2+n2]; 84125c28e83SPiotr Jasiukajtis t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0); 84225c28e83SPiotr Jasiukajtis t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1); 84325c28e83SPiotr Jasiukajtis t2 = (__vlibm_TBL_sincos_hi[j2+((n2+xsb2)&3)] * w2 + a2 * t2) + __vlibm_TBL_sincos_lo[j2+n2]; 84425c28e83SPiotr Jasiukajtis *py0 = t0; 84525c28e83SPiotr Jasiukajtis *py1 = t1; 84625c28e83SPiotr Jasiukajtis *py2 = ( a2 + t2 ); 84725c28e83SPiotr Jasiukajtis break; 84825c28e83SPiotr Jasiukajtis 84925c28e83SPiotr Jasiukajtis case 4: 85025c28e83SPiotr Jasiukajtis j0 = (xsb0 + 0x4000) & 0xffff8000; 85125c28e83SPiotr Jasiukajtis j1 = (xsb1 + 0x4000) & 0xffff8000; 85225c28e83SPiotr Jasiukajtis j2 = n2 & 1; 85325c28e83SPiotr Jasiukajtis HI(&t0) = j0; 85425c28e83SPiotr Jasiukajtis HI(&t1) = j1; 85525c28e83SPiotr Jasiukajtis LO(&t0) = 0; 85625c28e83SPiotr Jasiukajtis LO(&t1) = 0; 85725c28e83SPiotr Jasiukajtis x2_or_one[0] = x2; 85825c28e83SPiotr Jasiukajtis x2_or_one[2] = -x2; 85925c28e83SPiotr Jasiukajtis x0 = (x0 - t0) + y0; 86025c28e83SPiotr Jasiukajtis x1 = (x1 - t1) + y1; 86125c28e83SPiotr Jasiukajtis y2_or_zero[0] = y2; 86225c28e83SPiotr Jasiukajtis y2_or_zero[2] = -y2; 86325c28e83SPiotr Jasiukajtis z0 = x0 * x0; 86425c28e83SPiotr Jasiukajtis z1 = x1 * x1; 86525c28e83SPiotr Jasiukajtis z2 = x2 * x2; 86625c28e83SPiotr Jasiukajtis t0 = z0 * (qq1 + z0 * qq2); 86725c28e83SPiotr Jasiukajtis t1 = z1 * (qq1 + z1 * qq2); 86825c28e83SPiotr Jasiukajtis t2 = z2 * (poly3[j2] + z2 * poly4[j2]); 86925c28e83SPiotr Jasiukajtis w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 87025c28e83SPiotr Jasiukajtis w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 87125c28e83SPiotr Jasiukajtis t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2)); 87225c28e83SPiotr Jasiukajtis j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 87325c28e83SPiotr Jasiukajtis j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 87425c28e83SPiotr Jasiukajtis xsb0 = (xsb0 >> 30) & 2; 87525c28e83SPiotr Jasiukajtis xsb1 = (xsb1 >> 30) & 2; 87625c28e83SPiotr Jasiukajtis n0 ^= (xsb0 & ~(n0 << 1)); 87725c28e83SPiotr Jasiukajtis n1 ^= (xsb1 & ~(n1 << 1)); 87825c28e83SPiotr Jasiukajtis xsb0 |= 1; 87925c28e83SPiotr Jasiukajtis xsb1 |= 1; 88025c28e83SPiotr Jasiukajtis a0 = __vlibm_TBL_sincos_hi[j0+n0]; 88125c28e83SPiotr Jasiukajtis a1 = __vlibm_TBL_sincos_hi[j1+n1]; 88225c28e83SPiotr Jasiukajtis t0 = (__vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+n0]; 88325c28e83SPiotr Jasiukajtis t1 = (__vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+n1]; 88425c28e83SPiotr Jasiukajtis t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2); 88525c28e83SPiotr Jasiukajtis *py0 = ( a0 + t0 ); 88625c28e83SPiotr Jasiukajtis *py1 = ( a1 + t1 ); 88725c28e83SPiotr Jasiukajtis *py2 = t2; 88825c28e83SPiotr Jasiukajtis break; 88925c28e83SPiotr Jasiukajtis 89025c28e83SPiotr Jasiukajtis case 5: 89125c28e83SPiotr Jasiukajtis j0 = n0 & 1; 89225c28e83SPiotr Jasiukajtis j1 = (xsb1 + 0x4000) & 0xffff8000; 89325c28e83SPiotr Jasiukajtis j2 = n2 & 1; 89425c28e83SPiotr Jasiukajtis HI(&t1) = j1; 89525c28e83SPiotr Jasiukajtis LO(&t1) = 0; 89625c28e83SPiotr Jasiukajtis x0_or_one[0] = x0; 89725c28e83SPiotr Jasiukajtis x0_or_one[2] = -x0; 89825c28e83SPiotr Jasiukajtis x2_or_one[0] = x2; 89925c28e83SPiotr Jasiukajtis x2_or_one[2] = -x2; 90025c28e83SPiotr Jasiukajtis y0_or_zero[0] = y0; 90125c28e83SPiotr Jasiukajtis y0_or_zero[2] = -y0; 90225c28e83SPiotr Jasiukajtis x1 = (x1 - t1) + y1; 90325c28e83SPiotr Jasiukajtis y2_or_zero[0] = y2; 90425c28e83SPiotr Jasiukajtis y2_or_zero[2] = -y2; 90525c28e83SPiotr Jasiukajtis z0 = x0 * x0; 90625c28e83SPiotr Jasiukajtis z1 = x1 * x1; 90725c28e83SPiotr Jasiukajtis z2 = x2 * x2; 90825c28e83SPiotr Jasiukajtis t0 = z0 * (poly3[j0] + z0 * poly4[j0]); 90925c28e83SPiotr Jasiukajtis t1 = z1 * (qq1 + z1 * qq2); 91025c28e83SPiotr Jasiukajtis t2 = z2 * (poly3[j2] + z2 * poly4[j2]); 91125c28e83SPiotr Jasiukajtis t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0)); 91225c28e83SPiotr Jasiukajtis w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 91325c28e83SPiotr Jasiukajtis t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2)); 91425c28e83SPiotr Jasiukajtis j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 91525c28e83SPiotr Jasiukajtis xsb1 = (xsb1 >> 30) & 2; 91625c28e83SPiotr Jasiukajtis n1 ^= (xsb1 & ~(n1 << 1)); 91725c28e83SPiotr Jasiukajtis xsb1 |= 1; 91825c28e83SPiotr Jasiukajtis a1 = __vlibm_TBL_sincos_hi[j1+n1]; 91925c28e83SPiotr Jasiukajtis t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0); 92025c28e83SPiotr Jasiukajtis t1 = (__vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+n1]; 92125c28e83SPiotr Jasiukajtis t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2); 92225c28e83SPiotr Jasiukajtis *py0 = t0; 92325c28e83SPiotr Jasiukajtis *py1 = ( a1 + t1 ); 92425c28e83SPiotr Jasiukajtis *py2 = t2; 92525c28e83SPiotr Jasiukajtis break; 92625c28e83SPiotr Jasiukajtis 92725c28e83SPiotr Jasiukajtis case 6: 92825c28e83SPiotr Jasiukajtis j0 = (xsb0 + 0x4000) & 0xffff8000; 92925c28e83SPiotr Jasiukajtis j1 = n1 & 1; 93025c28e83SPiotr Jasiukajtis j2 = n2 & 1; 93125c28e83SPiotr Jasiukajtis HI(&t0) = j0; 93225c28e83SPiotr Jasiukajtis LO(&t0) = 0; 93325c28e83SPiotr Jasiukajtis x1_or_one[0] = x1; 93425c28e83SPiotr Jasiukajtis x1_or_one[2] = -x1; 93525c28e83SPiotr Jasiukajtis x2_or_one[0] = x2; 93625c28e83SPiotr Jasiukajtis x2_or_one[2] = -x2; 93725c28e83SPiotr Jasiukajtis x0 = (x0 - t0) + y0; 93825c28e83SPiotr Jasiukajtis y1_or_zero[0] = y1; 93925c28e83SPiotr Jasiukajtis y1_or_zero[2] = -y1; 94025c28e83SPiotr Jasiukajtis y2_or_zero[0] = y2; 94125c28e83SPiotr Jasiukajtis y2_or_zero[2] = -y2; 94225c28e83SPiotr Jasiukajtis z0 = x0 * x0; 94325c28e83SPiotr Jasiukajtis z1 = x1 * x1; 94425c28e83SPiotr Jasiukajtis z2 = x2 * x2; 94525c28e83SPiotr Jasiukajtis t0 = z0 * (qq1 + z0 * qq2); 94625c28e83SPiotr Jasiukajtis t1 = z1 * (poly3[j1] + z1 * poly4[j1]); 94725c28e83SPiotr Jasiukajtis t2 = z2 * (poly3[j2] + z2 * poly4[j2]); 94825c28e83SPiotr Jasiukajtis w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 94925c28e83SPiotr Jasiukajtis t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1)); 95025c28e83SPiotr Jasiukajtis t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2)); 95125c28e83SPiotr Jasiukajtis j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 95225c28e83SPiotr Jasiukajtis xsb0 = (xsb0 >> 30) & 2; 95325c28e83SPiotr Jasiukajtis n0 ^= (xsb0 & ~(n0 << 1)); 95425c28e83SPiotr Jasiukajtis xsb0 |= 1; 95525c28e83SPiotr Jasiukajtis a0 = __vlibm_TBL_sincos_hi[j0+n0]; 95625c28e83SPiotr Jasiukajtis t0 = (__vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+n0]; 95725c28e83SPiotr Jasiukajtis t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1); 95825c28e83SPiotr Jasiukajtis t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2); 95925c28e83SPiotr Jasiukajtis *py0 = ( a0 + t0 ); 96025c28e83SPiotr Jasiukajtis *py1 = t1; 96125c28e83SPiotr Jasiukajtis *py2 = t2; 96225c28e83SPiotr Jasiukajtis break; 96325c28e83SPiotr Jasiukajtis 96425c28e83SPiotr Jasiukajtis case 7: 96525c28e83SPiotr Jasiukajtis j0 = n0 & 1; 96625c28e83SPiotr Jasiukajtis j1 = n1 & 1; 96725c28e83SPiotr Jasiukajtis j2 = n2 & 1; 96825c28e83SPiotr Jasiukajtis x0_or_one[0] = x0; 96925c28e83SPiotr Jasiukajtis x0_or_one[2] = -x0; 97025c28e83SPiotr Jasiukajtis x1_or_one[0] = x1; 97125c28e83SPiotr Jasiukajtis x1_or_one[2] = -x1; 97225c28e83SPiotr Jasiukajtis x2_or_one[0] = x2; 97325c28e83SPiotr Jasiukajtis x2_or_one[2] = -x2; 97425c28e83SPiotr Jasiukajtis y0_or_zero[0] = y0; 97525c28e83SPiotr Jasiukajtis y0_or_zero[2] = -y0; 97625c28e83SPiotr Jasiukajtis y1_or_zero[0] = y1; 97725c28e83SPiotr Jasiukajtis y1_or_zero[2] = -y1; 97825c28e83SPiotr Jasiukajtis y2_or_zero[0] = y2; 97925c28e83SPiotr Jasiukajtis y2_or_zero[2] = -y2; 98025c28e83SPiotr Jasiukajtis z0 = x0 * x0; 98125c28e83SPiotr Jasiukajtis z1 = x1 * x1; 98225c28e83SPiotr Jasiukajtis z2 = x2 * x2; 98325c28e83SPiotr Jasiukajtis t0 = z0 * (poly3[j0] + z0 * poly4[j0]); 98425c28e83SPiotr Jasiukajtis t1 = z1 * (poly3[j1] + z1 * poly4[j1]); 98525c28e83SPiotr Jasiukajtis t2 = z2 * (poly3[j2] + z2 * poly4[j2]); 98625c28e83SPiotr Jasiukajtis t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0)); 98725c28e83SPiotr Jasiukajtis t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1)); 98825c28e83SPiotr Jasiukajtis t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2)); 98925c28e83SPiotr Jasiukajtis t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0); 99025c28e83SPiotr Jasiukajtis t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1); 99125c28e83SPiotr Jasiukajtis t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2); 99225c28e83SPiotr Jasiukajtis *py0 = t0; 99325c28e83SPiotr Jasiukajtis *py1 = t1; 99425c28e83SPiotr Jasiukajtis *py2 = t2; 99525c28e83SPiotr Jasiukajtis break; 99625c28e83SPiotr Jasiukajtis } 99725c28e83SPiotr Jasiukajtis 99825c28e83SPiotr Jasiukajtis x += stridex; 99925c28e83SPiotr Jasiukajtis y += stridey; 100025c28e83SPiotr Jasiukajtis i = 0; 100125c28e83SPiotr Jasiukajtis } while (--n > 0); 100225c28e83SPiotr Jasiukajtis 100325c28e83SPiotr Jasiukajtis if (i > 0) 100425c28e83SPiotr Jasiukajtis { 100525c28e83SPiotr Jasiukajtis double fn0, fn1, a0, a1, w0, w1, y0, y1; 100625c28e83SPiotr Jasiukajtis double t0, t1, z0, z1; 100725c28e83SPiotr Jasiukajtis unsigned j0, j1; 100825c28e83SPiotr Jasiukajtis int n0, n1; 100925c28e83SPiotr Jasiukajtis 101025c28e83SPiotr Jasiukajtis if (i > 1) 101125c28e83SPiotr Jasiukajtis { 101225c28e83SPiotr Jasiukajtis n1 = (int) (x1 * invpio2 + half[xsb1]); 101325c28e83SPiotr Jasiukajtis fn1 = (double) n1; 101425c28e83SPiotr Jasiukajtis n1 &= 3; 101525c28e83SPiotr Jasiukajtis a1 = x1 - fn1 * pio2_1; 101625c28e83SPiotr Jasiukajtis w1 = fn1 * pio2_2; 101725c28e83SPiotr Jasiukajtis x1 = a1 - w1; 101825c28e83SPiotr Jasiukajtis y1 = (a1 - x1) - w1; 101925c28e83SPiotr Jasiukajtis a1 = x1; 102025c28e83SPiotr Jasiukajtis w1 = fn1 * pio2_3 - y1; 102125c28e83SPiotr Jasiukajtis x1 = a1 - w1; 102225c28e83SPiotr Jasiukajtis y1 = (a1 - x1) - w1; 102325c28e83SPiotr Jasiukajtis a1 = x1; 102425c28e83SPiotr Jasiukajtis w1 = fn1 * pio2_3t - y1; 102525c28e83SPiotr Jasiukajtis x1 = a1 - w1; 102625c28e83SPiotr Jasiukajtis y1 = (a1 - x1) - w1; 102725c28e83SPiotr Jasiukajtis xsb1 = HI(&x1); 102825c28e83SPiotr Jasiukajtis if ((xsb1 & ~0x80000000) < thresh[n1&1]) 102925c28e83SPiotr Jasiukajtis { 103025c28e83SPiotr Jasiukajtis j1 = n1 & 1; 103125c28e83SPiotr Jasiukajtis x1_or_one[0] = x1; 103225c28e83SPiotr Jasiukajtis x1_or_one[2] = -x1; 103325c28e83SPiotr Jasiukajtis y1_or_zero[0] = y1; 103425c28e83SPiotr Jasiukajtis y1_or_zero[2] = -y1; 103525c28e83SPiotr Jasiukajtis z1 = x1 * x1; 103625c28e83SPiotr Jasiukajtis t1 = z1 * (poly3[j1] + z1 * poly4[j1]); 103725c28e83SPiotr Jasiukajtis t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1)); 103825c28e83SPiotr Jasiukajtis t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1); 103925c28e83SPiotr Jasiukajtis *py1 = t1; 104025c28e83SPiotr Jasiukajtis } 104125c28e83SPiotr Jasiukajtis else 104225c28e83SPiotr Jasiukajtis { 104325c28e83SPiotr Jasiukajtis j1 = (xsb1 + 0x4000) & 0xffff8000; 104425c28e83SPiotr Jasiukajtis HI(&t1) = j1; 104525c28e83SPiotr Jasiukajtis LO(&t1) = 0; 104625c28e83SPiotr Jasiukajtis x1 = (x1 - t1) + y1; 104725c28e83SPiotr Jasiukajtis z1 = x1 * x1; 104825c28e83SPiotr Jasiukajtis t1 = z1 * (qq1 + z1 * qq2); 104925c28e83SPiotr Jasiukajtis w1 = x1 * (one + z1 * (pp1 + z1 * pp2)); 105025c28e83SPiotr Jasiukajtis j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 105125c28e83SPiotr Jasiukajtis xsb1 = (xsb1 >> 30) & 2; 105225c28e83SPiotr Jasiukajtis n1 ^= (xsb1 & ~(n1 << 1)); 105325c28e83SPiotr Jasiukajtis xsb1 |= 1; 105425c28e83SPiotr Jasiukajtis a1 = __vlibm_TBL_sincos_hi[j1+n1]; 105525c28e83SPiotr Jasiukajtis t1 = (__vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)] * w1 + a1 * t1) + __vlibm_TBL_sincos_lo[j1+n1]; 105625c28e83SPiotr Jasiukajtis *py1 = ( a1 + t1 ); 105725c28e83SPiotr Jasiukajtis } 105825c28e83SPiotr Jasiukajtis } 105925c28e83SPiotr Jasiukajtis n0 = (int) (x0 * invpio2 + half[xsb0]); 106025c28e83SPiotr Jasiukajtis fn0 = (double) n0; 106125c28e83SPiotr Jasiukajtis n0 &= 3; 106225c28e83SPiotr Jasiukajtis a0 = x0 - fn0 * pio2_1; 106325c28e83SPiotr Jasiukajtis w0 = fn0 * pio2_2; 106425c28e83SPiotr Jasiukajtis x0 = a0 - w0; 106525c28e83SPiotr Jasiukajtis y0 = (a0 - x0) - w0; 106625c28e83SPiotr Jasiukajtis a0 = x0; 106725c28e83SPiotr Jasiukajtis w0 = fn0 * pio2_3 - y0; 106825c28e83SPiotr Jasiukajtis x0 = a0 - w0; 106925c28e83SPiotr Jasiukajtis y0 = (a0 - x0) - w0; 107025c28e83SPiotr Jasiukajtis a0 = x0; 107125c28e83SPiotr Jasiukajtis w0 = fn0 * pio2_3t - y0; 107225c28e83SPiotr Jasiukajtis x0 = a0 - w0; 107325c28e83SPiotr Jasiukajtis y0 = (a0 - x0) - w0; 107425c28e83SPiotr Jasiukajtis xsb0 = HI(&x0); 107525c28e83SPiotr Jasiukajtis if ((xsb0 & ~0x80000000) < thresh[n0&1]) 107625c28e83SPiotr Jasiukajtis { 107725c28e83SPiotr Jasiukajtis j0 = n0 & 1; 107825c28e83SPiotr Jasiukajtis x0_or_one[0] = x0; 107925c28e83SPiotr Jasiukajtis x0_or_one[2] = -x0; 108025c28e83SPiotr Jasiukajtis y0_or_zero[0] = y0; 108125c28e83SPiotr Jasiukajtis y0_or_zero[2] = -y0; 108225c28e83SPiotr Jasiukajtis z0 = x0 * x0; 108325c28e83SPiotr Jasiukajtis t0 = z0 * (poly3[j0] + z0 * poly4[j0]); 108425c28e83SPiotr Jasiukajtis t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0)); 108525c28e83SPiotr Jasiukajtis t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0); 108625c28e83SPiotr Jasiukajtis *py0 = t0; 108725c28e83SPiotr Jasiukajtis } 108825c28e83SPiotr Jasiukajtis else 108925c28e83SPiotr Jasiukajtis { 109025c28e83SPiotr Jasiukajtis j0 = (xsb0 + 0x4000) & 0xffff8000; 109125c28e83SPiotr Jasiukajtis HI(&t0) = j0; 109225c28e83SPiotr Jasiukajtis LO(&t0) = 0; 109325c28e83SPiotr Jasiukajtis x0 = (x0 - t0) + y0; 109425c28e83SPiotr Jasiukajtis z0 = x0 * x0; 109525c28e83SPiotr Jasiukajtis t0 = z0 * (qq1 + z0 * qq2); 109625c28e83SPiotr Jasiukajtis w0 = x0 * (one + z0 * (pp1 + z0 * pp2)); 109725c28e83SPiotr Jasiukajtis j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3; 109825c28e83SPiotr Jasiukajtis xsb0 = (xsb0 >> 30) & 2; 109925c28e83SPiotr Jasiukajtis n0 ^= (xsb0 & ~(n0 << 1)); 110025c28e83SPiotr Jasiukajtis xsb0 |= 1; 110125c28e83SPiotr Jasiukajtis a0 = __vlibm_TBL_sincos_hi[j0+n0]; 110225c28e83SPiotr Jasiukajtis t0 = (__vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)] * w0 + a0 * t0) + __vlibm_TBL_sincos_lo[j0+n0]; 110325c28e83SPiotr Jasiukajtis *py0 = ( a0 + t0 ); 110425c28e83SPiotr Jasiukajtis } 110525c28e83SPiotr Jasiukajtis } 110625c28e83SPiotr Jasiukajtis 110725c28e83SPiotr Jasiukajtis if (biguns) 110825c28e83SPiotr Jasiukajtis __vlibm_vsin_big(nsave, xsave, sxsave, ysave, sysave, 0x413921fb); 110925c28e83SPiotr Jasiukajtis } 1110