From 71f7dcd95cb9d8d05f270e51fff5e95722a10716 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Tue, 23 Nov 2021 15:18:54 -0700
Subject: [PATCH 01/20] Inlined code from variance functions

---
 Lib/statistics.py | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/Lib/statistics.py b/Lib/statistics.py
index 5c3f77df1549ddd..eec4afc0cacca42 100644
--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -840,7 +840,13 @@ def stdev(data, xbar=None):
     # Fixme: Despite the exact sum of squared deviations, some inaccuracy
     # remain because there are two rounding steps.  The first occurs in
     # the _convert() step for variance(), the second occurs in math.sqrt().
-    var = variance(data, xbar)
+    if iter(data) is data:
+        data = list(data)
+    n = len(data)
+    if n < 2:
+        raise StatisticsError('stdev requires at least two data points')
+    T, ss = _ss(data, xbar)
+    var = _convert(ss / (n - 1), T)
     try:
         return var.sqrt()
     except AttributeError:
@@ -859,7 +865,13 @@ def pstdev(data, mu=None):
     # Fixme: Despite the exact sum of squared deviations, some inaccuracy
     # remain because there are two rounding steps.  The first occurs in
     # the _convert() step for pvariance(), the second occurs in math.sqrt().
-    var = pvariance(data, mu)
+    if iter(data) is data:
+        data = list(data)
+    n = len(data)
+    if n < 1:
+        raise StatisticsError('pstdev requires at least one data point')
+    T, ss = _ss(data, mu)
+    var = _convert(ss / n, T)
     try:
         return var.sqrt()
     except AttributeError:

From 139e19f7dfc65649a0b80033aec2de054f342323 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Tue, 23 Nov 2021 15:21:47 -0700
Subject: [PATCH 02/20] Added helper functions for the float square root of a
 fraction

---
 Lib/statistics.py | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/Lib/statistics.py b/Lib/statistics.py
index eec4afc0cacca42..d194d5459202d06 100644
--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -304,6 +304,19 @@ def _fail_neg(values, errmsg='negative value'):
             raise StatisticsError(errmsg)
         yield x
 
+def _isqrt_frac_rto(n, m):
+    'Square root of n/m, rounded to the nearest integer using round-to-odd.'
+    a = math.isqrt(n*m) // m
+    return a | (a*a*m != n)
+
+def _sqrt_frac(n, m):
+    'Square root of n/m as a float, correctly rounded.'
+    q = (n.bit_length() - m.bit_length() - 109) // 2
+    if q >= 0:
+        return float(_isqrt_frac_rto(n, m << 2 * q) << q)
+    else:
+        return _isqrt_frac_rto(n << -2 * q, m) / (1 << -q)
+
 
 # === Measures of central tendency (averages) ===
 

From 06c2080835ac82d1675ebdb2ffd709b0f8451ae4 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Tue, 23 Nov 2021 15:31:08 -0700
Subject: [PATCH 03/20] Call helper functions

---
 Lib/statistics.py | 22 ++++++++--------------
 1 file changed, 8 insertions(+), 14 deletions(-)

diff --git a/Lib/statistics.py b/Lib/statistics.py
index d194d5459202d06..68c94e292f8092c 100644
--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -850,20 +850,17 @@ def stdev(data, xbar=None):
     1.0810874155219827
 
     """
-    # Fixme: Despite the exact sum of squared deviations, some inaccuracy
-    # remain because there are two rounding steps.  The first occurs in
-    # the _convert() step for variance(), the second occurs in math.sqrt().
     if iter(data) is data:
         data = list(data)
     n = len(data)
     if n < 2:
         raise StatisticsError('stdev requires at least two data points')
     T, ss = _ss(data, xbar)
-    var = _convert(ss / (n - 1), T)
-    try:
+    mss = ss / (n - 1)
+    if hasattr(T, 'sqrt'):
+        var = _convert(mss, T)
         return var.sqrt()
-    except AttributeError:
-        return math.sqrt(var)
+    return _sqrt_frac(mss.numerator, mss.denominator)
 
 
 def pstdev(data, mu=None):
@@ -875,20 +872,17 @@ def pstdev(data, mu=None):
     0.986893273527251
 
     """
