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.