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.
_______________________________________________
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"
