This is the mail archive of the
mailing list for the GSL project.
Re: GSL K0/K1
- From: Brian Gladman <brg at gladman dot plus dot com>
- To: gsl-discuss at sourceware dot org
- Date: Fri, 25 Mar 2016 22:19:16 +0000
- Subject: Re: GSL K0/K1
- Authentication-results: sourceware.org; auth=none
- References: <56F42FED dot 8060500 at colorado dot edu> <56F4304C dot 3000907 at colorado dot edu> <56F5998F dot 40902 at lanl dot gov>
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
> 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
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
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.