From: Brian Gladman <brg@gladman.plus.com>
To: gsl-discuss@sourceware.org
Subject: Re: GSL K0/K1
Date: Fri, 25 Mar 2016 22:20:00 -0000 [thread overview]
Message-ID: <56F5B964.7050401@gladman.plus.com> (raw)
In-Reply-To: <56F5998F.40902@lanl.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
next prev parent reply other threads:[~2016-03-25 22:20 UTC|newest]
Thread overview: 7+ messages / expand[flat|nested] mbox.gz Atom feed top
[not found] <56F42FED.8060500@colorado.edu>
2016-03-24 18:22 ` Patrick Alken
2016-03-25 20:03 ` Gerard Jungman
2016-03-25 21:18 ` Patrick Alken
2016-03-26 5:43 ` Pavel Holoborodko
2016-03-25 22:20 ` Brian Gladman [this message]
2016-03-25 21:24 ` error estimates (was: GSL K0/K1) Gerard Jungman
2016-03-25 21:37 ` error estimates Gerard Jungman
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=56F5B964.7050401@gladman.plus.com \
--to=brg@gladman.plus.com \
--cc=gsl-discuss@sourceware.org \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).