diff options
Diffstat (limited to 'sysdeps/ieee754/flt-32/e_j0f.c')
| -rw-r--r-- | sysdeps/ieee754/flt-32/e_j0f.c | 200 |
1 files changed, 47 insertions, 153 deletions
diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c index 8c499e614e..d2da43f929 100644 --- a/sysdeps/ieee754/flt-32/e_j0f.c +++ b/sysdeps/ieee754/flt-32/e_j0f.c @@ -13,29 +13,17 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: e_j0f.c,v 1.4 1995/05/10 20:45:25 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static float pzerof(float), qzerof(float); -#else -static float pzerof(), qzerof(); -#endif -#ifdef __STDC__ static const float -#else -static float -#endif -huge = 1e30, +huge = 1e30, one = 1.0, invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */ tpi = 6.3661974669e-01, /* 0x3f22f983 */ - /* R0/S0 on [0, 2.00] */ + /* R0/S0 on [0, 2.00] */ R02 = 1.5625000000e-02, /* 0x3c800000 */ R03 = -1.8997929874e-04, /* 0xb947352e */ R04 = 1.8295404516e-06, /* 0x35f58e88 */ @@ -45,18 +33,10 @@ S02 = 1.1692678527e-04, /* 0x38f53697 */ S03 = 5.1354652442e-07, /* 0x3509daa6 */ S04 = 1.1661400734e-09; /* 0x30a045e8 */ -#ifdef __STDC__ static const float zero = 0.0; -#else -static float zero = 0.0; -#endif -#ifdef __STDC__ - float __ieee754_j0f(float x) -#else - float __ieee754_j0f(x) - float x; -#endif +float +__ieee754_j0f(float x) { float z, s,c,ss,cc,r,u,v; int32_t hx,ix; @@ -72,7 +52,7 @@ static float zero = 0.0; if(ix<0x7f000000) { /* make sure x+x not overflow */ z = -__cosf(x+x); if ((s*c)<zero) cc = z/ss; - else ss = z/cc; + else ss = z/cc; } /* * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) @@ -87,8 +67,8 @@ static float zero = 0.0; } if(ix<0x39000000) { /* |x| < 2**-13 */ if(huge+x>one) { /* raise inexact if x != 0 */ - if(ix<0x32000000) return one; /* |x|<2**-27 */ - else return one - (float)0.25*x*x; + if(ix<0x32000000) return one; /* |x|<2**-27 */ + else return one - (float)0.25*x*x; } } z = x*x; @@ -101,12 +81,9 @@ static float zero = 0.0; return((one+u)*(one-u)+z*(r/s)); } } +strong_alias (__ieee754_j0f, __j0f_finite) -#ifdef __STDC__ static const float -#else -static float -#endif u00 = -7.3804296553e-02, /* 0xbd9726b5 */ u01 = 1.7666645348e-01, /* 0x3e34e80d */ u02 = -1.3818567619e-02, /* 0xbc626746 */ @@ -119,52 +96,48 @@ v02 = 7.6006865129e-05, /* 0x389f65e0 */ v03 = 2.5915085189e-07, /* 0x348b216c */ v04 = 4.4111031494e-10; /* 0x2ff280c2 */ -#ifdef __STDC__ - float __ieee754_y0f(float x) -#else - float __ieee754_y0f(x) - float x; -#endif +float +__ieee754_y0f(float x) { float z, s,c,ss,cc,u,v; int32_t hx,ix; GET_FLOAT_WORD(hx,x); - ix = 0x7fffffff&hx; + ix = 0x7fffffff&hx; /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf. */ if(ix>=0x7f800000) return one/(x+x*x); - if(ix==0) return -HUGE_VALF+x; /* -inf and overflow exception. */ - if(hx<0) return zero/(zero*x); - if(ix >= 0x40000000) { /* |x| >= 2.0 */ - /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) - * where x0 = x-pi/4 - * Better formula: - * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) - * = 1/sqrt(2) * (sin(x) + cos(x)) - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one. - */ + if(ix==0) return -HUGE_VALF+x; /* -inf and overflow exception. */ + if(hx<0) return zero/(zero*x); + if(ix >= 0x40000000) { /* |x| >= 2.0 */ + /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) + * where x0 = x-pi/4 + * Better formula: + * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) + * = 1/sqrt(2) * (sin(x) + cos(x)) + * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) + * = 1/sqrt(2) * (sin(x) - cos(x)) + * To avoid cancellation, use + * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + * to compute the worse one. + */ __sincosf (x, &s, &c); - ss = s-c; - cc = s+c; + ss = s-c; + cc = s+c; /* * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) */ - if(ix<0x7f000000) { /* make sure x+x not overflow */ - z = -__cosf(x+x); - if ((s*c)<zero) cc = z/ss; - else ss = z/cc; - } - if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x); - else { - u = pzerof(x); v = qzerof(x); - z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x); - } - return z; + if(ix<0x7f000000) { /* make sure x+x not overflow */ + z = -__cosf(x+x); + if ((s*c)<zero) cc = z/ss; + else ss = z/cc; + } + if(ix>0x48000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x); + else { + u = pzerof(x); v = qzerof(x); + z = invsqrtpi*(u*ss+v*cc)/__ieee754_sqrtf(x); + } + return z; } if(ix<=0x32000000) { /* x < 2**-27 */ return(u00 + tpi*__ieee754_logf(x)); @@ -174,21 +147,18 @@ v04 = 4.4111031494e-10; /* 0x2ff280c2 */ v = one+z*(v01+z*(v02+z*(v03+z*v04))); return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x))); } +strong_alias (__ieee754_y0f, __y0f_finite) /* The asymptotic expansions of pzero is * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. * For x >= 2, We approximate pzero by - * pzero(x) = 1 + (R/S) + * pzero(x) = 1 + (R/S) * where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10 - * S = 1 + pS0*s^2 + ... + pS4*s^10 + * S = 1 + pS0*s^2 + ... + pS4*s^10 * and * | pzero(x)-1-R/S | <= 2 ** ( -60.26) */ -#ifdef __STDC__ static const float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#else -static float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#endif 0.0000000000e+00, /* 0x00000000 */ -7.0312500000e-02, /* 0xbd900000 */ -8.0816707611e+00, /* 0xc1014e86 */ @@ -196,22 +166,14 @@ static float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -2.4852163086e+03, /* 0xc51b5376 */ -5.2530439453e+03, /* 0xc5a4285a */ }; -#ifdef __STDC__ static const float pS8[5] = { -#else -static float pS8[5] = { -#endif 1.1653436279e+02, /* 0x42e91198 */ 3.8337448730e+03, /* 0x456f9beb */ 4.0597855469e+04, /* 0x471e95db */ 1.1675296875e+05, /* 0x47e4087c */ 4.7627726562e+04, /* 0x473a0bba */ }; -#ifdef __STDC__ static const float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#else -static float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#endif -1.1412546255e-11, /* 0xad48c58a */ -7.0312492549e-02, /* 0xbd8fffff */ -4.1596107483e+00, /* 0xc0851b88 */ @@ -219,11 +181,7 @@ static float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -3.3123129272e+02, /* 0xc3a59d9b */ -3.4643338013e+02, /* 0xc3ad3779 */ }; -#ifdef __STDC__ static const float pS5[5] = { -#else -static float pS5[5] = { -#endif 6.0753936768e+01, /* 0x42730408 */ 1.0512523193e+03, /* 0x44836813 */ 5.9789707031e+03, /* 0x45bad7c4 */ @@ -231,11 +189,7 @@ static float pS5[5] = { 2.4060581055e+03, /* 0x451660ee */ }; -#ifdef __STDC__ static const float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#else -static float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#endif -2.5470459075e-09, /* 0xb12f081b */ -7.0311963558e-02, /* 0xbd8fffb8 */ -2.4090321064e+00, /* 0xc01a2d95 */ @@ -243,11 +197,7 @@ static float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -5.8079170227e+01, /* 0xc2685112 */ -3.1447946548e+01, /* 0xc1fb9565 */ }; -#ifdef __STDC__ static const float pS3[5] = { -#else -static float pS3[5] = { -#endif 3.5856033325e+01, /* 0x420f6c94 */ 3.6151397705e+02, /* 0x43b4c1ca */ 1.1936077881e+03, /* 0x44953373 */ @@ -255,11 +205,7 @@ static float pS3[5] = { 1.7358093262e+02, /* 0x432d94b8 */ }; -#ifdef __STDC__ static const float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#else -static float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#endif -8.8753431271e-08, /* 0xb3be98b7 */ -7.0303097367e-02, /* 0xbd8ffb12 */ -1.4507384300e+00, /* 0xbfb9b1cc */ @@ -267,11 +213,7 @@ static float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -1.1193166733e+01, /* 0xc1331736 */ -3.2336456776e+00, /* 0xc04ef40d */ }; -#ifdef __STDC__ static const float pS2[5] = { -#else -static float pS2[5] = { -#endif 2.2220300674e+01, /* 0x41b1c32d */ 1.3620678711e+02, /* 0x430834f0 */ 2.7047027588e+02, /* 0x43873c32 */ @@ -279,18 +221,10 @@ static float pS2[5] = { 1.4657617569e+01, /* 0x416a859a */ }; -#ifdef __STDC__ - static float pzerof(float x) -#else - static float pzerof(x) - float x; -#endif +static float +pzerof(float x) { -#ifdef __STDC__ const float *p,*q; -#else - float *p,*q; -#endif float z,r,s; int32_t ix; GET_FLOAT_WORD(ix,x); @@ -309,17 +243,13 @@ static float pS2[5] = { /* For x >= 8, the asymptotic expansions of qzero is * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. * We approximate pzero by - * qzero(x) = s*(-1.25 + (R/S)) + * qzero(x) = s*(-1.25 + (R/S)) * where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10 - * S = 1 + qS0*s^2 + ... + qS5*s^12 + * S = 1 + qS0*s^2 + ... + qS5*s^12 * and * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) */ -#ifdef __STDC__ static const float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#else -static float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ -#endif 0.0000000000e+00, /* 0x00000000 */ 7.3242187500e-02, /* 0x3d960000 */ 1.1768206596e+01, /* 0x413c4a93 */ @@ -327,11 +257,7 @@ static float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ 8.8591972656e+03, /* 0x460a6cca */ 3.7014625000e+04, /* 0x471096a0 */ }; -#ifdef __STDC__ static const float qS8[6] = { -#else -static float qS8[6] = { -#endif 1.6377603149e+02, /* 0x4323c6aa */ 8.0983447266e+03, /* 0x45fd12c2 */ 1.4253829688e+05, /* 0x480b3293 */ @@ -340,11 +266,7 @@ static float qS8[6] = { -3.4389928125e+05, /* 0xc8a7eb69 */ }; -#ifdef __STDC__ static const float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#else -static float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ -#endif 1.8408595828e-11, /* 0x2da1ec79 */ 7.3242180049e-02, /* 0x3d95ffff */ 5.8356351852e+00, /* 0x40babd86 */ @@ -352,11 +274,7 @@ static float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ 1.0272437744e+03, /* 0x448067cd */ 1.9899779053e+03, /* 0x44f8bf4b */ }; -#ifdef __STDC__ static const float qS5[6] = { -#else -static float qS5[6] = { -#endif 8.2776611328e+01, /* 0x42a58da0 */ 2.0778142090e+03, /* 0x4501dd07 */ 1.8847289062e+04, /* 0x46933e94 */ @@ -365,11 +283,7 @@ static float qS5[6] = { -5.3543427734e+03, /* 0xc5a752be */ }; -#ifdef __STDC__ static const float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#else -static float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ -#endif 4.3774099900e-09, /* 0x3196681b */ 7.3241114616e-02, /* 0x3d95ff70 */ 3.3442313671e+00, /* 0x405607e3 */ @@ -377,11 +291,7 @@ static float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ 1.7080809021e+02, /* 0x432acedf */ 1.6673394775e+02, /* 0x4326bbe4 */ }; -#ifdef __STDC__ static const float qS3[6] = { -#else -static float qS3[6] = { -#endif 4.8758872986e+01, /* 0x42430916 */ 7.0968920898e+02, /* 0x44316c1c */ 3.7041481934e+03, /* 0x4567825f */ @@ -390,11 +300,7 @@ static float qS3[6] = { -1.4924745178e+02, /* 0xc3153f59 */ }; -#ifdef __STDC__ static const float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#else -static float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ -#endif 1.5044444979e-07, /* 0x342189db */ 7.3223426938e-02, /* 0x3d95f62a */ 1.9981917143e+00, /* 0x3fffc4bf */ @@ -402,11 +308,7 @@ static float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ 3.1666231155e+01, /* 0x41fd5471 */ 1.6252708435e+01, /* 0x4182058c */ }; -#ifdef __STDC__ static const float qS2[6] = { -#else -static float qS2[6] = { -#endif 3.0365585327e+01, /* 0x41f2ecb8 */ 2.6934811401e+02, /* 0x4386ac8f */ 8.4478375244e+02, /* 0x44533229 */ @@ -415,18 +317,10 @@ static float qS2[6] = { -5.3109550476e+00, /* 0xc0a9f358 */ }; -#ifdef __STDC__ - static float qzerof(float x) -#else - static float qzerof(x) - float x; -#endif +static float +qzerof(float x) { -#ifdef __STDC__ const float *p,*q; -#else - float *p,*q; -#endif float s,r,z; int32_t ix; GET_FLOAT_WORD(ix,x); |