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
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 }