j0 (and y0) in the range 2 <= x < (p/2)*log(2)
Steve Kargl
sgk at troutmask.apl.washington.edu
Tue Sep 4 00:15:00 UTC 2018
Anyone know where the approximations for j0 (and y0) come from?
msun/src/e_j0.c states
* for x in (2,inf)
* j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
* where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
* as follow:
* cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
* = 1/sqrt(2) * (cos(x) + sin(x))
* sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
* = 1/sqrt(2) * (sin(x) - cos(x))
* (To avoid cancellation, use
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.)
p0(x) and q0(x) are divergent asymptotic series. If I extract
pzero() and qzero() from e_j0.c and compare the results against
summing truncated versions of p0(x) and q0(x), there are no obvious
connections.
Reading the documentation for the algorithms used in MPFR suggests
that x >= p/2*log(2), where p is precision of x, is required for use
of the large argument approximation for j0(x). In double precision,
p = 53, so we have x >= 18.368...
Consider x=18.4 and sum up to N = 31 in the asymptotic series:
% ./pq 30 18.4
p = 9.997932830701132e-01, q = -6.781826311540553e-03 <-- series
pp = 9.997932830701132e-01, qq = -6.781826311540509e-03 <-- pzero,qzero
ulp(p, pp) = 0.000000e+00
ulp(q, qq) = 2.550000e+01
This is almost reasonable if 25.5 ULP is acceptable in q0(x). Note
the series are computed in long double with 64 bits of precision.
Now, comparing x = 2 and summing N = 4 (best results).
% ./pq 4 2
p = 9.894313812255859e-01, q = -5.334472656250000e-02
pp = 9.862158212188928e-01, qq = -5.647769967932505e-02
ulp(p, pp) = 1.448159e+13
ulp(q, qq) = 2.257545e+14
For values of N > 4, the series start to diverge! So, how does msun
use the large argument approximation for j0(x)?
--
Steve
More information about the freebsd-numerics
mailing list