Date: Wed, 5 Sep 2018 22:19:06 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> Cc: Steve Kargl <sgk@troutmask.apl.washington.edu>, freebsd-numerics@freebsd.org Subject: Re: j0 (and y0) in the range 2 <= x < (p/2)*log(2) Message-ID: <20180905221116.X1704@besplex.bde.org> In-Reply-To: <20180905201540.D1142@besplex.bde.org> References: <20180903235724.GA95333@troutmask.apl.washington.edu> <20180905201540.D1142@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On Wed, 5 Sep 2018, Bruce Evans wrote:
> On Mon, 3 Sep 2018, Steve Kargl wrote:
>
>> Anyone know where the approximations for j0 (and y0) come from?
>
> I think they are ordinary minimax rational approximations for related
> functions. As you noticed, the asymptotic expansion doesn't work below
> about x = 8 (it is off by about 10% for j0(2). But we want to use the
... j0(2)).
> single formula given by the asymptotic expansion for all the subintervals:
>
> XX /*
> XX * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
> XX * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
> XX */
> XX if(ix>0x48000000) z = (invsqrtpi*cc)/sqrt(x);
> XX else {
> XX u = pzero(x); v = qzero(x);
> XX z = invsqrtpi*(u*cc-v*ss)/sqrt(x);
> XX }
> XX return z;
> XX }
>
> where pzero(s) is nominally 1 - 9/128 s**2 + ... and qzero(s) is nominally
> -1/8 *s + ....
These polynomials are actually only part of the numerators of pzero and
qzero (s = 1/x already gives a non-polynomial, and even more divisions are
used in pzero = 1 + R/S...).
> To work, pzero and qzero must not actually be these nominal functions.
> ...
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20180905221116.X1704>