X-Git-Url: http://nsz.repo.hu/git/?p=musl;a=blobdiff_plain;f=src%2Fmath%2Fj1f.c;h=6abde34927da2998ac0dfa94a651182612b98582;hp=0323ec7828bb97244ea1d34fa35f23016fe0f960;hb=e216951f509b71da193da2fc63e25b998740d58b;hpb=b69f695acedd4ce2798ef9ea28d834ceccc789bd diff --git a/src/math/j1f.c b/src/math/j1f.c index 0323ec78..6abde349 100644 --- a/src/math/j1f.c +++ b/src/math/j1f.c @@ -18,11 +18,38 @@ static float ponef(float), qonef(float); static const float -huge = 1e30, -one = 1.0, invsqrtpi = 5.6418961287e-01, /* 0x3f106ebb */ -tpi = 6.3661974669e-01, /* 0x3f22f983 */ +tpi = 6.3661974669e-01; /* 0x3f22f983 */ + +static float common(uint32_t ix, float x, int y1, int sign) +{ + double z,s,c,ss,cc; + + s = sinf(x); + if (y1) + s = -s; + c = cosf(x); + cc = s-c; + if (ix < 0x7f000000) { + ss = -s-c; + z = cosf(2*x); + if (s*c > 0) + cc = z/ss; + else + ss = z/cc; + if (ix < 0x58800000) { + if (y1) + ss = -ss; + cc = ponef(x)*cc-qonef(x)*ss; + } + } + if (sign) + cc = -cc; + return invsqrtpi*cc/sqrtf(x); +} + /* R0/S0 on [0,2] */ +static const float r00 = -6.2500000000e-02, /* 0xbd800000 */ r01 = 1.4070566976e-03, /* 0x3ab86cfd */ r02 = -1.5995563444e-05, /* 0xb7862e36 */ @@ -33,55 +60,28 @@ s03 = 1.1771846857e-06, /* 0x359dffc2 */ s04 = 5.0463624390e-09, /* 0x31ad6446 */ s05 = 1.2354227016e-11; /* 0x2d59567e */ -static const float zero = 0.0; - float j1f(float x) { - float z,s,c,ss,cc,r,u,v,y; - int32_t hx,ix; + float z,r,s; + uint32_t ix; + int sign; - GET_FLOAT_WORD(hx, x); - ix = hx & 0x7fffffff; + GET_FLOAT_WORD(ix, x); + sign = ix>>31; + ix &= 0x7fffffff; if (ix >= 0x7f800000) - return one/x; - y = fabsf(x); - if (ix >= 0x40000000) { /* |x| >= 2.0 */ - s = sinf(y); - c = cosf(y); - ss = -s-c; - cc = s-c; - if (ix < 0x7f000000) { /* make sure y+y not overflow */ - z = cosf(y+y); - if (s*c > zero) - cc = z/ss; - else - ss = z/cc; - } - /* - * 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); - else { - u = ponef(y); - v = qonef(y); - z = invsqrtpi*(u*cc-v*ss)/sqrtf(y); - } - if (hx < 0) - return -z; - return z; - } - if (ix < 0x32000000) { /* |x| < 2**-27 */ + return 1/(x*x); + if (ix >= 0x40000000) /* |x| >= 2 */ + return common(ix, fabsf(x), 0, sign); + if (ix >= 0x32000000) { /* |x| >= 2**-27 */ + z = x*x; + r = z*(r00+z*(r01+z*(r02+z*r03))); + s = 1+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); + z = 0.5f + r/s; + } else /* raise inexact if x!=0 */ - if (huge+x > one) - return (float)0.5*x; - } - z = x*x; - r = z*(r00+z*(r01+z*(r02+z*r03))); - s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); - r *= x; - return x*(float)0.5 + r/s; + z = 0.5f + x; + return z*x; } static const float U0[5] = { @@ -101,56 +101,24 @@ static const float V0[5] = { float y1f(float x) { - float z,s,c,ss,cc,u,v; - int32_t hx,ix; + float z,u,v; + uint32_t ix; - GET_FLOAT_WORD(hx, x); - ix = 0x7fffffff & hx; - /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ + GET_FLOAT_WORD(ix, x); + if ((ix & 0x7fffffff) == 0) + return -1/0.0f; + if (ix>>31) + return 0/0.0f; if (ix >= 0x7f800000) - return one/(x+x*x); - if (ix == 0) - return -one/zero; - if (hx < 0) - return zero/zero; - if (ix >= 0x40000000) { /* |x| >= 2.0 */ - s = sinf(x); - c = cosf(x); - ss = -s-c; - cc = s-c; - if (ix < 0x7f000000) { /* make sure x+x not overflow */ - z = cosf(x+x); - if (s*c > zero) - cc = z/ss; - else - ss = z/cc; - } - /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0)) - * where x0 = x-3pi/4 - * Better formula: - * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = -1/sqrt(2) * (cos(x) + sin(x)) - * To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one. - */ - if (ix > 0x48000000) - z = (invsqrtpi*ss)/sqrtf(x); - else { - u = ponef(x); - v = qonef(x); - z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); - } - return z; - } - if (ix <= 0x24800000) /* x < 2**-54 */ + return 1/x; + if (ix >= 0x40000000) /* |x| >= 2.0 */ + return common(ix,x,1,0); + if (ix < 0x32000000) /* x < 2**-27 */ return -tpi/x; z = x*x; u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); - v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); - return x*(u/v) + tpi*(j1f(x)*logf(x)-one/x); + v = 1.0f+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); + return x*(u/v) + tpi*(j1f(x)*logf(x)-1.0f/x); } /* For x >= 8, the asymptotic expansions of pone is @@ -230,19 +198,19 @@ static const float ps2[5] = { static float ponef(float x) { const float *p,*q; - float z,r,s; - int32_t ix; + float_t z,r,s; + uint32_t ix; 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;} - z = one/(x*x); + else /*ix >= 0x40000000*/ {p = pr2; q = ps2;} + z = 1.0f/(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])))); - return one + r/s; + s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); + return 1.0f + r/s; } /* For x >= 8, the asymptotic expansions of qone is @@ -326,17 +294,17 @@ static const float qs2[6] = { static float qonef(float x) { const float *p,*q; - float s,r,z; - int32_t ix; + float_t s,r,z; + uint32_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;} - z = one/(x*x); + else /*ix >= 0x40000000*/ {p = qr2; q = qs2;} + z = 1.0f/(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]))))); - return ((float).375 + r/s)/x; + s = 1.0f+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); + return (.375f + r/s)/x; }