From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 44725 invoked by alias); 25 Mar 2016 21:24:51 -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 44629 invoked by uid 89); 25 Mar 2016 21:24:44 -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=despair, stick, terrible, blind X-HELO: proofpoint4.lanl.gov Received: from proofpoint4.lanl.gov (HELO proofpoint4.lanl.gov) (204.121.3.52) by sourceware.org (qpsmtpd/0.93/v0.84-503-g423c35a) with (AES256-SHA encrypted) ESMTPS; Fri, 25 Mar 2016 21:24:27 +0000 Received: from mailrelay2.lanl.gov (mailrelay2.lanl.gov [128.165.4.103]) by mailgate4.lanl.gov (8.15.0.59/8.15.0.59) with ESMTP id u2PLOP7x022243 for ; Fri, 25 Mar 2016 15:24:25 -0600 Received: from localhost (localhost.localdomain [127.0.0.1]) by mailrelay2.lanl.gov (Postfix) with ESMTP id 679DEEF9E5C for ; Fri, 25 Mar 2016 15:24:25 -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 4360CEFA9D0 for ; Fri, 25 Mar 2016 15:24:25 -0600 (MDT) Subject: error estimates (was: GSL K0/K1) To: gsl-discuss@sourceware.org References: <56F42FED.8060500@colorado.edu> <56F4304C.3000907@colorado.edu> From: Gerard Jungman Message-ID: <56F5AC89.6050708@lanl.gov> Date: Fri, 25 Mar 2016 21:24: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/msg00028.txt.bz2 On the question of error estimates. I have addressed some of this before on this list, although in more general contexts. Here is what I have to say about error estimates in particular. Warning: This is probably going to be a long one, like many of my cries of despair on this list... First, some agreement on the behavior of function evaluations. The best possible behavior for an evaluation is to conform to the following specification: The function f(x) must return a result, for a machine represented number x, which is the closest possible machine representable value to the mathematically correct value, correctly rounded. For something like the platform-provided sin(x) and cos(x), this is definitely the specification. Of course, some implementations fall short. There are examples, and associated discussions to be googled... https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/ I have not checked this one myself. But you get the idea. I hope this one is not really true... anybody want to check? If a function conforms to the best-possible specification, then there is no need for an error estimate. The estimate is implicit in the contract. But this is hard enough to achieve, even for sin(x). If we stuck to that, there might no functions at all. So we sacrifice consistency and quality for coverage. The overriding problem with the special functions is the wide range of quality. Some of the functions are pretty good (although maybe not as good as hoped, i.e. K0/K1). Some are kind of difficult and lose in corner cases. And some are just plain nuts. Amongst the "just plain nuts" examples, the confluent hypergeometrics stick out. The error estimates are a kind of apology for this state of affairs, but they are obviously not very desirable either. Here is a quote from a message that I sent to this list several years ago: "This turns out to be a poor apology, since it tends to gum-up the works for the whole sublibrary, eating performance and occupying a large piece of intellectual real estate. The error estimation code must, at the least, be factored out. More to the point, it should likely be discarded in the main line of development. Other notions of error control should be investigated." So I think we are in agreement about what is desirable. But I have only the vaguest notions of how to get there. Here is a further quote from the same message: "Similar to the notion of sublibrary dependence discussed elsewhere, the special functions should be hierarchically organized. The dependencies should be made clear, and foundational levels of ironclad functions should be made explicit. Such ironclad functions should occupy the same place in the users mind as platform-native implementations of sin(x); they should be beyond question for daily use. Other functions, which suffer from implementation problems or are simply too complicated to guarantee the same level of correctness should be explicitly identified, as part of the design (not just the documentation!)." For example, functions which are ironclad should not bother to do error estimates. Functions which are harder to control should do something to indicate how they are behaving, whether this means returning an error estimate, or just returning some kind of "loss of accuracy" flag. Of course, the heuristics needed to identify loss of accuracy might be just as expensive as or equivalent to an actual error estimate. Functions which are "just plain nuts" might be pushed into an "experimental" sublibrary. Along with the problems involving variable quality, there is some actual intellectual confusion in this area as well. Examples of this are the versions of trig functions which return error estimates. Not only should these be among the ironclad functions (and probably not even provided by the library... what is the point), but the notion of "error estimate" expressed there is just confused. What those functions are trying to do is detect and report some form of ill-conditioning in the arguments, related to argument reduction. This is completely different from the theology expressed in the best-possible specification above. In that specification, the argument x is a known machine represented number, so there is no question of its conditioning. I remember writing these trig functions and puzzling over what I was doing. I was trying to fix some instability problem in the testing of functions which depended on trig evaluation; the ill-conditioning for large arguments was making it hard to create stable tests. But this was simply a problem with the tests and not with the functions. That is one day that I should have just sat on my hands and typed nothing; we would have been better off in the long run. I think the problem of argument conditioning is fundamentally a client-side issue. We should probably do nothing about it. If the client wants to evaluate sin(1.0e300), then it is probably not our problem. We should just do our best to return f(x) for the given machine-represented x. Ok, back to the real problem of error estimation and variable quality. What constitutes an ironclad function? In essence any function of a single variable (with no other parameters) must be ironclad. At the very least, it can be expressed with some form of Chebyshev expansion, argument transformation, etc, which should meet the best-possible specification. In other cases, there is an efficient high-quality algorithm which makes the evaluation easy, like the Lanczos algorithm for gamma(z). This includes functions of a single complex variable which it should be possible to treat in the same manner. What about functions with one "parameter", like besselK(nu, x)? These are harder and tend to have problems with corner cases or large arguments. But many of these should also be ironclad. Most of them have well-known stable methods, and the implementation only needs to take some care to get everything right. What about besselJ(nu, x)? This should also be ironclad. Although it may have some problems with argument conditioning, which can arise for many oscillatory functions, we have decided that argument-conditioning is not a problem we can address. So it should be no harder than sin(x), in principle. Stable iterative methods are applicable in many cases. Continued fraction methods are good, and I use these in many places. But they can also fail mysteriously. Gautschi has written some classic papers on this subject. Asymptotic methods are needed to fill the gaps. The tradeoff between iteration and asymptotic expansion has not been studied well enough. Certainly not by me. It would be best if the tradeoff point were detected automatically in some way, by the implementation. Otherwise, it becomes harder to understand what the code is doing, with many hard-wired constants, etc, that control the switch on methods. These matching problems are really tedious to think about explicitly; the machine should be tasked to figure that out. What about functions with more parameters? First the good news. All the elliptic functions and related elementary functions should be ironclad. A great deal is known about the evaluation of these. At the very least, you can always fall back on the arithmetic-geometric mean, which is stable and correct, though maybe not always the most efficient. Now the bad news: hypergeometric functions, including the regular and confluent functions. These are terrible. Probably all must go to the "experimental" ghetto. Functions with more parameters are harder to test as well. Maybe even impossible to test. Certainly, blind testing seems useless. You have to know what errors to look for. Another area of doubtful intellectual content: the "mode" argument for some of the evaluations. At the time, I wanted to be able to decrease accuracy at runtime, for performance gain in return. The obvious thing to do was to control the order of Chebyshev expansions. Fewer coefficients should be faster right? Well, who knows. It was never really studied properly. And it was only applicable for functions which depend on underlying Chebyshev expansions. Probably just a bad idea. The main library type is 'double'; the functions should probably just run as close to double precision as possible. Architectures have changes a lot since then as well, so spending extra cycles on some small local computations is not a big deal.