The following reply was made to PR bin/170206; it has been noted by GNATS.
From: Bruce Evans <b...@optusnet.com.au>
To: Stephen Montgomery-Smith <step...@missouri.edu>
Cc: Bruce Evans <b...@optusnet.com.au>,
Stephen Montgomery-Smith <step...@freebsd.org>,
freebsd-gnats-sub...@freebsd.org, freebsd-bugs@freebsd.org
Subject: Re: bin/170206: complex arcsinh, log, etc.
Date: Sat, 28 Jul 2012 15:25:04 +1000 (EST)
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
_______________________________________________
freebsd-bugs@freebsd.org mailing list
http://lists.freebsd.org/mailman/listinfo/freebsd-bugs
To unsubscribe, send any mail to "freebsd-bugs-unsubscr...@freebsd.org"
