| 1 |
207 |
jeremybenn |
/* ef_j0.c -- float version of e_j0.c.
|
| 2 |
|
|
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
| 3 |
|
|
*/
|
| 4 |
|
|
|
| 5 |
|
|
/*
|
| 6 |
|
|
* ====================================================
|
| 7 |
|
|
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
| 8 |
|
|
*
|
| 9 |
|
|
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
| 10 |
|
|
* Permission to use, copy, modify, and distribute this
|
| 11 |
|
|
* software is freely granted, provided that this notice
|
| 12 |
|
|
* is preserved.
|
| 13 |
|
|
* ====================================================
|
| 14 |
|
|
*/
|
| 15 |
|
|
|
| 16 |
|
|
#include "fdlibm.h"
|
| 17 |
|
|
|
| 18 |
|
|
#ifdef __STDC__
|
| 19 |
|
|
static float pzerof(float), qzerof(float);
|
| 20 |
|
|
#else
|
| 21 |
|
|
static float pzerof(), qzerof();
|
| 22 |
|
|
#endif
|
| 23 |
|
|
|
| 24 |
|
|
#ifdef __STDC__
|
| 25 |
|
|
static const float
|
| 26 |
|
|
#else
|
| 27 |
|
|
static float
|
| 28 |
|
|
#endif
|
| 29 |
|
|
huge = 1e30,
|
| 30 |
|
|
one = 1.0,
|
| 31 |
|
|
invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */
|
| 32 |
|
|
tpi = 6.3661974669e-01, /* 0x3f22f983 */
|
| 33 |
|
|
/* R0/S0 on [0, 2.00] */
|
| 34 |
|
|
R02 = 1.5625000000e-02, /* 0x3c800000 */
|
| 35 |
|
|
R03 = -1.8997929874e-04, /* 0xb947352e */
|
| 36 |
|
|
R04 = 1.8295404516e-06, /* 0x35f58e88 */
|
| 37 |
|
|
R05 = -4.6183270541e-09, /* 0xb19eaf3c */
|
| 38 |
|
|
S01 = 1.5619102865e-02, /* 0x3c7fe744 */
|
| 39 |
|
|
S02 = 1.1692678527e-04, /* 0x38f53697 */
|
| 40 |
|
|
S03 = 5.1354652442e-07, /* 0x3509daa6 */
|
| 41 |
|
|
S04 = 1.1661400734e-09; /* 0x30a045e8 */
|
| 42 |
|
|
|
| 43 |
|
|
#ifdef __STDC__
|
| 44 |
|
|
static const float zero = 0.0;
|
| 45 |
|
|
#else
|
| 46 |
|
|
static float zero = 0.0;
|
| 47 |
|
|
#endif
|
| 48 |
|
|
|
| 49 |
|
|
#ifdef __STDC__
|
| 50 |
|
|
float __ieee754_j0f(float x)
|
| 51 |
|
|
#else
|
| 52 |
|
|
float __ieee754_j0f(x)
|
| 53 |
|
|
float x;
|
| 54 |
|
|
#endif
|
| 55 |
|
|
{
|
| 56 |
|
|
float z, s,c,ss,cc,r,u,v;
|
| 57 |
|
|
__int32_t hx,ix;
|
| 58 |
|
|
|
| 59 |
|
|
GET_FLOAT_WORD(hx,x);
|
| 60 |
|
|
ix = hx&0x7fffffff;
|
| 61 |
|
|
if(!FLT_UWORD_IS_FINITE(ix)) return one/(x*x);
|
| 62 |
|
|
x = fabsf(x);
|
| 63 |
|
|
if(ix >= 0x40000000) { /* |x| >= 2.0 */
|
| 64 |
|
|
s = sinf(x);
|
| 65 |
|
|
c = cosf(x);
|
| 66 |
|
|
ss = s-c;
|
| 67 |
|
|
cc = s+c;
|
| 68 |
|
|
if(ix<=FLT_UWORD_HALF_MAX) { /* make sure x+x not overflow */
|
| 69 |
|
|
z = -cosf(x+x);
|
| 70 |
|
|
if ((s*c)<zero) cc = z/ss;
|
| 71 |
|
|
else ss = z/cc;
|
| 72 |
|
|
}
|
| 73 |
|
|
/*
|
| 74 |
|
|
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
|
| 75 |
|
|
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
|
| 76 |
|
|
*/
|
| 77 |
|
|
if(ix>0x80000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(x);
|
| 78 |
|
|
else {
|
| 79 |
|
|
u = pzerof(x); v = qzerof(x);
|
| 80 |
|
|
z = invsqrtpi*(u*cc-v*ss)/__ieee754_sqrtf(x);
|
| 81 |
|
|
}
|
| 82 |
|
|
return z;
|
| 83 |
|
|
}
|
| 84 |
|
|
if(ix<0x39000000) { /* |x| < 2**-13 */
|
| 85 |
|
|
if(huge+x>one) { /* raise inexact if x != 0 */
|
| 86 |
|
|
if(ix<0x32000000) return one; /* |x|<2**-27 */
|
| 87 |
|
|
else return one - (float)0.25*x*x;
|
| 88 |
|
|
}
|
| 89 |
|
|
}
|
| 90 |
|
|
z = x*x;
|
| 91 |
|
|
r = z*(R02+z*(R03+z*(R04+z*R05)));
|
| 92 |
|
|
s = one+z*(S01+z*(S02+z*(S03+z*S04)));
|
| 93 |
|
|
if(ix < 0x3F800000) { /* |x| < 1.00 */
|
| 94 |
|
|
return one + z*((float)-0.25+(r/s));
|
| 95 |
|
|
} else {
|
| 96 |
|
|
u = (float)0.5*x;
|
| 97 |
|
|
return((one+u)*(one-u)+z*(r/s));
|
| 98 |
|
|
}
|
| 99 |
|
|
}
|
| 100 |
|
|
|
| 101 |
|
|
#ifdef __STDC__
|
| 102 |
|
|
static const float
|
| 103 |
|
|
#else
|
| 104 |
|
|
static float
|
| 105 |
|
|
#endif
|
| 106 |
|
|
u00 = -7.3804296553e-02, /* 0xbd9726b5 */
|
| 107 |
|
|
u01 = 1.7666645348e-01, /* 0x3e34e80d */
|
| 108 |
|
|
u02 = -1.3818567619e-02, /* 0xbc626746 */
|
| 109 |
|
|
u03 = 3.4745343146e-04, /* 0x39b62a69 */
|
| 110 |
|
|
u04 = -3.8140706238e-06, /* 0xb67ff53c */
|
| 111 |
|
|
u05 = 1.9559013964e-08, /* 0x32a802ba */
|
| 112 |
|
|
u06 = -3.9820518410e-11, /* 0xae2f21eb */
|
| 113 |
|
|
v01 = 1.2730483897e-02, /* 0x3c509385 */
|
| 114 |
|
|
v02 = 7.6006865129e-05, /* 0x389f65e0 */
|
| 115 |
|
|
v03 = 2.5915085189e-07, /* 0x348b216c */
|
| 116 |
|
|
v04 = 4.4111031494e-10; /* 0x2ff280c2 */
|
| 117 |
|
|
|
| 118 |
|
|
#ifdef __STDC__
|
| 119 |
|
|
float __ieee754_y0f(float x)
|
| 120 |
|
|
#else
|
| 121 |
|
|
float __ieee754_y0f(x)
|
| 122 |
|
|
float x;
|
| 123 |
|
|
#endif
|
| 124 |
|
|
{
|
| 125 |
|
|
float z, s,c,ss,cc,u,v;
|
| 126 |
|
|
__int32_t hx,ix;
|
| 127 |
|
|
|
| 128 |
|
|
GET_FLOAT_WORD(hx,x);
|
| 129 |
|
|
ix = 0x7fffffff&hx;
|
| 130 |
|
|
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
|
| 131 |
|
|
if(!FLT_UWORD_IS_FINITE(ix)) return one/(x+x*x);
|
| 132 |
|
|
if(FLT_UWORD_IS_ZERO(ix)) return -one/zero;
|
| 133 |
|
|
if(hx<0) return zero/zero;
|
| 134 |
|
|
if(ix >= 0x40000000) { /* |x| >= 2.0 */
|
| 135 |
|
|
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
|
| 136 |
|
|
* where x0 = x-pi/4
|
| 137 |
|
|
* Better formula:
|
| 138 |
|
|
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
|
| 139 |
|
|
* = 1/sqrt(2) * (sin(x) + cos(x))
|
| 140 |
|
|
* sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
|
| 141 |
|
|
* = 1/sqrt(2) * (sin(x) - cos(x))
|
| 142 |
|
|
* To avoid cancellation, use
|
| 143 |
|
|
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
|
| 144 |
|
|
* to compute the worse one.
|
| 145 |
|
|
*/
|
| 146 |
|
|
s = sinf(x);
|
| 147 |
|
|
c = cosf(x);
|
| 148 |
|
|
ss = s-c;
|
| 149 |
|
|
cc = s+c;
|
| 150 |
|
|
/*
|
| 151 |
|
|
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
|
| 152 |
|
|
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
|
| 153 |
|
|
*/
|
| 154 |
|
|
if(ix<=FLT_UWORD_HALF_MAX) { /* make sure x+x not overflow */
|
| 155 |
|
|
z = -cosf(x+x);
|
| 156 |
|
|
if ((s*c)<zero) cc = z/ss;
|
| 157 |
|
|
else ss = z/cc;
|
| 158 |
|
|
}
|
| 159 |
|
|
if(ix>0x80000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x);
|
| 160 |
|
|
else {
|
| 161 |
|
|
u = pzerof(x); v = qzerof(x);
|
| 162 |
|
|
z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x);
|
| 163 |
|
|
}
|
| 164 |
|
|
return z;
|
| 165 |
|
|
}
|
| 166 |
|
|
if(ix<=0x32000000) { /* x < 2**-27 */
|
| 167 |
|
|
return(u00 + tpi*__ieee754_logf(x));
|
| 168 |
|
|
}
|
| 169 |
|
|
z = x*x;
|
| 170 |
|
|
u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
|
| 171 |
|
|
v = one+z*(v01+z*(v02+z*(v03+z*v04)));
|
| 172 |
|
|
return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x)));
|
| 173 |
|
|
}
|
| 174 |
|
|
|
| 175 |
|
|
/* The asymptotic expansions of pzero is
|
| 176 |
|
|
* 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x.
|
| 177 |
|
|
* For x >= 2, We approximate pzero by
|
| 178 |
|
|
* pzero(x) = 1 + (R/S)
|
| 179 |
|
|
* where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10
|
| 180 |
|
|
* S = 1 + pS0*s^2 + ... + pS4*s^10
|
| 181 |
|
|
* and
|
| 182 |
|
|
* | pzero(x)-1-R/S | <= 2 ** ( -60.26)
|
| 183 |
|
|
*/
|
| 184 |
|
|
#ifdef __STDC__
|
| 185 |
|
|
static const float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
|
| 186 |
|
|
#else
|
| 187 |
|
|
static float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
|
| 188 |
|
|
#endif
|
| 189 |
|
|
0.0000000000e+00, /* 0x00000000 */
|
| 190 |
|
|
-7.0312500000e-02, /* 0xbd900000 */
|
| 191 |
|
|
-8.0816707611e+00, /* 0xc1014e86 */
|
| 192 |
|
|
-2.5706311035e+02, /* 0xc3808814 */
|
| 193 |
|
|
-2.4852163086e+03, /* 0xc51b5376 */
|
| 194 |
|
|
-5.2530439453e+03, /* 0xc5a4285a */
|
| 195 |
|
|
};
|
| 196 |
|
|
#ifdef __STDC__
|
| 197 |
|
|
static const float pS8[5] = {
|
| 198 |
|
|
#else
|
| 199 |
|
|
static float pS8[5] = {
|
| 200 |
|
|
#endif
|
| 201 |
|
|
1.1653436279e+02, /* 0x42e91198 */
|
| 202 |
|
|
3.8337448730e+03, /* 0x456f9beb */
|
| 203 |
|
|
4.0597855469e+04, /* 0x471e95db */
|
| 204 |
|
|
1.1675296875e+05, /* 0x47e4087c */
|
| 205 |
|
|
4.7627726562e+04, /* 0x473a0bba */
|
| 206 |
|
|
};
|
| 207 |
|
|
#ifdef __STDC__
|
| 208 |
|
|
static const float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
|
| 209 |
|
|
#else
|
| 210 |
|
|
static float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
|
| 211 |
|
|
#endif
|
| 212 |
|
|
-1.1412546255e-11, /* 0xad48c58a */
|
| 213 |
|
|
-7.0312492549e-02, /* 0xbd8fffff */
|
| 214 |
|
|
-4.1596107483e+00, /* 0xc0851b88 */
|
| 215 |
|
|
-6.7674766541e+01, /* 0xc287597b */
|
| 216 |
|
|
-3.3123129272e+02, /* 0xc3a59d9b */
|
| 217 |
|
|
-3.4643338013e+02, /* 0xc3ad3779 */
|
| 218 |
|
|
};
|
| 219 |
|
|
#ifdef __STDC__
|
| 220 |
|
|
static const float pS5[5] = {
|
| 221 |
|
|
#else
|
| 222 |
|
|
static float pS5[5] = {
|
| 223 |
|
|
#endif
|
| 224 |
|
|
6.0753936768e+01, /* 0x42730408 */
|
| 225 |
|
|
1.0512523193e+03, /* 0x44836813 */
|
| 226 |
|
|
5.9789707031e+03, /* 0x45bad7c4 */
|
| 227 |
|
|
9.6254453125e+03, /* 0x461665c8 */
|
| 228 |
|
|
2.4060581055e+03, /* 0x451660ee */
|
| 229 |
|
|
};
|
| 230 |
|
|
|
| 231 |
|
|
#ifdef __STDC__
|
| 232 |
|
|
static const float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
|
| 233 |
|
|
#else
|
| 234 |
|
|
static float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
|
| 235 |
|
|
#endif
|
| 236 |
|
|
-2.5470459075e-09, /* 0xb12f081b */
|
| 237 |
|
|
-7.0311963558e-02, /* 0xbd8fffb8 */
|
| 238 |
|
|
-2.4090321064e+00, /* 0xc01a2d95 */
|
| 239 |
|
|
-2.1965976715e+01, /* 0xc1afba52 */
|
| 240 |
|
|
-5.8079170227e+01, /* 0xc2685112 */
|
| 241 |
|
|
-3.1447946548e+01, /* 0xc1fb9565 */
|
| 242 |
|
|
};
|
| 243 |
|
|
#ifdef __STDC__
|
| 244 |
|
|
static const float pS3[5] = {
|
| 245 |
|
|
#else
|
| 246 |
|
|
static float pS3[5] = {
|
| 247 |
|
|
#endif
|
| 248 |
|
|
3.5856033325e+01, /* 0x420f6c94 */
|
| 249 |
|
|
3.6151397705e+02, /* 0x43b4c1ca */
|
| 250 |
|
|
1.1936077881e+03, /* 0x44953373 */
|
| 251 |
|
|
1.1279968262e+03, /* 0x448cffe6 */
|
| 252 |
|
|
1.7358093262e+02, /* 0x432d94b8 */
|
| 253 |
|
|
};
|
| 254 |
|
|
|
| 255 |
|
|
#ifdef __STDC__
|
| 256 |
|
|
static const float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
|
| 257 |
|
|
#else
|
| 258 |
|
|
static float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
|
| 259 |
|
|
#endif
|
| 260 |
|
|
-8.8753431271e-08, /* 0xb3be98b7 */
|
| 261 |
|
|
-7.0303097367e-02, /* 0xbd8ffb12 */
|
| 262 |
|
|
-1.4507384300e+00, /* 0xbfb9b1cc */
|
| 263 |
|
|
-7.6356959343e+00, /* 0xc0f4579f */
|
| 264 |
|
|
-1.1193166733e+01, /* 0xc1331736 */
|
| 265 |
|
|
-3.2336456776e+00, /* 0xc04ef40d */
|
| 266 |
|
|
};
|
| 267 |
|
|
#ifdef __STDC__
|
| 268 |
|
|
static const float pS2[5] = {
|
| 269 |
|
|
#else
|
| 270 |
|
|
static float pS2[5] = {
|
| 271 |
|
|
#endif
|
| 272 |
|
|
2.2220300674e+01, /* 0x41b1c32d */
|
| 273 |
|
|
1.3620678711e+02, /* 0x430834f0 */
|
| 274 |
|
|
2.7047027588e+02, /* 0x43873c32 */
|
| 275 |
|
|
1.5387539673e+02, /* 0x4319e01a */
|
| 276 |
|
|
1.4657617569e+01, /* 0x416a859a */
|
| 277 |
|
|
};
|
| 278 |
|
|
|
| 279 |
|
|
#ifdef __STDC__
|
| 280 |
|
|
static float pzerof(float x)
|
| 281 |
|
|
#else
|
| 282 |
|
|
static float pzerof(x)
|
| 283 |
|
|
float x;
|
| 284 |
|
|
#endif
|
| 285 |
|
|
{
|
| 286 |
|
|
#ifdef __STDC__
|
| 287 |
|
|
const float *p,*q;
|
| 288 |
|
|
#else
|
| 289 |
|
|
float *p,*q;
|
| 290 |
|
|
#endif
|
| 291 |
|
|
float z,r,s;
|
| 292 |
|
|
__int32_t ix;
|
| 293 |
|
|
GET_FLOAT_WORD(ix,x);
|
| 294 |
|
|
ix &= 0x7fffffff;
|
| 295 |
|
|
if(ix>=0x41000000) {p = pR8; q= pS8;}
|
| 296 |
|
|
else if(ix>=0x40f71c58){p = pR5; q= pS5;}
|
| 297 |
|
|
else if(ix>=0x4036db68){p = pR3; q= pS3;}
|
| 298 |
|
|
else {p = pR2; q= pS2;}
|
| 299 |
|
|
z = one/(x*x);
|
| 300 |
|
|
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
|
| 301 |
|
|
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
|
| 302 |
|
|
return one+ r/s;
|
| 303 |
|
|
}
|
| 304 |
|
|
|
| 305 |
|
|
|
| 306 |
|
|
/* For x >= 8, the asymptotic expansions of qzero is
|
| 307 |
|
|
* -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
|
| 308 |
|
|
* We approximate qzero by
|
| 309 |
|
|
* qzero(x) = s*(-1.25 + (R/S))
|
| 310 |
|
|
* where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10
|
| 311 |
|
|
* S = 1 + qS0*s^2 + ... + qS5*s^12
|
| 312 |
|
|
* and
|
| 313 |
|
|
* | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22)
|
| 314 |
|
|
*/
|
| 315 |
|
|
#ifdef __STDC__
|
| 316 |
|
|
static const float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
|
| 317 |
|
|
#else
|
| 318 |
|
|
static float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */
|
| 319 |
|
|
#endif
|
| 320 |
|
|
0.0000000000e+00, /* 0x00000000 */
|
| 321 |
|
|
7.3242187500e-02, /* 0x3d960000 */
|
| 322 |
|
|
1.1768206596e+01, /* 0x413c4a93 */
|
| 323 |
|
|
5.5767340088e+02, /* 0x440b6b19 */
|
| 324 |
|
|
8.8591972656e+03, /* 0x460a6cca */
|
| 325 |
|
|
3.7014625000e+04, /* 0x471096a0 */
|
| 326 |
|
|
};
|
| 327 |
|
|
#ifdef __STDC__
|
| 328 |
|
|
static const float qS8[6] = {
|
| 329 |
|
|
#else
|
| 330 |
|
|
static float qS8[6] = {
|
| 331 |
|
|
#endif
|
| 332 |
|
|
1.6377603149e+02, /* 0x4323c6aa */
|
| 333 |
|
|
8.0983447266e+03, /* 0x45fd12c2 */
|
| 334 |
|
|
1.4253829688e+05, /* 0x480b3293 */
|
| 335 |
|
|
8.0330925000e+05, /* 0x49441ed4 */
|
| 336 |
|
|
8.4050156250e+05, /* 0x494d3359 */
|
| 337 |
|
|
-3.4389928125e+05, /* 0xc8a7eb69 */
|
| 338 |
|
|
};
|
| 339 |
|
|
|
| 340 |
|
|
#ifdef __STDC__
|
| 341 |
|
|
static const float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
|
| 342 |
|
|
#else
|
| 343 |
|
|
static float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */
|
| 344 |
|
|
#endif
|
| 345 |
|
|
1.8408595828e-11, /* 0x2da1ec79 */
|
| 346 |
|
|
7.3242180049e-02, /* 0x3d95ffff */
|
| 347 |
|
|
5.8356351852e+00, /* 0x40babd86 */
|
| 348 |
|
|
1.3511157227e+02, /* 0x43071c90 */
|
| 349 |
|
|
1.0272437744e+03, /* 0x448067cd */
|
| 350 |
|
|
1.9899779053e+03, /* 0x44f8bf4b */
|
| 351 |
|
|
};
|
| 352 |
|
|
#ifdef __STDC__
|
| 353 |
|
|
static const float qS5[6] = {
|
| 354 |
|
|
#else
|
| 355 |
|
|
static float qS5[6] = {
|
| 356 |
|
|
#endif
|
| 357 |
|
|
8.2776611328e+01, /* 0x42a58da0 */
|
| 358 |
|
|
2.0778142090e+03, /* 0x4501dd07 */
|
| 359 |
|
|
1.8847289062e+04, /* 0x46933e94 */
|
| 360 |
|
|
5.6751113281e+04, /* 0x475daf1d */
|
| 361 |
|
|
3.5976753906e+04, /* 0x470c88c1 */
|
| 362 |
|
|
-5.3543427734e+03, /* 0xc5a752be */
|
| 363 |
|
|
};
|
| 364 |
|
|
|
| 365 |
|
|
#ifdef __STDC__
|
| 366 |
|
|
static const float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
|
| 367 |
|
|
#else
|
| 368 |
|
|
static float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
|
| 369 |
|
|
#endif
|
| 370 |
|
|
4.3774099900e-09, /* 0x3196681b */
|
| 371 |
|
|
7.3241114616e-02, /* 0x3d95ff70 */
|
| 372 |
|
|
3.3442313671e+00, /* 0x405607e3 */
|
| 373 |
|
|
4.2621845245e+01, /* 0x422a7cc5 */
|
| 374 |
|
|
1.7080809021e+02, /* 0x432acedf */
|
| 375 |
|
|
1.6673394775e+02, /* 0x4326bbe4 */
|
| 376 |
|
|
};
|
| 377 |
|
|
#ifdef __STDC__
|
| 378 |
|
|
static const float qS3[6] = {
|
| 379 |
|
|
#else
|
| 380 |
|
|
static float qS3[6] = {
|
| 381 |
|
|
#endif
|
| 382 |
|
|
4.8758872986e+01, /* 0x42430916 */
|
| 383 |
|
|
7.0968920898e+02, /* 0x44316c1c */
|
| 384 |
|
|
3.7041481934e+03, /* 0x4567825f */
|
| 385 |
|
|
6.4604252930e+03, /* 0x45c9e367 */
|
| 386 |
|
|
2.5163337402e+03, /* 0x451d4557 */
|
| 387 |
|
|
-1.4924745178e+02, /* 0xc3153f59 */
|
| 388 |
|
|
};
|
| 389 |
|
|
|
| 390 |
|
|
#ifdef __STDC__
|
| 391 |
|
|
static const float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
|
| 392 |
|
|
#else
|
| 393 |
|
|
static float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */
|
| 394 |
|
|
#endif
|
| 395 |
|
|
1.5044444979e-07, /* 0x342189db */
|
| 396 |
|
|
7.3223426938e-02, /* 0x3d95f62a */
|
| 397 |
|
|
1.9981917143e+00, /* 0x3fffc4bf */
|
| 398 |
|
|
1.4495602608e+01, /* 0x4167edfd */
|
| 399 |
|
|
3.1666231155e+01, /* 0x41fd5471 */
|
| 400 |
|
|
1.6252708435e+01, /* 0x4182058c */
|
| 401 |
|
|
};
|
| 402 |
|
|
#ifdef __STDC__
|
| 403 |
|
|
static const float qS2[6] = {
|
| 404 |
|
|
#else
|
| 405 |
|
|
static float qS2[6] = {
|
| 406 |
|
|
#endif
|
| 407 |
|
|
3.0365585327e+01, /* 0x41f2ecb8 */
|
| 408 |
|
|
2.6934811401e+02, /* 0x4386ac8f */
|
| 409 |
|
|
8.4478375244e+02, /* 0x44533229 */
|
| 410 |
|
|
8.8293585205e+02, /* 0x445cbbe5 */
|
| 411 |
|
|
2.1266638184e+02, /* 0x4354aa98 */
|
| 412 |
|
|
-5.3109550476e+00, /* 0xc0a9f358 */
|
| 413 |
|
|
};
|
| 414 |
|
|
|
| 415 |
|
|
#ifdef __STDC__
|
| 416 |
|
|
static float qzerof(float x)
|
| 417 |
|
|
#else
|
| 418 |
|
|
static float qzerof(x)
|
| 419 |
|
|
float x;
|
| 420 |
|
|
#endif
|
| 421 |
|
|
{
|
| 422 |
|
|
#ifdef __STDC__
|
| 423 |
|
|
const float *p,*q;
|
| 424 |
|
|
#else
|
| 425 |
|
|
float *p,*q;
|
| 426 |
|
|
#endif
|
| 427 |
|
|
float s,r,z;
|
| 428 |
|
|
__int32_t ix;
|
| 429 |
|
|
GET_FLOAT_WORD(ix,x);
|
| 430 |
|
|
ix &= 0x7fffffff;
|
| 431 |
|
|
if(ix>=0x41000000) {p = qR8; q= qS8;}
|
| 432 |
|
|
else if(ix>=0x40f71c58){p = qR5; q= qS5;}
|
| 433 |
|
|
else if(ix>=0x4036db68){p = qR3; q= qS3;}
|
| 434 |
|
|
else {p = qR2; q= qS2;}
|
| 435 |
|
|
z = one/(x*x);
|
| 436 |
|
|
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
|
| 437 |
|
|
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
|
| 438 |
|
|
return (-(float).125 + r/s)/x;
|
| 439 |
|
|
}
|
