xref: /illumos-gate/usr/src/lib/libmvec/common/__vsin.c (revision 8b8bfb2a)
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
__vsin(int n,double * restrict x,int stridex,double * restrict y,int stridey)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