Mailing List Archive

Change by Raymond Hettinger <raymond.hettinger@gmail.com>:


----------
keywords: +patch
pull_requests: +20937
stage: -> patch review
pull_request: https://github.com/python/cpython/pull/21803

_______________________________________
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

Change by Raymond Hettinger <raymond.hettinger@gmail.com>:


Added file: https://bugs.python.org/file49378/test_hypot.py

_______________________________________
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

Change by Raymond Hettinger <raymond.hettinger@gmail.com>:


Removed file: https://bugs.python.org/file49378/test_hypot.py

_______________________________________
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

Raymond Hettinger <raymond.hettinger@gmail.com> 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 <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

Change by Raymond Hettinger <raymond.hettinger@gmail.com>:


Removed file: https://bugs.python.org/file49379/test_hypot.py

_______________________________________
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

Change by Raymond Hettinger <raymond.hettinger@gmail.com>:


Added file: https://bugs.python.org/file49380/test_hypot.py

_______________________________________
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

Change by Raymond Hettinger <raymond.hettinger@gmail.com>:


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

_______________________________________
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

Raymond Hettinger <raymond.hettinger@gmail.com> 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: [(-2.0, 14), (-1.0, 1721), (0.0, 6533), (1.0, 1722), (2.0, 10)]
scaled-by-2: [(-1.0, 841), (0.0, 8339), (1.0, 820)]

======== 3 dimensions ========
div-by-max: [(-2.0, 3), (-1.0, 1525), (0.0, 6933), (1.0, 1535), (2.0, 4)]
scaled-by-2: [(-1.0, 739), (0.0, 8515), (1.0, 746)]

======== 5 dimensions ========
div-by-max: [(-2.0, 2), (-1.0, 1377), (0.0, 7263), (1.0, 1358)]
scaled-by-2: [(-1.0, 695), (0.0, 8607), (1.0, 698)]

======== 10 dimensions ========
div-by-max: [(-2.0, 13), (-1.0, 1520), (0.0, 6873), (1.0, 1580), (2.0, 14)]
scaled-by-2: [(-1.0, 692), (0.0, 8605), (1.0, 703)]

======== 20 dimensions ========
div-by-max: [(-1.0, 1285), (0.0, 7310), (1.0, 1405)]
scaled-by-2: [(-1.0, 555), (0.0, 8822), (1.0, 623)]

----------

_______________________________________
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

Mark Dickinson <dickinsm@gmail.com> 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 it. There was a noticeable accuracy regression for two-argument hypot when hypot gained the ability to support multiple arguments.

----------

_______________________________________
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

Raymond Hettinger <raymond.hettinger@gmail.com> 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 PR which uses lossless power-of-two-scaling and Neumaier summation.

In GCC's libquadmath implementation, the comments say that the error is less than 1 ulp, falling short of being correctly rounded within ±½ ulp.

If the platform hypots have the nearly the same accuracy as the current PR, it may make sense to skip the special case and opt for consistent cross-platform results.

==================================================================

$ git diff
diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index d0621f59df..3a42ea5318 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -2457,6 +2457,10 @@ vector_norm(Py_ssize_t n, double *vec, double max, int found_nan)
if (max == 0.0 || n <= 1) {
return max;
}
+ if (n == 2) {
+ /* C library hypot() implementations tend to be very accurate */
+ return hypot(vec[0], vec[1]);
+ }
frexp(max, &max_e);
if (max_e >= -1023) {
scale = ldexp(1.0, -max_e)

$ gcc --version
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/c++/4.2.1
Apple clang version 11.0.3 (clang-1103.0.32.62)
Target: x86_64-apple-darwin19.6.0
Thread model: posix
InstalledDir: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin

$ ./python.exe test_hypot.py

======== 2 dimensions ========
platform hypot(): [(-1.0, 8398), (0.0, 83152), (1.0, 8450)]
scaled-by-2 : [(-1.0, 8412), (0.0, 83166), (1.0, 8422)]

----------

_______________________________________
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

Raymond Hettinger <raymond.hettinger@gmail.com> 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 distribution, I get 14.08% for both the clang hypot() and for the current version of the PR, showing that the two are fundamentally the same.

----------

_______________________________________
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

Raymond Hettinger <raymond.hettinger@gmail.com> 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 handle the 2-dimensional cases.
# Has 1 ulp error 0.88% of the time (0.54% per the paper).
# Can be reduced to 0.00% with a reliable fused multiply-add.
a, b = map(fabs, coords)
h = sqrt(a*a + b*b)
if h <= 2*b:
delta = h - b
x = a*(2*delta - a) + (delta - 2*(a - b))*delta
else:
delta = h - a
x = 2*delta*(a - 2*b) + (4*delta - b)*b + delta*delta
return h - x/(2*h)

----------

_______________________________________
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

Raymond Hettinger <raymond.hettinger@gmail.com> 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 reduced to 0% with a fused multiply-add.
a, b = map(fabs, coords)
if a < b:
a, b = b, a
h = sqrt(a*a + b*b)
if h <= 2*b:
delta = h - b
x = a*(2*delta - a) + (delta - 2*(a - b))*delta
else:
delta = h - a
x = 2*delta*(a - 2*b) + (4*delta - b)*b + delta*delta
return h - x/(2*h)

----------

_______________________________________
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

Tim Peters <tim@python.org> 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 so effectively, since it's a cheap way of estimating (& correcting for) the discrepancy between sum(a**2) and h**2.

Note that "descending order" is important: it's trying to cancel as many leading bits as possible as early as possible, so that lower-order bits come into play.

Then again ... why bother? ;-)

