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 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 *y = *x;
105 x += stridex;
106 y += stridey;
107 i = 0;
108 if ( --n <= 0 )
109 break;
110 goto LOOP0;
111 }
112 x0 = *x;
113 py0 = y;
114 x += stridex;
115 y += stridey;
116 i = 1;
117 if ( --n <= 0 )
118 break;
119
120 LOOP1:
121 xsb1 = HI(x);
122 hx1 = xsb1 & ~0x80000000;
123 if ( hx1 > 0x3fe921fb )
124 {
125 biguns = 2;
126 goto MEDIUM;
127 }
128 if ( hx1 < 0x3e400000 )
129 {
130 *y = *x;
131 x += stridex;
132 y += stridey;
133 i = 1;
134 if ( --n <= 0 )
135 break;
136 goto LOOP1;
137 }
138 x1 = *x;
139 py1 = y;
140 x += stridex;
141 y += stridey;
142 i = 2;
143 if ( --n <= 0 )
144 break;
145
146 LOOP2:
147 xsb2 = HI(x);
148 hx2 = xsb2 & ~0x80000000;
149 if ( hx2 > 0x3fe921fb )
150 {
151 biguns = 3;
152 goto MEDIUM;
153 }
154 if ( hx2 < 0x3e400000 )
155 {
156 *y = *x;
157 x += stridex;
158 y += stridey;
159 i = 2;
160 if ( --n <= 0 )
161 break;
162 goto LOOP2;
163 }
164 x2 = *x;
165 py2 = y;
166
167 i = ( hx0 - 0x3fc90000 ) >> 31;
168 i |= ( ( hx1 - 0x3fc90000 ) >> 30 ) & 2;
169 i |= ( ( hx2 - 0x3fc90000 ) >> 29 ) & 4;
170 switch ( i )
171 {
172 double a0, a1, a2, w0, w1, w2;
173 double t0, t1, t2, z0, z1, z2;
174 unsigned j0, j1, j2;
175
176 case 0:
177 j0 = ( xsb0 + 0x4000 ) & 0xffff8000;
178 j1 = ( xsb1 + 0x4000 ) & 0xffff8000;
179 j2 = ( xsb2 + 0x4000 ) & 0xffff8000;
180 HI(&t0) = j0;
181 HI(&t1) = j1;
182 HI(&t2) = j2;
183 LO(&t0) = 0;
184 LO(&t1) = 0;
185 LO(&t2) = 0;
186 x0 -= t0;
187 x1 -= t1;
188 x2 -= t2;
189 z0 = x0 * x0;
190 z1 = x1 * x1;
191 z2 = x2 * x2;
192 t0 = z0 * ( qq1 + z0 * qq2 );
193 t1 = z1 * ( qq1 + z1 * qq2 );
194 t2 = z2 * ( qq1 + z2 * qq2 );
195 w0 = x0 * ( one + z0 * ( pp1 + z0 * pp2 ) );
196 w1 = x1 * ( one + z1 * ( pp1 + z1 * pp2 ) );
197 w2 = x2 * ( one + z2 * ( pp1 + z2 * pp2 ) );
198 j0 = ( ( ( j0 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
199 j1 = ( ( ( j1 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
200 j2 = ( ( ( j2 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
201 xsb0 = ( xsb0 >> 30 ) & 2;
202 xsb1 = ( xsb1 >> 30 ) & 2;
203 xsb2 = ( xsb2 >> 30 ) & 2;
204 a0 = __vlibm_TBL_sincos_hi[j0+xsb0];
205 a1 = __vlibm_TBL_sincos_hi[j1+xsb1];
206 a2 = __vlibm_TBL_sincos_hi[j2+xsb2];
207 t0 = ( __vlibm_TBL_sincos_hi[j0+1] * w0 + a0 * t0 ) + __vlibm_TBL_sincos_lo[j0+xsb0];
208 t1 = ( __vlibm_TBL_sincos_hi[j1+1] * w1 + a1 * t1 ) + __vlibm_TBL_sincos_lo[j1+xsb1];
209 t2 = ( __vlibm_TBL_sincos_hi[j2+1] * w2 + a2 * t2 ) + __vlibm_TBL_sincos_lo[j2+xsb2];
210 *py0 = a0 + t0;
211 *py1 = a1 + t1;
212 *py2 = a2 + t2;
213 break;
214
215 case 1:
216 j1 = ( xsb1 + 0x4000 ) & 0xffff8000;
217 j2 = ( xsb2 + 0x4000 ) & 0xffff8000;
218 HI(&t1) = j1;
219 HI(&t2) = j2;
220 LO(&t1) = 0;
221 LO(&t2) = 0;
222 x1 -= t1;
223 x2 -= t2;
224 z0 = x0 * x0;
225 z1 = x1 * x1;
226 z2 = x2 * x2;
227 t0 = z0 * ( poly3[0] + z0 * poly4[0] );
228 t1 = z1 * ( qq1 + z1 * qq2 );
229 t2 = z2 * ( qq1 + z2 * qq2 );
230 t0 = z0 * ( poly1[0] + z0 * ( poly2[0] + t0 ) );
231 w1 = x1 * ( one + z1 * ( pp1 + z1 * pp2 ) );
232 w2 = x2 * ( one + z2 * ( pp1 + z2 * pp2 ) );
233 j1 = ( ( ( j1 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
234 j2 = ( ( ( j2 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
235 xsb1 = ( xsb1 >> 30 ) & 2;
236 xsb2 = ( xsb2 >> 30 ) & 2;
237 a1 = __vlibm_TBL_sincos_hi[j1+xsb1];
238 a2 = __vlibm_TBL_sincos_hi[j2+xsb2];
239 t0 = x0 + x0 * t0;
240 t1 = ( __vlibm_TBL_sincos_hi[j1+1] * w1 + a1 * t1 ) + __vlibm_TBL_sincos_lo[j1+xsb1];
241 t2 = ( __vlibm_TBL_sincos_hi[j2+1] * w2 + a2 * t2 ) + __vlibm_TBL_sincos_lo[j2+xsb2];
242 *py0 = t0;
243 *py1 = a1 + t1;
244 *py2 = a2 + t2;
245 break;
246
247 case 2:
248 j0 = ( xsb0 + 0x4000 ) & 0xffff8000;
249 j2 = ( xsb2 + 0x4000 ) & 0xffff8000;
250 HI(&t0) = j0;
251 HI(&t2) = j2;
252 LO(&t0) = 0;
253 LO(&t2) = 0;
254 x0 -= t0;
255 x2 -= t2;
256 z0 = x0 * x0;
257 z1 = x1 * x1;
258 z2 = x2 * x2;
259 t0 = z0 * ( qq1 + z0 * qq2 );
260 t1 = z1 * ( poly3[0] + z1 * poly4[0] );
261 t2 = z2 * ( qq1 + z2 * qq2 );
262 w0 = x0 * ( one + z0 * ( pp1 + z0 * pp2 ) );
263 t1 = z1 * ( poly1[0] + z1 * ( poly2[0] + t1 ) );
264 w2 = x2 * ( one + z2 * ( pp1 + z2 * pp2 ) );
265 j0 = ( ( ( j0 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
266 j2 = ( ( ( j2 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
267 xsb0 = ( xsb0 >> 30 ) & 2;
268 xsb2 = ( xsb2 >> 30 ) & 2;
269 a0 = __vlibm_TBL_sincos_hi[j0+xsb0];
270 a2 = __vlibm_TBL_sincos_hi[j2+xsb2];
271 t0 = ( __vlibm_TBL_sincos_hi[j0+1] * w0 + a0 * t0 ) + __vlibm_TBL_sincos_lo[j0+xsb0];
272 t1 = x1 + x1 * t1;
273 t2 = ( __vlibm_TBL_sincos_hi[j2+1] * w2 + a2 * t2 ) + __vlibm_TBL_sincos_lo[j2+xsb2];
274 *py0 = a0 + t0;
275 *py1 = t1;
276 *py2 = a2 + t2;
277 break;
278
279 case 3:
280 j2 = ( xsb2 + 0x4000 ) & 0xffff8000;
281 HI(&t2) = j2;
282 LO(&t2) = 0;
283 x2 -= t2;
284 z0 = x0 * x0;
285 z1 = x1 * x1;
286 z2 = x2 * x2;
287 t0 = z0 * ( poly3[0] + z0 * poly4[0] );
288 t1 = z1 * ( poly3[0] + z1 * poly4[0] );
289 t2 = z2 * ( qq1 + z2 * qq2 );
290 t0 = z0 * ( poly1[0] + z0 * ( poly2[0] + t0 ) );
291 t1 = z1 * ( poly1[0] + z1 * ( poly2[0] + t1 ) );
292 w2 = x2 * ( one + z2 * ( pp1 + z2 * pp2 ) );
293 j2 = ( ( ( j2 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
294 xsb2 = ( xsb2 >> 30 ) & 2;
295 a2 = __vlibm_TBL_sincos_hi[j2+xsb2];
296 t0 = x0 + x0 * t0;
297 t1 = x1 + x1 * t1;
298 t2 = ( __vlibm_TBL_sincos_hi[j2+1] * w2 + a2 * t2 ) + __vlibm_TBL_sincos_lo[j2+xsb2];
299 *py0 = t0;
300 *py1 = t1;
301 *py2 = a2 + t2;
302 break;
303
304 case 4:
305 j0 = ( xsb0 + 0x4000 ) & 0xffff8000;
306 j1 = ( xsb1 + 0x4000 ) & 0xffff8000;
307 HI(&t0) = j0;
308 HI(&t1) = j1;
309 LO(&t0) = 0;
310 LO(&t1) = 0;
311 x0 -= t0;
312 x1 -= t1;
313 z0 = x0 * x0;
314 z1 = x1 * x1;
315 z2 = x2 * x2;
316 t0 = z0 * ( qq1 + z0 * qq2 );
317 t1 = z1 * ( qq1 + z1 * qq2 );
318 t2 = z2 * ( poly3[0] + z2 * poly4[0] );
319 w0 = x0 * ( one + z0 * ( pp1 + z0 * pp2 ) );
320 w1 = x1 * ( one + z1 * ( pp1 + z1 * pp2 ) );
321 t2 = z2 * ( poly1[0] + z2 * ( poly2[0] + t2 ) );
322 j0 = ( ( ( j0 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
323 j1 = ( ( ( j1 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
324 xsb0 = ( xsb0 >> 30 ) & 2;
325 xsb1 = ( xsb1 >> 30 ) & 2;
326 a0 = __vlibm_TBL_sincos_hi[j0+xsb0];
327 a1 = __vlibm_TBL_sincos_hi[j1+xsb1];
328 t0 = ( __vlibm_TBL_sincos_hi[j0+1] * w0 + a0 * t0 ) + __vlibm_TBL_sincos_lo[j0+xsb0];
329 t1 = ( __vlibm_TBL_sincos_hi[j1+1] * w1 + a1 * t1 ) + __vlibm_TBL_sincos_lo[j1+xsb1];
330 t2 = x2 + x2 * t2;
331 *py0 = a0 + t0;
332 *py1 = a1 + t1;
333 *py2 = t2;
334 break;
335
336 case 5:
337 j1 = ( xsb1 + 0x4000 ) & 0xffff8000;
338 HI(&t1) = j1;
339 LO(&t1) = 0;
340 x1 -= t1;
341 z0 = x0 * x0;
342 z1 = x1 * x1;
343 z2 = x2 * x2;
344 t0 = z0 * ( poly3[0] + z0 * poly4[0] );
345 t1 = z1 * ( qq1 + z1 * qq2 );
346 t2 = z2 * ( poly3[0] + z2 * poly4[0] );
347 t0 = z0 * ( poly1[0] + z0 * ( poly2[0] + t0 ) );
348 w1 = x1 * ( one + z1 * ( pp1 + z1 * pp2 ) );
349 t2 = z2 * ( poly1[0] + z2 * ( poly2[0] + t2 ) );
350 j1 = ( ( ( j1 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
351 xsb1 = ( xsb1 >> 30 ) & 2;
352 a1 = __vlibm_TBL_sincos_hi[j1+xsb1];
353 t0 = x0 + x0 * t0;
354 t1 = ( __vlibm_TBL_sincos_hi[j1+1] * w1 + a1 * t1 ) + __vlibm_TBL_sincos_lo[j1+xsb1];
355 t2 = x2 + x2 * t2;
356 *py0 = t0;
357 *py1 = a1 + t1;
358 *py2 = t2;
359 break;
360
361 case 6:
362 j0 = ( xsb0 + 0x4000 ) & 0xffff8000;
363 HI(&t0) = j0;
364 LO(&t0) = 0;
365 x0 -= t0;
366 z0 = x0 * x0;
367 z1 = x1 * x1;
368 z2 = x2 * x2;
369 t0 = z0 * ( qq1 + z0 * qq2 );
370 t1 = z1 * ( poly3[0] + z1 * poly4[0] );
371 t2 = z2 * ( poly3[0] + z2 * poly4[0] );
372 w0 = x0 * ( one + z0 * ( pp1 + z0 * pp2 ) );
373 t1 = z1 * ( poly1[0] + z1 * ( poly2[0] + t1 ) );
374 t2 = z2 * ( poly1[0] + z2 * ( poly2[0] + t2 ) );
375 j0 = ( ( ( j0 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
376 xsb0 = ( xsb0 >> 30 ) & 2;
377 a0 = __vlibm_TBL_sincos_hi[j0+xsb0];
378 t0 = ( __vlibm_TBL_sincos_hi[j0+1] * w0 + a0 * t0 ) + __vlibm_TBL_sincos_lo[j0+xsb0];
379 t1 = x1 + x1 * t1;
380 t2 = x2 + x2 * t2;
381 *py0 = a0 + t0;
382 *py1 = t1;
383 *py2 = t2;
384 break;
385
386 case 7:
387 z0 = x0 * x0;
388 z1 = x1 * x1;
389 z2 = x2 * x2;
390 t0 = z0 * ( poly3[0] + z0 * poly4[0] );
391 t1 = z1 * ( poly3[0] + z1 * poly4[0] );
392 t2 = z2 * ( poly3[0] + z2 * poly4[0] );
393 t0 = z0 * ( poly1[0] + z0 * ( poly2[0] + t0 ) );
394 t1 = z1 * ( poly1[0] + z1 * ( poly2[0] + t1 ) );
395 t2 = z2 * ( poly1[0] + z2 * ( poly2[0] + t2 ) );
396 t0 = x0 + x0 * t0;
397 t1 = x1 + x1 * t1;
398 t2 = x2 + x2 * t2;
399 *py0 = t0;
400 *py1 = t1;
401 *py2 = t2;
402 break;
403 }
404
405 x += stridex;
406 y += stridey;
407 i = 0;
408 } while ( --n > 0 );
409
410 if ( i > 0 )
411 {
412 double a0, a1, w0, w1;
413 double t0, t1, z0, z1;
414 unsigned j0, j1;
415
416 if ( i > 1 )
417 {
418 if ( hx1 < 0x3fc90000 )
419 {
420 z1 = x1 * x1;
421 t1 = z1 * ( poly3[0] + z1 * poly4[0] );
422 t1 = z1 * ( poly1[0] + z1 * ( poly2[0] + t1 ) );
423 t1 = x1 + x1 * t1;
424 *py1 = t1;
425 }
426 else
427 {
428 j1 = ( xsb1 + 0x4000 ) & 0xffff8000;
429 HI(&t1) = j1;
430 LO(&t1) = 0;
431 x1 -= t1;
432 z1 = x1 * x1;
433 t1 = z1 * ( qq1 + z1 * qq2 );
434 w1 = x1 * ( one + z1 * ( pp1 + z1 * pp2 ) );
435 j1 = ( ( ( j1 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
436 xsb1 = ( xsb1 >> 30 ) & 2;
437 a1 = __vlibm_TBL_sincos_hi[j1+xsb1];
438 t1 = ( __vlibm_TBL_sincos_hi[j1+1] * w1 + a1 * t1 ) + __vlibm_TBL_sincos_lo[j1+xsb1];
439 *py1 = a1 + t1;
440 }
441 }
442 if ( hx0 < 0x3fc90000 )
443 {
444 z0 = x0 * x0;
445 t0 = z0 * ( poly3[0] + z0 * poly4[0] );
446 t0 = z0 * ( poly1[0] + z0 * ( poly2[0] + t0 ) );
447 t0 = x0 + x0 * t0;
448 *py0 = t0;
449 }
450 else
451 {
452 j0 = ( xsb0 + 0x4000 ) & 0xffff8000;
453 HI(&t0) = j0;
454 LO(&t0) = 0;
455 x0 -= t0;
456 z0 = x0 * x0;
457 t0 = z0 * ( qq1 + z0 * qq2 );
458 w0 = x0 * ( one + z0 * ( pp1 + z0 * pp2 ) );
459 j0 = ( ( ( j0 & ~0x80000000 ) - 0x3fc40000 ) >> 13 ) & ~0x3;
460 xsb0 = ( xsb0 >> 30 ) & 2;
461 a0 = __vlibm_TBL_sincos_hi[j0+xsb0];
462 t0 = ( __vlibm_TBL_sincos_hi[j0+1] * w0 + a0 * t0 ) + __vlibm_TBL_sincos_lo[j0+xsb0];
463 *py0 = a0 + t0;
464 }
465 }
466
467 return;
468
469 /*
470 * MEDIUM RANGE PROCESSING
471 * Jump here at first sign of medium range argument. We are a bit
472 * confused due to the jump.. fix up several variables and jump into
473 * the nth loop, same as was being processed above.
474 */
475
476 MEDIUM:
477
478 x0_or_one[1] = 1.0;
479 x1_or_one[1] = 1.0;
480 x2_or_one[1] = 1.0;
481 x0_or_one[3] = -1.0;
482 x1_or_one[3] = -1.0;
483 x2_or_one[3] = -1.0;
484 y0_or_zero[1] = 0.0;
485 y1_or_zero[1] = 0.0;
486 y2_or_zero[1] = 0.0;
487 y0_or_zero[3] = 0.0;
488 y1_or_zero[3] = 0.0;
489 y2_or_zero[3] = 0.0;
490
491 if ( biguns == 3 )
492 {
493 biguns = 0;
494 xsb0 = xsb0 >> 31;
495 xsb1 = xsb1 >> 31;
496 goto loop2;
497 }
498 else if ( biguns == 2 )
499 {
500 xsb0 = xsb0 >> 31;
501 biguns = 0;
502 goto loop1;
503 }
504 biguns = 0;
505
506 do
507 {
508 double fn0, fn1, fn2, a0, a1, a2, w0, w1, w2, y0, y1, y2;
509 unsigned hx;
510 int n0, n1, n2;
511
512 loop0:
513 hx = HI(x);
514 xsb0 = hx >> 31;
515 hx &= ~0x80000000;
516 if ( hx < 0x3e400000 )
517 {
518 *y = *x;
519 x += stridex;
520 y += stridey;
521 i = 0;
522 if ( --n <= 0 )
523 break;
524 goto loop0;
525 }
526 if ( hx > 0x413921fb )
527 {
528 if ( hx >= 0x7ff00000 )
529 {
530 x0 = *x;
531 *y = x0 - x0;
532 }
533 else
534 biguns = 1;
535 x += stridex;
536 y += stridey;
537 i = 0;
538 if ( --n <= 0 )
539 break;
540 goto loop0;
541 }
542 x0 = *x;
543 py0 = y;
544 x += stridex;
545 y += stridey;
546 i = 1;
547 if ( --n <= 0 )
548 break;
549
550 loop1:
551 hx = HI(x);
552 xsb1 = hx >> 31;
553 hx &= ~0x80000000;
554 if ( hx < 0x3e400000 )
555 {
556 *y = *x;
557 x += stridex;
558 y += stridey;
559 i = 1;
560 if ( --n <= 0 )
561 break;
562 goto loop1;
563 }
564 if ( hx > 0x413921fb )
565 {
566 if ( hx >= 0x7ff00000 )
567 {
568 x1 = *x;
569 *y = x1 - x1;
570 }
571 else
572 biguns = 1;
573 x += stridex;
574 y += stridey;
575 i = 1;
576 if ( --n <= 0 )
577 break;
578 goto loop1;
579 }
580 x1 = *x;
581 py1 = y;
582 x += stridex;
583 y += stridey;
584 i = 2;
585 if ( --n <= 0 )
586 break;
587
588 loop2:
589 hx = HI(x);
590 xsb2 = hx >> 31;
591 hx &= ~0x80000000;
592 if ( hx < 0x3e400000 )
593 {
594 *y = *x;
595 x += stridex;
596 y += stridey;
597 i = 2;
598 if ( --n <= 0 )
599 break;
600 goto loop2;
601 }
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 &= 3;
628 n1 &= 3;
629 n2 &= 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 &= 3;
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 &= 3;
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_vsin_big( nsave, xsave, sxsave, ysave, sysave, 0x413921fb );
1101 }