Re: GSL K0/K1

[Brian Gladman <brg@gladman.plus.com>]

On 25/03/2016 20:03, Gerard Jungman wrote:
> 
> As you say, this is the first I have seen of this.
> The results look very good.
> 
> I briefly compared the code to the old version.
> It looks like you changed the domains for some
> of the expansions as well as the coefficients
> and some of the fractional transformations.
> So that's a lot of differences to track.
> 
> Just glancing at the differences, I cannot tell
> why it is so much better. Do you know exactly
> what happened?
> 
> The main problem seemed to be for x in (1,2).
> The big break in the 1/x plots at x=2 suggests
> that there was a problem in the coefficients
> in that domain. Chebyshev expansions should
> not behave that way; they should provide
> essentially uniform error over the domain.
> 
> Was it some kind of typo in the tables?
> That's what it looks like to me. But
> maybe not... in fact I don't see how
> that would have happened.

The coefficients were not correctly rounded and this caused larger than
expected relative errors.  There is also a logarithmic singularity at x
= 0 that was computed in a way that led to some loss of precision
through cancellation.

> As for the method and error estimates for K0 and K1, I am
> trying to understand myself what is happening there. Most
> of it is straightforward; the only problem is the domain
> x <= 2.0. I wrote this code about 20 years ago and
> haven't looked at it since. The assumption was that
> the original SLATEC was basically best possible, so
> it did not get a lot of scrutiny.
> 
> I looked at the SLATEC source. In the code I refer to "bk0.f",
> although it seems to actually be called "besk0.f", so maybe
> there have been further changes to SLATEC since then?
> 
> Anyway, SLATEC uses this method where it evaluates I0 for
> some reason. So I do the same. I have no idea why they do
> this. 

For low x the modified Bessel function of the second kind for order 0 is
evaluated using the definition:

  K0(x) = -(log(x/2) + gamma)I0(x) + polynomial(x^2)

where I0(x) is the modified function of the first kind. In GSL the
polynomial component is developed as the chebyshev series expansion
under discussion.

Maybe there is an original paper somewhere that
> explains it. Clearly it is not necessary, and it's
> much cleaner to avoid this dependency as well.
> If the expansion was correct in the old code,
> then the error must have been creeping in
> through the evaluation of I0.

The dependence of Kn(x) on In(x) for low x is not avoidable because the
functions of the second kind have a logarithmic sigularity at x = 0.

And, as you surmise, additional loss of precision was creeping in
because the expansion coefficients for I0(x) were also incorrectly rounded.

   Brian Gladman