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