High accuracy math.hypot()

@tim-one

BPO	41513
Nosy	@tim-one, @rhettinger, @terryjreedy, @mdickinson, @vstinner, @serhiy-storchaka
PRs	bpo-41513: Improve speed and accuracy of math.hypot() #21803 bpo-41513: More accurate hypot() #21916 bpo-41513: Save unnecessary steps in the hypot() calculation #21994 Further improve accuracy of math.hypot() #22013 bpo-41513: Improve hypot() accuracy with three separate accumulators #22032 bpo-41513: Expand comments #22123 bpo-41513: Add docs and tests for hypot() #22238 bpo-41513: Remove broken tests that fail on Gentoo #22249 bpo-41513: Make steps after the loop consistent with the loop body #22315 bpo-41513: Second attempt to add accuracy tests for math.hypot() #22327 bpo-41513: Improve order of adding fractional values. Improve variable names. #22368
Files	test_hypot.py: Script for measuring accuracy test_hypot_accuracy.py: Cross version accuracy tester hypot_flow.pdf: Data flow visualization hypot.png: Updated data flow visualization test_hypot_commutativity.py: Test commutativity best_frac.py: Measure accuracy of hypot() variants

^{Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.}

Show more details

GitHub fields:

assignee = None
closed_at = <Date 2020-09-15.00:15:26.344>
created_at = <Date 2020-08-10.01:14:14.750>
labels = ['type-feature', 'library', '3.10']
title = 'High accuracy math.hypot()'
updated_at = <Date 2020-10-02.00:19:15.829>
user = 'https://github.com/rhettinger'

bugs.python.org fields:

activity = <Date 2020-10-02.00:19:15.829>
actor = 'rhettinger'
assignee = 'none'
closed = True
closed_date = <Date 2020-09-15.00:15:26.344>
closer = 'rhettinger'
components = ['Library (Lib)']
creation = <Date 2020-08-10.01:14:14.750>
creator = 'rhettinger'
dependencies = []
files = ['49380', '49399', '49424', '49439', '49448', '49484']
hgrepos = []
issue_num = 41513
keywords = ['patch']
message_count = 39.0
messages = ['375089', '375096', '375166', '375437', '375440', '375442', '375447', '375448', '375483', '375489', '375498', '375500', '375501', '375591', '375600', '375623', '375634', '375635', '375651', '375685', '375721', '375762', '375763', '375770', '375829', '375865', '375866', '375867', '375868', '375896', '376068', '376098', '376470', '376869', '376910', '376919', '376920', '377242', '377354']
nosy_count = 6.0
nosy_names = ['tim.peters', 'rhettinger', 'terry.reedy', 'mark.dickinson', 'vstinner', 'serhiy.storchaka']
pr_nums = ['21803', '21916', '21994', '22013', '22032', '22123', '22238', '22249', '22315', '22327', '22368']
priority = 'normal'
resolution = 'fixed'
stage = 'resolved'
status = 'closed'
superseder = None
type = 'enhancement'
url = 'https://bugs.python.org/issue41513'
versions = ['Python 3.10']

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) more accurate.

For accuracy, the current approach has the advantage that for the largest coordinate, (x/max)**2 is exactly 1.0. If that element is much larger than the others, it overpowers the 1/2 ulp error in each of the other divisions.

In contrast, scaling by a power of two is always exact, but the squaring of the largest scaled coordinate is no longer exact.

The two approaches each have some cases that are more accurate than the other. I generated 10,000 samples. In 66% of the cases, the two results agree. In 26% of the cases, scaling by a power of two was more accurate. In the remaining 8%, the division by max was more accurate. The standard deviation of relative errors was better using power-of-two scaling. These results were consistent whether using 2, 3, or 20 dimensions.

As expected, multiplying by a power of two was a bit faster than dividing by an arbitrary float.

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)]

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.

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)]

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.

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)

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)

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? ;-)

...
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)

One more implication: since the quality of the initial square root doesn't really much matter, instead of

    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)
    result += x / (2.0 * result)

