Date: Wed, 3 Dec 2014 07:49:44 -0800 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: freebsd-numerics@freebsd.org Subject: Re: bug in j0f() Message-ID: <20141203154944.GA1976@troutmask.apl.washington.edu> In-Reply-To: <20141203233540.Q44095@besplex.bde.org> References: <20141202214325.GA94909@troutmask.apl.washington.edu> <20141203000941.GA98467@troutmask.apl.washington.edu> <20141203002908.GA98589@troutmask.apl.washington.edu> <20141203023207.GA99054@troutmask.apl.washington.edu> <20141203233540.Q44095@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20141203154944.GA1976>