From owner-freebsd-bugs@FreeBSD.ORG Sat Jul 28 05:44:05 2012 Return-Path: Delivered-To: freebsd-bugs@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52]) by hub.freebsd.org (Postfix) with ESMTP id 9528B106564A; Sat, 28 Jul 2012 05:44:05 +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 3F0C88FC15; Sat, 28 Jul 2012 05:44:05 +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 q6S5i4w4048692; Sat, 28 Jul 2012 00:44:04 -0500 (CDT) (envelope-from stephen@missouri.edu) Message-ID: <50137C24.1060004@missouri.edu> Date: Sat, 28 Jul 2012 00:44:04 -0500 From: Stephen Montgomery-Smith User-Agent: Mozilla/5.0 (X11; Linux i686; rv:14.0) Gecko/20120714 Thunderbird/14.0 MIME-Version: 1.0 To: Bruce Evans References: <201207270247.q6R2lkeR021134@wilberforce.math.missouri.edu> <20120727233939.A7820@besplex.bde.org> <5012A96E.9090400@missouri.edu> <20120728142915.K909@besplex.bde.org> In-Reply-To: <20120728142915.K909@besplex.bde.org> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Cc: freebsd-bugs@freebsd.org, FreeBSD-gnats-submit@freebsd.org, Stephen Montgomery-Smith 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:44:05 -0000 On 07/28/2012 12:25 AM, Bruce Evans wrote: > 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. Yes, but in my situation, I am rather sure that when I am adding highest to lowest that this won't occur. I am starting with -1, then adding something close to 1, then adding lots of smaller terms. And I find it very plausible that the kind of situation you describe won't happen. x0*x0 is close to 1. x0*x1 is at most sqrt(DE) times smaller. And so on. So I think the kind of situation you describe should never happen. As I said, I don't have a mathematical proof that the kind of thing you describe can NEVER happen. I just have never observed it happen.