summaryrefslogtreecommitdiffstats
path: root/all_pairs/source/susan/wcclibm.c
diff options
context:
space:
mode:
Diffstat (limited to 'all_pairs/source/susan/wcclibm.c')
-rw-r--r--all_pairs/source/susan/wcclibm.c444
1 files changed, 444 insertions, 0 deletions
diff --git a/all_pairs/source/susan/wcclibm.c b/all_pairs/source/susan/wcclibm.c
new file mode 100644
index 0000000..bca8907
--- /dev/null
+++ b/all_pairs/source/susan/wcclibm.c
@@ -0,0 +1,444 @@
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. */
33static const int32_t susan_npio2_hw[] = {
340x3fc90f00, 0x40490f00, 0x4096cb00, 0x40c90f00, 0x40fb5300, 0x4116cb00,
350x412fed00, 0x41490f00, 0x41623100, 0x417b5300, 0x418a3a00, 0x4196cb00,
360x41a35c00, 0x41afed00, 0x41bc7e00, 0x41c90f00, 0x41d5a000, 0x41e23100,
370x41eec200, 0x41fb5300, 0x4203f200, 0x420a3a00, 0x42108300, 0x4216cb00,
380x421d1400, 0x42235c00, 0x4229a500, 0x422fed00, 0x42363600, 0x423c7e00,
390x4242c700, 0x42490f00
40};
41
42/*
43 * invpio2: 24 bits of 2/pi
44 * pio2_1: first 17 bit of pi/2
45 * pio2_1t: pi/2 - pio2_1
46 * pio2_2: second 17 bit of pi/2
47 * pio2_2t: pi/2 - (pio2_1+pio2_2)
48 * pio2_3: third 17 bit of pi/2
49 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
50 */
51
52static const float
53/* zero = 0.0000000000e+00f, /\* 0x00000000 *\/ */
54/* half = 5.0000000000e-01f, /\* 0x3f000000 *\/ */
55/* two8 = 2.5600000000e+02f, /\* 0x43800000 *\/ */
56susan_invpio2 = 6.3661980629e-01f, /* 0x3f22f984 */
57susan_pio2_1 = 1.5707855225e+00f, /* 0x3fc90f80 */
58susan_pio2_1t = 1.0804334124e-05f, /* 0x37354443 */
59susan_pio2_2 = 1.0804273188e-05f, /* 0x37354400 */
60susan_pio2_2t = 6.0770999344e-11f, /* 0x2e85a308 */
61susan_pio2_3 = 6.0770943833e-11f, /* 0x2e85a300 */
62susan_pio2_3t = 6.1232342629e-17f; /* 0x248d3132 */
63
64int32_t susan___ieee754_rem_pio2f(float x, float *y)
65{
66 float z,w,t,r,fn;
67 int32_t i,j,n,ix,hx;
68
69 GET_FLOAT_WORD(hx,x);
70 ix = hx&0x7fffffff;
71 if(ix<=0x3f490fd8) /* |x| ~<= pi/4 , no need for reduction */
72 {y[0] = x; y[1] = 0; return 0;}
73 if(ix<0x4016cbe4) { /* |x| < 3pi/4, special case with n=+-1 */
74 if(hx>0) {
75 z = x - susan_pio2_1;
76 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
77 y[0] = z - susan_pio2_1t;
78 y[1] = (z-y[0])-susan_pio2_1t;
79 } else { /* near pi/2, use 24+24+24 bit pi */
80 z -= susan_pio2_2;
81 y[0] = z - susan_pio2_2t;
82 y[1] = (z-y[0])-susan_pio2_2t;
83 }
84 return 1;
85 } else { /* negative x */
86 z = x + susan_pio2_1;
87 if((ix&0xfffffff0)!=0x3fc90fd0) { /* 24+24 bit pi OK */
88 y[0] = z + susan_pio2_1t;
89 y[1] = (z-y[0])+susan_pio2_1t;
90 } else { /* near pi/2, use 24+24+24 bit pi */
91 z += susan_pio2_2;
92 y[0] = z + susan_pio2_2t;
93 y[1] = (z-y[0])+susan_pio2_2t;
94 }
95 return -1;
96 }
97 }
98 if(ix<=0x43490f80) { /* |x| ~<= 2^7*(pi/2), medium size */
99 t = fabsf(x);
100 n = (int32_t) (t*susan_invpio2+susan_half);
101 fn = (float)n;
102 r = t-fn*susan_pio2_1;
103 w = fn*susan_pio2_1t; /* 1st round good to 40 bit */
104 if(n<32&&(int32_t)(ix&0xffffff00)!=susan_npio2_hw[n-1]) {
105 y[0] = r-w; /* quick check no cancellation */
106 } else {
107 u_int32_t high;
108 j = ix>>23;
109 y[0] = r-w;
110 GET_FLOAT_WORD(high,y[0]);
111 i = j-((high>>23)&0xff);
112 if(i>8) { /* 2nd iteration needed, good to 57 */
113 t = r;
114 w = fn*susan_pio2_2;
115 r = t-w;
116 w = fn*susan_pio2_2t-((t-r)-w);
117 y[0] = r-w;
118 GET_FLOAT_WORD(high,y[0]);
119 i = j-((high>>23)&0xff);
120 if(i>25) { /* 3rd iteration need, 74 bits acc */
121 t = r; /* will cover all possible cases */
122 w = fn*susan_pio2_3;
123 r = t-w;
124 w = fn*susan_pio2_3t-((t-r)-w);
125 y[0] = r-w;
126 }
127 }
128 }
129 y[1] = (r-y[0])-w;
130 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
131 else return n;
132 }
133 /*
134 * all other (large) arguments
135 */
136 if(ix>=0x7f800000) { /* x is inf or NaN */
137 y[0]=y[1]=x-x; return 0;
138 }
139
140 y[0]=y[1]=x-x; return 0; /* dummy initialisation */
141 return 0; /* doesn't happen for our input */
142}
143
144/* k_cosf.c -- float version of k_cos.c
145 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
146 */
147
148/*
149 * ====================================================
150 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
151 *
152 * Developed at SunPro, a Sun Microsystems, Inc. business.
153 * Permission to use, copy, modify, and distribute this
154 * software is freely granted, provided that this notice
155 * is preserved.
156 * ====================================================
157 */
158
159#if defined(LIBM_SCCS) && !defined(lint)
160static char rcsid[] = "$NetBSD: k_cosf.c,v 1.4 1995/05/10 20:46:23 jtc Exp $";
161#endif
162
163
164static const float
165/* one = 1.0000000000e+00, /\* 0x3f800000 *\/ */
166susan_C1 = 4.1666667908e-02f, /* 0x3d2aaaab */
167susan_C2 = -1.3888889225e-03f, /* 0xbab60b61 */
168susan_C3 = 2.4801587642e-05f, /* 0x37d00d01 */
169susan_C4 = -2.7557314297e-07f, /* 0xb493f27c */
170susan_C5 = 2.0875723372e-09f, /* 0x310f74f6 */
171susan_C6 = -1.1359647598e-11f; /* 0xad47d74e */
172
173float susan___kernel_cosf(float x, float y)
174{
175 float a,hz,z,r,qx;
176 int32_t ix;
177 GET_FLOAT_WORD(ix,x);
178 ix &= 0x7fffffff; /* ix = |x|'s high word*/
179 if(ix<0x32000000) { /* if x < 2**27 */
180 if(((int)x)==0) return susan_one; /* generate inexact */
181 }
182 z = x*x;
183 r = z*(susan_C1+z*(susan_C2+z*(susan_C3+z*(susan_C4+z*(susan_C5+z*susan_C6)))));
184 if(ix < 0x3e99999a) /* if |x| < 0.3 */
185 return susan_one - ((float)0.5f*z - (z*r - x*y));
186 else {
187 if(ix > 0x3f480000) { /* x > 0.78125 */
188 qx = (float)0.28125f;
189 } else {
190 SET_FLOAT_WORD(qx,ix-0x01000000); /* x/4 */
191 }
192 hz = (float)0.5f*z-qx;
193 a = susan_one-qx;
194 return a - (hz - (z*r-x*y));
195 }
196}
197
198/* k_sinf.c -- float version of k_sin.c
199 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.