bug in j0f()

[Steve Kargl]

Anyone object to the following patch?

Index: e_j0f.c
===================================================================
--- e_j0f.c	(revision 275211)
+++ e_j0f.c	(working copy)
@@ -62,7 +62,7 @@
 	 * 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>0x80000000) z = (invsqrtpi*cc)/sqrtf(x);
+		if(ix>0x4b800000) z = (invsqrtpi*cc)/sqrtf(x);
 		else {
 		    u = pzerof(x); v = qzerof(x);
 		    z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);


Clearly, the value 0x8000000 was a mistake in the original
translation from double to float.  In looking over other
codebases, one finds

OpenBSD
		if(ix>0x80000000U) z = (invsqrtpi*cc)/sqrtf(x);

NetBSD

#ifdef DEAD_CODE
		if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(x);
		else
#endif
		{

glibc
                if(ix>0x48000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(x);

OpenBSD is clearly wrong.  NetBSD suppresses a code optimization.
glibc simply copied the test from e_j0.c to e_j0f.c.  The glibc
choice is somewhat too small (although I haven't done any extensive
testing, yet).

prn(0x48000000); /* 0x1p17  = 1.31072000e+05 */
prn(0x4b800000); /* 0x1p24  = 1.67772160e+07 */

-- 
Steve