From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 23041 invoked by alias); 25 Mar 2016 22:20:03 -0000 Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org Received: (qmail 23029 invoked by uid 89); 25 Mar 2016 22:20:02 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-1.6 required=5.0 tests=BAYES_00,KAM_LAZY_DOMAIN_SECURITY,RCVD_IN_DNSWL_LOW autolearn=no version=3.3.2 spammy=cancellation, HTo:U*gsl-discuss, paper, domains X-HELO: avasout01.plus.net Received: from avasout01.plus.net (HELO avasout01.plus.net) (84.93.230.227) by sourceware.org (qpsmtpd/0.93/v0.84-503-g423c35a) with (AES128-SHA encrypted) ESMTPS; Fri, 25 Mar 2016 22:20:00 +0000 Received: from [192.168.1.122] ([80.229.7.237]) by avasout01 with smtp id aNKN1s00256qjnU01NKcky; Fri, 25 Mar 2016 22:19:57 +0000 X-CM-Score: 0.00 X-CNFS-Analysis: v=2.1 cv=bsGxfxui c=1 sm=1 tr=0 a=PIjflcjPJHybs4Jfvhhq9w==:117 a=PIjflcjPJHybs4Jfvhhq9w==:17 a=L9H7d07YOLsA:10 a=9cW_t1CCXrUA:10 a=s5jvgZ67dGcA:10 a=IkcTkHD0fZMA:10 a=XhmNM7dcSVINbxAkTD0A:9 a=QEXdDO2ut3YA:10 X-AUTH: gladman+brg@:2500 Subject: Re: GSL K0/K1 To: gsl-discuss@sourceware.org References: <56F42FED.8060500@colorado.edu> <56F4304C.3000907@colorado.edu> <56F5998F.40902@lanl.gov> From: Brian Gladman Message-ID: <56F5B964.7050401@gladman.plus.com> Date: Fri, 25 Mar 2016 22:20:00 -0000 User-Agent: Mozilla/5.0 (Windows NT 10.0; WOW64; rv:38.0) Gecko/20100101 Thunderbird/38.6.0 MIME-Version: 1.0 In-Reply-To: <56F5998F.40902@lanl.gov> Content-Type: text/plain; charset=utf-8 Content-Transfer-Encoding: 8bit X-SW-Source: 2016-q1/txt/msg00030.txt.bz2 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