Date: Tue, 2 Dec 2014 13:43:25 -0800 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-numerics@freebsd.org Subject: bug in j0f() Message-ID: <20141202214325.GA94909@troutmask.apl.washington.edu>
next in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20141202214325.GA94909>