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