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

[Bruce Evans]

On Fri, 7 Sep 2018, Steve Kargl wrote:

> On Fri, Sep 07, 2018 at 10:54:06AM -0700, Steve Kargl wrote:
>>
>> I'll probably try a different approach tomorrow.  j0(x),
>> cc, and ss are smoothly varying functions of x.  u and v
>> appear to be slowly varying (in small intervals), so we can
>> have
>>
>> j0(x)   = sqrt(1/pi) * (cc(x)  *u-ss(x)  *v) / sqrt(x).
>> j0(x+e) = sqrt(1/pi) * (cc(x+e)*u-ss(x+e)*v) / sqrt(x+e).
>>
>> where e is some small perturbution.  Two equation and two
>> unknowns is easily solved.  Simply need to do this across
>> the interval to generate u(x) and v(x).
>
> I may have have a better approach!
>
> j0(x) = sqrt(1/pi) * (cc(x) * u - ss(x) * v) / sqrt(x)
> y0(x) = sqrt(1/pi) * (cc(x) * u + ss(x) * v) / sqrt(x)
>
> j0(x) + y0(x) =   2 * sqrt(1/pi) * cc(x) * u
> j0(x) - y0(x) = - 2 * sqrt(1/pi) * ss(x) * v

Divided by sqrt(x) on the RHS.

> So, we have
>
> u = (j0 + y0) / (2 * sqrt(1/pi) * cc)
> v = (y0 - j0) / (s * sqrt(1/pi) * ss)

Multiplied by sqrt(x) on the RHS.  Also, change s back to 2.

> Thus, we can find u = r/s or u = 1 + r/s and v = r/s.

Good idea.

I checked that the u and v calculated by this work right in pari (they
give a precision almost as high as the precision used in pari when
they are used to implement j0 using libm's formula in pari) (57 decimal
digits with 100-digit precision in pari, even though I used the limit
formula A&S 9.1.2 for y0 because I couldn't find y0 in pari (use j_nu
and y_nu with nu = 1e-50 to get the 57 decimal digits of accuracy).

   (I just noticed that libm only supports integral nu, presumably because
   its recursion method only works for integers.  Current versions of
   pari document not very good accuracy for small args and/or small nu,
   and I'm using an 11 year old version, but nu = 1e-50 works perfectly.)

However, the u and v calculated by this have no resemblance to libm's
pzero and qzero (with the latter checked by pari to give a precision
of about 16 decimal digits when used to implement j0).  pzero / u - 1
is between -0.5 and -3.1.  qzero / v - 1 is between -1.2 and -1.0.  So
libm doesn't seem to be doing anything better than fixing a vague
approximation to the correct u and then solving for p.

Better use the simple product method.

Bruce