-    # Fixme: Despite the exact sum of squared deviations, some inaccuracy
-    # remain because there are two rounding steps.  The first occurs in
-    # the _convert() step for pvariance(), the second occurs in math.sqrt().
     if iter(data) is data:
         data = list(data)
     n = len(data)
     if n < 1:
         raise StatisticsError('pstdev requires at least one data point')
     T, ss = _ss(data, mu)
-    var = _convert(ss / n, T)
-    try:
+    mss = ss / n
+    if hasattr(T, 'sqrt'):
+        var = _convert(mss, T)
         return var.sqrt()
-    except AttributeError:
-        return math.sqrt(var)
+    return _sqrt_frac(mss.numerator, mss.denominator)
 
 
 # === Statistics for relations between two inputs ===

From e7ff885c29de53219f3402c5f6319e4c2693242d Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Tue, 23 Nov 2021 15:37:05 -0700
Subject: [PATCH 04/20] Add blurb

---
 .../next/Library/2021-11-23-15-36-56.bpo-45876.NO8Yaj.rst       | 2 ++
 1 file changed, 2 insertions(+)
 create mode 100644 Misc/NEWS.d/next/Library/2021-11-23-15-36-56.bpo-45876.NO8Yaj.rst

diff --git a/Misc/NEWS.d/next/Library/2021-11-23-15-36-56.bpo-45876.NO8Yaj.rst b/Misc/NEWS.d/next/Library/2021-11-23-15-36-56.bpo-45876.NO8Yaj.rst
new file mode 100644
index 000000000000000..889ed6ce3ffb2f7
--- /dev/null
+++ b/Misc/NEWS.d/next/Library/2021-11-23-15-36-56.bpo-45876.NO8Yaj.rst
@@ -0,0 +1,2 @@
+Improve the accuracy of stdev() and pstdev() in the statistics module.  When
+the inputs are floats or fractions, the output is a correctly rounded float

From 34d59a5ea2428c51101523032e20efefc4c257d6 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Tue, 23 Nov 2021 17:10:37 -0700
Subject: [PATCH 05/20] Fix over-specified test

---
 Lib/test/test_statistics.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py
index c0e427d9355f257..41fc5c3092a68ce 100644
--- a/Lib/test/test_statistics.py
+++ b/Lib/test/test_statistics.py
@@ -2175,7 +2175,7 @@ def test_compare_to_variance(self):
         # Test that stdev is, in fact, the square root of variance.
         data = [random.uniform(-2, 9) for _ in range(1000)]
         expected = math.sqrt(statistics.variance(data))
-        self.assertEqual(self.func(data), expected)
+        self.assertAlmostEqual(self.func(data), expected)
 
     def test_center_not_at_mean(self):
         data = (1.0, 2.0)

From f6c8a97a0b4bedbb8dcdc674efbad9e50e16360c Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Tue, 23 Nov 2021 20:30:14 -0700
Subject: [PATCH 06/20] Add a test for the _sqrt_frac() helper function

---
 Lib/test/test_statistics.py | 27 +++++++++++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py
index 41fc5c3092a68ce..f3d475c13634c93 100644
--- a/Lib/test/test_statistics.py
+++ b/Lib/test/test_statistics.py
@@ -2161,6 +2161,33 @@ def test_center_not_at_mean(self):
         self.assertEqual(self.func(data), 2.5)
         self.assertEqual(self.func(data, mu=0.5), 6.5)
 
+class TestSqrtHelper(unittest.TestCase):
+
+    def test_sqrt_frac(self):
+        # This helper function aspires to produce more accurate square roots
+        # of fractional inputs that can be had with math.sqrt() alone.
+        # For a inputs spanning large ranges, we test that
+        # 1) the result is close to math.sqrt()
+        # 2) the result is at least as good as the two adjacent float values
+
+        randrange = random.randrange
+        sqrt_frac = statistics._sqrt_frac
+
+        for i in range(10_000):
+            numerator: int = randrange(10 ** randrange(30))
+            denonimator: int = randrange(10 ** randrange(30)) + 1
+            with self.subTest(numerator=numerator, denonimator=denonimator):
+                x: Fraction = Fraction(numerator, denonimator)
+
+                root: float = sqrt_frac(numerator, denonimator)
+                self.assertTrue(math.isclose(root, math.sqrt(numerator / denonimator)))
+
+                r_up: float = math.nextafter(root, math.inf)
+                self.assertLessEqual(abs(Fraction(root)**2 - x), abs(Fraction(r_up)**2 - x))
+
+                r_down: float = math.nextafter(root, -math.inf)
+                self.assertLessEqual(abs(Fraction(root)**2 - x), abs(Fraction(r_down)**2 - x))
+
 class TestStdev(VarianceStdevMixin, NumericTestCase):
     # Tests for sample standard deviation.
     def setUp(self):

