Date: Thu, 4 Dec 2014 00:08:19 +1100 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Steve Kargl <sgk@troutmask.apl.washington.edu> Cc: freebsd-numerics@freebsd.org Subject: Re: bug in j0f() Message-ID: <20141203233540.Q44095@besplex.bde.org> In-Reply-To: <20141203023207.GA99054@troutmask.apl.washington.edu> 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>
next in thread | previous in thread | raw e-mail | index | archive | help
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. > 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. Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20141203233540.Q44095>