Skip site navigation (1)Skip section navigation (2)
Date:      Sat, 28 Jul 2012 00:44:04 -0500
From:      Stephen Montgomery-Smith <stephen@missouri.edu>
To:        Bruce Evans <brde@optusnet.com.au>
Cc:        freebsd-bugs@freebsd.org, FreeBSD-gnats-submit@freebsd.org, Stephen Montgomery-Smith <stephen@freebsd.org>
Subject:   Re: bin/170206: complex arcsinh, log, etc.
Message-ID:  <50137C24.1060004@missouri.edu>
In-Reply-To: <20120728142915.K909@besplex.bde.org>
References:  <201207270247.q6R2lkeR021134@wilberforce.math.missouri.edu> <20120727233939.A7820@besplex.bde.org> <5012A96E.9090400@missouri.edu> <20120728142915.K909@besplex.bde.org>

next in thread | previous in thread | raw e-mail | index | archive | help
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.



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?50137C24.1060004>