High accuracy math.hypot()

[rhettinger]

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

[rhettinger]

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.

[rhettinger]

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

[mdickinson]

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.

[rhettinger]

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

[rhettinger]

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.

[rhettinger]

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)

[rhettinger]

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)

[tim-one]

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

[tim-one]

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