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

[Steve Kargl]

On Wed, Sep 05, 2018 at 10:06:29PM +1000, 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
> single formula given by the asymptotic expansion for all the subintervals:

I've scoured the literature and web for methods of computing 
Bessel functions.  These functions are important to my real
work.  I have not found any paper, webpage, documentation, etc.
that describes what "the related functions" are.

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

If |x| > p/2*log(2) with p the precision of x, P0(x) and Q0(x) can be
used directly.  It is the 2 <= |x| < p/2*log(2) range that is the issue,
but your last paragraph below may be what I have been missing. 

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

Obviously, I had not thought about scaling P0(x) and Q0(x) with a
common factor to force agreement between j0(x) and its approximation.

-- 
Steve