j0 (and y0) in the range 2 <= x < (p/2)*log(2)
Bruce Evans
brde at optusnet.com.au
Wed Sep 5 12:19:11 UTC 2018
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
More information about the freebsd-numerics
mailing list