Repository navigation
High accuracy math.hypot() #85685
Description
Activity
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)
28 remaining items
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".
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.
Multiple tests failed on x86 Gentoo Non-Debug with X 3.x:
https://buildbot.python.org/all/#/builders/58/builds/73Example 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.)
I saw that the Gentoo bot reacted badly to the new tests, so I reverted them.
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:
bugs.python.org fields: