diff options
Diffstat (limited to 'lib/msun/src/e_j1f.c')
| -rw-r--r-- | lib/msun/src/e_j1f.c | 43 |
1 files changed, 25 insertions, 18 deletions
diff --git a/lib/msun/src/e_j1f.c b/lib/msun/src/e_j1f.c index 88e2d833d393..ec7f38101b53 100644 --- a/lib/msun/src/e_j1f.c +++ b/lib/msun/src/e_j1f.c @@ -16,10 +16,16 @@ #include <sys/cdefs.h> __FBSDID("$FreeBSD$"); +/* + * See e_j1.c for complete comments. + */ + #include "math.h" #include "math_private.h" -static float ponef(float), qonef(float); +static __inline float ponef(float), qonef(float); + +static const volatile float vone = 1, vzero = 0; static const float huge = 1e30, @@ -63,7 +69,7 @@ __ieee754_j1f(float x) * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) */ - if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(y); + if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(y); /* |x|>2**49 */ else { u = ponef(y); v = qonef(y); z = invsqrtpi*(u*cc-v*ss)/sqrtf(y); @@ -71,7 +77,7 @@ __ieee754_j1f(float x) if(hx<0) return -z; else return z; } - if(ix<0x32000000) { /* |x|<2**-27 */ + if(ix<0x39000000) { /* |x|<2**-13 */ if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */ } z = x*x; @@ -104,10 +110,9 @@ __ieee754_y1f(float x) GET_FLOAT_WORD(hx,x); ix = 0x7fffffff&hx; - /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ - if(ix>=0x7f800000) return one/(x+x*x); - if(ix==0) return -one/zero; - if(hx<0) return zero/zero; + if(ix>=0x7f800000) return vone/(x+x*x); + if(ix==0) return -one/vzero; + if(hx<0) return vzero/vzero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ s = sinf(x); c = cosf(x); @@ -129,14 +134,14 @@ __ieee754_y1f(float x) * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) * to compute the worse one. */ - if(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x); + if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */ else { u = ponef(x); v = qonef(x); z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); } return z; } - if(ix<=0x24800000) { /* x < 2**-54 */ + if(ix<=0x33000000) { /* x < 2**-25 */ return(-tpi/x); } z = x*x; @@ -219,7 +224,8 @@ static const float ps2[5] = { 8.3646392822e+00, /* 0x4105d590 */ }; - static float ponef(float x) +static __inline float +ponef(float x) { const float *p,*q; float z,r,s; @@ -227,9 +233,9 @@ static const float ps2[5] = { GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; if(ix>=0x41000000) {p = pr8; q= ps8;} - else if(ix>=0x40f71c58){p = pr5; q= ps5;} - else if(ix>=0x4036db68){p = pr3; q= ps3;} - else if(ix>=0x40000000){p = pr2; q= ps2;} + else if(ix>=0x409173eb){p = pr5; q= ps5;} + else if(ix>=0x4036d917){p = pr3; q= ps3;} + else {p = pr2; q= ps2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); @@ -315,17 +321,18 @@ static const float qs2[6] = { -4.9594988823e+00, /* 0xc09eb437 */ }; - static float qonef(float x) +static __inline float +qonef(float x) { const float *p,*q; float s,r,z; int32_t ix; GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; - if(ix>=0x40200000) {p = qr8; q= qs8;} - else if(ix>=0x40f71c58){p = qr5; q= qs5;} - else if(ix>=0x4036db68){p = qr3; q= qs3;} - else if(ix>=0x40000000){p = qr2; q= qs2;} + if(ix>=0x41000000) {p = qr8; q= qs8;} + else if(ix>=0x409173eb){p = qr5; q= qs5;} + else if(ix>=0x4036d917){p = qr3; q= qs3;} + else {p = qr2; q= qs2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); |