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