From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 115715 invoked by alias); 26 Mar 2016 05:43:40 -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 115694 invoked by uid 89); 26 Mar 2016 05:43:39 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=0.2 required=5.0 tests=AWL,BAYES_00,RCVD_IN_DNSWL_LOW autolearn=ham version=3.3.2 spammy=realtime, real-time, cancellation, disastrous X-HELO: mail-lb0-f179.google.com Received: from mail-lb0-f179.google.com (HELO mail-lb0-f179.google.com) (209.85.217.179) by sourceware.org (qpsmtpd/0.93/v0.84-503-g423c35a) with (AES128-GCM-SHA256 encrypted) ESMTPS; Sat, 26 Mar 2016 05:43:28 +0000 Received: by mail-lb0-f179.google.com with SMTP id vo2so4611416lbb.1 for ; Fri, 25 Mar 2016 22:43:28 -0700 (PDT) X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20130820; h=x-gm-message-state:mime-version:in-reply-to:references:date :message-id:subject:from:to:cc; bh=y5/vDGjnCctsZDX4qXxLpqcW0mWCl6pEabBzpxR8cuw=; b=jguqr+9mnOYp+Jh7MSD/Kd6fzy7o4SLmC3PiIHK/1XjOcLo3y62vlpSHvcKHJzNWJJ 90QqS0XfLSpP4W4Go+JKZ2iIJnzTrh2hiIcrOsCHxvqD2ZjTKyuh2Ev2cnmZrrm91IdR QJDinoC3sSchCvDXKd/BQs+pfecdzwTvxOS7+FOFaLbParYGphEeNZneDGNXZKK3/Fmr e5sbPLbMJcXwar0z5oUi20V6OV9IX7WBZQj4rbRU9/rcFXTS9poubHrQWMGm+nUL58hI xY+ZsEFRfLcQFCu3Bzd/XXcLnVcFu4gLYwM95cT95BNV8tsIkywrecIWd01XrYgy+Wf9 8VYg== X-Gm-Message-State: AD7BkJKkm3WhG5DYJ1T4uYVFlq2B7MlVihGbw9oE95Vj/Mm0HzJMFJ4TlVQzXxCfqmiMJDtCIAqlGbeGMJToXA== MIME-Version: 1.0 X-Received: by 10.112.135.39 with SMTP id pp7mr5379604lbb.43.1458971005410; Fri, 25 Mar 2016 22:43:25 -0700 (PDT) Received: by 10.112.74.196 with HTTP; Fri, 25 Mar 2016 22:43:25 -0700 (PDT) In-Reply-To: <56F5AB1F.8090006@colorado.edu> References: <56F42FED.8060500@colorado.edu> <56F4304C.3000907@colorado.edu> <56F5998F.40902@lanl.gov> <56F5AB1F.8090006@colorado.edu> Date: Sat, 26 Mar 2016 05:43:00 -0000 Message-ID: Subject: Re: GSL K0/K1 From: Pavel Holoborodko To: Patrick Alken Cc: gsl-discuss@sourceware.org Content-Type: text/plain; charset=UTF-8 X-SW-Source: 2016-q1/txt/msg00031.txt.bz2 Dear Patrick and Gerard, Let me comment on the situation with K0/1 as I see some confusion. A. The main issue was NOT the incorrectly rounded coefficients (it was also contributing to overall accuracy loss but wasn't the main issue). The issues with original SLATEC/GSL computations of K1/0 (in order of increasing importance): 1. (Minor) Incorrectly rounded coefficients for Chebyshev expansion. All Bessel routines suffer from the issue - the most affected is I0/1. 2. (Moderate) Cancellation near x=2 because -log(x/2) was computed as -lx+M_LN2 (subtraction of near-equal numbers near x=2). Fix for the issue improved accuracy by two times near x=2 - from 10.52eps to 5.11eps. 3. (Severe). In general, Chebyshev expansion for K0/1 on [0,2] is unstable when used in finite precision floating-point arithmetic. It suffers from cancellation error - tries to compute small values by subtracting large ones. Theoretically (when we ignore floating-point effects) it gives high accuracy. But it is quite disastrous when it comes to practical finite precision arithmetic. As a result, we had accuracy loss on first interval of approximation used in SLATEC/GSL for K0: [0,2] ~ 10.52 eps [2,8] ~ 1.25 eps (8,inf] ~ 1.25 eps In order to resolve the third issue, I derived rational approximation for [0,1] and new Chebyshev expansion for [1,8]. Rational approximation was built taking into account not only theoretical approximation order but also with minimization of effects of floating-point. I think this is very important design constrain which was overlooked in original SLATEC, etc. Theoretical bounds is informative but should be considered only together with practical limitations of finite precision arithmetic. Now we have for K0: [0,1] ~ 2.04 eps (rational approximation) [1,8] ~ 1.28 eps (new Chebyshev) (8,inf] ~ 1.25 eps (original Chebyshev from SLATEC) Results for K1 is similar. Interestingly, other software based on SLATEC code (including MATLAB) suffer from the same accuracy loss. I published some plots here: http://goo.gl/KUodAu B. Probably the real-time error analysis was needed when no uniform floating-point standard existed. I think we can safely drop this thanks to 100% widespread of IEEE 754. However I think error bounds should be computed based on uniform sampling in real conditions for the reasons above. Also it is very difficult to take into account all rounding effects coming from arithmetic operations, elementary functions, etc. included in expressions. Testing against extended precision combines everything automatically and gives meaningful error bound (for floating-point arithmetic of IEEE 754 standard). C. I0/1 have decisive influence of K0/1 near zero. For example: K0(x) = -(log(x/2)+g)*I0(x)+... In the same time I0/1 have power series with extremely fast convergency (+all terms are positive, no cancellation). Thus it has sense to use I0/1 for computing the K0/1 on [0,1]. In fact it is probably the optimal way. D. Accuracy of I0/1 in GSL. Chebyshev expansions are excellent for I0/1. The only issue here is incorrectly rounded coefficients of expansion - this increases relative error by ~2 times. I am still studying the code and there is a chance that rational approximations might give better accuracy when x->inf. Pavel. On Sat, Mar 26, 2016 at 6:18 AM, Patrick Alken wrote: > Gerard, > > Yes sorry for including you so late in the discussion! Pavel (cc'd) found > that some of the GSL Chebyshev coefficients copied over from SLATEC were > rounded incorrectly (ie: just the very last digit or two was off which > caused quite a large error near x=2 for K0 and K1). > > The original GSL / SLATEC implementations broke up the K0/K1 interval in 3 > pieces: > [0,2] > [2,8] > [8,inf] > > Pavel developed new approximations using the intervals: > > [0,1] (new rational approximation from Pavel) > [1,8] (new Chebyshev from Pavel) > [8,inf] (original GSL Chebyshev) > > which avoids the troublesome point x=2 and has an overall lower error when > compared with arbitrary precision libraries over many random points. > > In principle we could simply correct the SLATEC rounding errors, but Pavel's > new method has even lower error than the corrected Chebyshev series. > > During some exhaustive numerical testing with random values, Pavel found > fairly good upper limits to the error on the K0/K1 on each of the 3 > intervals, which is why we think we should change the error estimate to the > empirical values found through the numerical simulations. > > Pavel found similar problems in I0/I1 which he is now working on correcting. > > Patrick > > On 03/25/2016 02:03 PM, 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. >> >> 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 >> >> >