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

[Steve Kargl]

On Fri, Sep 07, 2018 at 02:02:20AM +1000, Bruce Evans wrote:
> On Wed, 5 Sep 2018, Steve Kargl wrote:
> 
> > On Thu, Sep 06, 2018 at 04:09:05AM +1000, Bruce Evans wrote:
> >> On Wed, 5 Sep 2018, Steve Kargl wrote:
> >>> 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.
> >>
> >> They are just the functions in the asymptotic expansion with errors corrected
> >> as I discussed.
> >
> > And as I noted, there is no documentation stating the approximations
> > pzero(x) and qzero(x) aren't approximations for the asymptotic series
> > P0(x) and Q0(x).  If you are correct, then pzero(x) and qzero(x) are
> > approximations to fudge*P0(x) and fudge*Q0(x).  What fudge is and how
> > it is determined is not documented.
> 
> The documentation (comments in the source code) is indeed deficient.  It
> is too verbose about routine general methods but says little about non-
> routine things.
> 
> I now think that no single fudge factor works for both P0 and Q0.  P0
> and Q0 have very different characteristics.  Merely choosing the
> truncation point for the asymptotic expansion makes them very different.
> One thing that can go wrong is that if you truncate to more than about
> 5 terms, a zero point for one or both of P0 and Q0 ends up in the
> subinterval being handled.  The fudge factor would need to be nearly
> infinite to compensate for the zero, and since P0 and Q0 won't have
> the zero at the same place, the compensation for one would destroy the
> other.
> 
> The infinities show up in attempts to calculate the fudge factor near the
> zero.
> 
> When P0 is truncated after the s^14 term, then on [2, 2.857] P0 (1/x)
> is strictly increasing from -5.61 to 0.95, with the zero in the middle
> giving singularities, but pzero is strictly increasing from 1 - 0.013 to
> 1 - 0.007.  However, when P0 is truncated after the s^6 term, it is
> similar to pzero (strictly increasing from 1 - 0.019 to 1 - 0.009).
> pzero is apparently based on the latter P0.
> 
> To find different fudge factors, I would try starting with a common one
> and then iteratively adjust to different ones.
> 
> j0 has a zero in the middle of at least the first subinterval, and the
> relative error should be minimized for this.  I think the choice of
> subintervals is related to this.  With only one zero per subinterval,
> the relative error near it can be minimized without messing up the
> relative error near other zeros.  The adjustment for this is absolutely
> tiny, so it is easier to arrange that it doesn't mess up the absolute
> error too much.
> 

I'm beginning to like my use of the product formula.  It
allows one to remove the zero from an interval.  I did some
code spelunking and found

#ifndef lint
static char sccsid[] = "@(#)j0.c 4.2 (Berkeley) 8/21/85";
#endif not lint

Found at http://www.retro11.de/ouxr/43bsd/usr/src/usr.lib/libm/j0.c.html

The algorithm is completely different than what we have in libm.
 
I won't have time to look at this until the weekend.

-- 
Steve