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