diff options
Diffstat (limited to 'sysdeps/ieee754')
| -rw-r--r-- | sysdeps/ieee754/flt-32/e_jnf.c | 19 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-96/e_j0l.c | 18 |
2 files changed, 17 insertions, 20 deletions
diff --git a/sysdeps/ieee754/flt-32/e_jnf.c b/sysdeps/ieee754/flt-32/e_jnf.c index 9e5279c..34c4d95 100644 --- a/sysdeps/ieee754/flt-32/e_jnf.c +++ b/sysdeps/ieee754/flt-32/e_jnf.c @@ -25,7 +25,6 @@ static const float #else static float #endif -invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */ two = 2.0000000000e+00, /* 0x40000000 */ one = 1.0000000000e+00; /* 0x3F800000 */ @@ -62,7 +61,7 @@ static float zero = 0.0000000000e+00; 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 */ + 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) */ @@ -90,11 +89,11 @@ static float zero = 0.0000000000e+00; } } else { /* use backward recurrence */ - /* x x^2 x^2 + /* x x^2 x^2 * J(n,x)/J(n-1,x) = ---- ------ ------ ..... * 2n - 2(n+1) - 2(n+2) * - * 1 1 1 + * 1 1 1 * (for large x) = ---- ------ ------ ..... * 2n 2(n+1) 2(n+2) * -- - ------ - ------ - @@ -107,7 +106,7 @@ static float zero = 0.0000000000e+00; * 1 * w - ----------------- * 1 - * w+h - --------- + * w+h - --------- * w+2h - ... * * To determine how many terms needed, let @@ -144,15 +143,15 @@ static float zero = 0.0000000000e+00; 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--){ + 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--){ + for(i=n-1,di=(float)(i+i);i>0;i--){ temp = b; b *= di; b = b/x - a; @@ -164,9 +163,9 @@ static float zero = 0.0000000000e+00; t /= b; b = one; } - } + } } - b = (t*__ieee754_j0f(x)/b); + b = (t*__ieee754_j0f(x)/b); } } if(sgn==1) return -b; else return b; diff --git a/sysdeps/ieee754/ldbl-96/e_j0l.c b/sysdeps/ieee754/ldbl-96/e_j0l.c index 79e13be..e8966f6 100644 --- a/sysdeps/ieee754/ldbl-96/e_j0l.c +++ b/sysdeps/ieee754/ldbl-96/e_j0l.c @@ -20,16 +20,16 @@ * for x in (0,2) * j0(x) = 1 - z/4 + z^2*R0/S0, where z = x*x; * for x in (2,inf) - * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) - * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) + * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) + * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) * as follow: * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) * = 1/sqrt(2) * (cos(x) + sin(x)) * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4) * = 1/sqrt(2) * (sin(x) - cos(x)) - * (To avoid cancellation, use + * (To avoid cancellation, use * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one.) + * to compute the worse one.) * * 3 Special cases * j0(nan)= nan @@ -47,8 +47,8 @@ * Note: For tiny x, U/V = u0 and j0(x)~1, hence * y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27) * 2. For x>=2. - * y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0)) - * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) + * y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0)) + * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) * by the method mentioned above. * 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0. */ @@ -71,8 +71,6 @@ static long double one = 1.0L, invsqrtpi = 5.6418958354775628694807945156077258584405e-1L, tpi = 6.3661977236758134307553505349005744813784e-1L, - j0z1 = 2.40482555769577276862163187932650662155139L, - j0z2 = 5.520078110286310649596604112813027425221865L, /* J0(x) = 1 - x^2 / 4 + x^4 R0(x^2) / S0(x^2) 0 <= x <= 2 @@ -275,7 +273,7 @@ __ieee754_y0l (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 + s^2 R(s^2) / S(s^2) + * pzero(x) = 1 + s^2 R(s^2) / S(s^2) */ #ifdef __STDC__ static const long double pR8[7] = { @@ -450,7 +448,7 @@ pzero (x) /* For x >= 8, the asymptotic expansions of qzero is * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. * We approximate qzero by - * qzero(x) = s*(-.125 + R(s^2) / S(s^2)) + * qzero(x) = s*(-.125 + R(s^2) / S(s^2)) */ #ifdef __STDC__ static const long double qR8[7] = { |