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