Date: Wed, 5 Sep 2018 08:21:04 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: freebsd-numerics@freebsd.org Subject: Re: j0 (and y0) in the range 2 <= x < (p/2)*log(2) Message-ID: <20180905152104.GA26453@troutmask.apl.washington.edu> In-Reply-To: <20180905201540.D1142@besplex.bde.org> References: <20180903235724.GA95333@troutmask.apl.washington.edu> <20180905201540.D1142@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20180905152104.GA26453>