Mailing List Archive

1 2  View All
[issue41513] Scale by power of two in math.hypot() [ In reply to ]
Tim Peters <tim@python.org> added the comment:

About speed, the fsum() version I posted ran about twice as fast as the all-Decimal approach, but the add_on() version runs a little slower than all-Decimal. I assume that's because fsum() is coded in C while the add_on() prototype makes mounds of additional Python-level function calls. Using released Windows 3.8.5 CPython and on argument vectors of length 50.

Is it worth it? Up to you ;-) There are good arguments to be made in favor of increased accuracy, and in favor of speed.

About "correctly rounded", such a claim can only be justified by proof, not by testing (unless testing exhaustively covers the entire space of inputs). Short of that, best you can claim is stuff like "max error < 0.51 ulp" (or whatever other bound can be proved).

----------

_______________________________________
Python tracker <report@bugs.python.org>
<https://bugs.python.org/issue41513>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/list-python-bugs%40lists.gossamer-threads.com
[issue41513] Scale by power of two in math.hypot() [ In reply to ]
Tim Peters <tim@python.org> added the comment:

Here's a "correct rounding" fail for the add_on approach:

xs = [16.000000000000004] * 9

decimal result = 48.00000000000001065814103642
which rounds to float 48.000000000000014

add_on result: 48.00000000000001

That's about 0.5000000000026 ulp too small - shocking ;-)

The fsum approach got this one right, but has similar failures on other inputs of the same form; i.e.,

[i + ulp(i)] * 9

for various integers `i`.

----------

_______________________________________
Python tracker <report@bugs.python.org>
<https://bugs.python.org/issue41513>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/list-python-bugs%40lists.gossamer-threads.com
[issue41513] Scale by power of two in math.hypot() [ In reply to ]
Tim Peters <tim@python.org> added the comment:

> That's about 0.5000000000026 ulp too small - shocking ;-)

Actually, that's an illusion due to the limited precision of Decimal. The rounded result is exactly 1/2 ulp away from the infinitely precise result - it's a nearest/even tie case.

To make analysis easy, just note that the mathematical hypot(*([x] * 9)) is exactly 3*x. For x = 16 + ulp(16), 3*x moves to a higher binade and has two trailing one bits. The last has to be rounded away, rounding up to leave the last retained bit even. Any source of error, no matter how tiny, can be enough so that doesn't happen.

----------

_______________________________________
Python tracker <report@bugs.python.org>
<https://bugs.python.org/issue41513>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/list-python-bugs%40lists.gossamer-threads.com
[issue41513] Scale by power of two in math.hypot() [ In reply to ]
Tim Peters <tim@python.org> added the comment:

Here's an amusing cautionary tale: when looking at correct-rounding failures for the fsum approach, I was baffled until I realized it was actually the _decimal_ method that was failing. Simplest example I have is 9 instances of b=4.999999999999999, which is 1 ulp less than 5. Of course the best-possible hypot is 3*b rounded, but that's not what default decimal precision computes (when rounded back). In fact, it bounces between "off by 1" and "correct" until decimal precision exceeds 96:

>>> pfailed = []
>>> for p in range(28, 300):
... with decimal.localcontext() as c:
... c.prec = p
... if float(sum([decimal.Decimal(b) ** 2] * 9).sqrt()) != 3*b:
... pfailed.append(p)
>>> pfailed
[.28, 29, 30, 31, 34, 36, 37, 39, 41, 42, 43, 44, 45, 57, 77, 96]

In a nutshell, it takes in the ballpark of 50 decimal digits to represent b exactly, then twice as many to represent its square exactly, then another to avoid losing a digit when summing 9 of those. So if decimal prec is less than about 100, it can lose info.

So the fsum method is actually more reliable (it doesn't lose any info until the fsum step).

----------

_______________________________________
Python tracker <report@bugs.python.org>
<https://bugs.python.org/issue41513>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/list-python-bugs%40lists.gossamer-threads.com
[issue41513] Scale by power of two in math.hypot() [ In reply to ]
Tim Peters <tim@python.org> added the comment:

Just FYI, if the "differential correction" step seems obscure to anyone, here's some insight, following a chain of mathematical equivalent respellings:

result + x / (2 * result) =
result + (sumsq - result**2) / (2 * result) =
result + (sumsq - result**2) / result / 2 =
result + (sumsq / result - result**2 / result) / 2 =
result + (sumsq / result - result) / 2 =
result + sumsq / result / 2 - result / 2 =
result / 2 + sumsq / result / 2 =
(result + sumsq / result) / 2

I hope the last line is an "aha!" for you then: it's one iteration of Newton's square root method, for improving a guess "result" for the square root of "sumsq". Which shouldn't be too surprising, since Newton's method is also derived from a first-order Taylor expansion around 0.

Note an implication: the quality of the initial square root is pretty much irrelevant, since a Newton iteration basically doubles the number of "good bits". Pretty much the whole trick relies on computing "sumsq - result**2" to greater than basic machine precision.

----------

_______________________________________
Python tracker <report@bugs.python.org>
<https://bugs.python.org/issue41513>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/list-python-bugs%40lists.gossamer-threads.com
[issue41513] Scale by power of two in math.hypot() [ In reply to ]
Tim Peters <tim@python.org> added the comment:

My apologies if nobody cares about this ;-) I just think it's helpful if we all understand what really helps here.