From e439945426324af092bc38341d15727df4882de3 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Tue, 23 Nov 2021 21:04:40 -0700
Subject: [PATCH 07/20] Increase the tested range

---
 Lib/test/test_statistics.py | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py
index f3d475c13634c93..a9304b5a9ab1a69 100644
--- a/Lib/test/test_statistics.py
+++ b/Lib/test/test_statistics.py
@@ -2164,22 +2164,21 @@ def test_center_not_at_mean(self):
 class TestSqrtHelper(unittest.TestCase):
 
     def test_sqrt_frac(self):
-        # This helper function aspires to produce more accurate square roots
-        # of fractional inputs that can be had with math.sqrt() alone.
+        # This helper function aspires to produce correctly rounded square roots of
+        # fractional inputs, more accurate than can be had with math.sqrt() alone.
         # For a inputs spanning large ranges, we test that
         # 1) the result is close to math.sqrt()
         # 2) the result is at least as good as the two adjacent float values
 
         randrange = random.randrange
-        sqrt_frac = statistics._sqrt_frac
 
-        for i in range(10_000):
-            numerator: int = randrange(10 ** randrange(30))
-            denonimator: int = randrange(10 ** randrange(30)) + 1
+        for i in range(50_000):
+            numerator: int = randrange(10 ** randrange(40))
+            denonimator: int = randrange(10 ** randrange(40)) + 1
             with self.subTest(numerator=numerator, denonimator=denonimator):
                 x: Fraction = Fraction(numerator, denonimator)
 
-                root: float = sqrt_frac(numerator, denonimator)
+                root: float = statistics._sqrt_frac(numerator, denonimator)
                 self.assertTrue(math.isclose(root, math.sqrt(numerator / denonimator)))
 
                 r_up: float = math.nextafter(root, math.inf)

From 520e216793e6b5ba2b02a58c4960037e90668ebc Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Tue, 23 Nov 2021 21:12:58 -0700
Subject: [PATCH 08/20] Add type hints to the internal function.

---
 Lib/statistics.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Lib/statistics.py b/Lib/statistics.py
index 68c94e292f8092c..c501961fb92973d 100644
--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -304,14 +304,14 @@ def _fail_neg(values, errmsg='negative value'):
             raise StatisticsError(errmsg)
         yield x
 
-def _isqrt_frac_rto(n, m):
+def _isqrt_frac_rto(n: int, m: int) -> float:
     'Square root of n/m, rounded to the nearest integer using round-to-odd.'
     a = math.isqrt(n*m) // m
     return a | (a*a*m != n)
 
-def _sqrt_frac(n, m):
+def _sqrt_frac(n: int, m: int) -> float:
     'Square root of n/m as a float, correctly rounded.'
-    q = (n.bit_length() - m.bit_length() - 109) // 2
+    q: int = (n.bit_length() - m.bit_length() - 109) // 2
     if q >= 0:
         return float(_isqrt_frac_rto(n, m << 2 * q) << q)
     else:

From 8c9d78e0054d269ff6b2b292a0196c65f105047f Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Thu, 25 Nov 2021 10:57:02 -0700
Subject: [PATCH 09/20] Fix test for correct rounding

---
 Lib/test/test_statistics.py | 27 +++++++++++++++------------
 1 file changed, 15 insertions(+), 12 deletions(-)

diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py
index a9304b5a9ab1a69..908a3d01fe0eb4a 100644
--- a/Lib/test/test_statistics.py
+++ b/Lib/test/test_statistics.py
@@ -2164,11 +2164,20 @@ def test_center_not_at_mean(self):
 class TestSqrtHelper(unittest.TestCase):
 
     def test_sqrt_frac(self):
