xref: /illumos-gate/usr/src/lib/libmvec/common/__vsincos.c (revision c81a25e9d3950dc5fab08d21f8be56d463b32c7a)
1 /*
2  * CDDL HEADER START
3  *
4  * The contents of this file are subject to the terms of the
5  * Common Development and Distribution License (the "License").
6  * You may not use this file except in compliance with the License.
7  *
8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9  * or http://www.opensolaris.org/os/licensing.
10  * See the License for the specific language governing permissions
11  * and limitations under the License.
12  *
13  * When distributing Covered Code, include this CDDL HEADER in each
14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15  * If applicable, add the following below this CDDL HEADER, with the
16  * fields enclosed by brackets "[]" replaced with your own identifying
17  * information: Portions Copyright [yyyy] [name of copyright owner]
18  *
19  * CDDL HEADER END
20  */
21 
22 /*
23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
24  */
25 /*
26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
27  * Use is subject to license terms.
28  */
29 
30 #include <sys/isa_defs.h>
31 #include <sys/ccompile.h>
32 
33 #ifdef _LITTLE_ENDIAN
34 #define HI(x)	*(1+(int*)x)
35 #define LO(x)	*(unsigned*)x
36 #else
37 #define HI(x)	*(int*)x
38 #define LO(x)	*(1+(unsigned*)x)
39 #endif
40 
41 #ifdef __RESTRICT
42 #define restrict _Restrict
43 #else
44 #define restrict
45 #endif
46 
47 /*
48  * vsincos.c
49  *
50  * Vector sine and cosine function.  Just slight modifications to vcos.c.
51  */
52 
53 extern const double __vlibm_TBL_sincos_hi[], __vlibm_TBL_sincos_lo[];
54 
55 static const double
56 	half[2]	= { 0.5, -0.5 },
57 	one		= 1.0,
58 	invpio2 = 0.636619772367581343075535,  /* 53 bits of pi/2 */
59 	pio2_1	= 1.570796326734125614166,  /* first 33 bits of pi/2 */
60 	pio2_2	= 6.077100506303965976596e-11, /* second 33 bits of pi/2 */
61 	pio2_3	= 2.022266248711166455796e-21, /* third 33 bits of pi/2 */
62 	pio2_3t	= 8.478427660368899643959e-32, /* pi/2 - pio2_3 */
63 	pp1		= -1.666666666605760465276263943134982554676e-0001,
64 	pp2		=  8.333261209690963126718376566146180944442e-0003,
65 	qq1		= -4.999999999977710986407023955908711557870e-0001,
66 	qq2		=  4.166654863857219350645055881018842089580e-0002,
67 	poly1[2]= { -1.666666666666629669805215138920301589656e-0001,
68 				-4.999999999999931701464060878888294524481e-0001 },
69 	poly2[2]= {  8.333333332390951295683993455280336376663e-0003,
70 				 4.166666666394861917535640593963708222319e-0002 },
71 	poly3[2]= { -1.984126237997976692791551778230098403960e-0004,
72 				-1.388888552656142867832756687736851681462e-0003 },
73 	poly4[2]= {  2.753403624854277237649987622848330351110e-0006,
74 				 2.478519423681460796618128289454530524759e-0005 };
75 
76 /* Don't __ the following; acomp will handle it */
77 extern double fabs(double);
78 extern void __vlibm_vsincos_big(int, double *, int, double *, int, double *, int, int);
79 
80 /*
81  * y[i*stridey] := sin( x[i*stridex] ), for i = 0..n.
82  * c[i*stridec] := cos( x[i*stridex] ), for i = 0..n.
83  *
84  * Calls __vlibm_vsincos_big to handle all elts which have abs >~ 1.647e+06.
85  * Argument reduction is done here for elts pi/4 < arg < 1.647e+06.
86  *
87  * elts < 2^-27 use the approximation 1.0 ~ cos(x).
88  */
89 void
__vsincos(int n,double * restrict x,int stridex,double * restrict y,int stridey,double * restrict c,int stridec)90 __vsincos(int n, double * restrict x, int stridex,
91 				double * restrict y, int stridey,
92 				double * restrict c, int stridec)
93 {
94 	double		x0_or_one[4], x1_or_one[4], x2_or_one[4];
95 	double		y0_or_zero[4], y1_or_zero[4], y2_or_zero[4];
96 	double		x0, x1, x2,
97 			*py0, *py1, *py2,
98 			*pc0, *pc1, *pc2,
99 			*xsave, *ysave, *csave;
100 	unsigned	hx0, hx1, hx2, xsb0, xsb1, xsb2;
101 	int		i, biguns, nsave, sxsave, sysave, scsave;
102 	volatile int	v __unused;
103 	nsave = n;
104 	xsave = x;
105 	sxsave = stridex;
106 	ysave = y;
107 	sysave = stridey;
108 	csave = c;
109 	scsave = stridec;
110 	biguns = 0;
111 
112 	do /* MAIN LOOP */
113 	{
114 
115 		/* Gotos here so _break_ exits MAIN LOOP. */
116 LOOP0:  /* Find first arg in right range. */
117 		xsb0 = HI(x); /* get most significant word */
118 		hx0 = xsb0 & ~0x80000000; /* mask off sign bit */
119 		if (hx0 > 0x3fe921fb) {
120 			/* Too big: arg reduction needed, so leave for second part */
121 			biguns = 1;
122 			x += stridex;
123 			y += stridey;
124 			c += stridec;
125 			i = 0;
126 			if (--n <= 0)
127 				break;
128 			goto LOOP0;
129 		}
130 		if (hx0 < 0x3e400000) {
131 			/* Too small.  cos x ~ 1, sin x ~ x. */
132 			v = *x;
133 			*c = 1.0;
134 			*y = *x;
135 			x += stridex;
136 			y += stridey;
137 			c += stridec;
138 			i = 0;
139 			if (--n <= 0)
140 				break;
141 			goto LOOP0;
142 		}
143 		x0 = *x;
144 		py0 = y;
145 		pc0 = c;
146 		x += stridex;
147 		y += stridey;
148 		c += stridec;
149 		i = 1;
150 		if (--n <= 0)
151 			break;
152 
153 LOOP1: /* Get second arg, same as above. */
154 		xsb1 = HI(x);
155 		hx1 = xsb1 & ~0x80000000;
156 		if (hx1 > 0x3fe921fb)
157 		{
158 			biguns = 1;
159 			x += stridex;
160 			y += stridey;
161 			c += stridec;
162 			i = 1;
163 			if (--n <= 0)
164 				break;
165 			goto LOOP1;
166 		}
167 		if (hx1 < 0x3e400000)
168 		{
169 			v = *x;
170 			*c = 1.0;
171 			*y = *x;
172 			x += stridex;
173 			y += stridey;
174 			c += stridec;
175 			i = 1;
176 			if (--n <= 0)
177 				break;
178 			goto LOOP1;
179 		}
180 		x1 = *x;
181 		py1 = y;
182 		pc1 = c;
183 		x += stridex;
184 		y += stridey;
185 		c += stridec;
186 		i = 2;
187 		if (--n <= 0)
188 			break;
189 
190 LOOP2: /* Get third arg, same as above. */
191 		xsb2 = HI(x);
192 		hx2 = xsb2 & ~0x80000000;
193 		if (hx2 > 0x3fe921fb)
194 		{
195 			biguns = 1;
196 			x += stridex;
197 			y += stridey;
198 			c += stridec;
199 			i = 2;
200 			if (--n <= 0)
201 				break;
202 			goto LOOP2;
203 		}
204 		if (hx2 < 0x3e400000)
205 		{
206 			v = *x;
207 			*c = 1.0;
208 			*y = *x;
209 			x += stridex;
210 			y += stridey;
211 			c += stridec;
212 			i = 2;
213 			if (--n <= 0)
214 				break;
215 			goto LOOP2;
216 		}
217 		x2 = *x;
218 		py2 = y;
219 		pc2 = c;
220 
221 		/*
222 		 * 0x3fc40000 = 5/32 ~ 0.15625
223 		 * Get msb after subtraction.  Will be 1 only if
224 		 * hx0 - 5/32 is negative.
225 		 */
226 		i = (hx2 - 0x3fc40000) >> 31;
227 		i |= ((hx1 - 0x3fc40000) >> 30) & 2;
228 		i |= ((hx0 - 0x3fc40000) >> 29) & 4;
229 		switch (i)
230 		{
231 			double		a1_0, a1_1, a1_2, a2_0, a2_1, a2_2;
232 			double		w0, w1, w2;
233 			double		t0, t1, t2, t1_0, t1_1, t1_2, t2_0, t2_1, t2_2;
234 			double		z0, z1, z2;
235 			unsigned	j0, j1, j2;
236 
237 		case 0: /* All are > 5/32 */
238 			j0 = (xsb0 + 0x4000) & 0xffff8000;
239 			j1 = (xsb1 + 0x4000) & 0xffff8000;
240 			j2 = (xsb2 + 0x4000) & 0xffff8000;
241 
242 			HI(&t0) = j0;
243 			HI(&t1) = j1;
244 			HI(&t2) = j2;
245 			LO(&t0) = 0;
246 			LO(&t1) = 0;
247 			LO(&t2) = 0;
248 
249 			x0 -= t0;
250 			x1 -= t1;
251 			x2 -= t2;
252 
253 			z0 = x0 * x0;
254 			z1 = x1 * x1;
255 			z2 = x2 * x2;
256 
257 			t0 = z0 * (qq1 + z0 * qq2);
258 			t1 = z1 * (qq1 + z1 * qq2);
259 			t2 = z2 * (qq1 + z2 * qq2);
260 
261 			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
262 			w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
263 			w2 = x2 * (one + z2 * (pp1 + z2 * pp2));
264 
265 			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
266 			j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
267 			j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
268 
269 			xsb0 = (xsb0 >> 30) & 2;
270 			xsb1 = (xsb1 >> 30) & 2;
271 			xsb2 = (xsb2 >> 30) & 2;
272 
273 			a1_0 = __vlibm_TBL_sincos_hi[j0+xsb0]; /* sin_hi(t) */
274 			a1_1 = __vlibm_TBL_sincos_hi[j1+xsb1];
275 			a1_2 = __vlibm_TBL_sincos_hi[j2+xsb2];
276 
277 			a2_0 = __vlibm_TBL_sincos_hi[j0+1];	/* cos_hi(t) */
278 			a2_1 = __vlibm_TBL_sincos_hi[j1+1];
279 			a2_2 = __vlibm_TBL_sincos_hi[j2+1];
280 				/* cos_lo(t) */
281 			t2_0 = __vlibm_TBL_sincos_lo[j0+1] - (a1_0*w0 - a2_0*t0);
282 			t2_1 = __vlibm_TBL_sincos_lo[j1+1] - (a1_1*w1 - a2_1*t1);
283 			t2_2 = __vlibm_TBL_sincos_lo[j2+1] - (a1_2*w2 - a2_2*t2);
284 
285 			*pc0 = a2_0 + t2_0;
286 			*pc1 = a2_1 + t2_1;
287 			*pc2 = a2_2 + t2_2;
288 
289 			t1_0 = a2_0*w0 + a1_0*t0;
290 			t1_1 = a2_1*w1 + a1_1*t1;
291 			t1_2 = a2_2*w2 + a1_2*t2;
292 
293 			t1_0 += __vlibm_TBL_sincos_lo[j0+xsb0]; /* sin_lo(t) */
294 			t1_1 += __vlibm_TBL_sincos_lo[j1+xsb1];
295 			t1_2 += __vlibm_TBL_sincos_lo[j2+xsb2];
296 
297 			*py0 = a1_0 + t1_0;
298 			*py1 = a1_1 + t1_1;
299 			*py2 = a1_2 + t1_2;
300 
301 			break;
302 
303 		case 1:
304 			j0 = (xsb0 + 0x4000) & 0xffff8000;
305 			j1 = (xsb1 + 0x4000) & 0xffff8000;
306 			HI(&t0) = j0;
307 			HI(&t1) = j1;
308 			LO(&t0) = 0;
309 			LO(&t1) = 0;
310 			x0 -= t0;
311 			x1 -= t1;
312 			z0 = x0 * x0;
313 			z1 = x1 * x1;
314 			z2 = x2 * x2;
315 			t0 = z0 * (qq1 + z0 * qq2);
316 			t1 = z1 * (qq1 + z1 * qq2);
317 			t2 = z2 * (poly3[1] + z2 * poly4[1]);
318 			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
319 			w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
320 			t2 = z2 * (poly1[1] + z2 * (poly2[1] + t2));
321 			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
322 			j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
323 			xsb0 = (xsb0 >> 30) & 2;
324 			xsb1 = (xsb1 >> 30) & 2;
325 
326 			a1_0 = __vlibm_TBL_sincos_hi[j0+xsb0]; /* sin_hi(t) */
327 			a1_1 = __vlibm_TBL_sincos_hi[j1+xsb1];
328 
329 			a2_0 = __vlibm_TBL_sincos_hi[j0+1];	/* cos_hi(t) */
330 			a2_1 = __vlibm_TBL_sincos_hi[j1+1];
331 				/* cos_lo(t) */
332 			t2_0 = __vlibm_TBL_sincos_lo[j0+1] - (a1_0*w0 - a2_0*t0);
333 			t2_1 = __vlibm_TBL_sincos_lo[j1+1] - (a1_1*w1 - a2_1*t1);
334 
335 			*pc0 = a2_0 + t2_0;
336 			*pc1 = a2_1 + t2_1;
337 			*pc2 = one + t2;
338 
339 			t1_0 = a2_0*w0 + a1_0*t0;
340 			t1_1 = a2_1*w1 + a1_1*t1;
341 			t2 = z2 * (poly3[0] + z2 * poly4[0]);
342 
343 			t1_0 += __vlibm_TBL_sincos_lo[j0+xsb0]; /* sin_lo(t) */
344 			t1_1 += __vlibm_TBL_sincos_lo[j1+xsb1];
345 			t2 = z2 * (poly1[0] + z2 * (poly2[0] + t2));
346 
347 			*py0 = a1_0 + t1_0;
348 			*py1 = a1_1 + t1_1;
349 			t2 = x2 + x2 * t2;
350 			*py2 = t2;
351 
352 			break;
353 
354 		case 2:
355 			j0 = (xsb0 + 0x4000) & 0xffff8000;
356 			j2 = (xsb2 + 0x4000) & 0xffff8000;
357 			HI(&t0) = j0;
358 			HI(&t2) = j2;
359 			LO(&t0) = 0;
360 			LO(&t2) = 0;
361 			x0 -= t0;
362 			x2 -= t2;
363 			z0 = x0 * x0;
364 			z1 = x1 * x1;
365 			z2 = x2 * x2;
366 			t0 = z0 * (qq1 + z0 * qq2);
367 			t1 = z1 * (poly3[1] + z1 * poly4[1]);
368 			t2 = z2 * (qq1 + z2 * qq2);
369 			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
370 			t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1));
371 			w2 = x2 * (one + z2 * (pp1 + z2 * pp2));
372 			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
373 			j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
374 			xsb0 = (xsb0 >> 30) & 2;
375 			xsb2 = (xsb2 >> 30) & 2;
376 
377 			a1_0 = __vlibm_TBL_sincos_hi[j0+xsb0]; /* sin_hi(t) */
378 			a1_2 = __vlibm_TBL_sincos_hi[j2+xsb2];
379 
380 			a2_0 = __vlibm_TBL_sincos_hi[j0+1];	/* cos_hi(t) */
381 			a2_2 = __vlibm_TBL_sincos_hi[j2+1];
382 				/* cos_lo(t) */
383 			t2_0 = __vlibm_TBL_sincos_lo[j0+1] - (a1_0*w0 - a2_0*t0);
384 			t2_2 = __vlibm_TBL_sincos_lo[j2+1] - (a1_2*w2 - a2_2*t2);
385 
386 			*pc0 = a2_0 + t2_0;
387 			*pc1 = one + t1;
388 			*pc2 = a2_2 + t2_2;
389 
390 			t1_0 = a2_0*w0 + a1_0*t0;
391 			t1 = z1 * (poly3[0] + z1 * poly4[0]);
392 			t1_2 = a2_2*w2 + a1_2*t2;
393 
394 			t1_0 += __vlibm_TBL_sincos_lo[j0+xsb0]; /* sin_lo(t) */
395 			t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1));
396 			t1_2 += __vlibm_TBL_sincos_lo[j2+xsb2];
397 
398 			*py0 = a1_0 + t1_0;
399 			t1 = x1 + x1 * t1;
400 			*py1 = t1;
401 			*py2 = a1_2 + t1_2;
402 
403 			break;
404 
405 		case 3:
406 			j0 = (xsb0 + 0x4000) & 0xffff8000;
407 			HI(&t0) = j0;
408 			LO(&t0) = 0;
409 			x0 -= t0;
410 			z0 = x0 * x0;
411 			z1 = x1 * x1;
412 			z2 = x2 * x2;
413 			t0 = z0 * (qq1 + z0 * qq2);
414 			t1 = z1 * (poly3[1] + z1 * poly4[1]);
415 			t2 = z2 * (poly3[1] + z2 * poly4[1]);
416 			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
417 			t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1));
418 			t2 = z2 * (poly1[1] + z2 * (poly2[1] + t2));
419 			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
420 			xsb0 = (xsb0 >> 30) & 2;
421 			a1_0 = __vlibm_TBL_sincos_hi[j0+xsb0]; /* sin_hi(t) */
422 
423 			a2_0 = __vlibm_TBL_sincos_hi[j0+1];	/* cos_hi(t) */
424 
425 			t2_0 = __vlibm_TBL_sincos_lo[j0+1] - (a1_0*w0 - a2_0*t0);
426 
427 			*pc0 = a2_0 + t2_0;
428 			*pc1 = one + t1;
429 			*pc2 = one + t2;
430 
431 			t1_0 = a2_0*w0 + a1_0*t0;
432 			t1 = z1 * (poly3[0] + z1 * poly4[0]);
433 			t2 = z2 * (poly3[0] + z2 * poly4[0]);
434 
435 			t1_0 += __vlibm_TBL_sincos_lo[j0+xsb0]; /* sin_lo(t) */
436 			t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1));
437 			t2 = z2 * (poly1[0] + z2 * (poly2[0] + t2));
438 
439 			*py0 = a1_0 + t1_0;
440 			t1 = x1 + x1 * t1;
441 			*py1 = t1;
442 			t2 = x2 + x2 * t2;
443 			*py2 = t2;
444 
445 			break;
446 
447 		case 4:
448 			j1 = (xsb1 + 0x4000) & 0xffff8000;
449 			j2 = (xsb2 + 0x4000) & 0xffff8000;
450 			HI(&t1) = j1;
451 			HI(&t2) = j2;
452 			LO(&t1) = 0;
453 			LO(&t2) = 0;
454 			x1 -= t1;
455 			x2 -= t2;
456 			z0 = x0 * x0;
457 			z1 = x1 * x1;
458 			z2 = x2 * x2;
459 			t0 = z0 * (poly3[1] + z0 * poly4[1]);
460 			t1 = z1 * (qq1 + z1 * qq2);
461 			t2 = z2 * (qq1 + z2 * qq2);
462 			t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0));
463 			w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
464 			w2 = x2 * (one + z2 * (pp1 + z2 * pp2));
465 			j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
466 			j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
467 			xsb1 = (xsb1 >> 30) & 2;
468 			xsb2 = (xsb2 >> 30) & 2;
469 
470 			a1_1 = __vlibm_TBL_sincos_hi[j1+xsb1];
471 			a1_2 = __vlibm_TBL_sincos_hi[j2+xsb2];
472 
473 			a2_1 = __vlibm_TBL_sincos_hi[j1+1];
474 			a2_2 = __vlibm_TBL_sincos_hi[j2+1];
475 				/* cos_lo(t) */
476 			t2_1 = __vlibm_TBL_sincos_lo[j1+1] - (a1_1*w1 - a2_1*t1);
477 			t2_2 = __vlibm_TBL_sincos_lo[j2+1] - (a1_2*w2 - a2_2*t2);
478 
479 			*pc0 = one + t0;
480 			*pc1 = a2_1 + t2_1;
481 			*pc2 = a2_2 + t2_2;
482 
483 			t0 = z0 * (poly3[0] + z0 * poly4[0]);
484 			t1_1 = a2_1*w1 + a1_1*t1;
485 			t1_2 = a2_2*w2 + a1_2*t2;
486 
487 			t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0));
488 			t1_1 += __vlibm_TBL_sincos_lo[j1+xsb1];
489 			t1_2 += __vlibm_TBL_sincos_lo[j2+xsb2];
490 
491 			t0 = x0 + x0 * t0;
492 			*py0 = t0;
493 			*py1 = a1_1 + t1_1;
494 			*py2 = a1_2 + t1_2;
495 
496 			break;
497 
498 		case 5:
499 			j1 = (xsb1 + 0x4000) & 0xffff8000;
500 			HI(&t1) = j1;
501 			LO(&t1) = 0;
502 			x1 -= t1;
503 			z0 = x0 * x0;
504 			z1 = x1 * x1;
505 			z2 = x2 * x2;
506 			t0 = z0 * (poly3[1] + z0 * poly4[1]);
507 			t1 = z1 * (qq1 + z1 * qq2);
508 			t2 = z2 * (poly3[1] + z2 * poly4[1]);
509 			t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0));
510 			w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
511 			t2 = z2 * (poly1[1] + z2 * (poly2[1] + t2));
512 			j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
513 			xsb1 = (xsb1 >> 30) & 2;
514 
515 			a1_1 = __vlibm_TBL_sincos_hi[j1+xsb1];
516 
517 			a2_1 = __vlibm_TBL_sincos_hi[j1+1];
518 
519 			t2_1 = __vlibm_TBL_sincos_lo[j1+1] - (a1_1*w1 - a2_1*t1);
520 
521 			*pc0 = one + t0;
522 			*pc1 = a2_1 + t2_1;
523 			*pc2 = one + t2;
524 
525 			t0 = z0 * (poly3[0] + z0 * poly4[0]);
526 			t1_1 = a2_1*w1 + a1_1*t1;
527 			t2 = z2 * (poly3[0] + z2 * poly4[0]);
528 
529 			t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0));
530 			t1_1 += __vlibm_TBL_sincos_lo[j1+xsb1];
531 			t2 = z2 * (poly1[0] + z2 * (poly2[0] + t2));
532 
533 			t0 = x0 + x0 * t0;
534 			*py0 = t0;
535 			*py1 = a1_1 + t1_1;
536 			t2 = x2 + x2 * t2;
537 			*py2 = t2;
538 
539 			break;
540 
541 		case 6:
542 			j2 = (xsb2 + 0x4000) & 0xffff8000;
543 			HI(&t2) = j2;
544 			LO(&t2) = 0;
545 			x2 -= t2;
546 			z0 = x0 * x0;
547 			z1 = x1 * x1;
548 			z2 = x2 * x2;
549 			t0 = z0 * (poly3[1] + z0 * poly4[1]);
550 			t1 = z1 * (poly3[1] + z1 * poly4[1]);
551 			t2 = z2 * (qq1 + z2 * qq2);
552 			t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0));
553 			t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1));
554 			w2 = x2 * (one + z2 * (pp1 + z2 * pp2));
555 			j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
556 			xsb2 = (xsb2 >> 30) & 2;
557 			a1_2 = __vlibm_TBL_sincos_hi[j2+xsb2];
558 
559 			a2_2 = __vlibm_TBL_sincos_hi[j2+1];
560 
561 			t2_2 = __vlibm_TBL_sincos_lo[j2+1] - (a1_2*w2 - a2_2*t2);
562 
563 			*pc0 = one + t0;
564 			*pc1 = one + t1;
565 			*pc2 = a2_2 + t2_2;
566 
567 			t0 = z0 * (poly3[0] + z0 * poly4[0]);
568 			t1 = z1 * (poly3[0] + z1 * poly4[0]);
569 			t1_2 = a2_2*w2 + a1_2*t2;
570 
571 			t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0));
572 			t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1));
573 			t1_2 += __vlibm_TBL_sincos_lo[j2+xsb2];
574 
575 			t0 = x0 + x0 * t0;
576 			*py0 = t0;
577 			t1 = x1 + x1 * t1;
578 			*py1 = t1;
579 			*py2 = a1_2 + t1_2;
580 
581 			break;
582 
583 		case 7: /* All are < 5/32 */
584 			z0 = x0 * x0;
585 			z1 = x1 * x1;
586 			z2 = x2 * x2;
587 			t0 = z0 * (poly3[1] + z0 * poly4[1]);
588 			t1 = z1 * (poly3[1] + z1 * poly4[1]);
589 			t2 = z2 * (poly3[1] + z2 * poly4[1]);
590 			t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0));
591 			t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1));
592 			t2 = z2 * (poly1[1] + z2 * (poly2[1] + t2));
593 			*pc0 = one + t0;
594 			*pc1 = one + t1;
595 			*pc2 = one + t2;
596 			t0 = z0 * (poly3[0] + z0 * poly4[0]);
597 			t1 = z1 * (poly3[0] + z1 * poly4[0]);
598 			t2 = z2 * (poly3[0] + z2 * poly4[0]);
599 			t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0));
600 			t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1));
601 			t2 = z2 * (poly1[0] + z2 * (poly2[0] + t2));
602 			t0 = x0 + x0 * t0;
603 			t1 = x1 + x1 * t1;
604 			t2 = x2 + x2 * t2;
605 			*py0 = t0;
606 			*py1 = t1;
607 			*py2 = t2;
608 			break;
609 		}
610 
611 		x += stridex;
612 		y += stridey;
613 		c += stridec;
614 		i = 0;
615 	} while (--n > 0); /* END MAIN LOOP */
616 
617 	/*
618 	 * CLEAN UP last 0, 1, or 2 elts.
619 	 */
620 	if (i > 0) /* Clean up elts at tail.  i < 3. */
621 	{
622 		double		a1_0, a1_1, a2_0, a2_1;
623 		double		w0, w1;
624 		double		t0, t1, t1_0, t1_1, t2_0, t2_1;
625 		double		z0, z1;
626 		unsigned	j0, j1;
627 
628 		if (i > 1)
629 		{
630 			if (hx1 < 0x3fc40000)
631 			{
632 				z1 = x1 * x1;
633 				t1 = z1 * (poly3[1] + z1 * poly4[1]);
634 				t1 = z1 * (poly1[1] + z1 * (poly2[1] + t1));
635 				t1 = one + t1;
636 				*pc1 = t1;
637 				t1 = z1 * (poly3[0] + z1 * poly4[0]);
638 				t1 = z1 * (poly1[0] + z1 * (poly2[0] + t1));
639 				t1 = x1 + x1 * t1;
640 				*py1 = t1;
641 			}
642 			else
643 			{
644 				j1 = (xsb1 + 0x4000) & 0xffff8000;
645 				HI(&t1) = j1;
646 				LO(&t1) = 0;
647 				x1 -= t1;
648 				z1 = x1 * x1;
649 				t1 = z1 * (qq1 + z1 * qq2);
650 				w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
651 				j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
652 				xsb1 = (xsb1 >> 30) & 2;
653 				a1_1 = __vlibm_TBL_sincos_hi[j1+xsb1];
654 				a2_1 = __vlibm_TBL_sincos_hi[j1+1];
655 				t2_1 = __vlibm_TBL_sincos_lo[j1+1] - (a1_1*w1 - a2_1*t1);
656 				*pc1 = a2_1 + t2_1;
657 				t1_1 = a2_1*w1 + a1_1*t1;
658 				t1_1 += __vlibm_TBL_sincos_lo[j1+xsb1];
659 				*py1 = a1_1 + t1_1;
660 			}
661 		}
662 		if (hx0 < 0x3fc40000)
663 		{
664 			z0 = x0 * x0;
665 			t0 = z0 * (poly3[1] + z0 * poly4[1]);
666 			t0 = z0 * (poly1[1] + z0 * (poly2[1] + t0));
667 			t0 = one + t0;
668 			*pc0 = t0;
669 			t0 = z0 * (poly3[0] + z0 * poly4[0]);
670 			t0 = z0 * (poly1[0] + z0 * (poly2[0] + t0));
671 			t0 = x0 + x0 * t0;
672 			*py0 = t0;
673 		}
674 		else
675 		{
676 			j0 = (xsb0 + 0x4000) & 0xffff8000;
677 			HI(&t0) = j0;
678 			LO(&t0) = 0;
679 			x0 -= t0;
680 			z0 = x0 * x0;
681 			t0 = z0 * (qq1 + z0 * qq2);
682 			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
683 			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
684 			xsb0 = (xsb0 >> 30) & 2;
685 			a1_0 = __vlibm_TBL_sincos_hi[j0+xsb0]; /* sin_hi(t) */
686 			a2_0 = __vlibm_TBL_sincos_hi[j0+1];	/* cos_hi(t) */
687 			t2_0 = __vlibm_TBL_sincos_lo[j0+1] - (a1_0*w0 - a2_0*t0);
688 			*pc0 = a2_0 + t2_0;
689 			t1_0 = a2_0*w0 + a1_0*t0;
690 			t1_0 += __vlibm_TBL_sincos_lo[j0+xsb0]; /* sin_lo(t) */
691 			*py0 = a1_0 + t1_0;
692 		}
693 	} /* END CLEAN UP */
694 
695 	if (!biguns)
696 		return;
697 
698 	/*
699 	 * Take care of BIGUNS.
700 	 */
701 	n = nsave;
702 	x = xsave;
703 	stridex = sxsave;
704 	y = ysave;
705 	stridey = sysave;
706 	c = csave;
707 	stridec = scsave;
708 	biguns = 0;
709 
710 	x0_or_one[1] = 1.0;
711 	x1_or_one[1] = 1.0;
712 	x2_or_one[1] = 1.0;
713 	x0_or_one[3] = -1.0;
714 	x1_or_one[3] = -1.0;
715 	x2_or_one[3] = -1.0;
716 	y0_or_zero[1] = 0.0;
717 	y1_or_zero[1] = 0.0;
718 	y2_or_zero[1] = 0.0;
719 	y0_or_zero[3] = 0.0;
720 	y1_or_zero[3] = 0.0;
721 	y2_or_zero[3] = 0.0;
722 
723 	do
724 	{
725 		double		fn0, fn1, fn2, a0, a1, a2, w0, w1, w2, y0, y1, y2;
726 		unsigned	hx;
727 		int			n0, n1, n2;
728 
729 		/*
730 		 * Find 3 more to work on: Not already done, not too big.
731 		 */
732 loop0:
733 		hx = HI(x);
734 		xsb0 = hx >> 31;
735 		hx &= ~0x80000000;
736 		if (hx <= 0x3fe921fb) /* Done above. */
737 		{
738 			x += stridex;
739 			y += stridey;
740 			c += stridec;
741 			i = 0;
742 			if (--n <= 0)
743 				break;
744 			goto loop0;
745 		}
746 		if (hx > 0x413921fb) /* (1.6471e+06) Too big: leave it. */
747 		{
748 			if (hx >= 0x7ff00000) /* Inf or NaN */
749 			{
750 				x0 = *x;
751 				*y = x0 - x0;
752 				*c = x0 - x0;
753 			}
754 			else {
755 				biguns = 1;
756 			}
757 			x += stridex;
758 			y += stridey;
759 			c += stridec;
760 			i = 0;
761 			if (--n <= 0)
762 				break;
763 			goto loop0;
764 		}
765 		x0 = *x;
766 		py0 = y;
767 		pc0 = c;
768 		x += stridex;
769 		y += stridey;
770 		c += stridec;
771 		i = 1;
772 		if (--n <= 0)
773 			break;
774 
775 loop1:
776 		hx = HI(x);
777 		xsb1 = hx >> 31;
778 		hx &= ~0x80000000;
779 		if (hx <= 0x3fe921fb)
780 		{
781 			x += stridex;
782 			y += stridey;
783 			c += stridec;
784 			i = 1;
785 			if (--n <= 0)
786 				break;
787 			goto loop1;
788 		}
789 		if (hx > 0x413921fb)
790 		{
791 			if (hx >= 0x7ff00000)
792 			{
793 				x1 = *x;
794 				*y = x1 - x1;
795 				*c = x1 - x1;
796 			}
797 			else {
798 				biguns = 1;
799 			}
800 			x += stridex;
801 			y += stridey;
802 			c += stridec;
803 			i = 1;
804 			if (--n <= 0)
805 				break;
806 			goto loop1;
807 		}
808 		x1 = *x;
809 		py1 = y;
810 		pc1 = c;
811 		x += stridex;
812 		y += stridey;
813 		c += stridec;
814 		i = 2;
815 		if (--n <= 0)
816 			break;
817 
818 loop2:
819 		hx = HI(x);
820 		xsb2 = hx >> 31;
821 		hx &= ~0x80000000;
822 		if (hx <= 0x3fe921fb)
823 		{
824 			x += stridex;
825 			y += stridey;
826 			c += stridec;
827 			i = 2;
828 			if (--n <= 0)
829 				break;
830 			goto loop2;
831 		}
832 		if (hx > 0x413921fb)
833 		{
834 			if (hx >= 0x7ff00000)
835 			{
836 				x2 = *x;
837 				*y = x2 - x2;
838 				*c = x2 - x2;
839 			}
840 			else {
841 				biguns = 1;
842 			}
843 			x += stridex;
844 			y += stridey;
845 			c += stridec;
846 			i = 2;
847 			if (--n <= 0)
848 				break;
849 			goto loop2;
850 		}
851 		x2 = *x;
852 		py2 = y;
853 		pc2 = c;
854 
855 		n0 = (int) (x0 * invpio2 + half[xsb0]);
856 		n1 = (int) (x1 * invpio2 + half[xsb1]);
857 		n2 = (int) (x2 * invpio2 + half[xsb2]);
858 		fn0 = (double) n0;
859 		fn1 = (double) n1;
860 		fn2 = (double) n2;
861 		n0 &= 3;
862 		n1 &= 3;
863 		n2 &= 3;
864 		a0 = x0 - fn0 * pio2_1;
865 		a1 = x1 - fn1 * pio2_1;
866 		a2 = x2 - fn2 * pio2_1;
867 		w0 = fn0 * pio2_2;
868 		w1 = fn1 * pio2_2;
869 		w2 = fn2 * pio2_2;
870 		x0 = a0 - w0;
871 		x1 = a1 - w1;
872 		x2 = a2 - w2;
873 		y0 = (a0 - x0) - w0;
874 		y1 = (a1 - x1) - w1;
875 		y2 = (a2 - x2) - w2;
876 		a0 = x0;
877 		a1 = x1;
878 		a2 = x2;
879 		w0 = fn0 * pio2_3 - y0;
880 		w1 = fn1 * pio2_3 - y1;
881 		w2 = fn2 * pio2_3 - y2;
882 		x0 = a0 - w0;
883 		x1 = a1 - w1;
884 		x2 = a2 - w2;
885 		y0 = (a0 - x0) - w0;
886 		y1 = (a1 - x1) - w1;
887 		y2 = (a2 - x2) - w2;
888 		a0 = x0;
889 		a1 = x1;
890 		a2 = x2;
891 		w0 = fn0 * pio2_3t - y0;
892 		w1 = fn1 * pio2_3t - y1;
893 		w2 = fn2 * pio2_3t - y2;
894 		x0 = a0 - w0;
895 		x1 = a1 - w1;
896 		x2 = a2 - w2;
897 		y0 = (a0 - x0) - w0;
898 		y1 = (a1 - x1) - w1;
899 		y2 = (a2 - x2) - w2;
900 		xsb2 = HI(&x2);
901 		i = ((xsb2 & ~0x80000000) - 0x3fc40000) >> 31;
902 		xsb1 = HI(&x1);
903 		i |= (((xsb1 & ~0x80000000) - 0x3fc40000) >> 30) & 2;
904 		xsb0 = HI(&x0);
905 		i |= (((xsb0 & ~0x80000000) - 0x3fc40000) >> 29) & 4;
906 		switch (i)
907 		{
908 			double		a1_0, a1_1, a1_2, a2_0, a2_1, a2_2;
909 			double		t0, t1, t2, t1_0, t1_1, t1_2, t2_0, t2_1, t2_2;
910 			double		z0, z1, z2;
911 			unsigned	j0, j1, j2;
912 
913 		case 0:
914 			j0 = (xsb0 + 0x4000) & 0xffff8000;
915 			j1 = (xsb1 + 0x4000) & 0xffff8000;
916 			j2 = (xsb2 + 0x4000) & 0xffff8000;
917 			HI(&t0) = j0;
918 			HI(&t1) = j1;
919 			HI(&t2) = j2;
920 			LO(&t0) = 0;
921 			LO(&t1) = 0;
922 			LO(&t2) = 0;
923 			x0 = (x0 - t0) + y0;
924 			x1 = (x1 - t1) + y1;
925 			x2 = (x2 - t2) + y2;
926 			z0 = x0 * x0;
927 			z1 = x1 * x1;
928 			z2 = x2 * x2;
929 			t0 = z0 * (qq1 + z0 * qq2);
930 			t1 = z1 * (qq1 + z1 * qq2);
931 			t2 = z2 * (qq1 + z2 * qq2);
932 			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
933 			w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
934 			w2 = x2 * (one + z2 * (pp1 + z2 * pp2));
935 			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
936 			j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
937 			j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
938 			xsb0 = (xsb0 >> 30) & 2;
939 			xsb1 = (xsb1 >> 30) & 2;
940 			xsb2 = (xsb2 >> 30) & 2;
941 			n0 ^= (xsb0 & ~(n0 << 1));
942 			n1 ^= (xsb1 & ~(n1 << 1));
943 			n2 ^= (xsb2 & ~(n2 << 1));
944 			xsb0 |= 1;
945 			xsb1 |= 1;
946 			xsb2 |= 1;
947 
948 			a1_0 = __vlibm_TBL_sincos_hi[j0+n0];
949 			a1_1 = __vlibm_TBL_sincos_hi[j1+n1];
950 			a1_2 = __vlibm_TBL_sincos_hi[j2+n2];
951 
952 			a2_0 = __vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)];
953 			a2_1 = __vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)];
954 			a2_2 = __vlibm_TBL_sincos_hi[j2+((n2+xsb2)&3)];
955 
956 			t2_0 = __vlibm_TBL_sincos_lo[j0+((n0+xsb0)&3)] - (a1_0*w0 - a2_0*t0);
957 			t2_1 = __vlibm_TBL_sincos_lo[j1+((n1+xsb1)&3)] - (a1_1*w1 - a2_1*t1);
958 			t2_2 = __vlibm_TBL_sincos_lo[j2+((n2+xsb2)&3)] - (a1_2*w2 - a2_2*t2);
959 
960 			w0 *= a2_0;
961 			w1 *= a2_1;
962 			w2 *= a2_2;
963 
964 			*pc0 = a2_0 + t2_0;
965 			*pc1 = a2_1 + t2_1;
966 			*pc2 = a2_2 + t2_2;
967 
968 			t1_0 = w0 + a1_0*t0;
969 			t1_1 = w1 + a1_1*t1;
970 			t1_2 = w2 + a1_2*t2;
971 
972 			t1_0 += __vlibm_TBL_sincos_lo[j0+n0];
973 			t1_1 += __vlibm_TBL_sincos_lo[j1+n1];
974 			t1_2 += __vlibm_TBL_sincos_lo[j2+n2];
975 
976 			*py0 = a1_0 + t1_0;
977 			*py1 = a1_1 + t1_1;
978 			*py2 = a1_2 + t1_2;
979 
980 			break;
981 
982 		case 1:
983 			j0 = (xsb0 + 0x4000) & 0xffff8000;
984 			j1 = (xsb1 + 0x4000) & 0xffff8000;
985 			j2 = n2 & 1;
986 			HI(&t0) = j0;
987 			HI(&t1) = j1;
988 			LO(&t0) = 0;
989 			LO(&t1) = 0;
990 			x2_or_one[0] = x2;
991 			x2_or_one[2] = -x2;
992 			x0 = (x0 - t0) + y0;
993 			x1 = (x1 - t1) + y1;
994 			y2_or_zero[0] = y2;
995 			y2_or_zero[2] = -y2;
996 			z0 = x0 * x0;
997 			z1 = x1 * x1;
998 			z2 = x2 * x2;
999 			t0 = z0 * (qq1 + z0 * qq2);
1000 			t1 = z1 * (qq1 + z1 * qq2);
1001 			t2 = z2 * (poly3[j2] + z2 * poly4[j2]);
1002 			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
1003 			w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
1004 			t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2));
1005 			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
1006 			j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
1007 			xsb0 = (xsb0 >> 30) & 2;
1008 			xsb1 = (xsb1 >> 30) & 2;
1009 			n0 ^= (xsb0 & ~(n0 << 1));
1010 			n1 ^= (xsb1 & ~(n1 << 1));
1011 			xsb0 |= 1;
1012 			xsb1 |= 1;
1013 			a1_0 = __vlibm_TBL_sincos_hi[j0+n0];
1014 			a1_1 = __vlibm_TBL_sincos_hi[j1+n1];
1015 
1016 			a2_0 = __vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)];
1017 			a2_1 = __vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)];
1018 
1019 			t2_0 = __vlibm_TBL_sincos_lo[j0+((n0+xsb0)&3)] - (a1_0*w0 - a2_0*t0);
1020 			t2_1 = __vlibm_TBL_sincos_lo[j1+((n1+xsb1)&3)] - (a1_1*w1 - a2_1*t1);
1021 			t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2);
1022 
1023 			*pc0 = a2_0 + t2_0;
1024 			*pc1 = a2_1 + t2_1;
1025 			*py2 = t2;
1026 
1027 			n2 = (n2 + 1) & 3;
1028 			j2 = (j2 + 1) & 1;
1029 			t2 = z2 * (poly3[j2] + z2 * poly4[j2]);
1030 
1031 			t1_0 = a2_0*w0 + a1_0*t0;
1032 			t1_1 = a2_1*w1 + a1_1*t1;
1033 			t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2));
1034 
1035 			t1_0 += __vlibm_TBL_sincos_lo[j0+n0];
1036 			t1_1 += __vlibm_TBL_sincos_lo[j1+n1];
1037 			t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2);
1038 
1039 			*py0 = a1_0 + t1_0;
1040 			*py1 = a1_1 + t1_1;
1041 			*pc2 = t2;
1042 
1043 			break;
1044 
1045 		case 2:
1046 			j0 = (xsb0 + 0x4000) & 0xffff8000;
1047 			j1 = n1 & 1;
1048 			j2 = (xsb2 + 0x4000) & 0xffff8000;
1049 			HI(&t0) = j0;
1050 			HI(&t2) = j2;
1051 			LO(&t0) = 0;
1052 			LO(&t2) = 0;
1053 			x1_or_one[0] = x1;
1054 			x1_or_one[2] = -x1;
1055 			x0 = (x0 - t0) + y0;
1056 			y1_or_zero[0] = y1;
1057 			y1_or_zero[2] = -y1;
1058 			x2 = (x2 - t2) + y2;
1059 			z0 = x0 * x0;
1060 			z1 = x1 * x1;
1061 			z2 = x2 * x2;
1062 			t0 = z0 * (qq1 + z0 * qq2);
1063 			t1 = z1 * (poly3[j1] + z1 * poly4[j1]);
1064 			t2 = z2 * (qq1 + z2 * qq2);
1065 			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
1066 			t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1));
1067 			w2 = x2 * (one + z2 * (pp1 + z2 * pp2));
1068 			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
1069 			j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
1070 			xsb0 = (xsb0 >> 30) & 2;
1071 			xsb2 = (xsb2 >> 30) & 2;
1072 			n0 ^= (xsb0 & ~(n0 << 1));
1073 			n2 ^= (xsb2 & ~(n2 << 1));
1074 			xsb0 |= 1;
1075 			xsb2 |= 1;
1076 
1077 			a1_0 = __vlibm_TBL_sincos_hi[j0+n0];
1078 			a1_2 = __vlibm_TBL_sincos_hi[j2+n2];
1079 
1080 			a2_0 = __vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)];
1081 			a2_2 = __vlibm_TBL_sincos_hi[j2+((n2+xsb2)&3)];
1082 
1083 			t2_0 = __vlibm_TBL_sincos_lo[j0+((n0+xsb0)&3)] - (a1_0*w0 - a2_0*t0);
1084 			t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1);
1085 			t2_2 = __vlibm_TBL_sincos_lo[j2+((n2+xsb2)&3)] - (a1_2*w2 - a2_2*t2);
1086 
1087 			*pc0 = a2_0 + t2_0;
1088 			*py1 = t1;
1089 			*pc2 = a2_2 + t2_2;
1090 
1091 			n1 = (n1 + 1) & 3;
1092 			j1 = (j1 + 1) & 1;
1093 			t1 = z1 * (poly3[j1] + z1 * poly4[j1]);
1094 
1095 			t1_0 = a2_0*w0 + a1_0*t0;
1096 			t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1));
1097 			t1_2 = a2_2*w2 + a1_2*t2;
1098 
1099 			t1_0 += __vlibm_TBL_sincos_lo[j0+n0];
1100 			t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1);
1101 			t1_2 += __vlibm_TBL_sincos_lo[j2+n2];
1102 
1103 			*py0 = a1_0 + t1_0;
1104 			*pc1 = t1;
1105 			*py2 = a1_2 + t1_2;
1106 
1107 			break;
1108 
1109 		case 3:
1110 			j0 = (xsb0 + 0x4000) & 0xffff8000;
1111 			j1 = n1 & 1;
1112 			j2 = n2 & 1;
1113 			HI(&t0) = j0;
1114 			LO(&t0) = 0;
1115 			x1_or_one[0] = x1;
1116 			x1_or_one[2] = -x1;
1117 			x2_or_one[0] = x2;
1118 			x2_or_one[2] = -x2;
1119 			x0 = (x0 - t0) + y0;
1120 			y1_or_zero[0] = y1;
1121 			y1_or_zero[2] = -y1;
1122 			y2_or_zero[0] = y2;
1123 			y2_or_zero[2] = -y2;
1124 			z0 = x0 * x0;
1125 			z1 = x1 * x1;
1126 			z2 = x2 * x2;
1127 			t0 = z0 * (qq1 + z0 * qq2);
1128 			t1 = z1 * (poly3[j1] + z1 * poly4[j1]);
1129 			t2 = z2 * (poly3[j2] + z2 * poly4[j2]);
1130 			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
1131 			t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1));
1132 			t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2));
1133 			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
1134 			xsb0 = (xsb0 >> 30) & 2;
1135 			n0 ^= (xsb0 & ~(n0 << 1));
1136 			xsb0 |= 1;
1137 
1138 			a1_0 = __vlibm_TBL_sincos_hi[j0+n0];
1139 			a2_0 = __vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)];
1140 
1141 			t2_0 = __vlibm_TBL_sincos_lo[j0+((n0+xsb0)&3)] - (a1_0*w0 - a2_0*t0);
1142 			t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1);
1143 			t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2);
1144 
1145 			*pc0 = a2_0 + t2_0;
1146 			*py1 = t1;
1147 			*py2 = t2;
1148 
1149 			n1 = (n1 + 1) & 3;
1150 			n2 = (n2 + 1) & 3;
1151 			j1 = (j1 + 1) & 1;
1152 			j2 = (j2 + 1) & 1;
1153 
1154 			t1_0 = a2_0*w0 + a1_0*t0;
1155 			t1 = z1 * (poly3[j1] + z1 * poly4[j1]);
1156 			t2 = z2 * (poly3[j2] + z2 * poly4[j2]);
1157 
1158 			t1_0 += __vlibm_TBL_sincos_lo[j0+n0];
1159 			t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1));
1160 			t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2));
1161 
1162 			t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1);
1163 			t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2);
1164 
1165 			*py0 = a1_0 + t1_0;
1166 			*pc1 = t1;
1167 			*pc2 = t2;
1168 
1169 			break;
1170 
1171 		case 4:
1172 			j0 = n0 & 1;
1173 			j1 = (xsb1 + 0x4000) & 0xffff8000;
1174 			j2 = (xsb2 + 0x4000) & 0xffff8000;
1175 			HI(&t1) = j1;
1176 			HI(&t2) = j2;
1177 			LO(&t1) = 0;
1178 			LO(&t2) = 0;
1179 			x0_or_one[0] = x0;
1180 			x0_or_one[2] = -x0;
1181 			y0_or_zero[0] = y0;
1182 			y0_or_zero[2] = -y0;
1183 			x1 = (x1 - t1) + y1;
1184 			x2 = (x2 - t2) + y2;
1185 			z0 = x0 * x0;
1186 			z1 = x1 * x1;
1187 			z2 = x2 * x2;
1188 			t0 = z0 * (poly3[j0] + z0 * poly4[j0]);
1189 			t1 = z1 * (qq1 + z1 * qq2);
1190 			t2 = z2 * (qq1 + z2 * qq2);
1191 			t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0));
1192 			w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
1193 			w2 = x2 * (one + z2 * (pp1 + z2 * pp2));
1194 			j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
1195 			j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
1196 			xsb1 = (xsb1 >> 30) & 2;
1197 			xsb2 = (xsb2 >> 30) & 2;
1198 			n1 ^= (xsb1 & ~(n1 << 1));
1199 			n2 ^= (xsb2 & ~(n2 << 1));
1200 			xsb1 |= 1;
1201 			xsb2 |= 1;
1202 
1203 			a1_1 = __vlibm_TBL_sincos_hi[j1+n1];
1204 			a1_2 = __vlibm_TBL_sincos_hi[j2+n2];
1205 
1206 			a2_1 = __vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)];
1207 			a2_2 = __vlibm_TBL_sincos_hi[j2+((n2+xsb2)&3)];
1208 
1209 			t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0);
1210 			t2_1 = __vlibm_TBL_sincos_lo[j1+((n1+xsb1)&3)] - (a1_1*w1 - a2_1*t1);
1211 			t2_2 = __vlibm_TBL_sincos_lo[j2+((n2+xsb2)&3)] - (a1_2*w2 - a2_2*t2);
1212 
1213 			*py0 = t0;
1214 			*pc1 = a2_1 + t2_1;
1215 			*pc2 = a2_2 + t2_2;
1216 
1217 			n0 = (n0 + 1) & 3;
1218 			j0 = (j0 + 1) & 1;
1219 			t0 = z0 * (poly3[j0] + z0 * poly4[j0]);
1220 
1221 			t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0));
1222 			t1_1 = a2_1*w1 + a1_1*t1;
1223 			t1_2 = a2_2*w2 + a1_2*t2;
1224 
1225 			t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0);
1226 			t1_1 += __vlibm_TBL_sincos_lo[j1+n1];
1227 			t1_2 += __vlibm_TBL_sincos_lo[j2+n2];
1228 
1229 			*py1 = a1_1 + t1_1;
1230 			*py2 = a1_2 + t1_2;
1231 			*pc0 = t0;
1232 
1233 			break;
1234 
1235 		case 5:
1236 			j0 = n0 & 1;
1237 			j1 = (xsb1 + 0x4000) & 0xffff8000;
1238 			j2 = n2 & 1;
1239 			HI(&t1) = j1;
1240 			LO(&t1) = 0;
1241 			x0_or_one[0] = x0;
1242 			x0_or_one[2] = -x0;
1243 			x2_or_one[0] = x2;
1244 			x2_or_one[2] = -x2;
1245 			y0_or_zero[0] = y0;
1246 			y0_or_zero[2] = -y0;
1247 			x1 = (x1 - t1) + y1;
1248 			y2_or_zero[0] = y2;
1249 			y2_or_zero[2] = -y2;
1250 			z0 = x0 * x0;
1251 			z1 = x1 * x1;
1252 			z2 = x2 * x2;
1253 			t0 = z0 * (poly3[j0] + z0 * poly4[j0]);
1254 			t1 = z1 * (qq1 + z1 * qq2);
1255 			t2 = z2 * (poly3[j2] + z2 * poly4[j2]);
1256 			t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0));
1257 			w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
1258 			t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2));
1259 			j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
1260 			xsb1 = (xsb1 >> 30) & 2;
1261 			n1 ^= (xsb1 & ~(n1 << 1));
1262 			xsb1 |= 1;
1263 
1264 			a1_1 = __vlibm_TBL_sincos_hi[j1+n1];
1265 			a2_1 = __vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)];
1266 
1267 			t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0);
1268 			t2_1 = __vlibm_TBL_sincos_lo[j1+((n1+xsb1)&3)] - (a1_1*w1 - a2_1*t1);
1269 			t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2);
1270 
1271 			*py0 = t0;
1272 			*pc1 = a2_1 + t2_1;
1273 			*py2 = t2;
1274 
1275 			n0 = (n0 + 1) & 3;
1276 			n2 = (n2 + 1) & 3;
1277 			j0 = (j0 + 1) & 1;
1278 			j2 = (j2 + 1) & 1;
1279 
1280 			t0 = z0 * (poly3[j0] + z0 * poly4[j0]);
1281 			t1_1 = a2_1*w1 + a1_1*t1;
1282 			t2 = z2 * (poly3[j2] + z2 * poly4[j2]);
1283 
1284 			t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0));
1285 			t1_1 += __vlibm_TBL_sincos_lo[j1+n1];
1286 			t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2));
1287 
1288 			t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0);
1289 			t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2);
1290 
1291 			*pc0 = t0;
1292 			*py1 = a1_1 + t1_1;
1293 			*pc2 = t2;
1294 
1295 			break;
1296 
1297 		case 6:
1298 			j0 = n0 & 1;
1299 			j1 = n1 & 1;
1300 			j2 = (xsb2 + 0x4000) & 0xffff8000;
1301 			HI(&t2) = j2;
1302 			LO(&t2) = 0;
1303 			x0_or_one[0] = x0;
1304 			x0_or_one[2] = -x0;
1305 			x1_or_one[0] = x1;
1306 			x1_or_one[2] = -x1;
1307 			y0_or_zero[0] = y0;
1308 			y0_or_zero[2] = -y0;
1309 			y1_or_zero[0] = y1;
1310 			y1_or_zero[2] = -y1;
1311 			x2 = (x2 - t2) + y2;
1312 			z0 = x0 * x0;
1313 			z1 = x1 * x1;
1314 			z2 = x2 * x2;
1315 			t0 = z0 * (poly3[j0] + z0 * poly4[j0]);
1316 			t1 = z1 * (poly3[j1] + z1 * poly4[j1]);
1317 			t2 = z2 * (qq1 + z2 * qq2);
1318 			t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0));
1319 			t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1));
1320 			w2 = x2 * (one + z2 * (pp1 + z2 * pp2));
1321 			j2 = (((j2 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
1322 			xsb2 = (xsb2 >> 30) & 2;
1323 			n2 ^= (xsb2 & ~(n2 << 1));
1324 			xsb2 |= 1;
1325 
1326 			a1_2 = __vlibm_TBL_sincos_hi[j2+n2];
1327 			a2_2 = __vlibm_TBL_sincos_hi[j2+((n2+xsb2)&3)];
1328 
1329 			t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0);
1330 			t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1);
1331 			t2_2 = __vlibm_TBL_sincos_lo[j2+((n2+xsb2)&3)] - (a1_2*w2 - a2_2*t2);
1332 
1333 			*py0 = t0;
1334 			*py1 = t1;
1335 			*pc2 = a2_2 + t2_2;
1336 
1337 			n0 = (n0 + 1) & 3;
1338 			n1 = (n1 + 1) & 3;
1339 			j0 = (j0 + 1) & 1;
1340 			j1 = (j1 + 1) & 1;
1341 
1342 			t0 = z0 * (poly3[j0] + z0 * poly4[j0]);
1343 			t1 = z1 * (poly3[j1] + z1 * poly4[j1]);
1344 			t1_2 = a2_2*w2 + a1_2*t2;
1345 
1346 			t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0));
1347 			t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1));
1348 			t1_2 += __vlibm_TBL_sincos_lo[j2+n2];
1349 
1350 			t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0);
1351 			t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1);
1352 
1353 			*pc0 = t0;
1354 			*pc1 = t1;
1355 			*py2 = a1_2 + t1_2;
1356 
1357 			break;
1358 
1359 		case 7:
1360 			j0 = n0 & 1;
1361 			j1 = n1 & 1;
1362 			j2 = n2 & 1;
1363 			x0_or_one[0] = x0;
1364 			x0_or_one[2] = -x0;
1365 			x1_or_one[0] = x1;
1366 			x1_or_one[2] = -x1;
1367 			x2_or_one[0] = x2;
1368 			x2_or_one[2] = -x2;
1369 			y0_or_zero[0] = y0;
1370 			y0_or_zero[2] = -y0;
1371 			y1_or_zero[0] = y1;
1372 			y1_or_zero[2] = -y1;
1373 			y2_or_zero[0] = y2;
1374 			y2_or_zero[2] = -y2;
1375 			z0 = x0 * x0;
1376 			z1 = x1 * x1;
1377 			z2 = x2 * x2;
1378 			t0 = z0 * (poly3[j0] + z0 * poly4[j0]);
1379 			t1 = z1 * (poly3[j1] + z1 * poly4[j1]);
1380 			t2 = z2 * (poly3[j2] + z2 * poly4[j2]);
1381 			t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0));
1382 			t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1));
1383 			t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2));
1384 			t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0);
1385 			t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1);
1386 			t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2);
1387 			*py0 = t0;
1388 			*py1 = t1;
1389 			*py2 = t2;
1390 
1391 			n0 = (n0 + 1) & 3;
1392 			n1 = (n1 + 1) & 3;
1393 			n2 = (n2 + 1) & 3;
1394 			j0 = (j0 + 1) & 1;
1395 			j1 = (j1 + 1) & 1;
1396 			j2 = (j2 + 1) & 1;
1397 			t0 = z0 * (poly3[j0] + z0 * poly4[j0]);
1398 			t1 = z1 * (poly3[j1] + z1 * poly4[j1]);
1399 			t2 = z2 * (poly3[j2] + z2 * poly4[j2]);
1400 			t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0));
1401 			t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1));
1402 			t2 = z2 * (poly1[j2] + z2 * (poly2[j2] + t2));
1403 			t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0);
1404 			t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1);
1405 			t2 = x2_or_one[n2] + (y2_or_zero[n2] + x2_or_one[n2] * t2);
1406 			*pc0 = t0;
1407 			*pc1 = t1;
1408 			*pc2 = t2;
1409 			break;
1410 		}
1411 
1412 		x += stridex;
1413 		y += stridey;
1414 		c += stridec;
1415 		i = 0;
1416 	} while (--n > 0);
1417 
1418 	if (i > 0)
1419 	{
1420 		double		a1_0, a1_1, a2_0, a2_1;
1421 		double		t0, t1, t1_0, t1_1, t2_0, t2_1;
1422 		double		fn0, fn1, a0, a1, w0, w1, y0, y1;
1423 		double		z0, z1;
1424 		unsigned	j0, j1;
1425 		int		n0, n1;
1426 
1427 		if (i > 1)
1428 		{
1429 			n1 = (int) (x1 * invpio2 + half[xsb1]);
1430 			fn1 = (double) n1;
1431 			n1 &= 3;
1432 			a1 = x1 - fn1 * pio2_1;
1433 			w1 = fn1 * pio2_2;
1434 			x1 = a1 - w1;
1435 			y1 = (a1 - x1) - w1;
1436 			a1 = x1;
1437 			w1 = fn1 * pio2_3 - y1;
1438 			x1 = a1 - w1;
1439 			y1 = (a1 - x1) - w1;
1440 			a1 = x1;
1441 			w1 = fn1 * pio2_3t - y1;
1442 			x1 = a1 - w1;
1443 			y1 = (a1 - x1) - w1;
1444 			xsb1 = HI(&x1);
1445 			if ((xsb1 & ~0x80000000) < 0x3fc40000)
1446 			{
1447 				j1 = n1 & 1;
1448 				x1_or_one[0] = x1;
1449 				x1_or_one[2] = -x1;
1450 				y1_or_zero[0] = y1;
1451 				y1_or_zero[2] = -y1;
1452 				z1 = x1 * x1;
1453 				t1 = z1 * (poly3[j1] + z1 * poly4[j1]);
1454 				t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1));
1455 				t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1);
1456 				*py1 = t1;
1457 				n1 = (n1 + 1) & 3;
1458 				j1 = (j1 + 1) & 1;
1459 				t1 = z1 * (poly3[j1] + z1 * poly4[j1]);
1460 				t1 = z1 * (poly1[j1] + z1 * (poly2[j1] + t1));
1461 				t1 = x1_or_one[n1] + (y1_or_zero[n1] + x1_or_one[n1] * t1);
1462 				*pc1 = t1;
1463 			}
1464 			else
1465 			{
1466 				j1 = (xsb1 + 0x4000) & 0xffff8000;
1467 				HI(&t1) = j1;
1468 				LO(&t1) = 0;
1469 				x1 = (x1 - t1) + y1;
1470 				z1 = x1 * x1;
1471 				t1 = z1 * (qq1 + z1 * qq2);
1472 				w1 = x1 * (one + z1 * (pp1 + z1 * pp2));
1473 				j1 = (((j1 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
1474 				xsb1 = (xsb1 >> 30) & 2;
1475 				n1 ^= (xsb1 & ~(n1 << 1));
1476 				xsb1 |= 1;
1477 				a1_1 = __vlibm_TBL_sincos_hi[j1+n1];
1478 				a2_1 = __vlibm_TBL_sincos_hi[j1+((n1+xsb1)&3)];
1479 				t2_1 = __vlibm_TBL_sincos_lo[j1+((n1+xsb1)&3)] - (a1_1*w1 - a2_1*t1);
1480 				*pc1 = a2_1 + t2_1;
1481 				t1_1 = a2_1*w1 + a1_1*t1;
1482 				t1_1 += __vlibm_TBL_sincos_lo[j1+n1];
1483 				*py1 = a1_1 + t1_1;
1484 			}
1485 		}
1486 		n0 = (int) (x0 * invpio2 + half[xsb0]);
1487 		fn0 = (double) n0;
1488 		n0 &= 3;
1489 		a0 = x0 - fn0 * pio2_1;
1490 		w0 = fn0 * pio2_2;
1491 		x0 = a0 - w0;
1492 		y0 = (a0 - x0) - w0;
1493 		a0 = x0;
1494 		w0 = fn0 * pio2_3 - y0;
1495 		x0 = a0 - w0;
1496 		y0 = (a0 - x0) - w0;
1497 		a0 = x0;
1498 		w0 = fn0 * pio2_3t - y0;
1499 		x0 = a0 - w0;
1500 		y0 = (a0 - x0) - w0;
1501 		xsb0 = HI(&x0);
1502 		if ((xsb0 & ~0x80000000) < 0x3fc40000)
1503 		{
1504 			j0 = n0 & 1;
1505 			x0_or_one[0] = x0;
1506 			x0_or_one[2] = -x0;
1507 			y0_or_zero[0] = y0;
1508 			y0_or_zero[2] = -y0;
1509 			z0 = x0 * x0;
1510 			t0 = z0 * (poly3[j0] + z0 * poly4[j0]);
1511 			t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0));
1512 			t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0);
1513 			*py0 = t0;
1514 			n0 = (n0 + 1) & 3;
1515 			j0 = (j0 + 1) & 1;
1516 			t0 = z0 * (poly3[j0] + z0 * poly4[j0]);
1517 			t0 = z0 * (poly1[j0] + z0 * (poly2[j0] + t0));
1518 			t0 = x0_or_one[n0] + (y0_or_zero[n0] + x0_or_one[n0] * t0);
1519 			*pc0 = t0;
1520 		}
1521 		else
1522 		{
1523 			j0 = (xsb0 + 0x4000) & 0xffff8000;
1524 			HI(&t0) = j0;
1525 			LO(&t0) = 0;
1526 			x0 = (x0 - t0) + y0;
1527 			z0 = x0 * x0;
1528 			t0 = z0 * (qq1 + z0 * qq2);
1529 			w0 = x0 * (one + z0 * (pp1 + z0 * pp2));
1530 			j0 = (((j0 & ~0x80000000) - 0x3fc40000) >> 13) & ~0x3;
1531 			xsb0 = (xsb0 >> 30) & 2;
1532 			n0 ^= (xsb0 & ~(n0 << 1));
1533 			xsb0 |= 1;
1534 			a1_0 = __vlibm_TBL_sincos_hi[j0+n0];
1535 			a2_0 = __vlibm_TBL_sincos_hi[j0+((n0+xsb0)&3)];
1536 			t2_0 = __vlibm_TBL_sincos_lo[j0+((n0+xsb0)&3)] - (a1_0*w0 - a2_0*t0);
1537 			*pc0 = a2_0 + t2_0;
1538 			t1_0 = a2_0*w0 + a1_0*t0;
1539 			t1_0 += __vlibm_TBL_sincos_lo[j0+n0];
1540 			*py0 = a1_0 + t1_0;
1541 		}
1542 	}
1543 
1544 	if (biguns) {
1545 		__vlibm_vsincos_big(nsave, xsave, sxsave, ysave, sysave, csave, scsave, 0x413921fb);
1546 	}
1547 }
1548