at the end, it should work just as well (in fact, probably a little better, due to getting more beneficial cancellation) to do:

    a = parts[0] - 1.0
    result = sqrt(a)
    x, y = split(result)
    result += (a - x*x - 2.0*x*y - y*y + parts[1]) / (2.0 * result)

at the end. Although this crucially relies on the doing the last-line chain of subtracts and adds "left to right".

New changeset 27de286 by Raymond Hettinger in branch 'master':
bpo-41513: Save unnecessary steps in the hypot() calculation (bpo-21994)
27de286

About test_frac.py, I changed the main loop like so:

            got = [float(expected)] # NEW
            for hypot in hypots:
                actual = hypot(*coords)
                got.append(float(actual)) # NEW
                err = (actual - expected) / expected
                bits = round(1 / err).bit_length()
                errs[hypot][bits] += 1
            if len(set(got)) > 1: # NEW
                print(got) # NEW

That is, to display every case where the four float results weren't identical.

Result: nothing was displayed (although it's still running the n=1000 chunk, there's no sign that will change). None of these variations made any difference to results users actually get.

Even the "worst" of these reliably develops dozens of "good bits" beyond IEEE double precision, but invisibly (under the covers, with no visible effect on delivered results).

So if there's something else that speeds the code, perhaps it's worth pursuing, but we're already long beyond the point of getting any payback for pursuing accuracy.

New changeset 67c998d by Raymond Hettinger in branch 'master':
bpo-41513: Expand comments and add references for a better understanding (GH-22123)
67c998d

New changeset 457d4e9 by Raymond Hettinger in branch 'master':
bpo-41513: Add docs and tests for hypot() (GH-22238)
457d4e9

Multiple tests failed on x86 Gentoo Non-Debug with X 3.x:
https://buildbot.python.org/all/#/builders/58/builds/73

Example of failure:

======================================================================
FAIL: testHypotAccuracy (test.test_math.MathTests) (hx='0x1.10106eb4b44a2p+29', hy='0x1.ef0596cdc97f8p+29')
----------------------------------------------------------------------

Traceback (most recent call last):
  File "/buildbot/buildarea/cpython/3.x.ware-gentoo-x86.nondebug/build/Lib/test/test_math.py", line 855, in testHypotAccuracy
    self.assertEqual(hypot(x, y), z)
AssertionError: 1184594898.793983 != 1184594898.7939832

(I didn't investigate the bug, I just report the failure.)

New changeset 95a8a0e by Raymond Hettinger in branch 'master':
bpo-41513: Remove broken tests that fail on Gentoo (GH-22249)
95a8a0e

I saw that the Gentoo bot reacted badly to the new tests, so I reverted them.

New changeset bc6b7fa by Raymond Hettinger in branch 'master':
bpo-41513: Add accuracy tests for math.hypot() (GH-22327)
bc6b7fa

New changeset 438e9fc by Raymond Hettinger in branch 'master':
bpo-41513: Improve order of adding fractional values. Improve variable names. (GH-22368)
438e9fc

Navigation Menu

Repository navigation

High accuracy math.hypot() #85685

28 remaining items

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

High accuracy math.hypot() #85685

Description

Activity

rhettinger commented on Aug 10, 2020

rhettinger commented on Aug 10, 2020

mdickinson commented on Aug 11, 2020

rhettinger commented on Aug 14, 2020

rhettinger commented on Aug 14, 2020

rhettinger commented on Aug 15, 2020

rhettinger commented on Aug 15, 2020

tim-one commented on Aug 15, 2020

tim-one commented on Aug 15, 2020

28 remaining items

tim-one commented on Aug 25, 2020

rhettinger commented on Aug 29, 2020

tim-one commented on Aug 30, 2020

rhettinger commented on Sep 6, 2020

rhettinger commented on Sep 14, 2020

vstinner commented on Sep 14, 2020

rhettinger commented on Sep 15, 2020

rhettinger commented on Sep 15, 2020

rhettinger commented on Sep 21, 2020

rhettinger commented on Sep 23, 2020

Metadata

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Issue actions