bug in j0f()

Steve Kargl sgk at troutmask.apl.washington.edu
Wed Dec 3 00:09:43 UTC 2014


On Tue, Dec 02, 2014 at 01:43:25PM -0800, Steve Kargl wrote:
> 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);

Exhaustive testing in the range 0x1p38 to 0x1p100 
indicated at the constant should be 0x54000000.

-- 
Steve


More information about the freebsd-numerics mailing list