----------

_______________________________________
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

Tim Peters <tim@python.org> 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 sum of squares in increasing order of magnitude.

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:

def split(x, T27=ldexp(1.0, 27)+1.0):
t = x * T27
hi = t - (t - x)
lo = x - hi
assert hi + lo == x
return hi, lo

# result := hypot(*xs)
parts = []
for x in xs:
a, b = split(x)
parts.append(a*a)
parts.append(2.0*a*b)
parts.append(b*b)
result = sqrt(fsum(parts))
a, b = split(result)
parts.append(-a*a)
parts.append(-2.0*a*b)
parts.append(-b*b)
x = fsum(parts) # best float approx to sum(x_i ** 2) - result**2
result += x / (2.0 * result)

----------

_______________________________________
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

Raymond Hettinger <raymond.hettinger@gmail.com> 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 implementing? Or should we declare victory with the current PR which is both faster and more accurate than what we have now?

Having 1-ulp error 17% of the time and correctly rounded 83% of the time is pretty darned good (and on par with C library code for the two-argument case).

Unless we go all-out with the technique you described, the paper shows that we're already near the limit of what can be done by trying to make the sum of squares more accurate, "... the correctly rounded square root of the correctly rounded a**2+b**2 can still be off by as much as one ulp. This hints at the possibility that working harder to compute a**2+b**2 accurately may not be the best path to a better answer".

FWIW, in my use cases, the properties that matter most are monotonicity, commutativity, cross-platform portability, and speed. Extra accuracy would nice to have but isn't essential and would likely never be noticed.

----------

_______________________________________
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

Tim Peters <tim@python.org> 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, and fsum's output should be unaffected by summand order).

----------

_______________________________________
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

Raymond Hettinger <raymond.hettinger@gmail.com> 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


----------

_______________________________________
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

Raymond Hettinger <raymond.hettinger@gmail.com> 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, I think we're done for now :-)

----------
resolution: -> fixed
stage: patch review -> resolved
status: open -> closed

_______________________________________
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

Change by Raymond Hettinger <raymond.hettinger@gmail.com>:


Added file: https://bugs.python.org/file49398/test_hypot_commutativity.py

_______________________________________
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

Change by Raymond Hettinger <raymond.hettinger@gmail.com>:


Added file: https://bugs.python.org/file49399/test_hypot_accuracy.py

_______________________________________
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

Raymond Hettinger <raymond.hettinger@gmail.com> 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. Together, these allow the algorithm to simulate 128-bit quad precision throughout.

# Veltkamp-Dekker splitting

def split(x, T27=ldexp(1.0, 27)+1.0):
t = x * T27
hi = t - (t - x)
lo = x - hi
assert hi + lo == x
return hi, lo

# Variant of Neumaier summation specialized
# for cases where |csum| >= |x| for every x.
# Establishes the invariant by setting *csum*
# to 1.0, only adding *x* values of smaller
# magnitude, and subtracting 1.0 at the end.

def zero():
return 1.0, 0.0 # csum, frac

def add_on(x, state):
csum, frac = state
assert fabs(csum) >= fabs(x)
oldcsum = csum
csum += x
frac += (oldcsum - csum) + x
return csum, frac

def to_float(state):
csum, frac = state
return csum - 1.0 + frac

def hypot(*xs):
max_value = max(xs, key=fabs)
_, max_e = frexp(max_value)
scalar = ldexp(1.0, -max_e)
parts = zero()
for x in xs:
x *= scalar
a, b = split(x)
parts = add_on(a*a, parts)
parts = add_on(2.0*a*b, parts)
parts = add_on(b*b, parts)
result = sqrt(to_float(parts))
a, b = split(result)
parts = add_on(-a*a, parts)
parts = add_on(-2.0*a*b, parts)
parts = add_on(-b*b, parts)
x = to_float(parts) # best float approx to sum(x_i ** 2) - result**2
result += x / (2.0 * result)
return result / scalar

----------

_______________________________________
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

Raymond Hettinger <raymond.hettinger@gmail.com> 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' 'hypot(1.7, 15.5, 0.3, 5.2)'
2000000 loops, best of 11: 101 nsec per loop
$ ./python.exe -m timeit -r11 -s 'from math import hypot' 'hypot(1.7, 15.81)'
5000000 loops, best of 11: 82.8 nsec per loop

# Current baseline for Python 3.10
$ ./python.exe -m timeit -r11 -s 'from math import hypot' 'hypot(1.7, 15.5, 0.3, 5.2)'
5000000 loops, best of 11: 65.6 nsec per loop
$ ./python.exe -m timeit -r11 -s 'from math import hypot' 'hypot(1.7, 15.81)'
5000000 loops, best of 11: 56.6 nsec per loop

# Current baseline for Python 3.9
$ ./python.exe -m timeit -r11 -s 'from math import hypot' 'hypot(1.7, 15.5, 0.3, 5.2)'
5000000 loops, best of 11: 72.2 nsec per loop
~/npython $ ./python.exe -m timeit -r11 -s 'from math import hypot' 'hypot(1.7, 15.81)'
5000000 loops, best of 11: 56.2 nsec per loop

----------
status: closed -> open

_______________________________________
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

Change by Raymond Hettinger <raymond.hettinger@gmail.com>:


----------
pull_requests: +21032
stage: resolved -> patch review
pull_request: https://github.com/python/cpython/pull/21916

_______________________________________
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