From owner-freebsd-bugs@FreeBSD.ORG Sat Jul 28 05:25:14 2012 Return-Path: Delivered-To: freebsd-bugs@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id BA365106566B; Sat, 28 Jul 2012 05:25:14 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail06.syd.optusnet.com.au (mail06.syd.optusnet.com.au [211.29.132.187]) by mx1.freebsd.org (Postfix) with ESMTP id 4DA018FC15; Sat, 28 Jul 2012 05:25:13 +0000 (UTC) Received: from c122-106-171-246.carlnfd1.nsw.optusnet.com.au (c122-106-171-246.carlnfd1.nsw.optusnet.com.au [122.106.171.246]) by mail06.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q6S5P4BL010909 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Sat, 28 Jul 2012 15:25:05 +1000 Date: Sat, 28 Jul 2012 15:25:04 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <5012A96E.9090400@missouri.edu> Message-ID: <20120728142915.K909@besplex.bde.org> References: <201207270247.q6R2lkeR021134@wilberforce.math.missouri.edu> <20120727233939.A7820@besplex.bde.org> <5012A96E.9090400@missouri.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-bugs@freebsd.org, FreeBSD-gnats-submit@freebsd.org, Stephen Montgomery-Smith , Bruce Evans Subject: Re: bin/170206: complex arcsinh, log, etc. X-BeenThere: freebsd-bugs@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Bug reports List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sat, 28 Jul 2012 05:25:14 -0000 On Fri, 27 Jul 2012, Stephen Montgomery-Smith wrote: > On 07/27/2012 09:26 AM, Bruce Evans wrote: > >> % hm1 = -1; >> % for (i=0;i<12;i++) hm1 += val[i]; >> % return (cpack(0.5 * log1p(hm1), atan2(y, x))); >> >> It is the trailing terms that I think don't work right here. You sort >> them and add from high to low, but normally it is necessary to add >> from low to high (consider terms [1, DBL_EPSILON/2, DBL_EPSILON/4]). >> Adding from high to low cancels with the -1 term, but then only >> particular values work right. Also, I don't see how adding the low >> terms without extra precision preserves enough precision. > > I understand what you are saying. But in this case adding in order of > smallest to largest (adding -1 last) won't work. If all the signs in the > same direction, it would work. But -1 has the wrong sign. No, even if all the signs are the same, adding from the highest to lowest can lose precision. Normally at most 1 ulp, while cancelation can lose almost 2**MANT_DIG ulps. Example: #define DE DBL_EPSILON // for clarity (1) 1 + DE/2 = 1 (half way case rounded down to even) (2) 1 + DE/2 + DE/2 = 1 (double rounding) (3) DE/2 + DE/2 + 1 = 1 + DE (null rounding) We want to add -1 to a value near 1 like the above. Now a leading 1 in the above will cancel with the -1, and the the order in (3) becomes the inaccurate one. Modify the above by shifting the epsilons and adding another 1 to get an example for our context: (2') -1 + 1 + DE + DE*DE/2 + DE*DE/2 = DE (The leading +-1 are intentionally grouped and cancel. The other terms are (2) multiplied by DE, and suffer the same double rounding.) (3') -1 + 1 + DE*DE/2 + DE*DE/2 + DE = DE + DE*DE (The leading +-1 are intentionally grouped and cancel as before. The lower terms must be added from low to high, as in (3).) The right order is perhaps more transparently described as always from low to high, with suitable and explicit grouping of terms using parentheses to reduce cancelation errors. But I don't like parentheses and prefer to depend on left to right evaluation. With some parentheses, the above becomes: (3'') (DE**2/2 + DE**2/2 + DE + (1 + -1) (Full parentheses for the left to right order would be unreadable, so although the order is critical, they shouldn't be used to emphasize it.) Here the cancelation is exact, but in general it gives a nonzero term which might need to be sorted into the other terms. Strictly ordering all the terms is usually unnecessary and slow, and is usually not done. Neither is the analysis needed to prove that it is unnecessary. Even the above examples (3'*) are sloppy about this. They only work because the cancelation is exact. In (3'), (-1 + 1) is added first. That is correct since it is lowest (0). In (3'') (1 + -1) is added last. That is also correct, although the term is not the highest, since it is 0. Usually the cancelation won't be exact and gives a term that is far from lowest. Assuming that it is still highest is the best sloppy sorting. Efficient evaluation of polynomials usually requires regrouping them in numerically dangerous ways. I do some analysis for this. Bruce