From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 113208 invoked by alias); 25 Mar 2016 20:03:33 -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 113192 invoked by uid 89); 25 Mar 2016 20:03:32 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-0.2 required=5.0 tests=BAYES_50,RP_MATCHES_RCVD,SPF_PASS autolearn=ham version=3.3.2 spammy=gsldiscuss, gsl-discuss, 7.2, uniform X-HELO: proofpoint5.lanl.gov Received: from proofpoint5.lanl.gov (HELO proofpoint5.lanl.gov) (204.121.3.53) by sourceware.org (qpsmtpd/0.93/v0.84-503-g423c35a) with (AES256-SHA encrypted) ESMTPS; Fri, 25 Mar 2016 20:03:30 +0000 Received: from mailrelay2.lanl.gov (mailrelay2.lanl.gov [128.165.4.103]) by mailgate5.lanl.gov (8.15.0.59/8.15.0.59) with ESMTP id u2PK3SRQ028786 for ; Fri, 25 Mar 2016 14:03:28 -0600 Received: from localhost (localhost.localdomain [127.0.0.1]) by mailrelay2.lanl.gov (Postfix) with ESMTP id 39D5DEFA97D for ; Fri, 25 Mar 2016 14:03:28 -0600 (MDT) X-NIE-2-Virus-Scanner: amavisd-new at mailrelay2.lanl.gov Received: from nostromo.lanl.gov (nostromo.lanl.gov [130.55.124.72]) by mailrelay2.lanl.gov (Postfix) with ESMTP id 0AFA2EFA978 for ; Fri, 25 Mar 2016 14:03:28 -0600 (MDT) Subject: Re: GSL K0/K1 To: gsl-discuss@sourceware.org References: <56F42FED.8060500@colorado.edu> <56F4304C.3000907@colorado.edu> From: Gerard Jungman Message-ID: <56F5998F.40902@lanl.gov> Date: Fri, 25 Mar 2016 20:03:00 -0000 User-Agent: Mozilla/5.0 (X11; Linux x86_64; rv:38.0) Gecko/20100101 Thunderbird/38.4.0 MIME-Version: 1.0 In-Reply-To: <56F4304C.3000907@colorado.edu> Content-Type: text/plain; charset=utf-8; format=flowed Content-Transfer-Encoding: 7bit X-Proofpoint-Virus-Version: vendor=fsecure engine=2.50.10432:5.15.96,1.0.3,0.0.0000 definitions=2016-03-25_06:2016-03-25,2016-03-25,1970-01-01 signatures=0 X-SW-Source: 2016-q1/txt/msg00026.txt.bz2 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. Maybe it was something more subtle. The basic method, of course, goes back to SLATEC, so I suppose that version has the same problems...? We can and should have a discussion about error estimates in general. But I will create another thread for that, so as not to expand or hijack this topic. 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. 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. Has anybody checked I0? If it is also messed up, then maybe we have found the culprit. The error estimate in that section is just an application of some derivative estimates for that particular expression. But why that expression is used is entirely a mystery to me. A brief comment on the random testing: This is a good way to estimate the coefficient in the uniform error bound. But in principle it should not be necessary. There must be rigorous bounds for these expansions based on the uniformity of the Chebyshev expansions. Can the necessary coefficients be evaluated directly as part of the computation of the expansions? That would be best of all. These error bounds could be documented somewhere for posterity, at the least. -- G. Jungman On 03/24/2016 12:22 PM, Patrick Alken wrote: > Pavel, > > The Bessel function results look good. If everything is tested you can > merge it all into the master branch. I'm cc'ing gsl-discuss so that this > discussion will start to be archived. > > Regarding the error estimates, an attempt was made early on in GSL to > provide error estimates for all the special functions. From section 7.2 > of the manual: > > ---- > The field val contains the value and the field err contains an estimate > of the absolute error in the value. > ---- > > I confess I don't fully understand all the error estimation, or even how > accurate the estimates are. I don't know if any studies were done to > verify the accuracy of these error estimates. One of the original GSL > developers (Gerard Jungman) wrote a lot of this code (I've cc'd him on > this email) - maybe he can help shed light on this. > > Since you've already done exhaustive comparisons of your new Bessel > functions against arbitrary precision libraries, maybe you could simply > change the error estimate to be: > > err = factor * GSL_DBL_EPSILON > > where factor are the numbers you've found in your error analysis. This > way the calculation would be very fast (and probably more accurate than > the current error estimation). > > We can wait and see if Gerard has any suggestions on this too. > > Gerard - since you haven't been following this email chain, Pavel has > found more accurate Chebyshev expansions for the Bessel K0/K1 functions > which have lower errors than the original GSL implementations when > compared with arbitrary precision libraries. Through lots of random > testing he has found upper limits on the errors of each function given > by some positive constant times GSL_DBL_EPSILON. We are wondering if its > safe to change the error estimation of these functions (I don't know the > origin of the current calculations for these error estimates). > > Patrick