j0 (and y0) in the range 2 <= x < (p/2)*log(2)

[Bruce Evans]

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