Fix a 20+ year bug by using an appropriate constant for

the transition from one asymptotic approximation to another
for the zeroth order Bessel and Neumann functions.

Reviewed by:	bde
This commit is contained in:
Steve Kargl 2014-12-04 15:57:58 +00:00
parent 26f0f92fa2
commit 668986107d

View File

@ -62,7 +62,7 @@ __ieee754_j0f(float x)
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) * 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) * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/ */
if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(x); if(ix>0x54000000) z = (invsqrtpi*cc)/sqrtf(x);
else { else {
u = pzerof(x); v = qzerof(x); u = pzerof(x); v = qzerof(x);
z = invsqrtpi*(u*cc-v*ss)/sqrtf(x); z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);
@ -136,7 +136,7 @@ __ieee754_y0f(float x)
if ((s*c)<zero) cc = z/ss; if ((s*c)<zero) cc = z/ss;
else ss = z/cc; else ss = z/cc;
} }
if(ix>0x80000000) z = (invsqrtpi*ss)/sqrtf(x); if(ix>0x54800000) z = (invsqrtpi*ss)/sqrtf(x);
else { else {
u = pzerof(x); v = qzerof(x); u = pzerof(x); v = qzerof(x);
z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);