-        # This helper function aspires to produce correctly rounded square roots of
-        # fractional inputs, more accurate than can be had with math.sqrt() alone.
-        # For a inputs spanning large ranges, we test that
-        # 1) the result is close to math.sqrt()
-        # 2) the result is at least as good as the two adjacent float values
+
+        def is_root_correctly_rounded(x: Fraction, root: float) -> bool:
+            if not x:
+                return root == 0.0
+
+            r_up: float = math.nextafter(root, math.inf)
+            r_down: float = math.nextafter(root, -math.inf)
+            assert r_down < root < r_up
+
+            frac_root: Fraction = Fraction(root)
+            half_way_up: Fraction = (frac_root + Fraction(r_up)) / 2
+            half_way_down: Fraction = (frac_root + Fraction(r_down)) / 2
+
+            return half_way_down ** 2 <= x <= half_way_up ** 2
 
         randrange = random.randrange
 
@@ -2177,15 +2186,9 @@ def test_sqrt_frac(self):
             denonimator: int = randrange(10 ** randrange(40)) + 1
             with self.subTest(numerator=numerator, denonimator=denonimator):
                 x: Fraction = Fraction(numerator, denonimator)
-
                 root: float = statistics._sqrt_frac(numerator, denonimator)
-                self.assertTrue(math.isclose(root, math.sqrt(numerator / denonimator)))
-
-                r_up: float = math.nextafter(root, math.inf)
-                self.assertLessEqual(abs(Fraction(root)**2 - x), abs(Fraction(r_up)**2 - x))
+                self.assertTrue(is_root_correctly_rounded(x, root))
 
-                r_down: float = math.nextafter(root, -math.inf)
-                self.assertLessEqual(abs(Fraction(root)**2 - x), abs(Fraction(r_down)**2 - x))
 
 class TestStdev(VarianceStdevMixin, NumericTestCase):
     # Tests for sample standard deviation.

From d36ea9bd0ce49d868244fd9f083863fcb70b8107 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <rhettinger@users.noreply.github.com>
Date: Thu, 25 Nov 2021 11:04:32 -0700
Subject: [PATCH 10/20] =?UTF-8?q?Simplify=20=E2=8C=8A=E2=88=9A(n/m)?=
 =?UTF-8?q?=E2=8C=8B=20calculation?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Co-authored-by: Mark Dickinson <dickinsm@gmail.com>
---
 Lib/statistics.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Lib/statistics.py b/Lib/statistics.py
index c501961fb92973d..4e80d131777291d 100644
--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -306,7 +306,7 @@ def _fail_neg(values, errmsg='negative value'):
 
 def _isqrt_frac_rto(n: int, m: int) -> float:
     'Square root of n/m, rounded to the nearest integer using round-to-odd.'
