summaryrefslogtreecommitdiffstats
path: root/baseline/source/fmref/wcclibm.c
diff options
context:
space:
mode:
Diffstat (limited to 'baseline/source/fmref/wcclibm.c')
-rw-r--r--baseline/source/fmref/wcclibm.c522
1 files changed, 522 insertions, 0 deletions
diff --git a/baseline/source/fmref/wcclibm.c b/baseline/source/fmref/wcclibm.c
new file mode 100644
index 0000000..39b8cbe
--- /dev/null
+++ b/baseline/source/fmref/wcclibm.c
@@ -0,0 +1,522 @@
1#include "math_private.h"
2#include "wcclibm.h"
3
4
5
6/* e_rem_pio2f.c -- float version of e_rem_pio2.c
7 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
8 */
9
10/*
11 * ====================================================
12 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
13 *
14 * Developed at SunPro, a Sun Microsystems, Inc. business.
15 * Permission to use, copy, modify, and distribute this
16 * software is freely granted, provided that this notice
17 * is preserved.
18 * ====================================================
19 */
20
21#if defined(LIBM_SCCS) && !defined(lint)
22static char rcsid[] = "$NetBSD: e_rem_pio2f.c,v 1.5 1995/05/10 20:46:03 jtc Exp $";
23#endif
24
25/* __ieee754_rem_pio2f(x,y)
26 *
27 * return the remainder of x rem pi/2 in y[0]+y[1]
28 * use __kernel_rem_pio2f()
29 */
30
31/* This array is like the one in e_rem_pio2.c, but the numbers are
32 single precision and the last 8 bits are forced to 0. */
33#ifdef __STDC__
34static const int32_t fmref_npio2_hw[] = {
35#else
36static int32_t fmref_npio2_hw[] = {
37#endif
380x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
390x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
400x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
410x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
420x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
430x4242c700, 0x42490f00
44};
45
46/*
47 * invpio2: 24 bits of 2/pi
48 * pio2_1: first 17 bit of pi/2
49 * pio2_1t: pi/2 - pio2_1
50 * pio2_2: second 17 bit of pi/2
51 * pio2_2t: pi/2 - (pio2_1+pio2_2)
52 * pio2_3: third 17 bit of pi/2
53 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
54 */
55
56#ifdef __STDC__
57static const float
58#else
59static float
60#endif
61/* zero = 0.0000000000e+00f, /\* 0x00000000 *\/ */
62/* half = 5.0000000000e-01f, /\* 0x3f000000 *\/ */
63/* two8 = 2.5600000000e+02f, /\* 0x43800000 *\/ */
64fmref_invpio2 = 6.3661980629e-01f, /* 0x3f22f984 */
65fmref_pio2_1 = 1.5707855225e+00f, /* 0x3fc90f80 */
66fmref_pio2_1t = 1.0804334124e-05f, /* 0x37354443 */
67fmref_pio2_2 = 1.0804273188e-05f, /* 0x37354400 */
68fmref_pio2_2t = 6.0770999344e-11f, /* 0x2e85a308 */
69fmref_pio2_3 = 6.0770943833e-11f, /* 0x2e85a300 */
70fmref_pio2_3t = 6.1232342629e-17f; /* 0x248d3132 */
71
72#ifdef __STDC__
73 int32_t fmref___ieee754_rem_pio2f(float x, float *y)
74#else
75 int32_t fmref___ieee754_rem_pio2f(x,y)
76 float x,y[];
77#endif
78{
79 float z,w,t,r,fn;
80 int32_t i,j,n,ix,hx;
81
82 GET_FLOAT_WORD(hx,x);
83 ix = hx&0x7fffffff;
84 if(ix<=0x3f490fd8) /* |x| ~<= pi/4 , no need for reduction */
85 {y[0] = x; y[1] = 0; return 0;}
86 if(ix<0x4016cbe4) { /* |x| < 3pi/4, special case with n=+-1 */
87 if(hx>0) {
88 z = x - fmref_pio2_1;
89 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
90 y[0] = z - fmref_pio2_1t;
91 y[1] = (z-y[0])-fmref_pio2_1t;
92 } else { /* near pi/2, use 24+24+24 bit pi */
93 z -= fmref_pio2_2;
94 y[0] = z - fmref_pio2_2t;
95 y[1] = (z-y[0])-fmref_pio2_2t;
96 }
97 return 1;
98 } else { /* negative x */
99 z = x + fmref_pio2_1;
100 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
101 y[0] = z + fmref_pio2_1t;
102 y[1] = (z-y[0])+fmref_pio2_1t;
103 } else { /* near pi/2, use 24+24+24 bit pi */
104 z += fmref_pio2_2;
105 y[0] = z + fmref_pio2_2t;
106 y[1] = (z-y[0])+fmref_pio2_2t;
107 }
108 return -1;
109 }
110 }
111 if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
112 t = fabsf(x);
113 n = (int32_t) (t*fmref_invpio2+fmref_half);
114 fn = (float)n;
115 r = t-fn*fmref_pio2_1;
116 w = fn*fmref_pio2_1t; /* 1st round good to 40 bit */
117 if(n<32&&(int32_t)(ix&0xffffff00)!=fmref_npio2_hw[n-1]) {
118 y[0] = r-w; /* quick check no cancellation */
119 } else {
120 u_int32_t high;
121 j = ix>>23;
122 y[0] = r-w;
123 GET_FLOAT_WORD(high,y[0]);
124 i = j-((high>>23)&0xff);
125 if(i>8) { /* 2nd iteration needed, good to 57 */
126 t = r;
127 w = fn*fmref_pio2_2;
128 r = t-w;
129 w = fn*fmref_pio2_2t-((t-r)-w);
130 y[0] = r-w;
131 GET_FLOAT_WORD(high,y[0]);
132 i = j-((high>>23)&0xff);
133 if(i>25) { /* 3rd iteration need, 74 bits acc */
134 t = r; /* will cover all possible cases */
135 w = fn*fmref_pio2_3;
136 r = t-w;
137 w = fn*fmref_pio2_3t-((t-r)-w);
138 y[0] = r-w;
139 }
140 }
141 }
142 y[1] = (r-y[0])-w;
143 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
144 else return n;
145 }
146 /*
147 * all other (large) arguments
148 */
149 if(ix>=0x7f800000) { /* x is inf or NaN */
150 y[0]=y[1]=x-x; return 0;
151 }
152
153 y[0]=y[1]=x-x; /* dummy initialization */
154 return 0; /* doesn't happen for our input */
155}
156
157/* k_cosf.c -- float version of k_cos.c
158 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
159 */
160
161/*
162 * ====================================================
163 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
164 *
165 * Developed at SunPro, a Sun Microsystems, Inc. business.
166 * Permission to use, copy, modify, and distribute this
167 * software is freely granted, provided that this notice
168 * is preserved.
169 * ====================================================
170 */
171
172#if defined(LIBM_SCCS) && !defined(lint)
173static char rcsid[] = "$NetBSD: k_cosf.c,v 1.4 1995/05/10 20:46:23 jtc Exp $";
174#endif
175
176
177#ifdef __STDC__
178static const float
179#else
180static float
181#endif
182/* one = 1.0000000000e+00, /\* 0x3f800000 *\/ */
183fmref_C1 = 4.1666667908e-02f, /* 0x3d2aaaab */
184fmref_C2 = -1.3888889225e-03f, /* 0xbab60b61 */
185fmref_C3 = 2.4801587642e-05f, /* 0x37d00d01 */
186fmref_C4 = -2.7557314297e-07f, /* 0xb493f27c */
187fmref_C5 = 2.0875723372e-09f, /* 0x310f74f6 */
188fmref_C6 = -1.1359647598e-11f; /* 0xad47d74e */
189
190#ifdef __STDC__
191 float fmref___kernel_cosf(float x, float y)
192#else
193 float fmref___kernel_cosf(x, y)
194 float x,y;
195#endif
196{
197 float a,hz,z,r,qx;
198 int32_t ix;
199 GET_FLOAT_WORD(ix,x);
200 ix &= 0x7fffffff; /* ix = |x|'s high word*/
201 if(ix<0x32000000) { /* if x < 2**27 */
202 if(((int)x)==0) return fmref_one; /* generate inexact */
203 }
204 z = x*x;
205 r = z*(fmref_C1+z*(fmref_C2+z*(fmref_C3+z*(fmref_C4+z*(fmref_C5+z*fmref_C6)))));
206 if(ix < 0x3e99999a) /* if |x| < 0.3 */
207 return fmref_one - ((float)0.5f*z - (z*r - x*y));
208 else {
209 if(ix > 0x3f480000) { /* x > 0.78125 */
210 qx = (float)0.28125f;
211 } else {
212 SET_FLOAT_WORD(qx,ix-0x01000000); /* x/4 */
213 }
214 hz = (float)0.5f*z-qx;
215 a = fmref_one-qx;
216 return a - (hz - (z*r-x*y));
217 }
218}
219
220/* k_sinf.c -- float version of k_sin.c
221 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
222 */
223
224/*
225 * ====================================================
226 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
227 *
228 * Developed at SunPro, a Sun Microsystems, Inc. business.
229 * Permission to use, copy, modify, and distribute this
230 * software is freely granted, provided that this notice
231 * is preserved.
232 * ====================================================
233 */
234
235#if defined(LIBM_SCCS) && !defined(lint)
236static char rcsid[] = "$NetBSD: k_sinf.c,v 1.4 1995/05/10 20:46:33 jtc Exp $";
237#endif
238
239
240#ifdef __STDC__
241static const float
242#else
243static float
244#endif
245/* half = 5.0000000000e-01f,/\* 0x3f000000 *\/ */
246fmref_S1 = -1.6666667163e-01f, /* 0xbe2aaaab */
247fmref_S2 = 8.3333337680e-03f, /* 0x3c088889 */
248fmref_S3 = -1.9841270114e-04f, /* 0xb9500d01 */
249fmref_S4 = 2.7557314297e-06f, /* 0x3638ef1b */
250fmref_S5 = -2.5050759689e-08f, /* 0xb2d72f34 */
251fmref_S6 = 1.5896910177e-10f; /* 0x2f2ec9d3 */
252
253#ifdef __STDC__
254 float fmref___kernel_sinf(float x, float y, int iy)
255#else
256 float fmref___kernel_sinf(x, y, iy)
257 float x,y; int iy; /* iy=0 if y is zero */
258#endif
259{
260 float z,r,v;
261 int32_t ix;
262 GET_FLOAT_WORD(ix,x);
263 ix &= 0x7fffffff; /* high word of x */
264 if(ix<0x32000000) /* |x| < 2**-27 */
265 {if((int)x==0) return x;} /* generate inexact */
266 z = x*x;
267 v = z*x;
268 r = fmref_S2+z*(fmref_S3+z*(fmref_S4+z*(fmref_S5+z*fmref_S6)));
269 if(iy==0) return x+v*(fmref_S1+z*r);
270 else return x-((z*(fmref_half*y-v*r)-y)-v*fmref_S1);
271}
272/* s_atanf.c -- float version of s_atan.c.
273 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
274 */
275
276/*
277 * ====================================================
278 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
279 *
280 * Developed at SunPro, a Sun Microsystems, Inc. business.
281 * Permission to use, copy, modify, and distribute this
282 * software is freely granted, provided that this notice
283 * is preserved.
284 * ====================================================
285 */
286
287#if defined(LIBM_SCCS) && !defined(lint)
288static char rcsid[] = "$NetBSD: s_atanf.c,v 1.4 1995/05/10 20:46:47 jtc Exp $";
289#endif
290
291
292#ifdef __STDC__
293static const float fmref_atanhi[] = {
294#else
295static float fmref_atanhi[] = {
296#endif
297 4.6364760399e-01f, /* atan(0.5)hi 0x3eed6338 */
298 7.8539812565e-01f, /* atan(1.0)hi 0x3f490fda */
299 9.8279368877e-01f, /* atan(1.5)hi 0x3f7b985e */
300 1.5707962513e+00f, /* atan(inf)hi 0x3fc90fda */
301};
302
303#ifdef __STDC__
304static const float fmref_atanlo[] = {
305#else
306static float fmref_atanlo[] = {
307#endif
308 5.0121582440e-09f, /* atan(0.5)lo 0x31ac3769 */
309 3.7748947079e-08f, /* atan(1.0)lo 0x33222168 */
310 3.4473217170e-08f, /* atan(1.5)lo 0x33140fb4 */
311 7.5497894159e-08f, /* atan(inf)lo 0x33a22168 */
312};
313
314#ifdef __STDC__
315static const float fmref_aT[] = {
316#else
317static float fmref_aT[] = {
318#endif
319 3.3333334327e-01f, /* 0x3eaaaaaa */
320 -2.0000000298e-01f, /* 0xbe4ccccd */
321 1.4285714924e-01f, /* 0x3e124925 */
322 -1.1111110449e-01f, /* 0xbde38e38 */
323 9.0908870101e-02f, /* 0x3dba2e6e */
324 -7.6918758452e-02f, /* 0xbd9d8795 */
325 6.6610731184e-02f, /* 0x3d886b35 */
326 -5.8335702866e-02f, /* 0xbd6ef16b */
327 4.9768779427e-02f, /* 0x3d4bda59 */
328 -3.6531571299e-02f, /* 0xbd15a221 */
329 1.6285819933e-02f, /* 0x3c8569d7 */
330};
331
332/* #ifdef __STDC__ */
333/* static const float */
334/* #else */
335/* static float */
336/* #endif */
337/* one = 1.0, */
338/* huge = 1.0e30; */
339
340#ifdef __STDC__
341 float fmref___atanf(float x)
342#else
343 float fmref___atanf(x)
344 float x;
345#endif
346{
347 float w,s1,s2,z;
348 int32_t ix,hx,id;
349
350 GET_FLOAT_WORD(hx,x);
351 ix = hx&0x7fffffff;
352 if(ix>=0x50800000) { /* if |x| >= 2^34 */
353 if(ix>0x7f800000)
354 return x+x; /* NaN */
355 if(hx>0) return fmref_atanhi[3]+fmref_atanlo[3];
356 else return -fmref_atanhi[3]-fmref_atanlo[3];
357 } if (ix < 0x3ee00000) { /* |x| < 0.4375 */
358 if (ix < 0x31000000) { /* |x| < 2^-29 */
359 if(fmref_huge+x>fmref_one) return x; /* raise inexact */
360 }
361 id = -1;
362 } else {
363 x = fabsf(x);
364 if (ix < 0x3f980000) { /* |x| < 1.1875 */
365 if (ix < 0x3f300000) { /* 7/16 <=|x|<11/16 */
366 id = 0; x = ((float)2.0f*x-fmref_one)/((float)2.0f+x);
367 } else { /* 11/16<=|x|< 19/16 */
368 id = 1; x = (x-fmref_one)/(x+fmref_one);
369 }
370 } else {
371 if (ix < 0x401c0000) { /* |x| < 2.4375 */
372 id = 2; x = (x-(float)1.5f)/(fmref_one+(float)1.5f*x);
373 } else { /* 2.4375 <= |x| < 2^66 */
374 id = 3; x = -(float)1.0f/x;
375 }
376 }}
377 /* end of argument reduction */
378 z = x*x;
379 w = z*z;
380 /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
381 s1 = z*(fmref_aT[0]+w*(fmref_aT[2]+w*(fmref_aT[4]+w*(fmref_aT[6]+w*(fmref_aT[8]+w*fmref_aT[10])))));
382 s2 = w*(fmref_aT[1]+w*(fmref_aT[3]+w*(fmref_aT[5]+w*(fmref_aT[7]+w*fmref_aT[9]))));
383 if (id<0) return x - x*(s1+s2);
384 else {
385 z = fmref_atanhi[id] - ((x*(s1+s2) - fmref_atanlo[id]) - x);
386 return (hx<0)? -z:z;
387 }
388}
389//weak_alias (__atanf, atanf)
390
391/* s_cosf.c -- float version of s_cos.c.
392 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
393 */
394
395/*
396 * ====================================================
397 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
398 *
399 * Developed at SunPro, a Sun Microsystems, Inc. business.
400 * Permission to use, copy, modify, and distribute this
401 * software is freely granted, provided that this notice
402 * is preserved.
403 * ====================================================
404 */
405
406/* #ifdef __STDC__ */
407/* static const float one=1.0; */
408/* #else */
409/* static float one=1.0; */
410/* #endif */
411
412#ifdef __STDC__
413 float fmref___cosf(float x)
414#else
415 float fmref___cosf(x)
416 float x;
417#endif
418{
419 float y[2],z=0.0f;
420 int32_t n,ix;
421
422 GET_FLOAT_WORD(ix,x);
423
424 /* |x| ~< pi/4 */
425 ix &= 0x7fffffff;
426 if(ix <= 0x3f490fd8) return fmref___kernel_cosf(x,z);
427
428 /* cos(Inf or NaN) is NaN */
429 else if (ix>=0x7f800000) return x-x;
430
431 /* argument reduction needed */
432 else {
433 n = fmref___ieee754_rem_pio2f(x,y);
434 switch(n&3) {
435 case 0: return fmref___kernel_cosf(y[0],y[1]);
436 case 1: return -fmref___kernel_sinf(y[0],y[1],1);
437 case 2: return -fmref___kernel_cosf(y[0],y[1]);
438 default:
439 return fmref___kernel_sinf(y[0],y[1],1);
440 }
441 }
442}
443
444/* s_sinf.c -- float version of s_sin.c.
445 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
446 */
447
448/*
449 * ====================================================
450 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
451 *
452 * Developed at SunPro, a Sun Microsystems, Inc. business.
453 * Permission to use, copy, modify, and distribute this
454 * software is freely granted, provided that this notice
455 * is preserved.
456 * ====================================================
457 */
458
459 #ifdef __STDC__
460 float fmref___sinf(float x)
461#else
462 float fmref___sinf(x)
463 float x;
464#endif
465{
466 float y[2],z=0.0;
467 int32_t n, ix;
468
469 GET_FLOAT_WORD(ix,x);
470
471 /* |x| ~< pi/4 */
472 ix &= 0x7fffffff;
473 if(ix <= 0x3f490fd8) return fmref___kernel_sinf(x,z,0);
474
475 /* sin(Inf or NaN) is NaN */
476 else if (ix>=0x7f800000) return x-x;
477
478 /* argument reduction needed */
479 else {
480 n = fmref___ieee754_rem_pio2f(x,y);
481 switch(n&3) {
482 case 0: return fmref___kernel_sinf(y[0],y[1],1);
483 case 1: return fmref___kernel_cosf(y[0],y[1]);
484 case 2: return -fmref___kernel_sinf(y[0],y[1],1);
485 default:
486 return -fmref___kernel_cosf(y[0],y[1]);
487 }
488 }
489}
490
491/* s_fabsf.c -- float version of s_fabs.c.
492 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
493 */
494
495/*
496 * ====================================================
497 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
498 *
499 * Developed at SunPro, a Sun Microsystems, Inc. business.
500 * Permission to use, copy, modify, and distribute this
501 * software is freely granted, provided that this notice
502 * is preserved.
503 * ====================================================
504 */
505
506/*
507 * fabsf(x) returns the absolute value of x.
508 */
509
510
511#ifdef __STDC__
512 float fmref___fabsf(float x)
513#else
514 float fmref___fabsf(x)
515 float x;
516#endif
517{
518 u_int32_t ix;
519 GET_FLOAT_WORD(ix,x);
520 SET_FLOAT_WORD(x,ix&0x7fffffff);
521 return x;
522}
generated by cgit v1.2.2 (git 2.25.0) at 2025-05-02 18:52:33 -0400