This is the mail archive of the
`gsl-discuss@sourceware.org`
mailing list for the GSL project.

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |

Other format: | [Raw text] |

*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 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

**References**:**GSL K0/K1***From:*Patrick Alken

**Re: GSL K0/K1***From:*Gerard Jungman

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |