bug in j0f()
Steve Kargl
sgk at troutmask.apl.washington.edu
Wed Dec 3 02:32:08 UTC 2014
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?
> > >
> > > 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.
> >
>
> Note, a similar wrong condition exists within y0f(). I have
> not tested y0f(), but propose making a similar change in y0f()
> as well.
>
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.
--
Steve
More information about the freebsd-numerics
mailing list