From 185c365109b7b1b54a21416b512ca5e1547d2f5d Mon Sep 17 00:00:00 2001 From: Michael Date: Tue, 3 Nov 2026 18:56:53 +0100 Subject: [PATCH 1/4] gh-161259: Fix statistics._sqrtprod() for subnormal products --- Lib/statistics.py | 5 ++-- Lib/test/test_statistics.py | 23 +++++++++++++++++++ ...-11-03-18-56-28.gh-issue-161259.iyZz6y.rst | 3 +++ 3 files changed, 29 insertions(+), 2 deletions(-) create mode 100644 Misc/NEWS.d/next/Library/2026-11-03-18-56-28.gh-issue-161259.iyZz6y.rst diff --git a/Lib/statistics.py b/Lib/statistics.py index eee3e6eb94e990..0308e2578e2c88 100644 --- a/Lib/statistics.py +++ b/Lib/statistics.py @@ -1776,9 +1776,10 @@ def _sqrtprod(x: float, y: float) -> float: return _sqrtprod(scale * x, scale * y) / scale return h - if not h: + if h < 2.0 ** -485: + # x * y is subnormal or underflowed, or the correction term + # below would be. Rescale non-zero inputs. if x and y: - # Non-zero inputs underflowed, so scale up, and recompute. # Scale: 1 / sqrt(sys.float_info.min * sys.float_info.epsilon) scale = 2.0 ** 537 return _sqrtprod(scale * x, scale * y) / scale diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py index 700c5ac304f717..8c35cb9c8f7010 100644 --- a/Lib/test/test_statistics.py +++ b/Lib/test/test_statistics.py @@ -2790,6 +2790,29 @@ def test_sqrtprod_helper_function_fundamentals(self): self.assertEqual(actual, expected) self.assertEqual(sign(actual), sign(expected)) + def test_sqrtprod_helper_function_subnormal_products(self): + # gh-161259: accuracy was lost when x * y is subnormal. + # The reference scales the inputs by an exact power of two so that + # the product is a normal float. + scale = 2.0**600 + for x, y in [ + (2e-200, 2e-124), + (1e-160, 3e-150), + (3e-170, 7e-140), + (1e-160, 1e-135), + ]: + expected = math.sqrt((scale * x) * (scale * y)) / scale + actual = statistics._sqrtprod(x, y) + with self.subTest(x=x, y=y, expected=expected, actual=actual): + self.assertAlmostEqual(actual / expected, 1.0, delta=2**-50) + + def test_correlation_tiny_inputs(self): + # gh-161259: correlation() is scale-invariant, even when the + # intermediate sums of squares are tiny. + x = [1e-100, 2e-100, 3e-100] + y = [1e-62, 2e-62, 3e-62] + self.assertAlmostEqual(statistics.correlation(x, y), 1.0) + @requires_IEEE_754 @skip_if_double_rounding @support.cpython_only # Allow for a weaker sumprod() implementation diff --git a/Misc/NEWS.d/next/Library/2026-11-03-18-56-28.gh-issue-161259.iyZz6y.rst b/Misc/NEWS.d/next/Library/2026-11-03-18-56-28.gh-issue-161259.iyZz6y.rst new file mode 100644 index 00000000000000..c59173bf3bae8f --- /dev/null +++ b/Misc/NEWS.d/next/Library/2026-11-03-18-56-28.gh-issue-161259.iyZz6y.rst @@ -0,0 +1,3 @@ +Fix :func:`statistics.correlation` returning inaccurate results for very +small inputs, where the product of the sums of squares is subnormal rather +than zero. From 80afcf0d61aa73d7e09cbb8e0e91a23abd75fab0 Mon Sep 17 00:00:00 2001 From: Michael Date: Fri, 6 Nov 2026 19:24:10 +0100 Subject: [PATCH 2/4] refactor(statistics): move ambigious constant to a variable + note --- Lib/statistics.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Lib/statistics.py b/Lib/statistics.py index 0308e2578e2c88..aaf3799716d8fb 100644 --- a/Lib/statistics.py +++ b/Lib/statistics.py @@ -1776,7 +1776,8 @@ def _sqrtprod(x: float, y: float) -> float: return _sqrtprod(scale * x, scale * y) / scale return h - if h < 2.0 ** -485: + subnormal_cutoff = 2.0 ** -485 # sqrt(sys.float_info.min / sys.float_info.epsilon) + if h < subnormal_cutoff: # x * y is subnormal or underflowed, or the correction term # below would be. Rescale non-zero inputs. if x and y: From fc16f40abb457dbf7fe4f56f6b618ee0a53752c0 Mon Sep 17 00:00:00 2001 From: Michael Date: Fri, 6 Nov 2026 19:25:07 +0100 Subject: [PATCH 3/4] chore: reword changelog entry --- .../Library/2026-11-03-18-56-28.gh-issue-161259.iyZz6y.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Misc/NEWS.d/next/Library/2026-11-03-18-56-28.gh-issue-161259.iyZz6y.rst b/Misc/NEWS.d/next/Library/2026-11-03-18-56-28.gh-issue-161259.iyZz6y.rst index c59173bf3bae8f..400f6e5b716941 100644 --- a/Misc/NEWS.d/next/Library/2026-11-03-18-56-28.gh-issue-161259.iyZz6y.rst +++ b/Misc/NEWS.d/next/Library/2026-11-03-18-56-28.gh-issue-161259.iyZz6y.rst @@ -1,3 +1,3 @@ -Fix :func:`statistics.correlation` returning inaccurate results for very -small inputs, where the product of the sums of squares is subnormal rather +Improved the accuracy of :func:`statistics.correlation` for very small +inputs, where the product of the sums of squares is subnormal rather than zero. From e82fa8c3e6af3826bd7afe6d11325c43315327ca Mon Sep 17 00:00:00 2001 From: Michael Date: Fri, 6 Nov 2026 19:45:26 +0100 Subject: [PATCH 4/4] test(statistics): add subnormal correlation edge-case tests --- Lib/test/test_statistics.py | 45 +++++++++++++++++++++++++++++++++---- 1 file changed, 41 insertions(+), 4 deletions(-) diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py index 8c35cb9c8f7010..266d92828fe154 100644 --- a/Lib/test/test_statistics.py +++ b/Lib/test/test_statistics.py @@ -2790,11 +2790,12 @@ def test_sqrtprod_helper_function_fundamentals(self): self.assertEqual(actual, expected) self.assertEqual(sign(actual), sign(expected)) + @requires_IEEE_754 def test_sqrtprod_helper_function_subnormal_products(self): - # gh-161259: accuracy was lost when x * y is subnormal. - # The reference scales the inputs by an exact power of two so that - # the product is a normal float. - scale = 2.0**600 + # gh-161259: accuracy was lost when the product or the correction + # term was subnormal. The reference scales the inputs by an exact + # power of two so that the product is a normal float. + scale = 2.0 ** 600 for x, y in [ (2e-200, 2e-124), (1e-160, 3e-150), @@ -2806,6 +2807,23 @@ def test_sqrtprod_helper_function_subnormal_products(self): with self.subTest(x=x, y=y, expected=expected, actual=actual): self.assertAlmostEqual(actual / expected, 1.0, delta=2**-50) + @requires_IEEE_754 + def test_sqrtprod_helper_function_subnormal_boundary(self): + # sqrt(3 * 7) * 2**e straddles the 2**-485 cutoff where the + # correction term can become subnormal. + scale = 2.0 ** 600 + for e in (-488, -487): + x = 3.0 * 2.0 ** e + y = 7.0 * 2.0 ** e + expected = math.sqrt((scale * x) * (scale * y)) / scale + actual = statistics._sqrtprod(x, y) + with self.subTest(e=e, expected=expected, actual=actual): + self.assertAlmostEqual(actual / expected, 1.0, delta=2**-50) + + cutoff = 2.0 ** -485 + self.assertEqual(statistics._sqrtprod(cutoff, cutoff), cutoff) + + @requires_IEEE_754 def test_correlation_tiny_inputs(self): # gh-161259: correlation() is scale-invariant, even when the # intermediate sums of squares are tiny. @@ -2813,6 +2831,25 @@ def test_correlation_tiny_inputs(self): y = [1e-62, 2e-62, 3e-62] self.assertAlmostEqual(statistics.correlation(x, y), 1.0) + @requires_IEEE_754 + def test_correlation_asymmetric_scales(self): + for a, b in [ + (2.0 ** 50, 2.0 ** -537), + (1.3 * 2.0 ** -260, 2.7 * 2.0 ** -270), + ]: + x = [-a, 0.0, a] + y = [0.0, -b, b] + with self.subTest(a=a, b=b): + self.assertAlmostEqual(statistics.correlation(x, y), 0.5, + delta=2**-50) + + @requires_IEEE_754 + def test_correlation_tiny_products_max_recursion(self): + # Make both sums of squares round to the smallest subnormal. + a = 2.0 ** -537 + x = [0.0, 0.0, 0.0, a] + self.assertEqual(statistics.correlation(x, x), 1.0) + @requires_IEEE_754 @skip_if_double_rounding @support.cpython_only # Allow for a weaker sumprod() implementation