OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_j1f.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_j1f.c */ |
2 /* | 2 /* |
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. | 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
4 */ | 4 */ |
5 /* | 5 /* |
6 * ==================================================== | 6 * ==================================================== |
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
8 * | 8 * |
9 * Developed at SunPro, a Sun Microsystems, Inc. business. | 9 * Developed at SunPro, a Sun Microsystems, Inc. business. |
10 * Permission to use, copy, modify, and distribute this | 10 * Permission to use, copy, modify, and distribute this |
11 * software is freely granted, provided that this notice | 11 * software is freely granted, provided that this notice |
12 * is preserved. | 12 * is preserved. |
13 * ==================================================== | 13 * ==================================================== |
14 */ | 14 */ |
15 | 15 |
16 #define _GNU_SOURCE | 16 #define _GNU_SOURCE |
17 #include "libm.h" | 17 #include "libm.h" |
18 | 18 |
19 static float ponef(float), qonef(float); | 19 static float ponef(float), qonef(float); |
20 | 20 |
21 static const float | 21 static const float invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */ |
22 invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */ | 22 tpi = 6.3661974669e-01; /* 0x3f22f983 */ |
23 tpi = 6.3661974669e-01; /* 0x3f22f983 */ | |
24 | 23 |
25 static float common(uint32_t ix, float x, int y1, int sign) | 24 static float common(uint32_t ix, float x, int y1, int sign) { |
26 { | 25 double z, s, c, ss, cc; |
27 » double z,s,c,ss,cc; | |
28 | 26 |
29 » s = sinf(x); | 27 s = sinf(x); |
30 » if (y1) | 28 if (y1) |
31 » » s = -s; | 29 s = -s; |
32 » c = cosf(x); | 30 c = cosf(x); |
33 » cc = s-c; | 31 cc = s - c; |
34 » if (ix < 0x7f000000) { | 32 if (ix < 0x7f000000) { |
35 » » ss = -s-c; | 33 ss = -s - c; |
36 » » z = cosf(2*x); | 34 z = cosf(2 * x); |
37 » » if (s*c > 0) | 35 if (s * c > 0) |
38 » » » cc = z/ss; | 36 cc = z / ss; |
39 » » else | 37 else |
40 » » » ss = z/cc; | 38 ss = z / cc; |
41 » » if (ix < 0x58800000) { | 39 if (ix < 0x58800000) { |
42 » » » if (y1) | 40 if (y1) |
43 » » » » ss = -ss; | 41 ss = -ss; |
44 » » » cc = ponef(x)*cc-qonef(x)*ss; | 42 cc = ponef(x) * cc - qonef(x) * ss; |
45 » » } | 43 } |
46 » } | 44 } |
47 » if (sign) | 45 if (sign) |
48 » » cc = -cc; | 46 cc = -cc; |
49 » return invsqrtpi*cc/sqrtf(x); | 47 return invsqrtpi * cc / sqrtf(x); |
50 } | 48 } |
51 | 49 |
52 /* R0/S0 on [0,2] */ | 50 /* R0/S0 on [0,2] */ |
53 static const float | 51 static const float r00 = -6.2500000000e-02, /* 0xbd800000 */ |
54 r00 = -6.2500000000e-02, /* 0xbd800000 */ | 52 r01 = 1.4070566976e-03, /* 0x3ab86cfd */ |
55 r01 = 1.4070566976e-03, /* 0x3ab86cfd */ | 53 r02 = -1.5995563444e-05, /* 0xb7862e36 */ |
56 r02 = -1.5995563444e-05, /* 0xb7862e36 */ | 54 r03 = 4.9672799207e-08, /* 0x335557d2 */ |
57 r03 = 4.9672799207e-08, /* 0x335557d2 */ | 55 s01 = 1.9153760746e-02, /* 0x3c9ce859 */ |
58 s01 = 1.9153760746e-02, /* 0x3c9ce859 */ | 56 s02 = 1.8594678841e-04, /* 0x3942fab6 */ |
59 s02 = 1.8594678841e-04, /* 0x3942fab6 */ | 57 s03 = 1.1771846857e-06, /* 0x359dffc2 */ |
60 s03 = 1.1771846857e-06, /* 0x359dffc2 */ | 58 s04 = 5.0463624390e-09, /* 0x31ad6446 */ |
61 s04 = 5.0463624390e-09, /* 0x31ad6446 */ | 59 s05 = 1.2354227016e-11; /* 0x2d59567e */ |
62 s05 = 1.2354227016e-11; /* 0x2d59567e */ | |
63 | 60 |
64 float j1f(float x) | 61 float j1f(float x) { |
65 { | 62 float z, r, s; |
66 » float z,r,s; | 63 uint32_t ix; |
67 » uint32_t ix; | 64 int sign; |
68 » int sign; | |
69 | 65 |
70 » GET_FLOAT_WORD(ix, x); | 66 GET_FLOAT_WORD(ix, x); |
71 » sign = ix>>31; | 67 sign = ix >> 31; |
72 » ix &= 0x7fffffff; | 68 ix &= 0x7fffffff; |
73 » if (ix >= 0x7f800000) | 69 if (ix >= 0x7f800000) |
74 » » return 1/(x*x); | 70 return 1 / (x * x); |
75 » if (ix >= 0x40000000) /* |x| >= 2 */ | 71 if (ix >= 0x40000000) /* |x| >= 2 */ |
76 » » return common(ix, fabsf(x), 0, sign); | 72 return common(ix, fabsf(x), 0, sign); |
77 » if (ix >= 0x32000000) { /* |x| >= 2**-27 */ | 73 if (ix >= 0x32000000) { /* |x| >= 2**-27 */ |
78 » » z = x*x; | 74 z = x * x; |
79 » » r = z*(r00+z*(r01+z*(r02+z*r03))); | 75 r = z * (r00 + z * (r01 + z * (r02 + z * r03))); |
80 » » s = 1+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); | 76 s = 1 + z * (s01 + z * (s02 + z * (s03 + z * (s04 + z * s05)))); |
81 » » z = 0.5f + r/s; | 77 z = 0.5f + r / s; |
82 » } else | 78 } else |
83 » » /* raise inexact if x!=0 */ | 79 /* raise inexact if x!=0 */ |
84 » » z = 0.5f + x; | 80 z = 0.5f + x; |
85 » return z*x; | 81 return z * x; |
86 } | 82 } |
87 | 83 |
88 static const float U0[5] = { | 84 static const float U0[5] = { |
89 -1.9605709612e-01, /* 0xbe48c331 */ | 85 -1.9605709612e-01, /* 0xbe48c331 */ |
90 5.0443872809e-02, /* 0x3d4e9e3c */ | 86 5.0443872809e-02, /* 0x3d4e9e3c */ |
91 -1.9125689287e-03, /* 0xbafaaf2a */ | 87 -1.9125689287e-03, /* 0xbafaaf2a */ |
92 2.3525259166e-05, /* 0x37c5581c */ | 88 2.3525259166e-05, /* 0x37c5581c */ |
93 -9.1909917899e-08, /* 0xb3c56003 */ | 89 -9.1909917899e-08, /* 0xb3c56003 */ |
94 }; | 90 }; |
95 static const float V0[5] = { | 91 static const float V0[5] = { |
96 1.9916731864e-02, /* 0x3ca3286a */ | 92 1.9916731864e-02, /* 0x3ca3286a */ |
97 2.0255257550e-04, /* 0x3954644b */ | 93 2.0255257550e-04, /* 0x3954644b */ |
98 1.3560879779e-06, /* 0x35b602d4 */ | 94 1.3560879779e-06, /* 0x35b602d4 */ |
99 6.2274145840e-09, /* 0x31d5f8eb */ | 95 6.2274145840e-09, /* 0x31d5f8eb */ |
100 1.6655924903e-11, /* 0x2d9281cf */ | 96 1.6655924903e-11, /* 0x2d9281cf */ |
101 }; | 97 }; |
102 | 98 |
103 float y1f(float x) | 99 float y1f(float x) { |
104 { | 100 float z, u, v; |
105 » float z,u,v; | 101 uint32_t ix; |
106 » uint32_t ix; | |
107 | 102 |
108 » GET_FLOAT_WORD(ix, x); | 103 GET_FLOAT_WORD(ix, x); |
109 » if ((ix & 0x7fffffff) == 0) | 104 if ((ix & 0x7fffffff) == 0) |
110 » » return -1/0.0f; | 105 return -1 / 0.0f; |
111 » if (ix>>31) | 106 if (ix >> 31) |
112 » » return 0/0.0f; | 107 return 0 / 0.0f; |
113 » if (ix >= 0x7f800000) | 108 if (ix >= 0x7f800000) |
114 » » return 1/x; | 109 return 1 / x; |
115 » if (ix >= 0x40000000) /* |x| >= 2.0 */ | 110 if (ix >= 0x40000000) /* |x| >= 2.0 */ |
116 » » return common(ix,x,1,0); | 111 return common(ix, x, 1, 0); |
117 » if (ix < 0x32000000) /* x < 2**-27 */ | 112 if (ix < 0x32000000) /* x < 2**-27 */ |
118 » » return -tpi/x; | 113 return -tpi / x; |
119 » z = x*x; | 114 z = x * x; |
120 » u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); | 115 u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * U0[4]))); |
121 » v = 1.0f+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); | 116 v = 1.0f + z * (V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * V0[4])))); |
122 » return x*(u/v) + tpi*(j1f(x)*logf(x)-1.0f/x); | 117 return x * (u / v) + tpi * (j1f(x) * logf(x) - 1.0f / x); |
123 } | 118 } |
124 | 119 |
125 /* For x >= 8, the asymptotic expansions of pone is | 120 /* For x >= 8, the asymptotic expansions of pone is |
126 * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. | 121 * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. |
127 * We approximate pone by | 122 * We approximate pone by |
128 * pone(x) = 1 + (R/S) | 123 * pone(x) = 1 + (R/S) |
129 * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 | 124 * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 |
130 * S = 1 + ps0*s^2 + ... + ps4*s^10 | 125 * S = 1 + ps0*s^2 + ... + ps4*s^10 |
131 * and | 126 * and |
132 * | pone(x)-1-R/S | <= 2 ** ( -60.06) | 127 * | pone(x)-1-R/S | <= 2 ** ( -60.06) |
133 */ | 128 */ |
134 | 129 |
135 static const float pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ | 130 static const float pr8[6] = { |
136 0.0000000000e+00, /* 0x00000000 */ | 131 /* for x in [inf, 8]=1/[0,0.125] */ |
137 1.1718750000e-01, /* 0x3df00000 */ | 132 0.0000000000e+00, /* 0x00000000 */ |
138 1.3239480972e+01, /* 0x4153d4ea */ | 133 1.1718750000e-01, /* 0x3df00000 */ |
139 4.1205184937e+02, /* 0x43ce06a3 */ | 134 1.3239480972e+01, /* 0x4153d4ea */ |
140 3.8747453613e+03, /* 0x45722bed */ | 135 4.1205184937e+02, /* 0x43ce06a3 */ |
141 7.9144794922e+03, /* 0x45f753d6 */ | 136 3.8747453613e+03, /* 0x45722bed */ |
| 137 7.9144794922e+03, /* 0x45f753d6 */ |
142 }; | 138 }; |
143 static const float ps8[5] = { | 139 static const float ps8[5] = { |
144 1.1420736694e+02, /* 0x42e46a2c */ | 140 1.1420736694e+02, /* 0x42e46a2c */ |
145 3.6509309082e+03, /* 0x45642ee5 */ | 141 3.6509309082e+03, /* 0x45642ee5 */ |
146 3.6956207031e+04, /* 0x47105c35 */ | 142 3.6956207031e+04, /* 0x47105c35 */ |
147 9.7602796875e+04, /* 0x47bea166 */ | 143 9.7602796875e+04, /* 0x47bea166 */ |
148 3.0804271484e+04, /* 0x46f0a88b */ | 144 3.0804271484e+04, /* 0x46f0a88b */ |
149 }; | 145 }; |
150 | 146 |
151 static const float pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ | 147 static const float pr5[6] = { |
152 1.3199052094e-11, /* 0x2d68333f */ | 148 /* for x in [8,4.5454]=1/[0.125,0.22001] */ |
153 1.1718749255e-01, /* 0x3defffff */ | 149 1.3199052094e-11, /* 0x2d68333f */ |
154 6.8027510643e+00, /* 0x40d9b023 */ | 150 1.1718749255e-01, /* 0x3defffff */ |
155 1.0830818176e+02, /* 0x42d89dca */ | 151 6.8027510643e+00, /* 0x40d9b023 */ |
156 5.1763616943e+02, /* 0x440168b7 */ | 152 1.0830818176e+02, /* 0x42d89dca */ |
157 5.2871520996e+02, /* 0x44042dc6 */ | 153 5.1763616943e+02, /* 0x440168b7 */ |
| 154 5.2871520996e+02, /* 0x44042dc6 */ |
158 }; | 155 }; |
159 static const float ps5[5] = { | 156 static const float ps5[5] = { |
160 5.9280597687e+01, /* 0x426d1f55 */ | 157 5.9280597687e+01, /* 0x426d1f55 */ |
161 9.9140142822e+02, /* 0x4477d9b1 */ | 158 9.9140142822e+02, /* 0x4477d9b1 */ |
162 5.3532670898e+03, /* 0x45a74a23 */ | 159 5.3532670898e+03, /* 0x45a74a23 */ |
163 7.8446904297e+03, /* 0x45f52586 */ | 160 7.8446904297e+03, /* 0x45f52586 */ |
164 1.5040468750e+03, /* 0x44bc0180 */ | 161 1.5040468750e+03, /* 0x44bc0180 */ |
165 }; | 162 }; |
166 | 163 |
167 static const float pr3[6] = { | 164 static const float pr3[6] = { |
168 3.0250391081e-09, /* 0x314fe10d */ | 165 3.0250391081e-09, /* 0x314fe10d */ |
169 1.1718686670e-01, /* 0x3defffab */ | 166 1.1718686670e-01, /* 0x3defffab */ |
170 3.9329774380e+00, /* 0x407bb5e7 */ | 167 3.9329774380e+00, /* 0x407bb5e7 */ |
171 3.5119403839e+01, /* 0x420c7a45 */ | 168 3.5119403839e+01, /* 0x420c7a45 */ |
172 9.1055007935e+01, /* 0x42b61c2a */ | 169 9.1055007935e+01, /* 0x42b61c2a */ |
173 4.8559066772e+01, /* 0x42423c7c */ | 170 4.8559066772e+01, /* 0x42423c7c */ |
174 }; | 171 }; |
175 static const float ps3[5] = { | 172 static const float ps3[5] = { |
176 3.4791309357e+01, /* 0x420b2a4d */ | 173 3.4791309357e+01, /* 0x420b2a4d */ |
177 3.3676245117e+02, /* 0x43a86198 */ | 174 3.3676245117e+02, /* 0x43a86198 */ |
178 1.0468714600e+03, /* 0x4482dbe3 */ | 175 1.0468714600e+03, /* 0x4482dbe3 */ |
179 8.9081134033e+02, /* 0x445eb3ed */ | 176 8.9081134033e+02, /* 0x445eb3ed */ |
180 1.0378793335e+02, /* 0x42cf936c */ | 177 1.0378793335e+02, /* 0x42cf936c */ |
181 }; | 178 }; |
182 | 179 |
183 static const float pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ | 180 static const float pr2[6] = { |
184 1.0771083225e-07, /* 0x33e74ea8 */ | 181 /* for x in [2.8570,2]=1/[0.3499,0.5] */ |
185 1.1717621982e-01, /* 0x3deffa16 */ | 182 1.0771083225e-07, /* 0x33e74ea8 */ |
186 2.3685150146e+00, /* 0x401795c0 */ | 183 1.1717621982e-01, /* 0x3deffa16 */ |
187 1.2242610931e+01, /* 0x4143e1bc */ | 184 2.3685150146e+00, /* 0x401795c0 */ |
188 1.7693971634e+01, /* 0x418d8d41 */ | 185 1.2242610931e+01, /* 0x4143e1bc */ |
189 5.0735230446e+00, /* 0x40a25a4d */ | 186 1.7693971634e+01, /* 0x418d8d41 */ |
| 187 5.0735230446e+00, /* 0x40a25a4d */ |
190 }; | 188 }; |
191 static const float ps2[5] = { | 189 static const float ps2[5] = { |
192 2.1436485291e+01, /* 0x41ab7dec */ | 190 2.1436485291e+01, /* 0x41ab7dec */ |
193 1.2529022980e+02, /* 0x42fa9499 */ | 191 1.2529022980e+02, /* 0x42fa9499 */ |
194 2.3227647400e+02, /* 0x436846c7 */ | 192 2.3227647400e+02, /* 0x436846c7 */ |
195 1.1767937469e+02, /* 0x42eb5bd7 */ | 193 1.1767937469e+02, /* 0x42eb5bd7 */ |
196 8.3646392822e+00, /* 0x4105d590 */ | 194 8.3646392822e+00, /* 0x4105d590 */ |
197 }; | 195 }; |
198 | 196 |
199 static float ponef(float x) | 197 static float ponef(float x) { |
200 { | 198 const float *p, *q; |
201 » const float *p,*q; | 199 float_t z, r, s; |
202 » float_t z,r,s; | 200 uint32_t ix; |
203 » uint32_t ix; | |
204 | 201 |
205 » GET_FLOAT_WORD(ix, x); | 202 GET_FLOAT_WORD(ix, x); |
206 » ix &= 0x7fffffff; | 203 ix &= 0x7fffffff; |
207 » if (ix >= 0x41000000){p = pr8; q = ps8;} | 204 if (ix >= 0x41000000) { |
208 » else if (ix >= 0x40f71c58){p = pr5; q = ps5;} | 205 p = pr8; |
209 » else if (ix >= 0x4036db68){p = pr3; q = ps3;} | 206 q = ps8; |
210 » else /*ix >= 0x40000000*/ {p = pr2; q = ps2;} | 207 } else if (ix >= 0x40f71c58) { |
211 » z = 1.0f/(x*x); | 208 p = pr5; |
212 » r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); | 209 q = ps5; |
213 » s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); | 210 } else if (ix >= 0x4036db68) { |
214 » return 1.0f + r/s; | 211 p = pr3; |
| 212 q = ps3; |
| 213 } else /*ix >= 0x40000000*/ { |
| 214 p = pr2; |
| 215 q = ps2; |
| 216 } |
| 217 z = 1.0f / (x * x); |
| 218 r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5])))); |
| 219 s = 1.0f + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4])))); |
| 220 return 1.0f + r / s; |
215 } | 221 } |
216 | 222 |
217 /* For x >= 8, the asymptotic expansions of qone is | 223 /* For x >= 8, the asymptotic expansions of qone is |
218 * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. | 224 * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. |
219 * We approximate pone by | 225 * We approximate pone by |
220 * qone(x) = s*(0.375 + (R/S)) | 226 * qone(x) = s*(0.375 + (R/S)) |
221 * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10 | 227 * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10 |
222 * S = 1 + qs1*s^2 + ... + qs6*s^12 | 228 * S = 1 + qs1*s^2 + ... + qs6*s^12 |
223 * and | 229 * and |
224 * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13) | 230 * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13) |
225 */ | 231 */ |
226 | 232 |
227 static const float qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ | 233 static const float qr8[6] = { |
228 0.0000000000e+00, /* 0x00000000 */ | 234 /* for x in [inf, 8]=1/[0,0.125] */ |
229 -1.0253906250e-01, /* 0xbdd20000 */ | 235 0.0000000000e+00, /* 0x00000000 */ |
230 -1.6271753311e+01, /* 0xc1822c8d */ | 236 -1.0253906250e-01, /* 0xbdd20000 */ |
231 -7.5960174561e+02, /* 0xc43de683 */ | 237 -1.6271753311e+01, /* 0xc1822c8d */ |
232 -1.1849806641e+04, /* 0xc639273a */ | 238 -7.5960174561e+02, /* 0xc43de683 */ |
233 -4.8438511719e+04, /* 0xc73d3683 */ | 239 -1.1849806641e+04, /* 0xc639273a */ |
| 240 -4.8438511719e+04, /* 0xc73d3683 */ |
234 }; | 241 }; |
235 static const float qs8[6] = { | 242 static const float qs8[6] = { |
236 1.6139537048e+02, /* 0x43216537 */ | 243 1.6139537048e+02, /* 0x43216537 */ |
237 7.8253862305e+03, /* 0x45f48b17 */ | 244 7.8253862305e+03, /* 0x45f48b17 */ |
238 1.3387534375e+05, /* 0x4802bcd6 */ | 245 1.3387534375e+05, /* 0x4802bcd6 */ |
239 7.1965775000e+05, /* 0x492fb29c */ | 246 7.1965775000e+05, /* 0x492fb29c */ |
240 6.6660125000e+05, /* 0x4922be94 */ | 247 6.6660125000e+05, /* 0x4922be94 */ |
241 -2.9449025000e+05, /* 0xc88fcb48 */ | 248 -2.9449025000e+05, /* 0xc88fcb48 */ |
242 }; | 249 }; |
243 | 250 |
244 static const float qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ | 251 static const float qr5[6] = { |
245 -2.0897993405e-11, /* 0xadb7d219 */ | 252 /* for x in [8,4.5454]=1/[0.125,0.22001] */ |
246 -1.0253904760e-01, /* 0xbdd1fffe */ | 253 -2.0897993405e-11, /* 0xadb7d219 */ |
247 -8.0564479828e+00, /* 0xc100e736 */ | 254 -1.0253904760e-01, /* 0xbdd1fffe */ |
248 -1.8366960144e+02, /* 0xc337ab6b */ | 255 -8.0564479828e+00, /* 0xc100e736 */ |
249 -1.3731937256e+03, /* 0xc4aba633 */ | 256 -1.8366960144e+02, /* 0xc337ab6b */ |
250 -2.6124443359e+03, /* 0xc523471c */ | 257 -1.3731937256e+03, /* 0xc4aba633 */ |
| 258 -2.6124443359e+03, /* 0xc523471c */ |
251 }; | 259 }; |
252 static const float qs5[6] = { | 260 static const float qs5[6] = { |
253 8.1276550293e+01, /* 0x42a28d98 */ | 261 8.1276550293e+01, /* 0x42a28d98 */ |
254 1.9917987061e+03, /* 0x44f8f98f */ | 262 1.9917987061e+03, /* 0x44f8f98f */ |
255 1.7468484375e+04, /* 0x468878f8 */ | 263 1.7468484375e+04, /* 0x468878f8 */ |
256 4.9851425781e+04, /* 0x4742bb6d */ | 264 4.9851425781e+04, /* 0x4742bb6d */ |
257 2.7948074219e+04, /* 0x46da5826 */ | 265 2.7948074219e+04, /* 0x46da5826 */ |
258 -4.7191835938e+03, /* 0xc5937978 */ | 266 -4.7191835938e+03, /* 0xc5937978 */ |
259 }; | 267 }; |
260 | 268 |
261 static const float qr3[6] = { | 269 static const float qr3[6] = { |
262 -5.0783124372e-09, /* 0xb1ae7d4f */ | 270 -5.0783124372e-09, /* 0xb1ae7d4f */ |
263 -1.0253783315e-01, /* 0xbdd1ff5b */ | 271 -1.0253783315e-01, /* 0xbdd1ff5b */ |
264 -4.6101160049e+00, /* 0xc0938612 */ | 272 -4.6101160049e+00, /* 0xc0938612 */ |
265 -5.7847221375e+01, /* 0xc267638e */ | 273 -5.7847221375e+01, /* 0xc267638e */ |
266 -2.2824453735e+02, /* 0xc3643e9a */ | 274 -2.2824453735e+02, /* 0xc3643e9a */ |
267 -2.1921012878e+02, /* 0xc35b35cb */ | 275 -2.1921012878e+02, /* 0xc35b35cb */ |
268 }; | 276 }; |
269 static const float qs3[6] = { | 277 static const float qs3[6] = { |
270 4.7665153503e+01, /* 0x423ea91e */ | 278 4.7665153503e+01, /* 0x423ea91e */ |
271 6.7386511230e+02, /* 0x4428775e */ | 279 6.7386511230e+02, /* 0x4428775e */ |
272 3.3801528320e+03, /* 0x45534272 */ | 280 3.3801528320e+03, /* 0x45534272 */ |
273 5.5477290039e+03, /* 0x45ad5dd5 */ | 281 5.5477290039e+03, /* 0x45ad5dd5 */ |
274 1.9031191406e+03, /* 0x44ede3d0 */ | 282 1.9031191406e+03, /* 0x44ede3d0 */ |
275 -1.3520118713e+02, /* 0xc3073381 */ | 283 -1.3520118713e+02, /* 0xc3073381 */ |
276 }; | 284 }; |
277 | 285 |
278 static const float qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ | 286 static const float qr2[6] = { |
279 -1.7838172539e-07, /* 0xb43f8932 */ | 287 /* for x in [2.8570,2]=1/[0.3499,0.5] */ |
280 -1.0251704603e-01, /* 0xbdd1f475 */ | 288 -1.7838172539e-07, /* 0xb43f8932 */ |
281 -2.7522056103e+00, /* 0xc0302423 */ | 289 -1.0251704603e-01, /* 0xbdd1f475 */ |
282 -1.9663616180e+01, /* 0xc19d4f16 */ | 290 -2.7522056103e+00, /* 0xc0302423 */ |
283 -4.2325313568e+01, /* 0xc2294d1f */ | 291 -1.9663616180e+01, /* 0xc19d4f16 */ |
284 -2.1371921539e+01, /* 0xc1aaf9b2 */ | 292 -4.2325313568e+01, /* 0xc2294d1f */ |
| 293 -2.1371921539e+01, /* 0xc1aaf9b2 */ |
285 }; | 294 }; |
286 static const float qs2[6] = { | 295 static const float qs2[6] = { |
287 2.9533363342e+01, /* 0x41ec4454 */ | 296 2.9533363342e+01, /* 0x41ec4454 */ |
288 2.5298155212e+02, /* 0x437cfb47 */ | 297 2.5298155212e+02, /* 0x437cfb47 */ |
289 7.5750280762e+02, /* 0x443d602e */ | 298 7.5750280762e+02, /* 0x443d602e */ |
290 7.3939318848e+02, /* 0x4438d92a */ | 299 7.3939318848e+02, /* 0x4438d92a */ |
291 1.5594900513e+02, /* 0x431bf2f2 */ | 300 1.5594900513e+02, /* 0x431bf2f2 */ |
292 -4.9594988823e+00, /* 0xc09eb437 */ | 301 -4.9594988823e+00, /* 0xc09eb437 */ |
293 }; | 302 }; |
294 | 303 |
295 static float qonef(float x) | 304 static float qonef(float x) { |
296 { | 305 const float *p, *q; |
297 » const float *p,*q; | 306 float_t s, r, z; |
298 » float_t s,r,z; | 307 uint32_t ix; |
299 » uint32_t ix; | |
300 | 308 |
301 » GET_FLOAT_WORD(ix, x); | 309 GET_FLOAT_WORD(ix, x); |
302 » ix &= 0x7fffffff; | 310 ix &= 0x7fffffff; |
303 » if (ix >= 0x40200000){p = qr8; q = qs8;} | 311 if (ix >= 0x40200000) { |
304 » else if (ix >= 0x40f71c58){p = qr5; q = qs5;} | 312 p = qr8; |
305 » else if (ix >= 0x4036db68){p = qr3; q = qs3;} | 313 q = qs8; |
306 » else /*ix >= 0x40000000*/ {p = qr2; q = qs2;} | 314 } else if (ix >= 0x40f71c58) { |
307 » z = 1.0f/(x*x); | 315 p = qr5; |
308 » r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); | 316 q = qs5; |
309 » s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); | 317 } else if (ix >= 0x4036db68) { |
310 » return (.375f + r/s)/x; | 318 p = qr3; |
| 319 q = qs3; |
| 320 } else /*ix >= 0x40000000*/ { |
| 321 p = qr2; |
| 322 q = qs2; |
| 323 } |
| 324 z = 1.0f / (x * x); |
| 325 r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5])))); |
| 326 s = 1.0f + |
| 327 z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5]))))); |
| 328 return (.375f + r / s) / x; |
311 } | 329 } |
OLD | NEW |