From 9a8e2ffbf55e3ef53f2581af32112431f5eb32ef Mon Sep 17 00:00:00 2001 From: Raymond Hettinger Date: Sat, 28 Jan 2023 07:50:18 -0600 Subject: [PATCH 1/6] Fix hardwired constant --- Lib/test/test_math.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index fcdec93484e000..2665a5ed3f9a44 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -1441,7 +1441,7 @@ def RelativeError(res, ex): return fabs(n / target_sumprod) def Trial(dotfunc, c, n): - ex = GenDot(10, c) + ex = GenDot(n, c) res = dotfunc(ex.x, ex.y) return RelativeError(res, ex) From de81730be8ed0c1a9c1711de7b04f12fdd47e686 Mon Sep 17 00:00:00 2001 From: Raymond Hettinger Date: Sat, 28 Jan 2023 08:25:52 -0600 Subject: [PATCH 2/6] Instead of median, test all examples better than the target condition number. --- Lib/test/test_math.py | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index 2665a5ed3f9a44..28999e51b88b3c 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -1440,17 +1440,16 @@ def RelativeError(res, ex): n = DotExact(list(x) + [-res], list(y) + [1]) return fabs(n / target_sumprod) - def Trial(dotfunc, c, n): - ex = GenDot(n, c) - res = dotfunc(ex.x, ex.y) - return RelativeError(res, ex) - - times = 1000 # Number of trials - n = 20 # Length of vectors - c = 1e30 # Target condition number - - relative_err = median(Trial(math.sumprod, c, n) for i in range(times)) - self.assertLess(relative_err, 1e-16) + vector_length = 20 + target_condition_number = 1e30 + for i in range(2000): + ex = GenDot(n=vector_length, c=target_condition_number) + if ex.condition > target_condition_number: + continue + res = math.sumprod(ex.x, ex.y) + relative_error = RelativeError(res, ex) + self.assertLess(relative_error, 1e-15, + (ex, res, relative_error)) def testModf(self): self.assertRaises(TypeError, math.modf) From 078c1ecad488ce38433d271dbfb60543105bd9de Mon Sep 17 00:00:00 2001 From: Raymond Hettinger Date: Sat, 28 Jan 2023 09:00:24 -0600 Subject: [PATCH 3/6] Stronger test --- Lib/test/test_math.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index 28999e51b88b3c..8c5410bad40644 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -1407,6 +1407,7 @@ def GenDot(n, c): and calculating yi such that some cancellation occurs. Finally, we permute the vectors x, y randomly and calculate the achieved condition number. """ + # Test generator from: https://www.tuhh.de/ti3/paper/rump/OgRuOi05.pdf assert n >= 6 n2 = n // 2 @@ -1442,13 +1443,13 @@ def RelativeError(res, ex): vector_length = 20 target_condition_number = 1e30 - for i in range(2000): + for i in range(10_000): ex = GenDot(n=vector_length, c=target_condition_number) if ex.condition > target_condition_number: continue res = math.sumprod(ex.x, ex.y) relative_error = RelativeError(res, ex) - self.assertLess(relative_error, 1e-15, + self.assertLessEqual(relative_error, 2.0 ** -53, (ex, res, relative_error)) def testModf(self): From e3dbe9bca1aace8a6bcbf222e29aa2b9168c12c9 Mon Sep 17 00:00:00 2001 From: Raymond Hettinger Date: Sat, 28 Jan 2023 09:47:22 -0600 Subject: [PATCH 4/6] Replace relative error test with absolute test. --- Lib/test/test_math.py | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index 8c5410bad40644..46c55b0c55fe10 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -1379,7 +1379,7 @@ def test_sumprod_extended_precision_accuracy(self): from fractions import Fraction from itertools import starmap from collections import namedtuple - from math import log2, exp2, fabs + from math import log2, exp2, fabs, ulp from random import choices, uniform, shuffle from statistics import median @@ -1436,21 +1436,18 @@ def GenDot(n, c): return DotExample(x, y, DotExact(x, y), Condition(x, y)) - def RelativeError(res, ex): - x, y, target_sumprod, condition = ex - n = DotExact(list(x) + [-res], list(y) + [1]) - return fabs(n / target_sumprod) + def AbsoluteError(res, ex): + return DotExact(list(ex.x) + [-res], list(ex.y) + [1]) vector_length = 20 target_condition_number = 1e30 - for i in range(10_000): - ex = GenDot(n=vector_length, c=target_condition_number) + for i in range(1_000_000): + ex = GenDot(n=vector_length, c=target_condition_number / 2) if ex.condition > target_condition_number: continue - res = math.sumprod(ex.x, ex.y) - relative_error = RelativeError(res, ex) - self.assertLessEqual(relative_error, 2.0 ** -53, - (ex, res, relative_error)) + result = math.sumprod(ex.x, ex.y) + error = AbsoluteError(result, ex) + self.assertLess(fabs(error), ulp(result), (ex, result, error)) def testModf(self): self.assertRaises(TypeError, math.modf) From df42d9f60905681c9d1b43586df39c9dd8f43562 Mon Sep 17 00:00:00 2001 From: Raymond Hettinger Date: Sat, 28 Jan 2023 09:54:43 -0600 Subject: [PATCH 5/6] Add comments and todo --- Lib/test/test_math.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index 46c55b0c55fe10..fa02b185d38559 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -1442,12 +1442,19 @@ def AbsoluteError(res, ex): vector_length = 20 target_condition_number = 1e30 for i in range(1_000_000): + # GenDot creates examples whether the condition numbers + # are centered about c but can be two orders of magnitude + # above or below. Here, we generate examples centered + # around half our target and then reject examples higher + # than our target. ex = GenDot(n=vector_length, c=target_condition_number / 2) if ex.condition > target_condition_number: continue result = math.sumprod(ex.x, ex.y) error = AbsoluteError(result, ex) self.assertLess(fabs(error), ulp(result), (ex, result, error)) + # XXX Compute the theoretical bound for n=20, K=3, cond=1e30. + # ??? Is 1e30 too aggressive (within practical but not theoretical bounds. def testModf(self): self.assertRaises(TypeError, math.modf) From b7d72debcaa712ef10302e3b903bd06c1717af43 Mon Sep 17 00:00:00 2001 From: Raymond Hettinger Date: Sat, 28 Jan 2023 12:53:51 -0600 Subject: [PATCH 6/6] Relate the test parameters back to the paper --- Lib/test/test_math.py | 58 +++++++++++++++++++++++++++++++++---------- 1 file changed, 45 insertions(+), 13 deletions(-) diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index fa02b185d38559..b90409c70e1c2b 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -1406,6 +1406,12 @@ def GenDot(n, c): and y is then constructed choosing xi randomly with decreasing exponent, and calculating yi such that some cancellation occurs. Finally, we permute the vectors x, y randomly and calculate the achieved condition number. + + In general the actual condition number is a little larger than the + anticipated one with quite some variation. However, this is unimportant + because in the following figures we plot the relative error against the + actually achieved condition number. + """ # Test generator from: https://www.tuhh.de/ti3/paper/rump/OgRuOi05.pdf @@ -1436,25 +1442,51 @@ def GenDot(n, c): return DotExample(x, y, DotExact(x, y), Condition(x, y)) - def AbsoluteError(res, ex): - return DotExact(list(ex.x) + [-res], list(ex.y) + [1]) + def AbsoluteError(x, y, result): + return DotExact(list(x) + [-result], list(y) + [1]) + + def err_gamma(n): # Definition (2.4) + return n * ulp(1.0) / (1 - n * ulp(1.0)) + + def dotk_bound(cond, n, K=3): # Proposition (5.11) + "Relative error bound for DotK" + return ulp(1.0) + 2 * err_gamma(4*n-2)**2 + 1/2 * err_gamma(4*n-2)**K * cond + + # For a vector length of 20, the triple length sumprod calculation's + # theoretical error bound starts to degrade above a condition + # number of 1e25 where we get just under 52 bits of accuracy: + # + # >>> -log2(dotk_bound(1e25, n=20, K=3)) + # 51.84038876713239 + # >>> -log2(dotk_bound(1e26, n=20, K=3)) + # 50.8823973719395 + # >>> -log2(dotk_bound(1e27, n=20, K=3)) + # 48.3334013169794 + # + # The paper notes that the theoretical error bound is pessimistic + # by several orders of magnitude, so our test can use a higher + # condition number. (See Table 6.2) + # + # Here, we test a vector length of 20 and condition number of 1e28 + # for an absolute error difference within 1 ulp. This should be + # sufficiently loose to alway pass, but tight enough to fail in case + # of an implementation error, an incorrect compiler optimization, + # or an inadequate libmath implementation of fma(). If the latter + # ends up being a problem, we can revise dl_mul() to use Dekker's + # mul12() algorithm which slows sumprod() by 20% but would not be + # dependent on libmath() meeting the C99 specification for fma(). vector_length = 20 - target_condition_number = 1e30 - for i in range(1_000_000): - # GenDot creates examples whether the condition numbers - # are centered about c but can be two orders of magnitude - # above or below. Here, we generate examples centered - # around half our target and then reject examples higher - # than our target. - ex = GenDot(n=vector_length, c=target_condition_number / 2) + target_condition_number = 1e28 + for i in range(5_000): + # Generate examples centered around a quarter of the target condition + # number. Reject actual condition numbers higher than our target. + ex = GenDot(n=vector_length, c=target_condition_number / 4) if ex.condition > target_condition_number: continue result = math.sumprod(ex.x, ex.y) - error = AbsoluteError(result, ex) + error = AbsoluteError(ex.x, ex.y, result) self.assertLess(fabs(error), ulp(result), (ex, result, error)) - # XXX Compute the theoretical bound for n=20, K=3, cond=1e30. - # ??? Is 1e30 too aggressive (within practical but not theoretical bounds. def testModf(self): self.assertRaises(TypeError, math.modf)