bug in j0f()

Steve Kargl sgk at troutmask.apl.washington.edu
Wed Dec 3 15:49:46 UTC 2014


On Thu, Dec 04, 2014 at 12:08:19AM +1100, Bruce Evans wrote:
> On Tue, 2 Dec 2014, Steve Kargl wrote:
> 
> > On Tue, Dec 02, 2014 at 04:29:08PM -0800, Steve Kargl wrote:
> >> On Tue, Dec 02, 2014 at 04:09:41PM -0800, Steve Kargl wrote:
> >>> On Tue, Dec 02, 2014 at 01:43:25PM -0800, Steve Kargl wrote:
> >>>> Anyone object to the following patch?
> 
> OK (with 0x54000000).
> 
> >>>> 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.
> 
> My tests agree.  Tested on amd64 and i386.
> 
> >> Note, a similar wrong condition exists within y0f().  I have
> >> not tested y0f(), but propose making a similar change in y0f()
> >> as well.
> 
> Not so exhaustive testing gave 0x54800000 on amd64.
> 

Thanks for confirming the values and suggestion the
above for y0f().

> > While I'm monologuing, I might as well point out that the
> > rational approximations in j0f (and y0f and most likely
> > j1f and y1f) are over-specified.  I suspect that the
> > polynomials in the rational approximation can be reduced
> > by one or two terms.
> 
> Also, the cutoffs of 2**-13 and 2**-27 are the same for both precisions,
> thus likely to be wrong for float precision.
> 

I haven't gotten that far into analyzing the code. I'll
probably get there in a few months. :-)

-- 
Steve


More information about the freebsd-numerics mailing list