OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_j0f.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_j0f.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 pzerof(float), qzerof(float); | 19 static float pzerof(float), qzerof(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 y0) | 24 static float common(uint32_t ix, float x, int y0) { |
26 { | 25 float z, s, c, ss, cc; |
27 » float z,s,c,ss,cc; | 26 /* |
28 » /* | 27 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) |
29 » * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) | 28 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) |
30 » * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) | 29 */ |
31 » */ | 30 s = sinf(x); |
32 » s = sinf(x); | 31 c = cosf(x); |
33 » c = cosf(x); | 32 if (y0) |
34 » if (y0) | 33 c = -c; |
35 » » c = -c; | 34 cc = s + c; |
36 » cc = s+c; | 35 if (ix < 0x7f000000) { |
37 » if (ix < 0x7f000000) { | 36 ss = s - c; |
38 » » ss = s-c; | 37 z = -cosf(2 * x); |
39 » » z = -cosf(2*x); | 38 if (s * c < 0) |
40 » » if (s*c < 0) | 39 cc = z / ss; |
41 » » » cc = z/ss; | 40 else |
42 » » else | 41 ss = z / cc; |
43 » » » ss = z/cc; | 42 if (ix < 0x58800000) { |
44 » » if (ix < 0x58800000) { | 43 if (y0) |
45 » » » if (y0) | 44 ss = -ss; |
46 » » » » ss = -ss; | 45 cc = pzerof(x) * cc - qzerof(x) * ss; |
47 » » » cc = pzerof(x)*cc-qzerof(x)*ss; | 46 } |
48 » » } | 47 } |
49 » } | 48 return invsqrtpi * cc / sqrtf(x); |
50 » return invsqrtpi*cc/sqrtf(x); | |
51 } | 49 } |
52 | 50 |
53 /* R0/S0 on [0, 2.00] */ | 51 /* R0/S0 on [0, 2.00] */ |
54 static const float | 52 static const float R02 = 1.5625000000e-02, /* 0x3c800000 */ |
55 R02 = 1.5625000000e-02, /* 0x3c800000 */ | 53 R03 = -1.8997929874e-04, /* 0xb947352e */ |
56 R03 = -1.8997929874e-04, /* 0xb947352e */ | 54 R04 = 1.8295404516e-06, /* 0x35f58e88 */ |
57 R04 = 1.8295404516e-06, /* 0x35f58e88 */ | 55 R05 = -4.6183270541e-09, /* 0xb19eaf3c */ |
58 R05 = -4.6183270541e-09, /* 0xb19eaf3c */ | 56 S01 = 1.5619102865e-02, /* 0x3c7fe744 */ |
59 S01 = 1.5619102865e-02, /* 0x3c7fe744 */ | 57 S02 = 1.1692678527e-04, /* 0x38f53697 */ |
60 S02 = 1.1692678527e-04, /* 0x38f53697 */ | 58 S03 = 5.1354652442e-07, /* 0x3509daa6 */ |
61 S03 = 5.1354652442e-07, /* 0x3509daa6 */ | 59 S04 = 1.1661400734e-09; /* 0x30a045e8 */ |
62 S04 = 1.1661400734e-09; /* 0x30a045e8 */ | |
63 | 60 |
64 float j0f(float x) | 61 float j0f(float x) { |
65 { | 62 float z, r, s; |
66 » float z,r,s; | 63 uint32_t ix; |
67 » uint32_t ix; | |
68 | 64 |
69 » GET_FLOAT_WORD(ix, x); | 65 GET_FLOAT_WORD(ix, x); |
70 » ix &= 0x7fffffff; | 66 ix &= 0x7fffffff; |
71 » if (ix >= 0x7f800000) | 67 if (ix >= 0x7f800000) |
72 » » return 1/(x*x); | 68 return 1 / (x * x); |
73 » x = fabsf(x); | 69 x = fabsf(x); |
74 | 70 |
75 » if (ix >= 0x40000000) { /* |x| >= 2 */ | 71 if (ix >= 0x40000000) { /* |x| >= 2 */ |
76 » » /* large ulp error near zeros */ | 72 /* large ulp error near zeros */ |
77 » » return common(ix, x, 0); | 73 return common(ix, x, 0); |
78 » } | 74 } |
79 » if (ix >= 0x3a000000) { /* |x| >= 2**-11 */ | 75 if (ix >= 0x3a000000) { /* |x| >= 2**-11 */ |
80 » » /* up to 4ulp error near 2 */ | 76 /* up to 4ulp error near 2 */ |
81 » » z = x*x; | 77 z = x * x; |
82 » » r = z*(R02+z*(R03+z*(R04+z*R05))); | 78 r = z * (R02 + z * (R03 + z * (R04 + z * R05))); |
83 » » s = 1+z*(S01+z*(S02+z*(S03+z*S04))); | 79 s = 1 + z * (S01 + z * (S02 + z * (S03 + z * S04))); |
84 » » return (1+x/2)*(1-x/2) + z*(r/s); | 80 return (1 + x / 2) * (1 - x / 2) + z * (r / s); |
85 » } | 81 } |
86 » if (ix >= 0x21800000) /* |x| >= 2**-60 */ | 82 if (ix >= 0x21800000) /* |x| >= 2**-60 */ |
87 » » x = 0.25f*x*x; | 83 x = 0.25f * x * x; |
88 » return 1 - x; | 84 return 1 - x; |
89 } | 85 } |
90 | 86 |
91 static const float | 87 static const float u00 = -7.3804296553e-02, /* 0xbd9726b5 */ |
92 u00 = -7.3804296553e-02, /* 0xbd9726b5 */ | 88 u01 = 1.7666645348e-01, /* 0x3e34e80d */ |
93 u01 = 1.7666645348e-01, /* 0x3e34e80d */ | 89 u02 = -1.3818567619e-02, /* 0xbc626746 */ |
94 u02 = -1.3818567619e-02, /* 0xbc626746 */ | 90 u03 = 3.4745343146e-04, /* 0x39b62a69 */ |
95 u03 = 3.4745343146e-04, /* 0x39b62a69 */ | 91 u04 = -3.8140706238e-06, /* 0xb67ff53c */ |
96 u04 = -3.8140706238e-06, /* 0xb67ff53c */ | 92 u05 = 1.9559013964e-08, /* 0x32a802ba */ |
97 u05 = 1.9559013964e-08, /* 0x32a802ba */ | 93 u06 = -3.9820518410e-11, /* 0xae2f21eb */ |
98 u06 = -3.9820518410e-11, /* 0xae2f21eb */ | 94 v01 = 1.2730483897e-02, /* 0x3c509385 */ |
99 v01 = 1.2730483897e-02, /* 0x3c509385 */ | 95 v02 = 7.6006865129e-05, /* 0x389f65e0 */ |
100 v02 = 7.6006865129e-05, /* 0x389f65e0 */ | 96 v03 = 2.5915085189e-07, /* 0x348b216c */ |
101 v03 = 2.5915085189e-07, /* 0x348b216c */ | 97 v04 = 4.4111031494e-10; /* 0x2ff280c2 */ |
102 v04 = 4.4111031494e-10; /* 0x2ff280c2 */ | |
103 | 98 |
104 float y0f(float x) | 99 float y0f(float x) { |
105 { | 100 float z, u, v; |
106 » float z,u,v; | 101 uint32_t ix; |
107 » uint32_t ix; | |
108 | 102 |
109 » GET_FLOAT_WORD(ix, x); | 103 GET_FLOAT_WORD(ix, x); |
110 » if ((ix & 0x7fffffff) == 0) | 104 if ((ix & 0x7fffffff) == 0) |
111 » » return -1/0.0f; | 105 return -1 / 0.0f; |
112 » if (ix>>31) | 106 if (ix >> 31) |
113 » » return 0/0.0f; | 107 return 0 / 0.0f; |
114 » if (ix >= 0x7f800000) | 108 if (ix >= 0x7f800000) |
115 » » return 1/x; | 109 return 1 / x; |
116 » if (ix >= 0x40000000) { /* |x| >= 2.0 */ | 110 if (ix >= 0x40000000) { /* |x| >= 2.0 */ |
117 » » /* large ulp error near zeros */ | 111 /* large ulp error near zeros */ |
118 » » return common(ix,x,1); | 112 return common(ix, x, 1); |
119 » } | 113 } |
120 » if (ix >= 0x39000000) { /* x >= 2**-13 */ | 114 if (ix >= 0x39000000) { /* x >= 2**-13 */ |
121 » » /* large ulp error at x ~= 0.89 */ | 115 /* large ulp error at x ~= 0.89 */ |
122 » » z = x*x; | 116 z = x * x; |
123 » » u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); | 117 u = u00 + |
124 » » v = 1+z*(v01+z*(v02+z*(v03+z*v04))); | 118 z * (u01 + z * (u02 + z * (u03 + z * (u04 + z * (u05 + z * u06))))); |
125 » » return u/v + tpi*(j0f(x)*logf(x)); | 119 v = 1 + z * (v01 + z * (v02 + z * (v03 + z * v04))); |
126 » } | 120 return u / v + tpi * (j0f(x) * logf(x)); |
127 » return u00 + tpi*logf(x); | 121 } |
| 122 return u00 + tpi * logf(x); |
128 } | 123 } |
129 | 124 |
130 /* The asymptotic expansions of pzero is | 125 /* The asymptotic expansions of pzero is |
131 * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. | 126 * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. |
132 * For x >= 2, We approximate pzero by | 127 * For x >= 2, We approximate pzero by |
133 * pzero(x) = 1 + (R/S) | 128 * pzero(x) = 1 + (R/S) |
134 * where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10 | 129 * where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10 |
135 * S = 1 + pS0*s^2 + ... + pS4*s^10 | 130 * S = 1 + pS0*s^2 + ... + pS4*s^10 |
136 * and | 131 * and |
137 * | pzero(x)-1-R/S | <= 2 ** ( -60.26) | 132 * | pzero(x)-1-R/S | <= 2 ** ( -60.26) |
138 */ | 133 */ |
139 static const float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ | 134 static const float pR8[6] = { |
140 0.0000000000e+00, /* 0x00000000 */ | 135 /* for x in [inf, 8]=1/[0,0.125] */ |
141 -7.0312500000e-02, /* 0xbd900000 */ | 136 0.0000000000e+00, /* 0x00000000 */ |
142 -8.0816707611e+00, /* 0xc1014e86 */ | 137 -7.0312500000e-02, /* 0xbd900000 */ |
143 -2.5706311035e+02, /* 0xc3808814 */ | 138 -8.0816707611e+00, /* 0xc1014e86 */ |
144 -2.4852163086e+03, /* 0xc51b5376 */ | 139 -2.5706311035e+02, /* 0xc3808814 */ |
145 -5.2530439453e+03, /* 0xc5a4285a */ | 140 -2.4852163086e+03, /* 0xc51b5376 */ |
| 141 -5.2530439453e+03, /* 0xc5a4285a */ |
146 }; | 142 }; |
147 static const float pS8[5] = { | 143 static const float pS8[5] = { |
148 1.1653436279e+02, /* 0x42e91198 */ | 144 1.1653436279e+02, /* 0x42e91198 */ |
149 3.8337448730e+03, /* 0x456f9beb */ | 145 3.8337448730e+03, /* 0x456f9beb */ |
150 4.0597855469e+04, /* 0x471e95db */ | 146 4.0597855469e+04, /* 0x471e95db */ |
151 1.1675296875e+05, /* 0x47e4087c */ | 147 1.1675296875e+05, /* 0x47e4087c */ |
152 4.7627726562e+04, /* 0x473a0bba */ | 148 4.7627726562e+04, /* 0x473a0bba */ |
153 }; | 149 }; |
154 static const float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ | 150 static const float pR5[6] = { |
155 -1.1412546255e-11, /* 0xad48c58a */ | 151 /* for x in [8,4.5454]=1/[0.125,0.22001] */ |
156 -7.0312492549e-02, /* 0xbd8fffff */ | 152 -1.1412546255e-11, /* 0xad48c58a */ |
157 -4.1596107483e+00, /* 0xc0851b88 */ | 153 -7.0312492549e-02, /* 0xbd8fffff */ |
158 -6.7674766541e+01, /* 0xc287597b */ | 154 -4.1596107483e+00, /* 0xc0851b88 */ |
159 -3.3123129272e+02, /* 0xc3a59d9b */ | 155 -6.7674766541e+01, /* 0xc287597b */ |
160 -3.4643338013e+02, /* 0xc3ad3779 */ | 156 -3.3123129272e+02, /* 0xc3a59d9b */ |
| 157 -3.4643338013e+02, /* 0xc3ad3779 */ |
161 }; | 158 }; |
162 static const float pS5[5] = { | 159 static const float pS5[5] = { |
163 6.0753936768e+01, /* 0x42730408 */ | 160 6.0753936768e+01, /* 0x42730408 */ |
164 1.0512523193e+03, /* 0x44836813 */ | 161 1.0512523193e+03, /* 0x44836813 */ |
165 5.9789707031e+03, /* 0x45bad7c4 */ | 162 5.9789707031e+03, /* 0x45bad7c4 */ |
166 9.6254453125e+03, /* 0x461665c8 */ | 163 9.6254453125e+03, /* 0x461665c8 */ |
167 2.4060581055e+03, /* 0x451660ee */ | 164 2.4060581055e+03, /* 0x451660ee */ |
168 }; | 165 }; |
169 | 166 |
170 static const float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ | 167 static const float pR3[6] = { |
171 -2.5470459075e-09, /* 0xb12f081b */ | 168 /* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ |
172 -7.0311963558e-02, /* 0xbd8fffb8 */ | 169 -2.5470459075e-09, /* 0xb12f081b */ |
173 -2.4090321064e+00, /* 0xc01a2d95 */ | 170 -7.0311963558e-02, /* 0xbd8fffb8 */ |
174 -2.1965976715e+01, /* 0xc1afba52 */ | 171 -2.4090321064e+00, /* 0xc01a2d95 */ |
175 -5.8079170227e+01, /* 0xc2685112 */ | 172 -2.1965976715e+01, /* 0xc1afba52 */ |
176 -3.1447946548e+01, /* 0xc1fb9565 */ | 173 -5.8079170227e+01, /* 0xc2685112 */ |
| 174 -3.1447946548e+01, /* 0xc1fb9565 */ |
177 }; | 175 }; |
178 static const float pS3[5] = { | 176 static const float pS3[5] = { |
179 3.5856033325e+01, /* 0x420f6c94 */ | 177 3.5856033325e+01, /* 0x420f6c94 */ |
180 3.6151397705e+02, /* 0x43b4c1ca */ | 178 3.6151397705e+02, /* 0x43b4c1ca */ |
181 1.1936077881e+03, /* 0x44953373 */ | 179 1.1936077881e+03, /* 0x44953373 */ |
182 1.1279968262e+03, /* 0x448cffe6 */ | 180 1.1279968262e+03, /* 0x448cffe6 */ |
183 1.7358093262e+02, /* 0x432d94b8 */ | 181 1.7358093262e+02, /* 0x432d94b8 */ |
184 }; | 182 }; |
185 | 183 |
186 static const float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ | 184 static const float pR2[6] = { |
187 -8.8753431271e-08, /* 0xb3be98b7 */ | 185 /* for x in [2.8570,2]=1/[0.3499,0.5] */ |
188 -7.0303097367e-02, /* 0xbd8ffb12 */ | 186 -8.8753431271e-08, /* 0xb3be98b7 */ |
189 -1.4507384300e+00, /* 0xbfb9b1cc */ | 187 -7.0303097367e-02, /* 0xbd8ffb12 */ |
190 -7.6356959343e+00, /* 0xc0f4579f */ | 188 -1.4507384300e+00, /* 0xbfb9b1cc */ |
191 -1.1193166733e+01, /* 0xc1331736 */ | 189 -7.6356959343e+00, /* 0xc0f4579f */ |
192 -3.2336456776e+00, /* 0xc04ef40d */ | 190 -1.1193166733e+01, /* 0xc1331736 */ |
| 191 -3.2336456776e+00, /* 0xc04ef40d */ |
193 }; | 192 }; |
194 static const float pS2[5] = { | 193 static const float pS2[5] = { |
195 2.2220300674e+01, /* 0x41b1c32d */ | 194 2.2220300674e+01, /* 0x41b1c32d */ |
196 1.3620678711e+02, /* 0x430834f0 */ | 195 1.3620678711e+02, /* 0x430834f0 */ |
197 2.7047027588e+02, /* 0x43873c32 */ | 196 2.7047027588e+02, /* 0x43873c32 */ |
198 1.5387539673e+02, /* 0x4319e01a */ | 197 1.5387539673e+02, /* 0x4319e01a */ |
199 1.4657617569e+01, /* 0x416a859a */ | 198 1.4657617569e+01, /* 0x416a859a */ |
200 }; | 199 }; |
201 | 200 |
202 static float pzerof(float x) | 201 static float pzerof(float x) { |
203 { | 202 const float *p, *q; |
204 » const float *p,*q; | 203 float_t z, r, s; |
205 » float_t z,r,s; | 204 uint32_t ix; |
206 » uint32_t ix; | |
207 | 205 |
208 » GET_FLOAT_WORD(ix, x); | 206 GET_FLOAT_WORD(ix, x); |
209 » ix &= 0x7fffffff; | 207 ix &= 0x7fffffff; |
210 » if (ix >= 0x41000000){p = pR8; q = pS8;} | 208 if (ix >= 0x41000000) { |
211 » else if (ix >= 0x40f71c58){p = pR5; q = pS5;} | 209 p = pR8; |
212 » else if (ix >= 0x4036db68){p = pR3; q = pS3;} | 210 q = pS8; |
213 » else /*ix >= 0x40000000*/ {p = pR2; q = pS2;} | 211 } else if (ix >= 0x40f71c58) { |
214 » z = 1.0f/(x*x); | 212 p = pR5; |
215 » r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); | 213 q = pS5; |
216 » s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); | 214 } else if (ix >= 0x4036db68) { |
217 » return 1.0f + r/s; | 215 p = pR3; |
| 216 q = pS3; |
| 217 } else /*ix >= 0x40000000*/ { |
| 218 p = pR2; |
| 219 q = pS2; |
| 220 } |
| 221 z = 1.0f / (x * x); |
| 222 r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5])))); |
| 223 s = 1.0f + z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * q[4])))); |
| 224 return 1.0f + r / s; |
218 } | 225 } |
219 | 226 |
220 | |
221 /* For x >= 8, the asymptotic expansions of qzero is | 227 /* For x >= 8, the asymptotic expansions of qzero is |
222 * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. | 228 * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. |
223 * We approximate pzero by | 229 * We approximate pzero by |
224 * qzero(x) = s*(-1.25 + (R/S)) | 230 * qzero(x) = s*(-1.25 + (R/S)) |
225 * where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10 | 231 * where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10 |
226 * S = 1 + qS0*s^2 + ... + qS5*s^12 | 232 * S = 1 + qS0*s^2 + ... + qS5*s^12 |
227 * and | 233 * and |
228 * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) | 234 * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) |
229 */ | 235 */ |
230 static const float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ | 236 static const float qR8[6] = { |
231 0.0000000000e+00, /* 0x00000000 */ | 237 /* for x in [inf, 8]=1/[0,0.125] */ |
232 7.3242187500e-02, /* 0x3d960000 */ | 238 0.0000000000e+00, /* 0x00000000 */ |
233 1.1768206596e+01, /* 0x413c4a93 */ | 239 7.3242187500e-02, /* 0x3d960000 */ |
234 5.5767340088e+02, /* 0x440b6b19 */ | 240 1.1768206596e+01, /* 0x413c4a93 */ |
235 8.8591972656e+03, /* 0x460a6cca */ | 241 5.5767340088e+02, /* 0x440b6b19 */ |
236 3.7014625000e+04, /* 0x471096a0 */ | 242 8.8591972656e+03, /* 0x460a6cca */ |
| 243 3.7014625000e+04, /* 0x471096a0 */ |
237 }; | 244 }; |
238 static const float qS8[6] = { | 245 static const float qS8[6] = { |
239 1.6377603149e+02, /* 0x4323c6aa */ | 246 1.6377603149e+02, /* 0x4323c6aa */ |
240 8.0983447266e+03, /* 0x45fd12c2 */ | 247 8.0983447266e+03, /* 0x45fd12c2 */ |
241 1.4253829688e+05, /* 0x480b3293 */ | 248 1.4253829688e+05, /* 0x480b3293 */ |
242 8.0330925000e+05, /* 0x49441ed4 */ | 249 8.0330925000e+05, /* 0x49441ed4 */ |
243 8.4050156250e+05, /* 0x494d3359 */ | 250 8.4050156250e+05, /* 0x494d3359 */ |
244 -3.4389928125e+05, /* 0xc8a7eb69 */ | 251 -3.4389928125e+05, /* 0xc8a7eb69 */ |
245 }; | 252 }; |
246 | 253 |
247 static const float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ | 254 static const float qR5[6] = { |
248 1.8408595828e-11, /* 0x2da1ec79 */ | 255 /* for x in [8,4.5454]=1/[0.125,0.22001] */ |
249 7.3242180049e-02, /* 0x3d95ffff */ | 256 1.8408595828e-11, /* 0x2da1ec79 */ |
250 5.8356351852e+00, /* 0x40babd86 */ | 257 7.3242180049e-02, /* 0x3d95ffff */ |
251 1.3511157227e+02, /* 0x43071c90 */ | 258 5.8356351852e+00, /* 0x40babd86 */ |
252 1.0272437744e+03, /* 0x448067cd */ | 259 1.3511157227e+02, /* 0x43071c90 */ |
253 1.9899779053e+03, /* 0x44f8bf4b */ | 260 1.0272437744e+03, /* 0x448067cd */ |
| 261 1.9899779053e+03, /* 0x44f8bf4b */ |
254 }; | 262 }; |
255 static const float qS5[6] = { | 263 static const float qS5[6] = { |
256 8.2776611328e+01, /* 0x42a58da0 */ | 264 8.2776611328e+01, /* 0x42a58da0 */ |
257 2.0778142090e+03, /* 0x4501dd07 */ | 265 2.0778142090e+03, /* 0x4501dd07 */ |
258 1.8847289062e+04, /* 0x46933e94 */ | 266 1.8847289062e+04, /* 0x46933e94 */ |
259 5.6751113281e+04, /* 0x475daf1d */ | 267 5.6751113281e+04, /* 0x475daf1d */ |
260 3.5976753906e+04, /* 0x470c88c1 */ | 268 3.5976753906e+04, /* 0x470c88c1 */ |
261 -5.3543427734e+03, /* 0xc5a752be */ | 269 -5.3543427734e+03, /* 0xc5a752be */ |
262 }; | 270 }; |
263 | 271 |
264 static const float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ | 272 static const float qR3[6] = { |
265 4.3774099900e-09, /* 0x3196681b */ | 273 /* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ |
266 7.3241114616e-02, /* 0x3d95ff70 */ | 274 4.3774099900e-09, /* 0x3196681b */ |
267 3.3442313671e+00, /* 0x405607e3 */ | 275 7.3241114616e-02, /* 0x3d95ff70 */ |
268 4.2621845245e+01, /* 0x422a7cc5 */ | 276 3.3442313671e+00, /* 0x405607e3 */ |
269 1.7080809021e+02, /* 0x432acedf */ | 277 4.2621845245e+01, /* 0x422a7cc5 */ |
270 1.6673394775e+02, /* 0x4326bbe4 */ | 278 1.7080809021e+02, /* 0x432acedf */ |
| 279 1.6673394775e+02, /* 0x4326bbe4 */ |
271 }; | 280 }; |
272 static const float qS3[6] = { | 281 static const float qS3[6] = { |
273 4.8758872986e+01, /* 0x42430916 */ | 282 4.8758872986e+01, /* 0x42430916 */ |
274 7.0968920898e+02, /* 0x44316c1c */ | 283 7.0968920898e+02, /* 0x44316c1c */ |
275 3.7041481934e+03, /* 0x4567825f */ | 284 3.7041481934e+03, /* 0x4567825f */ |
276 6.4604252930e+03, /* 0x45c9e367 */ | 285 6.4604252930e+03, /* 0x45c9e367 */ |
277 2.5163337402e+03, /* 0x451d4557 */ | 286 2.5163337402e+03, /* 0x451d4557 */ |
278 -1.4924745178e+02, /* 0xc3153f59 */ | 287 -1.4924745178e+02, /* 0xc3153f59 */ |
279 }; | 288 }; |
280 | 289 |
281 static const float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ | 290 static const float qR2[6] = { |
282 1.5044444979e-07, /* 0x342189db */ | 291 /* for x in [2.8570,2]=1/[0.3499,0.5] */ |
283 7.3223426938e-02, /* 0x3d95f62a */ | 292 1.5044444979e-07, /* 0x342189db */ |
284 1.9981917143e+00, /* 0x3fffc4bf */ | 293 7.3223426938e-02, /* 0x3d95f62a */ |
285 1.4495602608e+01, /* 0x4167edfd */ | 294 1.9981917143e+00, /* 0x3fffc4bf */ |
286 3.1666231155e+01, /* 0x41fd5471 */ | 295 1.4495602608e+01, /* 0x4167edfd */ |
287 1.6252708435e+01, /* 0x4182058c */ | 296 3.1666231155e+01, /* 0x41fd5471 */ |
| 297 1.6252708435e+01, /* 0x4182058c */ |
288 }; | 298 }; |
289 static const float qS2[6] = { | 299 static const float qS2[6] = { |
290 3.0365585327e+01, /* 0x41f2ecb8 */ | 300 3.0365585327e+01, /* 0x41f2ecb8 */ |
291 2.6934811401e+02, /* 0x4386ac8f */ | 301 2.6934811401e+02, /* 0x4386ac8f */ |
292 8.4478375244e+02, /* 0x44533229 */ | 302 8.4478375244e+02, /* 0x44533229 */ |
293 8.8293585205e+02, /* 0x445cbbe5 */ | 303 8.8293585205e+02, /* 0x445cbbe5 */ |
294 2.1266638184e+02, /* 0x4354aa98 */ | 304 2.1266638184e+02, /* 0x4354aa98 */ |
295 -5.3109550476e+00, /* 0xc0a9f358 */ | 305 -5.3109550476e+00, /* 0xc0a9f358 */ |
296 }; | 306 }; |
297 | 307 |
298 static float qzerof(float x) | 308 static float qzerof(float x) { |
299 { | 309 const float *p, *q; |
300 » const float *p,*q; | 310 float_t s, r, z; |
301 » float_t s,r,z; | 311 uint32_t ix; |
302 » uint32_t ix; | |
303 | 312 |
304 » GET_FLOAT_WORD(ix, x); | 313 GET_FLOAT_WORD(ix, x); |
305 » ix &= 0x7fffffff; | 314 ix &= 0x7fffffff; |
306 » if (ix >= 0x41000000){p = qR8; q = qS8;} | 315 if (ix >= 0x41000000) { |
307 » else if (ix >= 0x40f71c58){p = qR5; q = qS5;} | 316 p = qR8; |
308 » else if (ix >= 0x4036db68){p = qR3; q = qS3;} | 317 q = qS8; |
309 » else /*ix >= 0x40000000*/ {p = qR2; q = qS2;} | 318 } else if (ix >= 0x40f71c58) { |
310 » z = 1.0f/(x*x); | 319 p = qR5; |
311 » r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); | 320 q = qS5; |
312 » s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); | 321 } else if (ix >= 0x4036db68) { |
313 » return (-.125f + r/s)/x; | 322 p = qR3; |
| 323 q = qS3; |
| 324 } else /*ix >= 0x40000000*/ { |
| 325 p = qR2; |
| 326 q = qS2; |
| 327 } |
| 328 z = 1.0f / (x * x); |
| 329 r = p[0] + z * (p[1] + z * (p[2] + z * (p[3] + z * (p[4] + z * p[5])))); |
| 330 s = 1.0f + |
| 331 z * (q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * q[5]))))); |
| 332 return (-.125f + r / s) / x; |
314 } | 333 } |
OLD | NEW |