> Pretty much the whole trick relies on computing
> "sumsq - result**2" to greater than basic machine
> precision.

But why is that necessary at all? Because we didn't compute the sqrt of the sum of the squares to begin with - the input to sqrt was a _rounded_-to-double approximation to the sum of squares. The correction step tries to account for the sum-of-squares bits sqrt() didn't see to begin with.

So, e.g., instead of doing

result = sqrt(to_float(parts))
a, b = split(result)
parts = add_on(-a*a, parts)
parts = add_on(-b*b, parts)
parts = add_on(-2.0*a*b, parts)
x = to_float(parts)
result += x / (2.0 * result)

it works just as well to do just

result = float((D(parts[0] - 1.0) + D(parts[1])).sqrt())

where D is the default-precision decimal.Decimal constructor. Then sqrt sees all the sum-of-squares bits we have (although the "+" rounds away those following the 28th decimal digit).

----------

_______________________________________
Python tracker <report@bugs.python.org>
<https://bugs.python.org/issue41513>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/list-python-bugs%40lists.gossamer-threads.com
[issue41513] Scale by power of two in math.hypot() [ In reply to ]
Raymond Hettinger <raymond.hettinger@gmail.com> added the comment:

> My apologies if nobody cares about this ;-)

I care :-)

Am in crunch right now, so won't have a chance to work through it for a week or so.

----------

_______________________________________
Python tracker <report@bugs.python.org>
<https://bugs.python.org/issue41513>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/list-python-bugs%40lists.gossamer-threads.com
[issue41513] Scale by power of two in math.hypot() [ In reply to ]
Tim Peters <tim@python.org> added the comment:

> won't have a chance to work through it for a week or so

These have been much more in the way of FYI glosses. There's no "suggestion" here to be pursued - just trying to get a deeper understanding of code already written :-)

While I can concoct any number of cases where the add_on() method fails correct rounding, all the ones I dug into turned out to be "oops! went the wrong way" in nearest/even tie cases. By any non-anal measure, its accuracy is excellent.

----------

_______________________________________
Python tracker <report@bugs.python.org>
<https://bugs.python.org/issue41513>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/list-python-bugs%40lists.gossamer-threads.com
[issue41513] Scale by power of two in math.hypot() [ In reply to ]
Terry J. Reedy <tjreedy@udel.edu> added the comment:

Tim, I have been reading this, without examining all the details, to learn more about FP accuracy problems.

----------
nosy: +terry.reedy

_______________________________________
Python tracker <report@bugs.python.org>
<https://bugs.python.org/issue41513>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/list-python-bugs%40lists.gossamer-threads.com
[issue41513] Scale by power of two in math.hypot() [ In reply to ]
Raymond Hettinger <raymond.hettinger@gmail.com> added the comment:

I'm ready to more forward on this one. I've improved the assertions, explanations, and variable names. I also added references to show that the foundations are firm. After tens of millions of trials, I haven't found a single misrounding. I can't prove "correctly rounded" but am content with "highly accurate". The speed is comparable to what we had in 3.7, so I think it is a net win all around.

----------
priority: low -> normal
resolution: fixed ->

_______________________________________
Python tracker <report@bugs.python.org>
<https://bugs.python.org/issue41513>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/list-python-bugs%40lists.gossamer-threads.com
[issue41513] Scale by power of two in math.hypot() [ In reply to ]
Change by Raymond Hettinger <raymond.hettinger@gmail.com>:


----------
Removed message: https://bugs.python.org/msg375827

_______________________________________
Python tracker <report@bugs.python.org>
<https://bugs.python.org/issue41513>
_______________________________________
_______________________________________________
Python-bugs-list mailing list
Unsubscribe: https://mail.python.org/mailman/options/python-bugs-list/list-python-bugs%40lists.gossamer-threads.com

1 2  View All