From 8c8693cf7962be7b0e53c491287aeb2fb13ae6f2 Mon Sep 17 00:00:00 2001 From: Ed Schouten Date: Thu, 12 Feb 2015 13:48:13 +0100 Subject: [PATCH] Remove [jy][01n]f(). X/Open only standardizes the double versions. --- include/openlibm_math.h | 8 - src/Make.files | 5 +- src/e_j0f.c | 344 ---------------------------------------- src/e_j1f.c | 340 --------------------------------------- src/e_jnf.c | 200 ----------------------- src/math_private.h | 6 - test/libm-test.c | 6 + 7 files changed, 8 insertions(+), 901 deletions(-) delete mode 100644 src/e_j0f.c delete mode 100644 src/e_j1f.c delete mode 100644 src/e_jnf.c diff --git a/include/openlibm_math.h b/include/openlibm_math.h index d97475c..4d53605 100644 --- a/include/openlibm_math.h +++ b/include/openlibm_math.h @@ -384,14 +384,6 @@ float fminf(float, float) __pure2; * float versions of BSD math library entry points */ #if __BSD_VISIBLE -float dremf(float, float); -float j0f(float); -float j1f(float); -float jnf(int, float); -float y0f(float); -float y1f(float); -float ynf(int, float); - /* * Float versions of reentrant version of lgamma; passes signgam back by * reference as the second argument; user must allocate space for signgam. diff --git a/src/Make.files b/src/Make.files index 1489aa8..855c4d3 100644 --- a/src/Make.files +++ b/src/Make.files @@ -1,9 +1,8 @@ $(CUR_SRCS) = common.c \ e_acos.c e_acosf.c e_acosh.c e_acoshf.c e_asin.c e_asinf.c \ e_atan2.c e_atan2f.c e_atanh.c e_atanhf.c e_cosh.c e_coshf.c e_exp.c \ - e_expf.c e_fmod.c e_fmodf.c \ - e_hypot.c e_hypotf.c e_j0.c e_j0f.c e_j1.c e_j1f.c \ - e_jn.c e_jnf.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \ + e_expf.c e_fmod.c e_fmodf.c e_hypot.c e_hypotf.c e_j0.c e_j1.c \ + e_jn.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \ e_lgammal.c e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c \ e_pow.c e_powf.c e_remainder.c e_remainderf.c \ e_rem_pio2.c e_rem_pio2f.c \ diff --git a/src/e_j0f.c b/src/e_j0f.c deleted file mode 100644 index 74d753b..0000000 --- a/src/e_j0f.c +++ /dev/null @@ -1,344 +0,0 @@ -/* e_j0f.c -- float version of e_j0.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#include - -#include "cdefs-compat.h" -//__FBSDID("$FreeBSD: src/lib/msun/src/e_j0f.c,v 1.8 2008/02/22 02:30:35 das Exp $"); - -#include - -#include "math_private.h" - -static float pzerof(float), qzerof(float); - -static const float -huge = 1e30, -one = 1.0, -invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */ -tpi = 6.3661974669e-01, /* 0x3f22f983 */ - /* R0/S0 on [0, 2.00] */ -R02 = 1.5625000000e-02, /* 0x3c800000 */ -R03 = -1.8997929874e-04, /* 0xb947352e */ -R04 = 1.8295404516e-06, /* 0x35f58e88 */ -R05 = -4.6183270541e-09, /* 0xb19eaf3c */ -S01 = 1.5619102865e-02, /* 0x3c7fe744 */ -S02 = 1.1692678527e-04, /* 0x38f53697 */ -S03 = 5.1354652442e-07, /* 0x3509daa6 */ -S04 = 1.1661400734e-09; /* 0x30a045e8 */ - -static const float zero = 0.0; - -DLLEXPORT float -__ieee754_j0f(float x) -{ - float z, s,c,ss,cc,r,u,v; - int32_t hx,ix; - - GET_FLOAT_WORD(hx,x); - ix = hx&0x7fffffff; - if(ix>=0x7f800000) return one/(x*x); - x = fabsf(x); - 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)0x80000000) z = (invsqrtpi*cc)/sqrtf(x); - else { - u = pzerof(x); v = qzerof(x); - z = invsqrtpi*(u*cc-v*ss)/sqrtf(x); - } - return z; - } - 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; - } - } - z = x*x; - r = z*(R02+z*(R03+z*(R04+z*R05))); - s = one+z*(S01+z*(S02+z*(S03+z*S04))); - if(ix < 0x3F800000) { /* |x| < 1.00 */ - return one + z*((float)-0.25+(r/s)); - } else { - u = (float)0.5*x; - return((one+u)*(one-u)+z*(r/s)); - } -} - -static const float -u00 = -7.3804296553e-02, /* 0xbd9726b5 */ -u01 = 1.7666645348e-01, /* 0x3e34e80d */ -u02 = -1.3818567619e-02, /* 0xbc626746 */ -u03 = 3.4745343146e-04, /* 0x39b62a69 */ -u04 = -3.8140706238e-06, /* 0xb67ff53c */ -u05 = 1.9559013964e-08, /* 0x32a802ba */ -u06 = -3.9820518410e-11, /* 0xae2f21eb */ -v01 = 1.2730483897e-02, /* 0x3c509385 */ -v02 = 7.6006865129e-05, /* 0x389f65e0 */ -v03 = 2.5915085189e-07, /* 0x348b216c */ -v04 = 4.4111031494e-10; /* 0x2ff280c2 */ - -DLLEXPORT 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; - /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(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 >= 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. - */ - s = sinf(x); - c = cosf(x); - 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)0x80000000) z = (invsqrtpi*ss)/sqrtf(x); - else { - u = pzerof(x); v = qzerof(x); - z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); - } - return z; - } - if(ix<=0x32000000) { /* x < 2**-27 */ - return(u00 + tpi*__ieee754_logf(x)); - } - z = x*x; - u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); - v = one+z*(v01+z*(v02+z*(v03+z*v04))); - return(u/v + tpi*(__ieee754_j0f(x)*__ieee754_logf(x))); -} - -/* 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) - * where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10 - * S = 1 + pS0*s^2 + ... + pS4*s^10 - * and - * | pzero(x)-1-R/S | <= 2 ** ( -60.26) - */ -static const float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0000000000e+00, /* 0x00000000 */ - -7.0312500000e-02, /* 0xbd900000 */ - -8.0816707611e+00, /* 0xc1014e86 */ - -2.5706311035e+02, /* 0xc3808814 */ - -2.4852163086e+03, /* 0xc51b5376 */ - -5.2530439453e+03, /* 0xc5a4285a */ -}; -static const float pS8[5] = { - 1.1653436279e+02, /* 0x42e91198 */ - 3.8337448730e+03, /* 0x456f9beb */ - 4.0597855469e+04, /* 0x471e95db */ - 1.1675296875e+05, /* 0x47e4087c */ - 4.7627726562e+04, /* 0x473a0bba */ -}; -static const float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - -1.1412546255e-11, /* 0xad48c58a */ - -7.0312492549e-02, /* 0xbd8fffff */ - -4.1596107483e+00, /* 0xc0851b88 */ - -6.7674766541e+01, /* 0xc287597b */ - -3.3123129272e+02, /* 0xc3a59d9b */ - -3.4643338013e+02, /* 0xc3ad3779 */ -}; -static const float pS5[5] = { - 6.0753936768e+01, /* 0x42730408 */ - 1.0512523193e+03, /* 0x44836813 */ - 5.9789707031e+03, /* 0x45bad7c4 */ - 9.6254453125e+03, /* 0x461665c8 */ - 2.4060581055e+03, /* 0x451660ee */ -}; - -static const float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - -2.5470459075e-09, /* 0xb12f081b */ - -7.0311963558e-02, /* 0xbd8fffb8 */ - -2.4090321064e+00, /* 0xc01a2d95 */ - -2.1965976715e+01, /* 0xc1afba52 */ - -5.8079170227e+01, /* 0xc2685112 */ - -3.1447946548e+01, /* 0xc1fb9565 */ -}; -static const float pS3[5] = { - 3.5856033325e+01, /* 0x420f6c94 */ - 3.6151397705e+02, /* 0x43b4c1ca */ - 1.1936077881e+03, /* 0x44953373 */ - 1.1279968262e+03, /* 0x448cffe6 */ - 1.7358093262e+02, /* 0x432d94b8 */ -}; - -static const float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - -8.8753431271e-08, /* 0xb3be98b7 */ - -7.0303097367e-02, /* 0xbd8ffb12 */ - -1.4507384300e+00, /* 0xbfb9b1cc */ - -7.6356959343e+00, /* 0xc0f4579f */ - -1.1193166733e+01, /* 0xc1331736 */ - -3.2336456776e+00, /* 0xc04ef40d */ -}; -static const float pS2[5] = { - 2.2220300674e+01, /* 0x41b1c32d */ - 1.3620678711e+02, /* 0x430834f0 */ - 2.7047027588e+02, /* 0x43873c32 */ - 1.5387539673e+02, /* 0x4319e01a */ - 1.4657617569e+01, /* 0x416a859a */ -}; - - /* Note: This function is only called for ix>=0x40000000 (see above) */ - static float pzerof(float x) -{ - const float *p,*q; - float z,r,s; - int32_t ix; - GET_FLOAT_WORD(ix,x); - ix &= 0x7fffffff; - assert(ix>=0x40000000 && ix<=0x48000000); - if(ix>=0x41000000) {p = pR8; q= pS8;} - else if(ix>=0x40f71c58){p = pR5; q= pS5;} - else if(ix>=0x4036db68){p = pR3; q= pS3;} - else {p = pR2; q= pS2;} - 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])))); - return one+ r/s; -} - - -/* 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)) - * where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10 - * S = 1 + qS0*s^2 + ... + qS5*s^12 - * and - * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) - */ -static const float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0000000000e+00, /* 0x00000000 */ - 7.3242187500e-02, /* 0x3d960000 */ - 1.1768206596e+01, /* 0x413c4a93 */ - 5.5767340088e+02, /* 0x440b6b19 */ - 8.8591972656e+03, /* 0x460a6cca */ - 3.7014625000e+04, /* 0x471096a0 */ -}; -static const float qS8[6] = { - 1.6377603149e+02, /* 0x4323c6aa */ - 8.0983447266e+03, /* 0x45fd12c2 */ - 1.4253829688e+05, /* 0x480b3293 */ - 8.0330925000e+05, /* 0x49441ed4 */ - 8.4050156250e+05, /* 0x494d3359 */ - -3.4389928125e+05, /* 0xc8a7eb69 */ -}; - -static const float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - 1.8408595828e-11, /* 0x2da1ec79 */ - 7.3242180049e-02, /* 0x3d95ffff */ - 5.8356351852e+00, /* 0x40babd86 */ - 1.3511157227e+02, /* 0x43071c90 */ - 1.0272437744e+03, /* 0x448067cd */ - 1.9899779053e+03, /* 0x44f8bf4b */ -}; -static const float qS5[6] = { - 8.2776611328e+01, /* 0x42a58da0 */ - 2.0778142090e+03, /* 0x4501dd07 */ - 1.8847289062e+04, /* 0x46933e94 */ - 5.6751113281e+04, /* 0x475daf1d */ - 3.5976753906e+04, /* 0x470c88c1 */ - -5.3543427734e+03, /* 0xc5a752be */ -}; - -static const float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - 4.3774099900e-09, /* 0x3196681b */ - 7.3241114616e-02, /* 0x3d95ff70 */ - 3.3442313671e+00, /* 0x405607e3 */ - 4.2621845245e+01, /* 0x422a7cc5 */ - 1.7080809021e+02, /* 0x432acedf */ - 1.6673394775e+02, /* 0x4326bbe4 */ -}; -static const float qS3[6] = { - 4.8758872986e+01, /* 0x42430916 */ - 7.0968920898e+02, /* 0x44316c1c */ - 3.7041481934e+03, /* 0x4567825f */ - 6.4604252930e+03, /* 0x45c9e367 */ - 2.5163337402e+03, /* 0x451d4557 */ - -1.4924745178e+02, /* 0xc3153f59 */ -}; - -static const float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - 1.5044444979e-07, /* 0x342189db */ - 7.3223426938e-02, /* 0x3d95f62a */ - 1.9981917143e+00, /* 0x3fffc4bf */ - 1.4495602608e+01, /* 0x4167edfd */ - 3.1666231155e+01, /* 0x41fd5471 */ - 1.6252708435e+01, /* 0x4182058c */ -}; -static const float qS2[6] = { - 3.0365585327e+01, /* 0x41f2ecb8 */ - 2.6934811401e+02, /* 0x4386ac8f */ - 8.4478375244e+02, /* 0x44533229 */ - 8.8293585205e+02, /* 0x445cbbe5 */ - 2.1266638184e+02, /* 0x4354aa98 */ - -5.3109550476e+00, /* 0xc0a9f358 */ -}; - - /* Note: This function is only called for ix>=0x40000000 (see above) */ - static float qzerof(float x) -{ - const float *p,*q; - float s,r,z; - int32_t ix; - GET_FLOAT_WORD(ix,x); - ix &= 0x7fffffff; - assert(ix>=0x40000000 && ix<=0x48000000); - if(ix>=0x41000000) {p = qR8; q= qS8;} - else if(ix>=0x40f71c58){p = qR5; q= qS5;} - else if(ix>=0x4036db68){p = qR3; q= qS3;} - else {p = qR2; q= qS2;} - 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]))))); - return (-(float).125 + r/s)/x; -} diff --git a/src/e_j1f.c b/src/e_j1f.c deleted file mode 100644 index afed84f..0000000 --- a/src/e_j1f.c +++ /dev/null @@ -1,340 +0,0 @@ -/* e_j1f.c -- float version of e_j1.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#include - -#include "cdefs-compat.h" -//__FBSDID("$FreeBSD: src/lib/msun/src/e_j1f.c,v 1.8 2008/02/22 02:30:35 das Exp $"); - -#include - -#include "math_private.h" - -static float ponef(float), qonef(float); - -static const float -huge = 1e30, -one = 1.0, -invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */ -tpi = 6.3661974669e-01, /* 0x3f22f983 */ - /* R0/S0 on [0,2] */ -r00 = -6.2500000000e-02, /* 0xbd800000 */ -r01 = 1.4070566976e-03, /* 0x3ab86cfd */ -r02 = -1.5995563444e-05, /* 0xb7862e36 */ -r03 = 4.9672799207e-08, /* 0x335557d2 */ -s01 = 1.9153760746e-02, /* 0x3c9ce859 */ -s02 = 1.8594678841e-04, /* 0x3942fab6 */ -s03 = 1.1771846857e-06, /* 0x359dffc2 */ -s04 = 5.0463624390e-09, /* 0x31ad6446 */ -s05 = 1.2354227016e-11; /* 0x2d59567e */ - -static const float zero = 0.0; - -DLLEXPORT float -__ieee754_j1f(float x) -{ - float z, s,c,ss,cc,r,u,v,y; - int32_t hx,ix; - - GET_FLOAT_WORD(hx,x); - ix = hx&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; - else return z; - } - if(ix<0x32000000) { /* |x|<2**-27 */ - if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */ - } - 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); -} - -static const float U0[5] = { - -1.9605709612e-01, /* 0xbe48c331 */ - 5.0443872809e-02, /* 0x3d4e9e3c */ - -1.9125689287e-03, /* 0xbafaaf2a */ - 2.3525259166e-05, /* 0x37c5581c */ - -9.1909917899e-08, /* 0xb3c56003 */ -}; -static const float V0[5] = { - 1.9916731864e-02, /* 0x3ca3286a */ - 2.0255257550e-04, /* 0x3954644b */ - 1.3560879779e-06, /* 0x35b602d4 */ - 6.2274145840e-09, /* 0x31d5f8eb */ - 1.6655924903e-11, /* 0x2d9281cf */ -}; - -DLLEXPORT float -__ieee754_y1f(float x) -{ - float z, s,c,ss,cc,u,v; - int32_t hx,ix; - - 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 >= 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(-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*(__ieee754_j1f(x)*__ieee754_logf(x)-one/x)); -} - -/* For x >= 8, the asymptotic expansions of pone is - * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. - * We approximate pone by - * pone(x) = 1 + (R/S) - * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 - * S = 1 + ps0*s^2 + ... + ps4*s^10 - * and - * | pone(x)-1-R/S | <= 2 ** ( -60.06) - */ - -static const float pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0000000000e+00, /* 0x00000000 */ - 1.1718750000e-01, /* 0x3df00000 */ - 1.3239480972e+01, /* 0x4153d4ea */ - 4.1205184937e+02, /* 0x43ce06a3 */ - 3.8747453613e+03, /* 0x45722bed */ - 7.9144794922e+03, /* 0x45f753d6 */ -}; -static const float ps8[5] = { - 1.1420736694e+02, /* 0x42e46a2c */ - 3.6509309082e+03, /* 0x45642ee5 */ - 3.6956207031e+04, /* 0x47105c35 */ - 9.7602796875e+04, /* 0x47bea166 */ - 3.0804271484e+04, /* 0x46f0a88b */ -}; - -static const float pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - 1.3199052094e-11, /* 0x2d68333f */ - 1.1718749255e-01, /* 0x3defffff */ - 6.8027510643e+00, /* 0x40d9b023 */ - 1.0830818176e+02, /* 0x42d89dca */ - 5.1763616943e+02, /* 0x440168b7 */ - 5.2871520996e+02, /* 0x44042dc6 */ -}; -static const float ps5[5] = { - 5.9280597687e+01, /* 0x426d1f55 */ - 9.9140142822e+02, /* 0x4477d9b1 */ - 5.3532670898e+03, /* 0x45a74a23 */ - 7.8446904297e+03, /* 0x45f52586 */ - 1.5040468750e+03, /* 0x44bc0180 */ -}; - -static const float pr3[6] = { - 3.0250391081e-09, /* 0x314fe10d */ - 1.1718686670e-01, /* 0x3defffab */ - 3.9329774380e+00, /* 0x407bb5e7 */ - 3.5119403839e+01, /* 0x420c7a45 */ - 9.1055007935e+01, /* 0x42b61c2a */ - 4.8559066772e+01, /* 0x42423c7c */ -}; -static const float ps3[5] = { - 3.4791309357e+01, /* 0x420b2a4d */ - 3.3676245117e+02, /* 0x43a86198 */ - 1.0468714600e+03, /* 0x4482dbe3 */ - 8.9081134033e+02, /* 0x445eb3ed */ - 1.0378793335e+02, /* 0x42cf936c */ -}; - -static const float pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - 1.0771083225e-07, /* 0x33e74ea8 */ - 1.1717621982e-01, /* 0x3deffa16 */ - 2.3685150146e+00, /* 0x401795c0 */ - 1.2242610931e+01, /* 0x4143e1bc */ - 1.7693971634e+01, /* 0x418d8d41 */ - 5.0735230446e+00, /* 0x40a25a4d */ -}; -static const float ps2[5] = { - 2.1436485291e+01, /* 0x41ab7dec */ - 1.2529022980e+02, /* 0x42fa9499 */ - 2.3227647400e+02, /* 0x436846c7 */ - 1.1767937469e+02, /* 0x42eb5bd7 */ - 8.3646392822e+00, /* 0x4105d590 */ -}; - - /* Note: This function is only called for ix>=0x40000000 (see above) */ - static float ponef(float x) -{ - const float *p,*q; - float z,r,s; - int32_t ix; - GET_FLOAT_WORD(ix,x); - ix &= 0x7fffffff; - assert(ix>=0x40000000 && ix<=0x48000000); - if(ix>=0x41000000) {p = pr8; q= ps8;} - else if(ix>=0x40f71c58){p = pr5; q= ps5;} - else if(ix>=0x4036db68){p = pr3; q= ps3;} - else {p = pr2; q= ps2;} - 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])))); - return one+ r/s; -} - - -/* For x >= 8, the asymptotic expansions of qone is - * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. - * We approximate pone by - * qone(x) = s*(0.375 + (R/S)) - * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10 - * S = 1 + qs1*s^2 + ... + qs6*s^12 - * and - * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13) - */ - -static const float qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0000000000e+00, /* 0x00000000 */ - -1.0253906250e-01, /* 0xbdd20000 */ - -1.6271753311e+01, /* 0xc1822c8d */ - -7.5960174561e+02, /* 0xc43de683 */ - -1.1849806641e+04, /* 0xc639273a */ - -4.8438511719e+04, /* 0xc73d3683 */ -}; -static const float qs8[6] = { - 1.6139537048e+02, /* 0x43216537 */ - 7.8253862305e+03, /* 0x45f48b17 */ - 1.3387534375e+05, /* 0x4802bcd6 */ - 7.1965775000e+05, /* 0x492fb29c */ - 6.6660125000e+05, /* 0x4922be94 */ - -2.9449025000e+05, /* 0xc88fcb48 */ -}; - -static const float qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - -2.0897993405e-11, /* 0xadb7d219 */ - -1.0253904760e-01, /* 0xbdd1fffe */ - -8.0564479828e+00, /* 0xc100e736 */ - -1.8366960144e+02, /* 0xc337ab6b */ - -1.3731937256e+03, /* 0xc4aba633 */ - -2.6124443359e+03, /* 0xc523471c */ -}; -static const float qs5[6] = { - 8.1276550293e+01, /* 0x42a28d98 */ - 1.9917987061e+03, /* 0x44f8f98f */ - 1.7468484375e+04, /* 0x468878f8 */ - 4.9851425781e+04, /* 0x4742bb6d */ - 2.7948074219e+04, /* 0x46da5826 */ - -4.7191835938e+03, /* 0xc5937978 */ -}; - -static const float qr3[6] = { - -5.0783124372e-09, /* 0xb1ae7d4f */ - -1.0253783315e-01, /* 0xbdd1ff5b */ - -4.6101160049e+00, /* 0xc0938612 */ - -5.7847221375e+01, /* 0xc267638e */ - -2.2824453735e+02, /* 0xc3643e9a */ - -2.1921012878e+02, /* 0xc35b35cb */ -}; -static const float qs3[6] = { - 4.7665153503e+01, /* 0x423ea91e */ - 6.7386511230e+02, /* 0x4428775e */ - 3.3801528320e+03, /* 0x45534272 */ - 5.5477290039e+03, /* 0x45ad5dd5 */ - 1.9031191406e+03, /* 0x44ede3d0 */ - -1.3520118713e+02, /* 0xc3073381 */ -}; - -static const float qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - -1.7838172539e-07, /* 0xb43f8932 */ - -1.0251704603e-01, /* 0xbdd1f475 */ - -2.7522056103e+00, /* 0xc0302423 */ - -1.9663616180e+01, /* 0xc19d4f16 */ - -4.2325313568e+01, /* 0xc2294d1f */ - -2.1371921539e+01, /* 0xc1aaf9b2 */ -}; -static const float qs2[6] = { - 2.9533363342e+01, /* 0x41ec4454 */ - 2.5298155212e+02, /* 0x437cfb47 */ - 7.5750280762e+02, /* 0x443d602e */ - 7.3939318848e+02, /* 0x4438d92a */ - 1.5594900513e+02, /* 0x431bf2f2 */ - -4.9594988823e+00, /* 0xc09eb437 */ -}; - - /* Note: This function is only called for ix>=0x40000000 (see above) */ - static float qonef(float x) -{ - const float *p,*q; - float s,r,z; - int32_t ix; - GET_FLOAT_WORD(ix,x); - ix &= 0x7fffffff; - assert(ix>=0x40000000 && ix<=0x48000000); - if(ix>=0x40200000) {p = qr8; q= qs8;} - else if(ix>=0x40f71c58){p = qr5; q= qs5;} - else if(ix>=0x4036db68){p = qr3; q= qs3;} - else {p = qr2; q= qs2;} - 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]))))); - return ((float).375 + r/s)/x; -} diff --git a/src/e_jnf.c b/src/e_jnf.c deleted file mode 100644 index a2b83ee..0000000 --- a/src/e_jnf.c +++ /dev/null @@ -1,200 +0,0 @@ -/* e_jnf.c -- float version of e_jn.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#include "cdefs-compat.h" -//__FBSDID("$FreeBSD: src/lib/msun/src/e_jnf.c,v 1.11 2010/11/13 10:54:10 uqs Exp $"); - -#include - -#include "math_private.h" - -static const float -two = 2.0000000000e+00, /* 0x40000000 */ -one = 1.0000000000e+00; /* 0x3F800000 */ - -static const float zero = 0.0000000000e+00; - -DLLEXPORT float -__ieee754_jnf(int n, float x) -{ - int32_t i,hx,ix, sgn; - float a, b, temp, di; - float z, w; - - /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) - * Thus, J(-n,x) = J(n,-x) - */ - GET_FLOAT_WORD(hx,x); - ix = 0x7fffffff&hx; - /* if J(n,NaN) is NaN */ - if(ix>0x7f800000) return x+x; - if(n<0){ - n = -n; - x = -x; - hx ^= 0x80000000; - } - if(n==0) return(__ieee754_j0f(x)); - if(n==1) return(__ieee754_j1f(x)); - sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */ - x = fabsf(x); - if(ix==0||ix>=0x7f800000) /* if x is 0 or inf */ - b = zero; - else if((float)n<=x) { - /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ - a = __ieee754_j0f(x); - b = __ieee754_j1f(x); - for(i=1;i33) /* underflow */ - b = zero; - else { - temp = x*(float)0.5; b = temp; - for (a=one,i=2;i<=n;i++) { - a *= (float)i; /* a = n! */ - b *= temp; /* b = (x/2)^n */ - } - b = b/a; - } - } else { - /* use backward recurrence */ - /* x x^2 x^2 - * J(n,x)/J(n-1,x) = ---- ------ ------ ..... - * 2n - 2(n+1) - 2(n+2) - * - * 1 1 1 - * (for large x) = ---- ------ ------ ..... - * 2n 2(n+1) 2(n+2) - * -- - ------ - ------ - - * x x x - * - * Let w = 2n/x and h=2/x, then the above quotient - * is equal to the continued fraction: - * 1 - * = ----------------------- - * 1 - * w - ----------------- - * 1 - * w+h - --------- - * w+2h - ... - * - * To determine how many terms needed, let - * Q(0) = w, Q(1) = w(w+h) - 1, - * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), - * When Q(k) > 1e4 good for single - * When Q(k) > 1e9 good for double - * When Q(k) > 1e17 good for quadruple - */ - /* determine k */ - float t,v; - float q0,q1,h,tmp; int32_t k,m; - w = (n+n)/(float)x; h = (float)2.0/(float)x; - q0 = w; z = w+h; q1 = w*z - (float)1.0; k=1; - while(q1<(float)1.0e9) { - k += 1; z += h; - tmp = z*q1 - q0; - q0 = q1; - q1 = tmp; - } - m = n+n; - for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t); - a = t; - b = one; - /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) - * Hence, if n*(log(2n/x)) > ... - * single 8.8722839355e+01 - * double 7.09782712893383973096e+02 - * long double 1.1356523406294143949491931077970765006170e+04 - * then recurrent value may overflow and the result is - * likely underflow to zero - */ - tmp = n; - v = two/x; - tmp = tmp*__ieee754_logf(fabsf(v*tmp)); - if(tmp<(float)8.8721679688e+01) { - for(i=n-1,di=(float)(i+i);i>0;i--){ - temp = b; - b *= di; - b = b/x - a; - a = temp; - di -= two; - } - } else { - for(i=n-1,di=(float)(i+i);i>0;i--){ - temp = b; - b *= di; - b = b/x - a; - a = temp; - di -= two; - /* scale b to avoid spurious overflow */ - if(b>(float)1e10) { - a /= b; - t /= b; - b = one; - } - } - } - z = __ieee754_j0f(x); - w = __ieee754_j1f(x); - if (fabsf(z) >= fabsf(w)) - b = (t*z/b); - else - b = (t*w/a); - } - } - if(sgn==1) return -b; else return b; -} - -DLLEXPORT float -__ieee754_ynf(int n, float x) -{ - int32_t i,hx,ix,ib; - int32_t sign; - float a, b, temp; - - GET_FLOAT_WORD(hx,x); - ix = 0x7fffffff&hx; - /* if Y(n,NaN) is NaN */ - if(ix>0x7f800000) return x+x; - if(ix==0) return -one/zero; - if(hx<0) return zero/zero; - sign = 1; - if(n<0){ - n = -n; - sign = 1 - ((n&1)<<1); - } - if(n==0) return(__ieee754_y0f(x)); - if(n==1) return(sign*__ieee754_y1f(x)); - if(ix==0x7f800000) return zero; - - a = __ieee754_y0f(x); - b = __ieee754_y1f(x); - /* quit if b is -inf */ - GET_FLOAT_WORD(ib,b); - for(i=1;i0) return b; else return -b; -} diff --git a/src/math_private.h b/src/math_private.h index 7bfd76b..50745f8 100644 --- a/src/math_private.h +++ b/src/math_private.h @@ -308,12 +308,6 @@ irint(double x) #define __ieee754_log2f log2f #define __ieee754_sinhf sinhf #define __ieee754_hypotf hypotf -#define __ieee754_j0f j0f -#define __ieee754_j1f j1f -#define __ieee754_y0f y0f -#define __ieee754_y1f y1f -#define __ieee754_jnf jnf -#define __ieee754_ynf ynf #define __ieee754_remainderf remainderf /* fdlibm kernel function */ diff --git a/test/libm-test.c b/test/libm-test.c index b5903ce..99f4718 100644 --- a/test/libm-test.c +++ b/test/libm-test.c @@ -2918,6 +2918,7 @@ isnormal_test (void) print_max_error ("isnormal", 0, 0); } +#ifdef TEST_DOUBLE static void j0_test (void) { @@ -3055,6 +3056,7 @@ jn_test (void) print_max_error ("jn", DELTAjn, 0); } +#endif static void @@ -4121,6 +4123,7 @@ trunc_test (void) print_max_error ("trunc", 0, 0); } +#ifdef TEST_DOUBLE static void y0_test (void) { @@ -4256,6 +4259,7 @@ yn_test (void) print_max_error ("yn", DELTAyn, 0); } +#endif @@ -4534,6 +4538,7 @@ main (int argc, char **argv) ctanh_test (); #endif +#ifdef TEST_DOUBLE /* Bessel functions: */ j0_test (); j1_test (); @@ -4541,6 +4546,7 @@ main (int argc, char **argv) y0_test (); y1_test (); yn_test (); +#endif if (output_ulps) fclose (ulps_file);