125c28e8Piotr Jasiukajtis/*
225c28e8Piotr Jasiukajtis * CDDL HEADER START
325c28e8Piotr Jasiukajtis *
425c28e8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e8Piotr Jasiukajtis *
825c28e8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e8Piotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e8Piotr Jasiukajtis * and limitations under the License.
1225c28e8Piotr Jasiukajtis *
1325c28e8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e8Piotr Jasiukajtis *
1925c28e8Piotr Jasiukajtis * CDDL HEADER END
2025c28e8Piotr Jasiukajtis */
2125c28e8Piotr Jasiukajtis
2225c28e8Piotr Jasiukajtis/*
2325c28e8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
2425c28e8Piotr Jasiukajtis */
2525c28e8Piotr Jasiukajtis/*
2625c28e8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
2725c28e8Piotr Jasiukajtis * Use is subject to license terms.
2825c28e8Piotr Jasiukajtis */
2925c28e8Piotr Jasiukajtis
3025c28e8Piotr Jasiukajtis#include <sys/isa_defs.h>
3125c28e8Piotr Jasiukajtis#include <sys/ccompile.h>
3225c28e8Piotr Jasiukajtis
3325c28e8Piotr Jasiukajtis#ifdef _LITTLE_ENDIAN
3425c28e8Piotr Jasiukajtis#define HI(x)	*(1+(int*)x)
3525c28e8Piotr Jasiukajtis#define LO(x)	*(unsigned*)x
3625c28e8Piotr Jasiukajtis#else
3725c28e8Piotr Jasiukajtis#define HI(x)	*(int*)x
3825c28e8Piotr Jasiukajtis#define LO(x)	*(1+(unsigned*)x)
3925c28e8Piotr Jasiukajtis#endif
4025c28e8Piotr Jasiukajtis
4125c28e8Piotr Jasiukajtis#ifdef __RESTRICT
4225c28e8Piotr Jasiukajtis#define restrict _Restrict
4325c28e8Piotr Jasiukajtis#else
4425c28e8Piotr Jasiukajtis#define restrict
4525c28e8Piotr Jasiukajtis#endif
4625c28e8Piotr Jasiukajtis
4725c28e8Piotr Jasiukajtis/*
4825c28e8Piotr Jasiukajtis * vsincos.c
4925c28e8Piotr Jasiukajtis *
5025c28e8Piotr Jasiukajtis * Vector sine and cosine function.  Just slight modifications to vcos.c.
5125c28e8Piotr Jasiukajtis */
5225c28e8Piotr Jasiukajtis
5325c28e8Piotr Jasiukajtisextern const double __vlibm_TBL_sincos_hi[], __vlibm_TBL_sincos_lo[];
5425c28e8Piotr Jasiukajtis
5525c28e8Piotr Jasiukajtisstatic const double
5625c28e8Piotr Jasiukajtis	half[2]	= { 0.5, -0.5 },
5725c28e8Piotr Jasiukajtis	one		= 1.0,
5825c28e8Piotr Jasiukajtis	invpio2 = 0.636619772367581343075535,  /* 53 bits of pi/2 */
5925c28e8Piotr Jasiukajtis	pio2_1	= 1.570796326734125614166,  /* first 33 bits of pi/2 */
6025c28e8Piotr Jasiukajtis	pio2_2	= 6.077100506303965976596e-11, /* second 33 bits of pi/2 */
6125c28e8Piotr Jasiukajtis	pio2_3	= 2.022266248711166455796e-21, /* third 33 bits of pi/2 */
6225c28e8Piotr Jasiukajtis	pio2_3t	= 8.478427660368899643959e-32, /* pi/2 - pio2_3 */
6325c28e8Piotr Jasiukajtis	pp1		= -1.666666666605760465276263943134982554676e-0001,
6425c28e8Piotr Jasiukajtis	pp2		=  8.333261209690963126718376566146180944442e-0003,
6525c28e8Piotr Jasiukajtis	qq1		= -4.999999999977710986407023955908711557870e-0001,
6625c28e8Piotr Jasiukajtis	qq2		=  4.166654863857219350645055881018842089580e-0002,
6725c28e8Piotr Jasiukajtis	poly1[2]= { -1.666666666666629669805215138920301589656e-0001,
6825c28e8Piotr Jasiukajtis				-4.999999999999931701464060878888294524481e-0001 },
6925c28e8Piotr Jasiukajtis	poly2[2]= {  8.333333332390951295683993455280336376663e-0003,
7025c28e8Piotr Jasiukajtis				 4.166666666394861917535640593963708222319e-0002 },
7125c28e8Piotr Jasiukajtis	poly3[2]= { -1.984126237997976692791551778230098403960e-0004,
7225c28e8Piotr Jasiukajtis				-1.388888552656142867832756687736851681462e-0003 },
7325c28e8Piotr Jasiukajtis	poly4[2]= {  2.753403624854277237649987622848330351110e-0006,
7425c28e8Piotr Jasiukajtis				 2.478519423681460796618128289454530524759e-0005 };
7525c28e8Piotr Jasiukajtis
7625c28e8Piotr Jasiukajtis/* Don't __ the following; acomp will handle it */
7725c28e8Piotr Jasiukajtisextern double fabs(double);
7825c28e8Piotr Jasiukajtisextern void __vlibm_vsincos_big(int, double *, int, double *, int, double *, int, int);
7925c28e8Piotr Jasiukajtis
8025c28e8Piotr Jasiukajtis/*
8125c28e8Piotr Jasiukajtis * y[i*stridey] := sin( x[i*stridex] ), for i = 0..n.
8225c28e8Piotr Jasiukajtis * c[i*stridec] := cos( x[i*stridex] ), for i = 0..n.
8325c28e8Piotr Jasiukajtis *
8425c28e8Piotr Jasiukajtis * Calls __vlibm_vsincos_big to handle all elts which have abs >~ 1.647e+06.
8525c28e8Piotr Jasiukajtis * Argument reduction is done here for elts pi/4 < arg < 1.647e+06.
8625c28e8Piotr Jasiukajtis *
8725c28e8Piotr Jasiukajtis * elts < 2^-27 use the approximation 1.0 ~ cos(x).
8825c28e8Piotr Jasiukajtis */
8925c28e8Piotr Jasiukajtisvoid
9025c28e8Piotr Jasiukajtis__vsincos(int n, double * restrict x, int stridex,
9125c28e8Piotr Jasiukajtis				double * restrict y, int stridey,
9225c28e8Piotr Jasiukajtis				double * restrict c, int stridec)
9325c28e8Piotr Jasiukajtis{
9425c28e8Piotr Jasiukajtis	double		x0_or_one[4], x1_or_one[4], x2_or_one[4];
9525c28e8Piotr Jasiukajtis	double		y0_or_zero[4], y1_or_zero[4], y2_or_zero[4];
9625c28e8Piotr Jasiukajtis	double		x0, x1, x2,
9725c28e8Piotr Jasiukajtis			*py0, *py1, *py2,
9825c28e8Piotr Jasiukajtis			*pc0, *pc1, *pc2,
9925c28e8Piotr Jasiukajtis			*xsave, *ysave, *csave;
10025c28e8Piotr Jasiukajtis	unsigned	hx0, hx1, hx2, xsb0, xsb1, xsb2;
10125c28e8Piotr Jasiukajtis	int		i, biguns, nsave, sxsave, sysave, scsave;
102c81a25eToomas Soome	volatile int	v __unused;
10325c28e8Piotr Jasiukajtis	nsave = n;
10425c28e8Piotr Jasiukajtis	xsave = x;
10525c28e8Piotr Jasiukajtis	sxsave = stridex;
10625c28e8Piotr Jasiukajtis	ysave = y;
10725c28e8Piotr Jasiukajtis	sysave = stridey;
10825c28e8Piotr Jasiukajtis	csave = c;
10925c28e8Piotr Jasiukajtis	scsave = stridec;
11025c28e8Piotr Jasiukajtis	biguns = 0;
11125c28e8Piotr Jasiukajtis
11225c28e8Piotr Jasiukajtis	do /* MAIN LOOP */
11325c28e8Piotr Jasiukajtis	{
11425c28e8Piotr Jasiukajtis
11525c28e8Piotr Jasiukajtis		/* Gotos here so _break_ exits MAIN LOOP. */
11625c28e8Piotr JasiukajtisLOOP0:  /* Find first arg in right range. */
11725c28e8Piotr Jasiukajtis		xsb0 = HI(x); /* get most significant word */
11825c28e8Piotr Jasiukajtis		hx0 = xsb0 & ~0x80000000; /* mask off sign bit */
11925c28e8Piotr Jasiukajtis		if (hx0 > 0x3fe921fb) {
12025c28e8Piotr Jasiukajtis			/* Too big: arg reduction needed, so leave for second part */
12125c28e8Piotr Jasiukajtis			biguns = 1;
12225c28e8Piotr Jasiukajtis			x += stridex;
12325c28e8Piotr Jasiukajtis			y += stridey;
12425c28e8Piotr Jasiukajtis			c += stridec;
12525c28e8Piotr Jasiukajtis			i = 0;
12625c28e8Piotr Jasiukajtis			if (--n <= 0)
12725c28e8Piotr Jasiukajtis				break;
12825c28e8Piotr Jasiukajtis			goto LOOP0;
12925c28e8Piotr Jasiukajtis		}
13025c28e8Piotr Jasiukajtis		if (hx0 < 0x3e400000) {
13125c28e8Piotr Jasiukajtis			/* Too small.  cos x ~ 1, sin x ~ x. */
13225c28e8Piotr Jasiukajtis			v = *x;
13325c28e8Piotr Jasiukajtis			*c = 1.0;
13425c28e8Piotr Jasiukajtis			*y = *x;
13525c28e8Piotr Jasiukajtis			x += stridex;
13625c28e8Piotr Jasiukajtis			y += stridey;
13725c28e8Piotr Jasiukajtis			c += stridec;
13825c28e8Piotr Jasiukajtis			i = 0;
13925c28e8Piotr Jasiukajtis			if (--n <= 0)
14025c28e8Piotr Jasiukajtis				break;
14125c28e8Piotr Jasiukajtis			goto LOOP0;
14225c28e8Piotr Jasiukajtis		}
14325c28e8Piotr Jasiukajtis		x0 = *x;
14425c28e8Piotr Jasiukajtis		py0 = y;
14525c28e8Piotr Jasiukajtis		pc0 = c;
14625c28e8Piotr Jasiukajtis		x += stridex;
14725c28e8Piotr Jasiukajtis		y += stridey;
14825c28e8Piotr Jasiukajtis		c += stridec;
14925c28e8Piotr Jasiukajtis		i = 1;
15025c28e8Piotr Jasiukajtis		if (--n <= 0)
15125c28e8Piotr Jasiukajtis			break;
15225c28e8Piotr Jasiukajtis
15325c28e8Piotr JasiukajtisLOOP1: /* Get second arg, same as above. */
15425c28e8Piotr Jasiukajtis		xsb1 = HI(x);
15525c28e8Piotr Jasiukajtis		hx1 = xsb1 & ~0x80000000;
15625c28e8Piotr Jasiukajtis		if (hx1 > 0x3fe921fb)
15725c28e8Piotr Jasiukajtis		{
15825c28e8Piotr Jasiukajtis			biguns = 1;
15925c28e8Piotr Jasiukajtis			x += stridex;
16025c28e8Piotr Jasiukajtis			y += stridey;
16125c28e8Piotr Jasiukajtis			c += stridec;
16225c28e8Piotr Jasiukajtis			i = 1;
16325c28e8Piotr Jasiukajtis			if (--n <= 0)
16425c28e8Piotr Jasiukajtis				break;
16525c28e8Piotr Jasiukajtis			goto LOOP1;
16625c28e8Piotr Jasiukajtis		}
16725c28e8Piotr Jasiukajtis		if (hx1 < 0x3e400000)
16825c28e8Piotr Jasiukajtis		{
16925c28e8Piotr Jasiukajtis			v = *x;
17025c28e8Piotr Jasiukajtis			*c = 1.0;
17125c28e8Piotr Jasiukajtis			*y = *x;
17225c28e8Piotr Jasiukajtis			x += stridex;
17325c28e8Piotr Jasiukajtis			y += stridey;
17425c28e8Piotr Jasiukajtis			c += stridec;
17525c28e8Piotr Jasiukajtis			i = 1;
17625c28e8Piotr Jasiukajtis			if (--n <= 0)
17725c28e8Piotr Jasiukajtis				break;
17825c28e8Piotr Jasiukajtis			goto LOOP1;
17925c28e8Piotr Jasiukajtis		}
18025c28e8Piotr Jasiukajtis		x1 = *x;
18125c28e8Piotr Jasiukajtis		py1 = y;
18225c28e8Piotr Jasiukajtis		pc1 = c;
18325c28e8Piotr Jasiukajtis		x += stridex;
18425c28e8Piotr Jasiukajtis		y += stridey;
18525c28e8Piotr Jasiukajtis		c += stridec;
18625c28e8Piotr Jasiukajtis		i = 2;
18725c28e8Piotr Jasiukajtis		if (--n <= 0)
18825c28e8Piotr Jasiukajtis			break;
18925c28e8Piotr Jasiukajtis
19025c28e8Piotr JasiukajtisLOOP2: /* Get third arg, same as above. */
19125c28e8Piotr Jasiukajtis		xsb2 = HI(x);
19225c28e8Piotr Jasiukajtis		hx2 = xsb2 & ~0x80000000;
19325c28e8Piotr Jasiukajtis		if (hx2 > 0x3fe921fb)
19425c28e8Piotr Jasiukajtis		{
19525c28e8Piotr Jasiukajtis			biguns = 1;
19625c28e8Piotr Jasiukajtis			x += stridex;
19725c28e8Piotr Jasiukajtis			y += stridey;
19825c28e8Piotr Jasiukajtis			c += stridec;
19925c28e8Piotr Jasiukajtis			i = 2;
20025c28e8Piotr Jasiukajtis			if (--n <= 0)
20125c28e8Piotr Jasiukajtis				break;
20225c28e8Piotr Jasiukajtis			goto LOOP2;
20325c28e8Piotr Jasiukajtis		}
20425c28e8Piotr Jasiukajtis		if (hx2 < 0x3e400000)
20525c28e8Piotr Jasiukajtis		{
20625c28e8Piotr Jasiukajtis			v = *x;
20725c28e8Piotr Jasiukajtis			*c = 1.0;
20825c28e8Piotr Jasiukajtis			*y = *x;
20925c28e8Piotr Jasiukajtis			x += stridex;
21025c28e8Piotr Jasiukajtis			y += stridey;
21125c28e8Piotr Jasiukajtis			c += stridec;
21225c28e8Piotr Jasiukajtis			i = 2;
21325c28e8Piotr Jasiukajtis			if (--n <= 0)
21425c28e8Piotr Jasiukajtis				break;
21525c28e8Piotr Jasiukajtis			goto LOOP2;
21625c28e8Piotr Jasiukajtis		}
21725c28e8Piotr Jasiukajtis		x2 = *x;
21825c28e8Piotr Jasiukajtis		py2 = y;
21925c28e8Piotr Jasiukajtis		pc2 = c;
22025c28e8Piotr Jasiukajtis
22125c28e8Piotr Jasiukajtis		/*
22225c28e8Piotr Jasiukajtis		 * 0x3fc40000 = 5/32 ~ 0.15625
22325c28e8Piotr Jasiukajtis		 * Get msb after subtraction.  Will be 1 only if
22425c28e8Piotr Jasiukajtis		 * hx0 - 5/32 is negative.
22525c28e8Piotr Jasiukajtis		 */
22625c28e8Piotr Jasiukajtis		i = (hx2 - 0x3fc40000) >> 31;
22725c28e8Piotr Jasiukajtis		i |= ((hx1 - 0x3fc40000) >> 30) & 2;
22825c28e8Piotr Jasiukajtis		i |= ((hx0 - 0x3fc40000) >> 29) & 4;
22925c28e8Piotr Jasiukajtis		switch (i)
23025c28e8Piotr Jasiukajtis		{
23125c28e8Piotr Jasiukajtis			double		a1_0, a1_1, a1_2, a2_0, a2_1, a2_2;
23225c28e8Piotr Jasiukajtis			double		w0, w1, w2;
23325c28e8Piotr Jasiukajtis			double		t0, t1, t2, t1_0, t1_1, t1_2, t2_0, t2_1, t2_2;
23425c28e8Piotr Jasiukajtis			double		z0, z1, z2;
23525c28e8Piotr Jasiukajtis			unsigned	j0, j1, j2;
23625c28e8Piotr Jasiukajtis
23725c28e8Piotr Jasiukajtis		case 0: /* All are > 5/32 */
23825c28e8Piotr Jasiukajtis			j0 = (xsb0 + 0x4000) & 0xffff8000;
23925c28e8Piotr Jasiukajtis			j1 = (xsb1 + 0x4000) & 0xffff8000;
24025c28e8Piotr Jasiukajtis			j2 = (xsb2 + 0x4000) & 0xffff8000;
24125c28e8Piotr Jasiukajtis
24225c28e8Piotr Jasiukajtis			HI(&t0) = j0;
24325c28e8Piotr Jasiukajtis			HI(&t1) = j1;
24425c28e8Piotr Jasiukajtis			HI(&t2) = j2;
24525c28e8Piotr Jasiukajtis			LO(&t0) = 0;
24625c28e8Piotr Jasiukajtis			LO(&t1) = 0;
24725c28e8Piotr Jasiukajtis			LO(&t2) = 0;
24825c28e8Piotr Jasiukajtis
24925c28e8Piotr Jasiukajtis			x0 -= t0;
25025c28e8Piotr Jasiukajtis			x1 -= t1;
25125c28e8Piotr Jasiukajtis			x2 -= t2;
25225c28e8Piotr Jasiukajtis
25325c28e8Piotr Jasiukajtis			z0 = x0 * x0;
25425c28e8Piotr Jasiukajtis			z1 = x1 * x1;
25525c28e8Piotr Jasiukajtis			z2 = x2 * x2;
25625c28e8Piotr Jasiukajtis
25725c28e8Piotr Jasiukajtis			t0 = z0 * (qq1 + z0 * qq2);
25825c28e8Piotr Jasiukajtis			t1 = z1 * (qq1 + z1 * qq2);
25925c28e8Piotr Jasiukajtis			t2 = z2 * (qq1 + z2 * qq2);
26025c28e8Piotr Jasiukajtis
26125c28e8Piotr Jasiukajtis			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
26225c28e8Piotr Jasiukajtis			w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
26325c28e8Piotr Jasiukajtis			w2 = x2 * (one + z2 * (pp1 + z2 * pp2));
26425c28e8Piotr Jasiukajtis
26525c28e8Piotr Jasiukajtis			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
26625c28e8Piotr Jasiukajtis			j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
26725c28e8Piotr Jasiukajtis			j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
26825c28e8Piotr Jasiukajtis
26925c28e8Piotr Jasiukajtis			xsb0 = (xsb0 >> 30) & 2;
27025c28e8Piotr Jasiukajtis			xsb1 = (xsb1 >> 30) & 2;
27125c28e8Piotr Jasiukajtis			xsb2 = (xsb2 >> 30) & 2;
27225c28e8Piotr Jasiukajtis
27325c28e8Piotr Jasiukajtis			a1_0 = __vlibm_TBL_sincos_hi[j0+xsb0]; /* sin_hi(t) */
27425c28e8Piotr Jasiukajtis			a1_1 = __vlibm_TBL_sincos_hi[j1+xsb1];
27525c28e8Piotr Jasiukajtis			a1_2 = __vlibm_TBL_sincos_hi[j2+xsb2];
27625c28e8Piotr Jasiukajtis
27725c28e8Piotr Jasiukajtis			a2_0 = __vlibm_TBL_sincos_hi[j0+1];	/* cos_hi(t) */
27825c28e8Piotr Jasiukajtis			a2_1 = __vlibm_TBL_sincos_hi[j1+1];
27925c28e8Piotr Jasiukajtis			a2_2 = __vlibm_TBL_sincos_hi[j2+1];
28025c28e8Piotr Jasiukajtis				/* cos_lo(t) */
28125c28e8Piotr Jasiukajtis			t2_0 = __vlibm_TBL_sincos_lo[j0+1] - (a1_0*w0 - a2_0*t0);
28225c28e8Piotr Jasiukajtis			t2_1 = __vlibm_TBL_sincos_lo[j1+1] - (a1_1*w1 - a2_1*t1);
28325c28e8Piotr Jasiukajtis			t2_2 = __vlibm_TBL_sincos_lo[j2+1] - (a1_2*w2 - a2_2*t2);
28425c28e8Piotr Jasiukajtis
28525c28e8Piotr Jasiukajtis			*pc0 = a2_0 + t2_0;
28625c28e8Piotr Jasiukajtis			*pc1 = a2_1 + t2_1;
28725c28e8Piotr Jasiukajtis			*pc2 = a2_2 + t2_2;
28825c28e8Piotr Jasiukajtis
28925c28e8Piotr Jasiukajtis			t1_0 = a2_0*w0 + a1_0*t0;
29025c28e8Piotr Jasiukajtis			t1_1 = a2_1*w1 + a1_1*t1;
29125c28e8Piotr Jasiukajtis			t1_2 = a2_2*w2 + a1_2*t2;
29225c28e8Piotr Jasiukajtis
29325c28e8Piotr Jasiukajtis			t1_0 += __vlibm_TBL_sincos_lo[j0+xsb0]; /* sin_lo(t) */
29425c28e8Piotr Jasiukajtis			t1_1 += __vlibm_TBL_sincos_lo[j1+xsb1];
29525c28e8Piotr Jasiukajtis			t1_2 += __vlibm_TBL_sincos_lo[j2+xsb2];
29625c28e8Piotr Jasiukajtis
29725c28e8Piotr Jasiukajtis			*py0 = a1_0 + t1_0;
29825c28e8Piotr Jasiukajtis			*py1 = a1_1 + t1_1;
29925c28e8Piotr Jasiukajtis			*py2 = a1_2 + t1_2;
30025c28e8Piotr Jasiukajtis
30125c28e8Piotr Jasiukajtis			break;
30225c28e8Piotr Jasiukajtis
30325c28e8Piotr Jasiukajtis		case 1:
30425c28e8Piotr Jasiukajtis			j0 = (xsb0 + 0x4000) & 0xffff8000;
30525c28e8Piotr Jasiukajtis			j1 = (xsb1 + 0x4000) & 0xffff8000;
30625c28e8Piotr Jasiukajtis			HI(&t0) = j0;
30725c28e8Piotr Jasiukajtis			HI(&t1) = j1;
30825c28e8Piotr Jasiukajtis			LO(&t0) = 0;
30925c28e8Piotr Jasiukajtis			LO(&t1) = 0;
31025c28e8Piotr Jasiukajtis			x0 -= t0;
31125c28e8Piotr Jasiukajtis			x1 -= t1;
31225c28e8Piotr Jasiukajtis			z0 = x0 * x0;
31325c28e8Piotr Jasiukajtis			z1 = x1 * x1;
31425c28e8Piotr Jasiukajtis			z2 = x2 * x2;
31525c28e8Piotr Jasiukajtis			t0 = z0 * (qq1 + z0 * qq2);
31625c28e8Piotr Jasiukajtis			t1 = z1 * (qq1 + z1 * qq2);
31725c28e8Piotr Jasiukajtis			t2 = z2 * (poly3[1] + z2 * poly4[1]);
31825c28e8Piotr Jasiukajtis			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
31925c28e8Piotr Jasiukajtis			w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
32025c28e8Piotr Jasiukajtis			t2 = z2 * (poly1[1] + z2 * (poly2[1] + t2));
32125c28e8Piotr Jasiukajtis			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
32225c28e8Piotr Jasiukajtis			j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
32325c28e8Piotr Jasiukajtis			xsb0 = (xsb0 >> 30) & 2;
32425c28e8Piotr Jasiukajtis			xsb1 = (xsb1 >> 30) & 2;
32525c28e8Piotr Jasiukajtis
32625c28e8Piotr Jasiukajtis			a1_0 = __vlibm_TBL_sincos_hi[j0+xsb0]; /* sin_hi(t) */
32725c28e8Piotr Jasiukajtis			a1_1 = __vlibm_TBL_sincos_hi[j1+xsb1];
32825c28e8Piotr Jasiukajtis
32925c28e8Piotr Jasiukajtis			a2_0 = __vlibm_TBL_sincos_hi[j0+1];	/* cos_hi(t) */
33025c28e8Piotr Jasiukajtis			a2_1 = __vlibm_TBL_sincos_hi[j1+1];
33125c28e8Piotr Jasiukajtis				/* cos_lo(t) */
33225c28e8Piotr Jasiukajtis			t2_0 = __vlibm_TBL_sincos_lo[j0+1] - (a1_0*w0 - a2_0*t0);
33325c28e8Piotr Jasiukajtis			t2_1 = __vlibm_TBL_sincos_lo[j1+1] - (a1_1*w1 - a2_1*t1);
33425c28e8Piotr Jasiukajtis
33525c28e8Piotr Jasiukajtis			*pc0 = a2_0 + t2_0;
33625c28e8Piotr Jasiukajtis			*pc1 = a2_1 + t2_1;
33725c28e8Piotr Jasiukajtis			*pc2 = one + t2;
33825c28e8Piotr Jasiukajtis
33925c28e8Piotr Jasiukajtis			t1_0 = a2_0*w0 + a1_0*t0;
34025c28e8Piotr Jasiukajtis			t1_1 = a2_1*w1 + a1_1*t1;
34125c28e8Piotr Jasiukajtis			t2 = z2 * (poly3[0] + z2 * poly4[0]);
34225c28e8Piotr Jasiukajtis
34325c28e8Piotr Jasiukajtis			t1_0 += __vlibm_TBL_sincos_lo[j0+xsb0]; /* sin_lo(t) */
34425c28e8Piotr Jasiukajtis			t1_1 += __vlibm_TBL_sincos_lo[j1+xsb1];
34525c28e8Piotr Jasiukajtis			t2 = z2 * (poly1[0] + z2 * (poly2[0] + t2));
34625c28e8Piotr Jasiukajtis
34725c28e8Piotr Jasiukajtis			*py0 = a1_0 + t1_0;
34825c28e8Piotr Jasiukajtis			*py1 = a1_1 + t1_1;
34925c28e8Piotr Jasiukajtis			t2 = x2 + x2 * t2;
35025c28e8Piotr Jasiukajtis			*py2 = t2;
35125c28e8Piotr Jasiukajtis
35225c28e8Piotr Jasiukajtis			break;
35325c28e8Piotr Jasiukajtis
35425c28e8Piotr Jasiukajtis		case 2:
35525c28e8Piotr Jasiukajtis			j0 = (xsb0 + 0x4000) & 0xffff8000;
35625c28e8Piotr Jasiukajtis			j2 = (xsb2 + 0x4000) & 0xffff8000;
35725c28e8Piotr Jasiukajtis			HI(&t0) = j0;
35825c28e8Piotr Jasiukajtis			HI(&t2) = j2;
35925c28e8Piotr Jasiukajtis			LO(&t0) = 0;
36025c28e8Piotr Jasiukajtis			LO(&t2) = 0;
36125c28e8Piotr Jasiukajtis			x0 -= t0;
36225c28e8Piotr Jasiukajtis			x2 -= t2;
36325c28e8Piotr Jasiukajtis			z0 = x0 * x0;
36425c28e8Piotr Jasiukajtis			z1 = x1 * x1;
36525c28e8Piotr Jasiukajtis			z2 = x2 * x2;
36625c28e8Piotr Jasiukajtis			t0 = z0 * (qq1 + z0 * qq2);
36725c28e8Piotr Jasiukajtis			t1 = z1 * (poly3[1] + z1 * poly4[1]);
36825c28e8Piotr Jasiukajtis			t2 = z2 * (qq1 + z2 * qq2);
36925c28e8Piotr Jasiukajtis			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
37025c28e8Piotr Jasiukajtis			t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1));
37125c28e8Piotr Jasiukajtis			w2 = x2 * (one + z2 * (pp1 + z2 * pp2));
37225c28e8Piotr Jasiukajtis			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
37325c28e8Piotr Jasiukajtis			j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
37425c28e8Piotr Jasiukajtis			xsb0 = (xsb0 >> 30) & 2;
37525c28e8Piotr Jasiukajtis			xsb2 = (xsb2 >> 30) & 2;
37625c28e8Piotr Jasiukajtis
37725c28e8Piotr Jasiukajtis			a1_0 = __vlibm_TBL_sincos_hi[j0+xsb0]; /* sin_hi(t) */
37825c28e8Piotr Jasiukajtis			a1_2 = __vlibm_TBL_sincos_hi[j2+xsb2];
37925c28e8Piotr Jasiukajtis
38025c28e8Piotr Jasiukajtis			a2_0 = __vlibm_TBL_sincos_hi[j0+1];	/* cos_hi(t) */
38125c28e8Piotr Jasiukajtis			a2_2 = __vlibm_TBL_sincos_hi[j2+1];
38225c28e8Piotr Jasiukajtis				/* cos_lo(t) */
38325c28e8Piotr Jasiukajtis			t2_0 = __vlibm_TBL_sincos_lo[j0+1] - (a1_0*w0 - a2_0*t0);
38425c28e8Piotr Jasiukajtis			t2_2 = __vlibm_TBL_sincos_lo[j2+1] - (a1_2*w2 - a2_2*t2);
38525c28e8Piotr Jasiukajtis
38625c28e8Piotr Jasiukajtis			*pc0 = a2_0 + t2_0;
38725c28e8Piotr Jasiukajtis			*pc1 = one + t1;
38825c28e8Piotr Jasiukajtis			*pc2 = a2_2 + t2_2;
38925c28e8Piotr Jasiukajtis
39025c28e8Piotr Jasiukajtis			t1_0 = a2_0*w0 + a1_0*t0;
39125c28e8Piotr Jasiukajtis			t1 = z1 * (poly3[0] + z1 * poly4[0]);
39225c28e8Piotr Jasiukajtis			t1_2 = a2_2*w2 + a1_2*t2;
39325c28e8Piotr Jasiukajtis
39425c28e8Piotr Jasiukajtis			t1_0 += __vlibm_TBL_sincos_lo[j0+xsb0]; /* sin_lo(t) */
39525c28e8Piotr Jasiukajtis			t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1));
39625c28e8Piotr Jasiukajtis			t1_2 += __vlibm_TBL_sincos_lo[j2+xsb2];
39725c28e8Piotr Jasiukajtis
39825c28e8Piotr Jasiukajtis			*py0 = a1_0 + t1_0;
39925c28e8Piotr Jasiukajtis			t1 = x1 + x1 * t1;
40025c28e8Piotr Jasiukajtis			*py1 = t1;
40125c28e8Piotr Jasiukajtis			*py2 = a1_2 + t1_2;
40225c28e8Piotr Jasiukajtis
40325c28e8Piotr Jasiukajtis			break;
40425c28e8Piotr Jasiukajtis
40525c28e8Piotr Jasiukajtis		case 3:
40625c28e8Piotr Jasiukajtis			j0 = (xsb0 + 0x4000) & 0xffff8000;
40725c28e8Piotr Jasiukajtis			HI(&t0) = j0;
40825c28e8Piotr Jasiukajtis			LO(&t0) = 0;
40925c28e8Piotr Jasiukajtis			x0 -= t0;
41025c28e8Piotr Jasiukajtis			z0 = x0 * x0;
41125c28e8Piotr Jasiukajtis			z1 = x1 * x1;
41225c28e8Piotr Jasiukajtis			z2 = x2 * x2;
41325c28e8Piotr Jasiukajtis			t0 = z0 * (qq1 + z0 * qq2);
41425c28e8Piotr Jasiukajtis			t1 = z1 * (poly3[1] + z1 * poly4[1]);
41525c28e8Piotr Jasiukajtis			t2 = z2 * (poly3[1] + z2 * poly4[1]);
41625c28e8Piotr Jasiukajtis			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
41725c28e8Piotr Jasiukajtis			t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1));
41825c28e8Piotr Jasiukajtis			t2 = z2 * (poly1[1] + z2 * (poly2[1] + t2));
41925c28e8Piotr Jasiukajtis			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
42025c28e8Piotr Jasiukajtis			xsb0 = (xsb0 >> 30) & 2;
42125c28e8Piotr Jasiukajtis			a1_0 = __vlibm_TBL_sincos_hi[j0+xsb0]; /* sin_hi(t) */
42225c28e8Piotr Jasiukajtis
42325c28e8Piotr Jasiukajtis			a2_0 = __vlibm_TBL_sincos_hi[j0+1];	/* cos_hi(t) */
42425c28e8Piotr Jasiukajtis
42525c28e8Piotr Jasiukajtis			t2_0 = __vlibm_TBL_sincos_lo[j0+1] - (a1_0*w0 - a2_0*t0);
42625c28e8Piotr Jasiukajtis
42725c28e8Piotr Jasiukajtis			*pc0 = a2_0 + t2_0;
42825c28e8Piotr Jasiukajtis			*pc1 = one + t1;
42925c28e8Piotr Jasiukajtis			*pc2 = one + t2;
43025c28e8Piotr Jasiukajtis
43125c28e8Piotr Jasiukajtis			t1_0 = a2_0*w0 + a1_0*t0;
43225c28e8Piotr Jasiukajtis			t1 = z1 * (poly3[0] + z1 * poly4[0]);
43325c28e8Piotr Jasiukajtis			t2 = z2 * (poly3[0] + z2 * poly4[0]);
43425c28e8Piotr Jasiukajtis
43525c28e8Piotr Jasiukajtis			t1_0 += __vlibm_TBL_sincos_lo[j0+xsb0]; /* sin_lo(t) */
43625c28e8Piotr Jasiukajtis			t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1));
43725c28e8Piotr Jasiukajtis			t2 = z2 * (poly1[0] + z2 * (poly2[0] + t2));
43825c28e8Piotr Jasiukajtis
43925c28e8Piotr Jasiukajtis			*py0 = a1_0 + t1_0;
44025c28e8Piotr Jasiukajtis			t1 = x1 + x1 * t1;
44125c28e8Piotr Jasiukajtis			*py1 = t1;
44225c28e8Piotr Jasiukajtis			t2 = x2 + x2 * t2;
44325c28e8Piotr Jasiukajtis			*py2 = t2;
44425c28e8Piotr Jasiukajtis
44525c28e8Piotr Jasiukajtis			break;
44625c28e8Piotr Jasiukajtis
44725c28e8Piotr Jasiukajtis		case 4:
44825c28e8Piotr Jasiukajtis			j1 = (xsb1 + 0x4000) & 0xffff8000;
44925c28e8Piotr Jasiukajtis			j2 = (xsb2 + 0x4000) & 0xffff8000;
45025c28e8Piotr Jasiukajtis			HI(&t1) = j1;
45125c28e8Piotr Jasiukajtis			HI(&t2) = j2;
45225c28e8Piotr Jasiukajtis			LO(&t1) = 0;
45325c28e8Piotr Jasiukajtis			LO(&t2) = 0;
45425c28e8Piotr Jasiukajtis			x1 -= t1;
45525c28e8Piotr Jasiukajtis			x2 -= t2;
45625c28e8Piotr Jasiukajtis			z0 = x0 * x0;
45725c28e8Piotr Jasiukajtis			z1 = x1 * x1;
45825c28e8Piotr Jasiukajtis			z2 = x2 * x2;
45925c28e8Piotr Jasiukajtis			t0 = z0 * (poly3[1] + z0 * poly4[1]);
46025c28e8Piotr Jasiukajtis			t1 = z1 * (qq1 + z1 * qq2);
46125c28e8Piotr Jasiukajtis			t2 = z2 * (qq1 + z2 * qq2);
46225c28e8Piotr Jasiukajtis			t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0));
46325c28e8Piotr Jasiukajtis			w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
46425c28e8Piotr Jasiukajtis			w2 = x2 * (one + z2 * (pp1 + z2 * pp2));
46525c28e8Piotr Jasiukajtis			j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
46625c28e8Piotr Jasiukajtis			j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
46725c28e8Piotr Jasiukajtis			xsb1 = (xsb1 >> 30) & 2;
46825c28e8Piotr Jasiukajtis			xsb2 = (xsb2 >> 30) & 2;
46925c28e8Piotr Jasiukajtis
47025c28e8Piotr Jasiukajtis			a1_1 = __vlibm_TBL_sincos_hi[j1+xsb1];
47125c28e8Piotr Jasiukajtis			a1_2 = __vlibm_TBL_sincos_hi[j2+xsb2];
47225c28e8Piotr Jasiukajtis
47325c28e8Piotr Jasiukajtis			a2_1 = __vlibm_TBL_sincos_hi[j1+1];
47425c28e8Piotr Jasiukajtis			a2_2 = __vlibm_TBL_sincos_hi[j2+1];
47525c28e8Piotr Jasiukajtis				/* cos_lo(t) */
47625c28e8Piotr Jasiukajtis			t2_1 = __vlibm_TBL_sincos_lo[j1+1] - (a1_1*w1 - a2_1*t1);
47725c28e8Piotr Jasiukajtis			t2_2 = __vlibm_TBL_sincos_lo[j2+1] - (a1_2*w2 - a2_2*t2);
47825c28e8Piotr Jasiukajtis
47925c28e8Piotr Jasiukajtis			*pc0 = one + t0;
48025c28e8Piotr Jasiukajtis			*pc1 = a2_1 + t2_1;
48125c28e8Piotr Jasiukajtis			*pc2 = a2_2 + t2_2;
48225c28e8Piotr Jasiukajtis
48325c28e8Piotr Jasiukajtis			t0 = z0 * (poly3[0] + z0 * poly4[0]);
48425c28e8Piotr Jasiukajtis			t1_1 = a2_1*w1 + a1_1*t1;
48525c28e8Piotr Jasiukajtis			t1_2 = a2_2*w2 + a1_2*t2;
48625c28e8Piotr Jasiukajtis
48725c28e8Piotr Jasiukajtis			t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0));
48825c28e8Piotr Jasiukajtis			t1_1 += __vlibm_TBL_sincos_lo[j1+xsb1];
48925c28e8Piotr Jasiukajtis			t1_2 += __vlibm_TBL_sincos_lo[j2+xsb2];
49025c28e8Piotr Jasiukajtis
49125c28e8Piotr Jasiukajtis			t0 = x0 + x0 * t0;
49225c28e8Piotr Jasiukajtis			*py0 = t0;
49325c28e8Piotr Jasiukajtis			*py1 = a1_1 + t1_1;
49425c28e8Piotr Jasiukajtis			*py2 = a1_2 + t1_2;
49525c28e8Piotr Jasiukajtis
49625c28e8Piotr Jasiukajtis			break;
49725c28e8Piotr Jasiukajtis
49825c28e8Piotr Jasiukajtis		case 5:
49925c28e8Piotr Jasiukajtis			j1 = (xsb1 + 0x4000) & 0xffff8000;
50025c28e8Piotr Jasiukajtis			HI(&t1) = j1;
50125c28e8Piotr Jasiukajtis			LO(&t1) = 0;
50225c28e8Piotr Jasiukajtis			x1 -= t1;
50325c28e8Piotr Jasiukajtis			z0 = x0 * x0;
50425c28e8Piotr Jasiukajtis			z1 = x1 * x1;
50525c28e8Piotr Jasiukajtis			z2 = x2 * x2;
50625c28e8Piotr Jasiukajtis			t0 = z0 * (poly3[1] + z0 * poly4[1]);
50725c28e8Piotr Jasiukajtis			t1 = z1 * (qq1 + z1 * qq2);
50825c28e8Piotr Jasiukajtis			t2 = z2 * (poly3[1] + z2 * poly4[1]);
50925c28e8Piotr Jasiukajtis			t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0));
51025c28e8Piotr Jasiukajtis			w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
51125c28e8Piotr Jasiukajtis			t2 = z2 * (poly1[1] + z2 * (poly2[1] + t2));
51225c28e8Piotr Jasiukajtis			j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
51325c28e8Piotr Jasiukajtis			xsb1 = (xsb1 >> 30) & 2;
51425c28e8Piotr Jasiukajtis
51525c28e8Piotr Jasiukajtis			a1_1 = __vlibm_TBL_sincos_hi[j1+xsb1];
51625c28e8Piotr Jasiukajtis
51725c28e8Piotr Jasiukajtis			a2_1 = __vlibm_TBL_sincos_hi[j1+1];
51825c28e8Piotr Jasiukajtis
51925c28e8Piotr Jasiukajtis			t2_1 = __vlibm_TBL_sincos_lo[j1+1] - (a1_1*w1 - a2_1*t1);
52025c28e8Piotr Jasiukajtis
52125c28e8Piotr Jasiukajtis			*pc0 = one + t0;
52225c28e8Piotr Jasiukajtis			*pc1 = a2_1 + t2_1;
52325c28e8Piotr Jasiukajtis			*pc2 = one + t2;
52425c28e8Piotr Jasiukajtis
52525c28e8Piotr Jasiukajtis			t0 = z0 * (poly3[0] + z0 * poly4[0]);
52625c28e8Piotr Jasiukajtis			t1_1 = a2_1*w1 + a1_1*t1;
52725c28e8Piotr Jasiukajtis			t2 = z2 * (poly3[0] + z2 * poly4[0]);
52825c28e8Piotr Jasiukajtis
52925c28e8Piotr Jasiukajtis			t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0));
53025c28e8Piotr Jasiukajtis			t1_1 += __vlibm_TBL_sincos_lo[j1+xsb1];
53125c28e8Piotr Jasiukajtis			t2 = z2 * (poly1[0] + z2 * (poly2[0] + t2));
53225c28e8Piotr Jasiukajtis
53325c28e8Piotr Jasiukajtis			t0 = x0 + x0 * t0;
53425c28e8Piotr Jasiukajtis			*py0 = t0;
53525c28e8Piotr Jasiukajtis			*py1 = a1_1 + t1_1;
53625c28e8Piotr Jasiukajtis			t2 = x2 + x2 * t2;
53725c28e8Piotr Jasiukajtis			*py2 = t2;
53825c28e8Piotr Jasiukajtis
53925c28e8Piotr Jasiukajtis			break;
54025c28e8Piotr Jasiukajtis
54125c28e8Piotr Jasiukajtis		case 6:
54225c28e8Piotr Jasiukajtis			j2 = (xsb2 + 0x4000) & 0xffff8000;
54325c28e8Piotr Jasiukajtis			HI(&t2) = j2;
54425c28e8Piotr Jasiukajtis			LO(&t2) = 0;
54525c28e8Piotr Jasiukajtis			x2 -= t2;
54625c28e8Piotr Jasiukajtis			z0 = x0 * x0;
54725c28e8Piotr Jasiukajtis			z1 = x1 * x1;
54825c28e8Piotr Jasiukajtis			z2 = x2 * x2;
54925c28e8Piotr Jasiukajtis			t0 = z0 * (poly3[1] + z0 * poly4[1]);
55025c28e8Piotr Jasiukajtis			t1 = z1 * (poly3[1] + z1 * poly4[1]);
55125c28e8Piotr Jasiukajtis			t2 = z2 * (qq1 + z2 * qq2);
55225c28e8Piotr Jasiukajtis			t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0));
55325c28e8Piotr Jasiukajtis			t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1));
55425c28e8Piotr Jasiukajtis			w2 = x2 * (one + z2 * (pp1 + z2 * pp2));
55525c28e8Piotr Jasiukajtis			j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
55625c28e8Piotr Jasiukajtis			xsb2 = (xsb2 >> 30) & 2;
55725c28e8Piotr Jasiukajtis			a1_2 = __vlibm_TBL_sincos_hi[j2+xsb2];
55825c28e8Piotr Jasiukajtis
55925c28e8Piotr Jasiukajtis			a2_2 = __vlibm_TBL_sincos_hi[j2+1];
56025c28e8Piotr Jasiukajtis
56125c28e8Piotr Jasiukajtis			t2_2 = __vlibm_TBL_sincos_lo[j2+1] - (a1_2*w2 - a2_2*t2);
56225c28e8Piotr Jasiukajtis
56325c28e8Piotr Jasiukajtis			*pc0 = one + t0;
56425c28e8Piotr Jasiukajtis			*pc1 = one + t1;
56525c28e8Piotr Jasiukajtis			*pc2 = a2_2 + t2_2;
56625c28e8Piotr Jasiukajtis
56725c28e8Piotr Jasiukajtis			t0 = z0 * (poly3[0] + z0 * poly4[0]);
56825c28e8Piotr Jasiukajtis			t1 = z1 * (poly3[0] + z1 * poly4[0]);
56925c28e8Piotr Jasiukajtis			t1_2 = a2_2*w2 + a1_2*t2;
57025c28e8Piotr Jasiukajtis
57125c28e8Piotr Jasiukajtis			t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0));
57225c28e8Piotr Jasiukajtis			t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1));
57325c28e8Piotr Jasiukajtis			t1_2 += __vlibm_TBL_sincos_lo[j2+xsb2];
57425c28e8Piotr Jasiukajtis
57525c28e8Piotr Jasiukajtis			t0 = x0 + x0 * t0;
57625c28e8Piotr Jasiukajtis			*py0 = t0;
57725c28e8Piotr Jasiukajtis			t1 = x1 + x1 * t1;
57825c28e8Piotr Jasiukajtis			*py1 = t1;
57925c28e8Piotr Jasiukajtis			*py2 = a1_2 + t1_2;
58025c28e8Piotr Jasiukajtis
58125c28e8Piotr Jasiukajtis			break;
58225c28e8Piotr Jasiukajtis
58325c28e8Piotr Jasiukajtis		case 7: /* All are < 5/32 */
58425c28e8Piotr Jasiukajtis			z0 = x0 * x0;
58525c28e8Piotr Jasiukajtis			z1 = x1 * x1;
58625c28e8Piotr Jasiukajtis			z2 = x2 * x2;
58725c28e8Piotr Jasiukajtis			t0 = z0 * (poly3[1] + z0 * poly4[1]);
58825c28e8Piotr Jasiukajtis			t1 = z1 * (poly3[1] + z1 * poly4[1]);
58925c28e8Piotr Jasiukajtis			t2 = z2 * (poly3[1] + z2 * poly4[1]);
59025c28e8Piotr Jasiukajtis			t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0));
59125c28e8Piotr Jasiukajtis			t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1));
59225c28e8Piotr Jasiukajtis			t2 = z2 * (poly1[1] + z2 * (poly2[1] + t2));
59325c28e8Piotr Jasiukajtis			*pc0 = one + t0;
59425c28e8Piotr Jasiukajtis			*pc1 = one + t1;
59525c28e8Piotr Jasiukajtis			*pc2 = one + t2;
59625c28e8Piotr Jasiukajtis			t0 = z0 * (poly3[0] + z0 * poly4[0]);
59725c28e8Piotr Jasiukajtis			t1 = z1 * (poly3[0] + z1 * poly4[0]);
59825c28e8Piotr Jasiukajtis			t2 = z2 * (poly3[0] + z2 * poly4[0]);
59925c28e8Piotr Jasiukajtis			t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0));
60025c28e8Piotr Jasiukajtis			t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1));
60125c28e8Piotr Jasiukajtis			t2 = z2 * (poly1[0] + z2 * (poly2[0] + t2));
60225c28e8Piotr Jasiukajtis			t0 = x0 + x0 * t0;
60325c28e8Piotr Jasiukajtis			t1 = x1 + x1 * t1;
60425c28e8Piotr Jasiukajtis			t2 = x2 + x2 * t2;
60525c28e8Piotr Jasiukajtis			*py0 = t0;
60625c28e8Piotr Jasiukajtis			*py1 = t1;
60725c28e8Piotr Jasiukajtis			*py2 = t2;
60825c28e8Piotr Jasiukajtis			break;
60925c28e8Piotr Jasiukajtis		}
61025c28e8Piotr Jasiukajtis
61125c28e8Piotr Jasiukajtis		x += stridex;
61225c28e8Piotr Jasiukajtis		y += stridey;
61325c28e8Piotr Jasiukajtis		c += stridec;
61425c28e8Piotr Jasiukajtis		i = 0;
61525c28e8Piotr Jasiukajtis	} while (--n > 0); /* END MAIN LOOP */
61625c28e8Piotr Jasiukajtis
61725c28e8Piotr Jasiukajtis	/*
61825c28e8Piotr Jasiukajtis	 * CLEAN UP last 0, 1, or 2 elts.
61925c28e8Piotr Jasiukajtis	 */
62025c28e8Piotr Jasiukajtis	if (i > 0) /* Clean up elts at tail.  i < 3. */
62125c28e8Piotr Jasiukajtis	{
62225c28e8Piotr Jasiukajtis		double		a1_0, a1_1, a2_0, a2_1;
62325c28e8Piotr Jasiukajtis		double		w0, w1;
62425c28e8Piotr Jasiukajtis		double		t0, t1, t1_0, t1_1, t2_0, t2_1;
62525c28e8Piotr Jasiukajtis		double		z0, z1;
62625c28e8Piotr Jasiukajtis		unsigned	j0, j1;
62725c28e8Piotr Jasiukajtis
62825c28e8Piotr Jasiukajtis		if (i > 1)
62925c28e8Piotr Jasiukajtis		{
63025c28e8Piotr Jasiukajtis			if (hx1 < 0x3fc40000)
63125c28e8Piotr Jasiukajtis			{
63225c28e8Piotr Jasiukajtis				z1 = x1 * x1;
63325c28e8Piotr Jasiukajtis				t1 = z1 * (poly3[1] + z1 * poly4[1]);
63425c28e8Piotr Jasiukajtis				t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1));
63525c28e8Piotr Jasiukajtis				t1 = one + t1;
63625c28e8Piotr Jasiukajtis				*pc1 = t1;
63725c28e8Piotr Jasiukajtis				t1 = z1 * (poly3[0] + z1 * poly4[0]);
63825c28e8Piotr Jasiukajtis				t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1));
63925c28e8Piotr Jasiukajtis				t1 = x1 + x1 * t1;
64025c28e8Piotr Jasiukajtis				*py1 = t1;
64125c28e8Piotr Jasiukajtis			}
64225c28e8Piotr Jasiukajtis			else
64325c28e8Piotr Jasiukajtis			{
64425c28e8Piotr Jasiukajtis				j1 = (xsb1 + 0x4000) & 0xffff8000;
64525c28e8Piotr Jasiukajtis				HI(&t1) = j1;
64625c28e8Piotr Jasiukajtis				LO(&t1) = 0;
64725c28e8Piotr Jasiukajtis				x1 -= t1;
64825c28e8Piotr Jasiukajtis				z1 = x1 * x1;
64925c28e8Piotr Jasiukajtis				t1 = z1 * (qq1 + z1 * qq2);
65025c28e8Piotr Jasiukajtis				w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
65125c28e8Piotr Jasiukajtis				j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
65225c28e8Piotr Jasiukajtis				xsb1 = (xsb1 >> 30) & 2;
65325c28e8Piotr Jasiukajtis				a1_1 = __vlibm_TBL_sincos_hi[j1+xsb1];
65425c28e8Piotr Jasiukajtis				a2_1 = __vlibm_TBL_sincos_hi[j1+1];
65525c28e8Piotr Jasiukajtis				t2_1 = __vlibm_TBL_sincos_lo[j1+1] - (a1_1*w1 - a2_1*t1);
65625c28e8Piotr Jasiukajtis				*pc1 = a2_1 + t2_1;
65725c28e8Piotr Jasiukajtis				t1_1 = a2_1*w1 + a1_1*t1;
65825c28e8Piotr Jasiukajtis				t1_1 += __vlibm_TBL_sincos_lo[j1+xsb1];
65925c28e8Piotr Jasiukajtis				*py1 = a1_1 + t1_1;
66025c28e8Piotr Jasiukajtis			}
66125c28e8Piotr Jasiukajtis		}
66225c28e8Piotr Jasiukajtis		if (hx0 < 0x3fc40000)
66325c28e8Piotr Jasiukajtis		{
66425c28e8Piotr Jasiukajtis			z0 = x0 * x0;
66525c28e8Piotr Jasiukajtis			t0 = z0 * (poly3[1] + z0 * poly4[1]);
66625c28e8Piotr Jasiukajtis			t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0));
66725c28e8Piotr Jasiukajtis			t0 = one + t0;
66825c28e8Piotr Jasiukajtis			*pc0 = t0;
66925c28e8Piotr Jasiukajtis			t0 = z0 * (poly3[0] + z0 * poly4[0]);
67025c28e8Piotr Jasiukajtis			t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0));
67125c28e8Piotr Jasiukajtis			t0 = x0 + x0 * t0;
67225c28e8Piotr Jasiukajtis			*py0 = t0;
67325c28e8Piotr Jasiukajtis		}
67425c28e8Piotr Jasiukajtis		else
67525c28e8Piotr Jasiukajtis		{
67625c28e8Piotr Jasiukajtis			j0 = (xsb0 + 0x4000) & 0xffff8000;
67725c28e8Piotr Jasiukajtis			HI(&t0) = j0;
67825c28e8Piotr Jasiukajtis			LO(&t0) = 0;
67925c28e8Piotr Jasiukajtis			x0 -= t0;
68025c28e8Piotr Jasiukajtis			z0 = x0 * x0;
68125c28e8Piotr Jasiukajtis			t0 = z0 * (qq1 + z0 * qq2);
68225c28e8Piotr Jasiukajtis			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
68325c28e8Piotr Jasiukajtis			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
68425c28e8Piotr Jasiukajtis			xsb0 = (xsb0 >> 30) & 2;
68525c28e8Piotr Jasiukajtis			a1_0 = __vlibm_TBL_sincos_hi[j0+xsb0]; /* sin_hi(t) */
68625c28e8Piotr Jasiukajtis			a2_0 = __vlibm_TBL_sincos_hi[j0+1];	/* cos_hi(t) */
68725c28e8Piotr Jasiukajtis			t2_0 = __vlibm_TBL_sincos_lo[j0+1] - (a1_0*w0 - a2_0*t0);
68825c28e8Piotr Jasiukajtis			*pc0 = a2_0 + t2_0;
68925c28e8Piotr Jasiukajtis			t1_0 = a2_0*w0 + a1_0*t0;
69025c28e8Piotr Jasiukajtis			t1_0 += __vlibm_TBL_sincos_lo[j0+xsb0]; /* sin_lo(t) */
69125c28e8Piotr Jasiukajtis			*py0 = a1_0 + t1_0;
69225c28e8Piotr Jasiukajtis		}
69325c28e8Piotr Jasiukajtis	} /* END CLEAN UP */
69425c28e8Piotr Jasiukajtis
69525c28e8Piotr Jasiukajtis	if (!biguns)
69625c28e8Piotr Jasiukajtis		return;
69725c28e8Piotr Jasiukajtis
69825c28e8Piotr Jasiukajtis	/*
69925c28e8Piotr Jasiukajtis	 * Take care of BIGUNS.
70025c28e8Piotr Jasiukajtis	 */
70125c28e8Piotr Jasiukajtis	n = nsave;
70225c28e8Piotr Jasiukajtis	x = xsave;
70325c28e8Piotr Jasiukajtis	stridex = sxsave;
70425c28e8Piotr Jasiukajtis	y = ysave;
70525c28e8Piotr Jasiukajtis	stridey = sysave;
70625c28e8Piotr Jasiukajtis	c = csave;
70725c28e8Piotr Jasiukajtis	stridec = scsave;
70825c28e8Piotr Jasiukajtis	biguns = 0;
70925c28e8Piotr Jasiukajtis
71025c28e8Piotr Jasiukajtis	x0_or_one[1] = 1.0;
71125c28e8Piotr Jasiukajtis	x1_or_one[1] = 1.0;
71225c28e8Piotr Jasiukajtis	x2_or_one[1] = 1.0;
71325c28e8Piotr Jasiukajtis	x0_or_one[3] = -1.0;
71425c28e8Piotr Jasiukajtis	x1_or_one[3] = -1.0;
71525c28e8Piotr Jasiukajtis	x2_or_one[3] = -1.0;
71625c28e8Piotr Jasiukajtis	y0_or_zero[1] = 0.0;
71725c28e8Piotr Jasiukajtis	y1_or_zero[1] = 0.0;
71825c28e8Piotr Jasiukajtis	y2_or_zero[1] = 0.0;
71925c28e8Piotr Jasiukajtis	y0_or_zero[3] = 0.0;
72025c28e8Piotr Jasiukajtis	y1_or_zero[3] = 0.0;
72125c28e8Piotr Jasiukajtis	y2_or_zero[3] = 0.0;
72225c28e8Piotr Jasiukajtis
72325c28e8Piotr Jasiukajtis	do
72425c28e8Piotr Jasiukajtis	{
72525c28e8Piotr Jasiukajtis		double		fn0, fn1, fn2, a0, a1, a2, w0, w1, w2, y0, y1, y2;
72625c28e8Piotr Jasiukajtis		unsigned	hx;
72725c28e8Piotr Jasiukajtis		int			n0, n1, n2;
72825c28e8Piotr Jasiukajtis
72925c28e8Piotr Jasiukajtis		/*
73025c28e8Piotr Jasiukajtis		 * Find 3 more to work on: Not already done, not too big.
73125c28e8Piotr Jasiukajtis		 */
73225c28e8Piotr Jasiukajtisloop0:
73325c28e8Piotr Jasiukajtis		hx = HI(x);
73425c28e8Piotr Jasiukajtis		xsb0 = hx >> 31;
73525c28e8Piotr Jasiukajtis		hx &= ~0x80000000;
73625c28e8Piotr Jasiukajtis		if (hx <= 0x3fe921fb) /* Done above. */
73725c28e8Piotr Jasiukajtis		{
73825c28e8Piotr Jasiukajtis			x += stridex;
73925c28e8Piotr Jasiukajtis			y += stridey;
74025c28e8Piotr Jasiukajtis			c += stridec;
74125c28e8Piotr Jasiukajtis			i = 0;
74225c28e8Piotr Jasiukajtis			if (--n <= 0)
74325c28e8Piotr Jasiukajtis				break;
74425c28e8Piotr Jasiukajtis			goto loop0;
74525c28e8Piotr Jasiukajtis		}
74625c28e8Piotr Jasiukajtis		if (hx > 0x413921fb) /* (1.6471e+06) Too big: leave it. */
74725c28e8Piotr Jasiukajtis		{
74825c28e8Piotr Jasiukajtis			if (hx >= 0x7ff00000) /* Inf or NaN */
74925c28e8Piotr Jasiukajtis			{
75025c28e8Piotr Jasiukajtis				x0 = *x;
75125c28e8Piotr Jasiukajtis				*y = x0 - x0;
75225c28e8Piotr Jasiukajtis				*c = x0 - x0;
75325c28e8Piotr Jasiukajtis			}
75425c28e8Piotr Jasiukajtis			else {
75525c28e8Piotr Jasiukajtis				biguns = 1;
75625c28e8Piotr Jasiukajtis			}
75725c28e8Piotr Jasiukajtis			x += stridex;
75825c28e8Piotr Jasiukajtis			y += stridey;
75925c28e8Piotr Jasiukajtis			c += stridec;
76025c28e8Piotr Jasiukajtis			i = 0;
76125c28e8Piotr Jasiukajtis			if (--n <= 0)
76225c28e8Piotr Jasiukajtis				break;
76325c28e8Piotr Jasiukajtis			goto loop0;
76425c28e8Piotr Jasiukajtis		}
76525c28e8Piotr Jasiukajtis		x0 = *x;
76625c28e8Piotr Jasiukajtis		py0 = y;
76725c28e8Piotr Jasiukajtis		pc0 = c;
76825c28e8Piotr Jasiukajtis		x += stridex;
76925c28e8Piotr Jasiukajtis		y += stridey;
77025c28e8Piotr Jasiukajtis		c += stridec;
77125c28e8Piotr Jasiukajtis		i = 1;
77225c28e8Piotr Jasiukajtis		if (--n <= 0)
77325c28e8Piotr Jasiukajtis			break;
77425c28e8Piotr Jasiukajtis
77525c28e8Piotr Jasiukajtisloop1:
77625c28e8Piotr Jasiukajtis		hx = HI(x);
77725c28e8Piotr Jasiukajtis		xsb1 = hx >> 31;
77825c28e8Piotr Jasiukajtis		hx &= ~0x80000000;
77925c28e8Piotr Jasiukajtis		if (hx <= 0x3fe921fb)
78025c28e8Piotr Jasiukajtis		{
78125c28e8Piotr Jasiukajtis			x += stridex;
78225c28e8Piotr Jasiukajtis			y += stridey;
78325c28e8Piotr Jasiukajtis			c += stridec;
78425c28e8Piotr Jasiukajtis			i = 1;
78525c28e8Piotr Jasiukajtis			if (--n <= 0)
78625c28e8Piotr Jasiukajtis				break;
78725c28e8Piotr Jasiukajtis			goto loop1;
78825c28e8Piotr Jasiukajtis		}
78925c28e8Piotr Jasiukajtis		if (hx > 0x413921fb)
79025c28e8Piotr Jasiukajtis		{
79125c28e8Piotr Jasiukajtis			if (hx >= 0x7ff00000)
79225c28e8Piotr Jasiukajtis			{
79325c28e8Piotr Jasiukajtis				x1 = *x;
79425c28e8Piotr Jasiukajtis				*y = x1 - x1;
79525c28e8Piotr Jasiukajtis				*c = x1 - x1;
79625c28e8Piotr Jasiukajtis			}
79725c28e8Piotr Jasiukajtis			else {
79825c28e8Piotr Jasiukajtis				biguns = 1;
79925c28e8Piotr Jasiukajtis			}
80025c28e8Piotr Jasiukajtis			x += stridex;
80125c28e8Piotr Jasiukajtis			y += stridey;
80225c28e8Piotr Jasiukajtis			c += stridec;
80325c28e8Piotr Jasiukajtis			i = 1;
80425c28e8Piotr Jasiukajtis			if (--n <= 0)
80525c28e8Piotr Jasiukajtis				break;
80625c28e8Piotr Jasiukajtis			goto loop1;
80725c28e8Piotr Jasiukajtis		}
80825c28e8Piotr Jasiukajtis		x1 = *x;
80925c28e8Piotr Jasiukajtis		py1 = y;
81025c28e8Piotr Jasiukajtis		pc1 = c;
81125c28e8Piotr Jasiukajtis		x += stridex;
81225c28e8Piotr Jasiukajtis		y += stridey;
81325c28e8Piotr Jasiukajtis		c += stridec;
81425c28e8Piotr Jasiukajtis		i = 2;
81525c28e8Piotr Jasiukajtis		if (--n <= 0)
81625c28e8Piotr Jasiukajtis			break;
81725c28e8Piotr Jasiukajtis
81825c28e8Piotr Jasiukajtisloop2:
81925c28e8Piotr Jasiukajtis		hx = HI(x);
82025c28e8Piotr Jasiukajtis		xsb2 = hx >> 31;
82125c28e8Piotr Jasiukajtis		hx &= ~0x80000000;
82225c28e8Piotr Jasiukajtis		if (hx <= 0x3fe921fb)
82325c28e8Piotr Jasiukajtis		{
82425c28e8Piotr Jasiukajtis			x += stridex;
82525c28e8Piotr Jasiukajtis			y += stridey;
82625c28e8Piotr Jasiukajtis			c += stridec;
82725c28e8Piotr Jasiukajtis			i = 2;
82825c28e8Piotr Jasiukajtis			if (--n <= 0)
82925c28e8Piotr Jasiukajtis				break;
83025c28e8Piotr Jasiukajtis			goto loop2;
83125c28e8Piotr Jasiukajtis		}
83225c28e8Piotr Jasiukajtis		if (hx > 0x413921fb)
83325c28e8Piotr Jasiukajtis		{
83425c28e8Piotr Jasiukajtis			if (hx >= 0x7ff00000)
83525c28e8Piotr Jasiukajtis			{
83625c28e8Piotr Jasiukajtis				x2 = *x;
83725c28e8Piotr Jasiukajtis				*y = x2 - x2;
83825c28e8Piotr Jasiukajtis				*c = x2 - x2;
83925c28e8Piotr Jasiukajtis			}
84025c28e8Piotr Jasiukajtis			else {
84125c28e8Piotr Jasiukajtis				biguns = 1;
84225c28e8Piotr Jasiukajtis			}
84325c28e8Piotr Jasiukajtis			x += stridex;
84425c28e8Piotr Jasiukajtis			y += stridey;
84525c28e8Piotr Jasiukajtis			c += stridec;
84625c28e8Piotr Jasiukajtis			i = 2;
84725c28e8Piotr Jasiukajtis			if (--n <= 0)
84825c28e8Piotr Jasiukajtis				break;
84925c28e8Piotr Jasiukajtis			goto loop2;
85025c28e8Piotr Jasiukajtis		}
85125c28e8Piotr Jasiukajtis		x2 = *x;
85225c28e8Piotr Jasiukajtis		py2 = y;
85325c28e8Piotr Jasiukajtis		pc2 = c;
85425c28e8Piotr Jasiukajtis
85525c28e8Piotr Jasiukajtis		n0 = (int) (x0 * invpio2 + half[xsb0]);
85625c28e8Piotr Jasiukajtis		n1 = (int) (x1 * invpio2 + half[xsb1]);
85725c28e8Piotr Jasiukajtis		n2 = (int) (x2 * invpio2 + half[xsb2]);
85825c28e8Piotr Jasiukajtis		fn0 = (double) n0;
85925c28e8Piotr Jasiukajtis		fn1 = (double) n1;
86025c28e8Piotr Jasiukajtis		fn2 = (double) n2;
86125c28e8Piotr Jasiukajtis		n0 &= 3;
86225c28e8Piotr Jasiukajtis		n1 &= 3;
86325c28e8Piotr Jasiukajtis		n2 &= 3;
86425c28e8Piotr Jasiukajtis		a0 = x0 - fn0 * pio2_1;
86525c28e8Piotr Jasiukajtis		a1 = x1 - fn1 * pio2_1;
86625c28e8Piotr Jasiukajtis		a2 = x2 - fn2 * pio2_1;
86725c28e8Piotr Jasiukajtis		w0 = fn0 * pio2_2;
86825c28e8Piotr Jasiukajtis		w1 = fn1 * pio2_2;
86925c28e8Piotr Jasiukajtis		w2 = fn2 * pio2_2;
87025c28e8Piotr Jasiukajtis		x0 = a0 - w0;
87125c28e8Piotr Jasiukajtis		x1 = a1 - w1;
87225c28e8Piotr Jasiukajtis		x2 = a2 - w2;
87325c28e8Piotr Jasiukajtis		y0 = (a0 - x0) - w0;
87425c28e8Piotr Jasiukajtis		y1 = (a1 - x1) - w1;
87525c28e8Piotr Jasiukajtis		y2 = (a2 - x2) - w2;
87625c28e8Piotr Jasiukajtis		a0 = x0;
87725c28e8Piotr Jasiukajtis		a1 = x1;
87825c28e8Piotr Jasiukajtis		a2 = x2;
87925c28e8Piotr Jasiukajtis		w0 = fn0 * pio2_3 - y0;
88025c28e8Piotr Jasiukajtis		w1 = fn1 * pio2_3 - y1;
88125c28e8Piotr Jasiukajtis		w2 = fn2 * pio2_3 - y2;
88225c28e8Piotr Jasiukajtis		x0 = a0 - w0;
88325c28e8Piotr Jasiukajtis		x1 = a1 - w1;
88425c28e8Piotr Jasiukajtis		x2 = a2 - w2;
88525c28e8Piotr Jasiukajtis		y0 = (a0 - x0) - w0;
88625c28e8Piotr Jasiukajtis		y1 = (a1 - x1) - w1;
88725c28e8Piotr Jasiukajtis		y2 = (a2 - x2) - w2;
88825c28e8Piotr Jasiukajtis		a0 = x0;
88925c28e8Piotr Jasiukajtis		a1 = x1;
89025c28e8Piotr Jasiukajtis		a2 = x2;
89125c28e8Piotr Jasiukajtis		w0 = fn0 * pio2_3t - y0;
89225c28e8Piotr Jasiukajtis		w1 = fn1 * pio2_3t - y1;
89325c28e8Piotr Jasiukajtis		w2 = fn2 * pio2_3t - y2;
89425c28e8Piotr Jasiukajtis		x0 = a0 - w0;
89525c28e8Piotr Jasiukajtis		x1 = a1 - w1;
89625c28e8Piotr Jasiukajtis		x2 = a2 - w2;
89725c28e8Piotr Jasiukajtis		y0 = (a0 - x0) - w0;
89825c28e8Piotr Jasiukajtis		y1 = (a1 - x1) - w1;
89925c28e8Piotr Jasiukajtis		y2 = (a2 - x2) - w2;
90025c28e8Piotr Jasiukajtis		xsb2 = HI(&x2);
90125c28e8Piotr Jasiukajtis		i = ((xsb2 & ~0x80000000) - 0x3fc40000) >> 31;
90225c28e8Piotr Jasiukajtis		xsb1 = HI(&x1);
90325c28e8Piotr Jasiukajtis		i |= (((xsb1 & ~0x80000000) - 0x3fc40000) >> 30) & 2;
90425c28e8Piotr Jasiukajtis		xsb0 = HI(&x0);
90525c28e8Piotr Jasiukajtis		i |= (((xsb0 & ~0x80000000) - 0x3fc40000) >> 29) & 4;
90625c28e8Piotr Jasiukajtis		switch (i)
90725c28e8Piotr Jasiukajtis		{
90825c28e8Piotr Jasiukajtis			double