[issue41513] Scale by power of two in math.hypot()

Change by Raymond Hettinger : -- Removed message: https://bugs.python.org/msg375827 ___ Python tracker ___ ___ Python-bugs-list

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger 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

[issue41513] Scale by power of two in math.hypot()

Terry J. Reedy 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

[issue41513] Scale by power of two in math.hypot()

Tim Peters 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

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger 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

[issue41513] Scale by power of two in math.hypot()

Tim Peters 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

[issue41513] Scale by power of two in math.hypot()

Tim Peters 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) /

[issue41513] Scale by power of two in math.hypot()

Tim Peters 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.999, which is 1

[issue41513] Scale by power of two in math.hypot()

Tim Peters added the comment: > That's about 0.50026 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

[issue41513] Scale by power of two in math.hypot()

Tim Peters added the comment: Here's a "correct rounding" fail for the add_on approach: xs = [16.004] * 9 decimal result = 48.01065814103642 which rounds to float 48.014 add_on result: 48.01 That's about 0.50026 ulp too small -

[issue41513] Scale by power of two in math.hypot()

Tim Peters 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

[issue41513] Scale by power of two in math.hypot()

Change by Raymond Hettinger : -- pull_requests: +21032 stage: resolved -> patch review pull_request: https://github.com/python/cpython/pull/21916 ___ Python tracker ___

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger added the comment: Hmm, I tried-out a C implementation and the timings weren't bad, 60% slower in exchange for always being correctly rounded. Do you all think that would be worth it? # New correctly rounded $ ./python.exe -m timeit -r11 -s 'from math import hypot'

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger added the comment: Here's a much cheaper way to get correctly rounded results almost all of the time. It uses Serhiy's idea for scaling by a power of two, Tim's idea for Veltkamp-Dekker splitting, my variant of Neumaier summation, and the paper's differential correction.

[issue41513] Scale by power of two in math.hypot()

Change by Raymond Hettinger : Added file: https://bugs.python.org/file49399/test_hypot_accuracy.py ___ Python tracker ___ ___

[issue41513] Scale by power of two in math.hypot()

Change by Raymond Hettinger : Added file: https://bugs.python.org/file49398/test_hypot_commutativity.py ___ Python tracker ___ ___

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger added the comment: If someone thinks there is a case for using the C library hypot() for the two-argument form, feel free to reopen this. Likewise, if someone thinks there is a case for doing the expensive but more accurate algorithm, go ahead and reopen this. Otherwise,

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger added the comment: New changeset fff3c28052e6b0750d6218e00acacd2fded4991a by Raymond Hettinger in branch 'master': bpo-41513: Improve speed and accuracy of math.hypot() (GH-21803) https://github.com/python/cpython/commit/fff3c28052e6b0750d6218e00acacd2fded4991a --

[issue41513] Scale by power of two in math.hypot()

Tim Peters added the comment: Oh no - I wouldn't use this as a default implementation. Too expensive. There is one aspect you may find especially attractive, though: unlike even the Decimal approach, it should be 100% insensitive to argument order (no info is lost before fsum() is called,

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger added the comment: > Cheapest way I know of that "seemingly always" reproduces > the Decimal result (when that's rounded back to float) > combines fsum(), Veltkamp splitting, and the correction > trick from the paper. That's impressive. Do you think this is worth

[issue41513] Scale by power of two in math.hypot()

Tim Peters added the comment: > ... > one at a time subtract a**2 (an argument squared) in descending > order of magnitude > ... But that doesn't really help unless the sum of squares was computed without care to begin with. Can do as well by skipping that but instead computing the original

[issue41513] Scale by power of two in math.hypot()

Tim Peters added the comment: Cute: for any number of arguments, try computing h**2, then one at a time subtract a**2 (an argument squared) in descending order of magnitude. Call that (h**2 - a1**2 - a2**2 - ...) x. Then h -= x/(2*h) That should reduce errors too, although not nearly

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger added the comment: Update: def hypot(*coords): # Corrected, unfused: https://arxiv.org/pdf/1904.09481.pdf # Simplified version omits scaling and handling of wide ranges # Has 1 ulp error 0.44% of the time (0.54% per the paper). # Can be

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger added the comment: Per the referenced paper, here's what is involved in further reducing error: def hypot(*coords): # Corrected, unfused: https://arxiv.org/pdf/1904.09481.pdf # Simplified version omits scaling and handling of wide ranges. # Only

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger added the comment: This paper addresses the topic directly: https://arxiv.org/pdf/1904.09481.pdf I tried to reproduce the results shown in Table 1 of the paper. For clib hypot(), they get 1 ulp errors 12.91% of the time. On my build, using the same gaussian

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger added the comment: > While we're doing this, any chance we could special-case > the two-argument hypot to use the libm hypot directly? Yes, that is an easy amendment (see below). I just tried it out on the macOs build but didn't see a change in accuracy from the current

[issue41513] Scale by power of two in math.hypot()

Mark Dickinson added the comment: Fine by me in principle; I haven't had a chance to look at the code yet. While we're doing this, any chance we could special-case the two-argument hypot to use the libm hypot directly? On many platforms the libm hypot will be correctly rounded, or close to

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger added the comment: Added a revised script for testing accuracy measured in ULPs. It shows an improvement for all dimensions tested, with the best for 2 dimensions. Also, the maximum error is now 1 ulp; formerly, it was 2 ulps. 2 dimensions div-by-max:

[issue41513] Scale by power of two in math.hypot()

Change by Raymond Hettinger : -- Removed message: https://bugs.python.org/msg375095 ___ Python tracker ___ ___ Python-bugs-list

[issue41513] Scale by power of two in math.hypot()

Change by Raymond Hettinger : Removed file: https://bugs.python.org/file49379/test_hypot.py ___ Python tracker ___ ___ Python-bugs-list

[issue41513] Scale by power of two in math.hypot()

Change by Raymond Hettinger : Added file: https://bugs.python.org/file49380/test_hypot.py ___ Python tracker ___ ___ Python-bugs-list

[issue41513] Scale by power of two in math.hypot()

Raymond Hettinger added the comment: Added a revised script for testing accuracy measured in ULPs. -- Added file: https://bugs.python.org/file49379/test_hypot.py ___ Python tracker

[issue41513] Scale by power of two in math.hypot()

Change by Raymond Hettinger : Removed file: https://bugs.python.org/file49378/test_hypot.py ___ Python tracker ___ ___ Python-bugs-list

[issue41513] Scale by power of two in math.hypot()

Change by Raymond Hettinger : Added file: https://bugs.python.org/file49378/test_hypot.py ___ Python tracker ___ ___ Python-bugs-list

[issue41513] Scale by power of two in math.hypot()

Change by Raymond Hettinger : -- keywords: +patch pull_requests: +20937 stage: -> patch review pull_request: https://github.com/python/cpython/pull/21803 ___ Python tracker

[issue41513] Scale by power of two in math.hypot()

New submission from Raymond Hettinger : I'd like to resurrect Serhiy's idea for using power-of-two scaling in math.hypot() rather than dividing each coordinate by the maximum value. I drafted a patch. It looks to be a little faster than what we have now and is generally (but not always)