From owner-freebsd-current@FreeBSD.ORG Mon Jul 9 03:41:46 2012 Return-Path: Delivered-To: freebsd-current@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 82C46106566C for ; Mon, 9 Jul 2012 03:41:46 +0000 (UTC) (envelope-from stephen@missouri.edu) Received: from wilberforce.math.missouri.edu (wilberforce.math.missouri.edu [128.206.184.213]) by mx1.freebsd.org (Postfix) with ESMTP id 534178FC14 for ; Mon, 9 Jul 2012 03:41:46 +0000 (UTC) Received: from [127.0.0.1] (wilberforce.math.missouri.edu [128.206.184.213]) by wilberforce.math.missouri.edu (8.14.5/8.14.5) with ESMTP id q693fhiA028963; Sun, 8 Jul 2012 22:41:45 -0500 (CDT) (envelope-from stephen@missouri.edu) Message-ID: <4FFA52F8.2080700@missouri.edu> Date: Sun, 08 Jul 2012 22:41:44 -0500 From: Stephen Montgomery-Smith User-Agent: Mozilla/5.0 (X11; Linux i686; rv:13.0) Gecko/20120615 Thunderbird/13.0.1 MIME-Version: 1.0 To: Steve Kargl References: <20120528233035.GA77157@troutmask.apl.washington.edu> <4FC40DEA.8030703@missouri.edu> <20120529000756.GA77386@troutmask.apl.washington.edu> <4FC43C8F.5090509@missouri.edu> <20120529045612.GB4445@server.rulingia.com> <20120708124047.GA44061@zim.MIT.EDU> <210816F0-7ED7-4481-ABFF-C94A700A3EA0@bsdimp.com> <4FF9DA46.2010502@missouri.edu> <20120708235848.GB53462@troutmask.apl.washington.edu> <4FFA25EA.5090705@missouri.edu> <20120709020107.GA53977@troutmask.apl.washington.edu> In-Reply-To: <20120709020107.GA53977@troutmask.apl.washington.edu> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Cc: freebsd-current@freebsd.org Subject: Re: Use of C99 extra long double math functions after r236148 X-BeenThere: freebsd-current@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Discussions about the use of FreeBSD-current List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 09 Jul 2012 03:41:46 -0000 On 07/08/2012 09:01 PM, Steve Kargl wrote: > Have you read Goldberg's paper? I must admit that I had not. I found it at: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html > Not to mention, I've seen way too many examples of 'x - y' > where cancellation of significant digits causes > problems. Throw in rather poor estimates of function > results with real poor ULP and you have problems. I think that if a program includes "x-y" where cancellation causes huge loss of digits, then the program should be considered highly flawed. The programmer might get lucky, and a few extra ULP might save his skin, but it would be better if the program catastrophically failed so that he would realize they need to do something smarter. I got my first scientific calculator around 1975, and I remember my surprise when I discovered that (1/3)*3-1 on my calculator did not produce zero. Years later I bought another calculator, and imagine my further surprise when (1/3)*3-1 did produce zero. They must have done something very smart, I thought. The problem is, these kinds of tricks can only help you so much. Calculations like arccos(arccos(arccos(arccos(cos(cos(cos(cos(x)))))))))-x are probably not going to be correct no matter how smart the programmer is. The paper by Goldberg has some really nice stuff. I teach a class on numerical methods. One example I present is the problem using equation (4) of Goldberg's paper to solve quadratic equations. I tell them that the smart thing to do is to use equation (4) or equation (5) depending on whether b has the same sign as sqrt(b^2-4ac). This is a very good illustration of how to overcome the "x-y" problem, and in my opinion it has to be understood by the programmer, not the writer of the square-root package. Trying to compensate by getting that lost drop of ULP out of square root is looking in the wrong direction. But there is other stuff in his paper I don't like, and it comes across as nit-picking rather than really thinking through why he wants the extra accuracy. I feel like he is saving the pennies, but losing the dollars. For example, computing the quadratic formula when b^2 and 4ac are roughly equal (discussed in his proof of Theorem 4). He says that roughly half the digits of the final answer will be contaminated. He recommends performing the calculation with double the precision, and then retaining what is left. The reason I don't like his solution is this. The contamination of half the digits is real, not some artifact of the computing method. Unless you know EXACTLY what a, b and c are, only half the digits of the quadratic formula are valid, no matter how you calculate it. For example, maybe a is calculated by some other part of the program as 1/3. Since 1/3 cannot be represented exactly in base 2, it will not be exactly 1/3. That part of the program doesn't know that this number is going to be used later on in the quadratic formula, so it computes it in single precision. Now the quadratic formula algorithm converts a to double precision. But a will still be 1/3 to single precision. The only way to avoid this error is to perform EVERY pertinent calculation prior to the use of quadratic formula in double precision. (And suppose instead the data comes from a scientific experiment - the programmer needs to be performing an error analysis, not hoping his implementation of quadratic formula is performed in double precision.) Similarly, I am not a fan of the Kahan Summation Formula (Theorem 8). There are two reasons why I might compute large sums. One is accounting, when I better be using high enough precision that the arithmetic is exact. The other is something like numerical integration. But in the latter case, unless the summands have a very special form, it is likely that each summand will only be accurate to within e. Thus the true sum is only known to within n*e, and so the naive method really hasn't lost any precision at all. And typically n*e will be so small compared to the true answer that it won't matter. If it does matter, the problem is that n is too big. The algorithm will take decades to run. Focus efforts on producing a different numerical integration method, instead of getting that lost drop of ULP. I cannot envision any non-contrived circumstance where the Kahan summation formula is going to give an answer that is significantly more reliable than naive summation. I can envision a circumstance when the reduction in speed of the Kahan summation formula over naive summation could be significant. I think the example I like best that illustrates floating point errors is inverting badly conditioned matrices. And any scientist who works with large matrices had better know what ill-conditioned means. The proper answer is that if your matrix is ill-conditioned, give up. You need to look at the problem where your matrix came from, and rethink how you concert it into a matrix problem, so that the matrix is not ill-conditioned. Chasing that last drop of ULP simply will not help one iota. Stephen