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

[Bruce Evans]

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
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 + ....

To work, pzero and qzero must not actually be these nominal functions.

A non-tricky implementation would use :

 	z = pzero(x) / qzero(x);

where pzero and qzero depend on the subinterval like now and have all the
other terms folded into them.  We already do this for |x| < 2 (except we
spell the rational function r/s and don't fold all of the terms directly
into r and s).  This might be the best implementation, and I would also
try using smaller intervals with polynomial approximations.

However, we use the asympotic formula for all cases above x = 2.  It is
off by at least 10% unless pzero and qzero are adjusted.  But 10% is not
much, and adjusting only pzero and qzero apparently works to compensate
for this.  To use the same formula for all cases, the complicated factors
cc, ss and sqrt(x) must not be folded into pzero and qzero; instead these
factors supply some of the errors which must be compensated for.

To find the adjustments, I would try applying the same scale factor to
the nominal pzero and qzero.  E.g., for j0, the error factor is
invsqrtpi*(u*cc-v*ss)/sqrt(x) / j0(x).  Divide the nominal pzero and qzero
by this.  Apply normal minimax methods to the modfided pzero and qzero.

Bruce