-    a = math.isqrt(n*m) // m
+    a = math.isqrt(n // m)
     return a | (a*a*m != n)
 
 def _sqrt_frac(n: int, m: int) -> float:

From 28fdccf215c7e65fa4cb3038d0b634bcbbb29134 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Thu, 25 Nov 2021 13:27:38 -0700
Subject: [PATCH 11/20] Add comment and beef-up tests

---
 Lib/statistics.py           |  1 +
 Lib/test/test_statistics.py | 13 ++++++++++++-
 2 files changed, 13 insertions(+), 1 deletion(-)

diff --git a/Lib/statistics.py b/Lib/statistics.py
index 4e80d131777291d..3d7f29e86299bbc 100644
--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -311,6 +311,7 @@ def _isqrt_frac_rto(n: int, m: int) -> float:
 
 def _sqrt_frac(n: int, m: int) -> float:
     'Square root of n/m as a float, correctly rounded.'
+    # The constant 109 is:  3 + 2 * sys.float_info.mant_dig
     q: int = (n.bit_length() - m.bit_length() - 109) // 2
     if q >= 0:
         return float(_isqrt_frac_rto(n, m << 2 * q) << q)
diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py
index 908a3d01fe0eb4a..8fa03556dcee318 100644
--- a/Lib/test/test_statistics.py
+++ b/Lib/test/test_statistics.py
@@ -15,7 +15,7 @@
 import sys
 import unittest
 from test import support
-from test.support import import_helper
+from test.support import import_helper, requires_IEEE_754
 
 from decimal import Decimal
 from fractions import Fraction
@@ -2163,6 +2163,7 @@ def test_center_not_at_mean(self):
 
 class TestSqrtHelper(unittest.TestCase):
 
+    @requires_IEEE_754
     def test_sqrt_frac(self):
 
         def is_root_correctly_rounded(x: Fraction, root: float) -> bool:
@@ -2189,6 +2190,16 @@ def is_root_correctly_rounded(x: Fraction, root: float) -> bool:
                 root: float = statistics._sqrt_frac(numerator, denonimator)
                 self.assertTrue(is_root_correctly_rounded(x, root))
 
+        # Verify that corner cases and error handling match math.sqrt()
+        self.assertEqual(statistics._sqrt_frac(0, 1), 0.0)
+        with self.assertRaises(ValueError):
+            statistics._sqrt_frac(-1, 1)
+        with self.assertRaises(ValueError):
+            statistics._sqrt_frac(1, -1)
+
+        # The result is well defined if both inputs are negative
+        self.assertAlmostEqual(statistics._sqrt_frac(-2, -1), math.sqrt(2.0))
+
 
 class TestStdev(VarianceStdevMixin, NumericTestCase):
     # Tests for sample standard deviation.

From a17f2d9246e4297b044199dfad4d7c39f0cfe2a4 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Thu, 25 Nov 2021 13:37:12 -0700
Subject: [PATCH 12/20] Test for zero denominator

---
 Lib/test/test_statistics.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py
index 8fa03556dcee318..0c4643eb4c8cd30 100644
--- a/Lib/test/test_statistics.py
+++ b/Lib/test/test_statistics.py
@@ -2197,6 +2197,10 @@ def is_root_correctly_rounded(x: Fraction, root: float) -> bool:
         with self.assertRaises(ValueError):
             statistics._sqrt_frac(1, -1)
 
+        # Error handling for zero denominator matches that for Fraction(1, 0)
+        with self.assertRaises(ZeroDivisionError):
+            statistics._sqrt_frac(1, 0)
+
         # The result is well defined if both inputs are negative
         self.assertAlmostEqual(statistics._sqrt_frac(-2, -1), math.sqrt(2.0))
 

From 2aa88b02f23ecf4f0d6585323e3715fcbb848b62 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Thu, 25 Nov 2021 13:43:22 -0700
Subject: [PATCH 13/20] Add algorithmic references

---
 Lib/statistics.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/Lib/statistics.py b/Lib/statistics.py
index 3d7f29e86299bbc..50c126164760827 100644
--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -305,12 +305,14 @@ def _fail_neg(values, errmsg='negative value'):
         yield x
 
 def _isqrt_frac_rto(n: int, m: int) -> float:
-    'Square root of n/m, rounded to the nearest integer using round-to-odd.'
+    """Square root of n/m, rounded to the nearest integer using round-to-odd."""
+    # Refernce: https://www.lri.fr/~melquion/doc/05-imacs17_1-expose.pdf
     a = math.isqrt(n // m)
     return a | (a*a*m != n)
 
 def _sqrt_frac(n: int, m: int) -> float:
-    'Square root of n/m as a float, correctly rounded.'
+    """Square root of n/m as a float, correctly rounded."""
+    # See algorithm sketch at: https://bugs.python.org/msg406911
     # The constant 109 is:  3 + 2 * sys.float_info.mant_dig
     q: int = (n.bit_length() - m.bit_length() - 109) // 2
     if q >= 0:

From 64bff139b647875a880974d564aa93f47b69d7f7 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Thu, 25 Nov 2021 17:20:21 -0700
Subject: [PATCH 14/20] Add test for the _isqrt_frac_rto() helper function.

---
 Lib/test/test_statistics.py | 15 ++++++++++++++-
 1 file changed, 14 insertions(+), 1 deletion(-)

diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py
index 0c4643eb4c8cd30..8cd1bf033043e1f 100644
--- a/Lib/test/test_statistics.py
+++ b/Lib/test/test_statistics.py
@@ -9,6 +9,7 @@
 import copy
 import decimal
 import doctest
+import itertools
 import math
 import pickle
 import random
@@ -2161,7 +2162,15 @@ def test_center_not_at_mean(self):
         self.assertEqual(self.func(data), 2.5)
         self.assertEqual(self.func(data, mu=0.5), 6.5)
 
-class TestSqrtHelper(unittest.TestCase):
+class TestSqrtHelpers(unittest.TestCase):
+
+    def test_isqrt_frac_rto(self):
+        # For all fractions n/m, the root should be an odd number
+        # or an exact root.
+        for n, m in itertools.product(range(100), range(1, 100)):
+            r = statistics._isqrt_frac_rto(n, m)
+            self.assertIsInstance(r, int)
+            self.assertTrue(r&1 or r*r*m == n, (n, m))
 
     @requires_IEEE_754
     def test_sqrt_frac(self):
@@ -2170,14 +2179,18 @@ def is_root_correctly_rounded(x: Fraction, root: float) -> bool:
             if not x:
                 return root == 0.0
 
+            # Extract adjacent representable floats
             r_up: float = math.nextafter(root, math.inf)
             r_down: float = math.nextafter(root, -math.inf)
             assert r_down < root < r_up
 
+            # Convert to fractions for exact arithmetic
             frac_root: Fraction = Fraction(root)
             half_way_up: Fraction = (frac_root + Fraction(r_up)) / 2
             half_way_down: Fraction = (frac_root + Fraction(r_down)) / 2
 
+            # Check a closed interval.
+            # Does not test for a midpoint rounding rule.
             return half_way_down ** 2 <= x <= half_way_up ** 2
 
         randrange = random.randrange

From 4fde15799272604960c3a48ac8a18f044df96125 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Fri, 26 Nov 2021 09:40:46 -0700
Subject: [PATCH 15/20] Compute the 109 instead of hard-wiring it

---
 Lib/statistics.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/Lib/statistics.py b/Lib/statistics.py
index 50c126164760827..e4be5b8a1b55a63 100644
--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -130,6 +130,7 @@
 import math
 import numbers
 import random
+import sys
 
 from fractions import Fraction
 from decimal import Decimal
@@ -306,15 +307,17 @@ def _fail_neg(values, errmsg='negative value'):
 
 def _isqrt_frac_rto(n: int, m: int) -> float:
     """Square root of n/m, rounded to the nearest integer using round-to-odd."""
-    # Refernce: https://www.lri.fr/~melquion/doc/05-imacs17_1-expose.pdf
+    # Reference: https://www.lri.fr/~melquion/doc/05-imacs17_1-expose.pdf
     a = math.isqrt(n // m)
     return a | (a*a*m != n)
 
+# For 53 bit precision floats, the _sqrt_frac() shift is 109.
+_sqrt_shift: int = 2 * sys.float_info.mant_dig + 3
+
 def _sqrt_frac(n: int, m: int) -> float:
     """Square root of n/m as a float, correctly rounded."""
     # See algorithm sketch at: https://bugs.python.org/msg406911
-    # The constant 109 is:  3 + 2 * sys.float_info.mant_dig
-    q: int = (n.bit_length() - m.bit_length() - 109) // 2
+    q: int = (n.bit_length() - m.bit_length() - _sqrt_shift) // 2
     if q >= 0:
         return float(_isqrt_frac_rto(n, m << 2 * q) << q)
     else:

From acc58ff1d8a091f8b621a0b4dae5647e5e0fd02f Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Fri, 26 Nov 2021 10:00:44 -0700
Subject: [PATCH 16/20] Stronger test for _isqrt_frac_rto()

---
 Lib/test/test_statistics.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py
index 8cd1bf033043e1f..bf288794897613d 100644
--- a/Lib/test/test_statistics.py
+++ b/Lib/test/test_statistics.py
@@ -2165,12 +2165,16 @@ def test_center_not_at_mean(self):
 class TestSqrtHelpers(unittest.TestCase):
 
     def test_isqrt_frac_rto(self):
-        # For all fractions n/m, the root should be an odd number
-        # or an exact root.
         for n, m in itertools.product(range(100), range(1, 100)):
             r = statistics._isqrt_frac_rto(n, m)
             self.assertIsInstance(r, int)
-            self.assertTrue(r&1 or r*r*m == n, (n, m))
+            if r*r*m == n:
+                # Root is exact
+                continue
+            # Inexact, so the root should be odd
+            self.assertEqual(r&1, 1)
+            # Verify correct rounding
+            self.assertTrue(m * (r - 1)**2 < n < m * (r + 1)**2)
 
     @requires_IEEE_754
     def test_sqrt_frac(self):

From eba6e8d013d1dfc64e8dcad42a9d211a326e875f Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Fri, 26 Nov 2021 10:02:19 -0700
Subject: [PATCH 17/20] Bigger range

---
 Lib/test/test_statistics.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py
index bf288794897613d..983e202a9f9f506 100644
--- a/Lib/test/test_statistics.py
+++ b/Lib/test/test_statistics.py
@@ -2165,7 +2165,7 @@ def test_center_not_at_mean(self):
 class TestSqrtHelpers(unittest.TestCase):
 
     def test_isqrt_frac_rto(self):
-        for n, m in itertools.product(range(100), range(1, 100)):
+        for n, m in itertools.product(range(100), range(1, 1000)):
             r = statistics._isqrt_frac_rto(n, m)
             self.assertIsInstance(r, int)
             if r*r*m == n:

From a4354c7a9ab1c55dd0dcdbf42bdf413e4f1be44d Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Fri, 26 Nov 2021 10:21:36 -0700
Subject: [PATCH 18/20] Bigger range

---
 Lib/test/test_statistics.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Lib/test/test_statistics.py b/Lib/test/test_statistics.py
index 983e202a9f9f506..771a03e707ee01d 100644
--- a/Lib/test/test_statistics.py
+++ b/Lib/test/test_statistics.py
@@ -2199,9 +2199,9 @@ def is_root_correctly_rounded(x: Fraction, root: float) -> bool:
 
         randrange = random.randrange
 
-        for i in range(50_000):
-            numerator: int = randrange(10 ** randrange(40))
-            denonimator: int = randrange(10 ** randrange(40)) + 1
+        for i in range(60_000):
+            numerator: int = randrange(10 ** randrange(50))
+            denonimator: int = randrange(10 ** randrange(50)) + 1
             with self.subTest(numerator=numerator, denonimator=denonimator):
                 x: Fraction = Fraction(numerator, denonimator)
                 root: float = statistics._sqrt_frac(numerator, denonimator)

From 95330ee6ce1e49812c2e373b483c24406cf8badc Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Fri, 26 Nov 2021 10:28:18 -0700
Subject: [PATCH 19/20] Replace float() call with int/int division to be
 parallel with the other code path.

---
 Lib/statistics.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Lib/statistics.py b/Lib/statistics.py
index e4be5b8a1b55a63..4065bdc433fdb88 100644
--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -319,7 +319,7 @@ def _sqrt_frac(n: int, m: int) -> float:
     # See algorithm sketch at: https://bugs.python.org/msg406911
     q: int = (n.bit_length() - m.bit_length() - _sqrt_shift) // 2
     if q >= 0:
-        return float(_isqrt_frac_rto(n, m << 2 * q) << q)
+        return (_isqrt_frac_rto(n, m << 2 * q) << q) / 1
     else:
         return _isqrt_frac_rto(n << -2 * q, m) / (1 << -q)
 

From 14ff0d1cfe3a9a1aa66efdccc23be48f63e4e9d1 Mon Sep 17 00:00:00 2001
From: Raymond Hettinger <raymond.hettinger@gmail.com>
Date: Fri, 26 Nov 2021 14:37:05 -0700
Subject: [PATCH 20/20] Factor out division. Update proof link. Remove internal
 type declaration

---
 Lib/statistics.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/Lib/statistics.py b/Lib/statistics.py
index 4065bdc433fdb88..cf8eaa0a61e624f 100644
--- a/Lib/statistics.py
+++ b/Lib/statistics.py
@@ -316,12 +316,15 @@ def _isqrt_frac_rto(n: int, m: int) -> float:
 
 def _sqrt_frac(n: int, m: int) -> float:
     """Square root of n/m as a float, correctly rounded."""
-    # See algorithm sketch at: https://bugs.python.org/msg406911
-    q: int = (n.bit_length() - m.bit_length() - _sqrt_shift) // 2
+    # See principle and proof sketch at: https://bugs.python.org/msg407078
+    q = (n.bit_length() - m.bit_length() - _sqrt_shift) // 2
     if q >= 0:
-        return (_isqrt_frac_rto(n, m << 2 * q) << q) / 1
+        numerator = _isqrt_frac_rto(n, m << 2 * q) << q
+        denominator = 1
     else:
-        return _isqrt_frac_rto(n << -2 * q, m) / (1 << -q)
+        numerator = _isqrt_frac_rto(n << -2 * q, m)
+        denominator = 1 << -q
+    return numerator / denominator   # Convert to float
 
 
 # === Measures of central tendency (averages) ===