diff -r dca18f0ec280 Lib/statistics.py
--- a/Lib/statistics.py	Sat Oct 01 03:11:04 2016 +0000
+++ b/Lib/statistics.py	Sat Oct 01 13:14:00 2016 +0100
@@ -478,6 +478,55 @@
 assert type(_nth_root) is type(lambda: None)
 
 
+def _frexp_gen(n):
+    """Version of math.frexp that works for floats and arbitrary integers.
+
+    This version of frexp avoids the OverflowError that occurs for math.frexp
+    applied to integers outside the range of a float.
+
+    Examples
+    --------
+    >>> _frexp_gen(10)
+    (0.625, 4)
+    >>> _frexp_gen(16)
+    (0.5, 5)
+    >>> _frexp_gen(-10)
+    (-0.625, 4)
+    >>> _frexp_gen(2**1023)
+    (0.5, 1024)
+    >>> _frexp_gen(2**1024)
+    (0.5, 1025)
+    >>> _frexp_gen(2**1024 - 1)
+    (0.5, 1025)
+    >>> _frexp_gen(2**1024 - 2**970)  # round-ties-to-even check
+    (0.5, 1025)
+    >>> _frexp_gen(2**1024 + 2**971)  # round-ties-to-even check
+    (0.5, 1025)
+    >>> _frexp_gen(2**1024 + 3 * 2**971)  # round-ties-to-even check
+    (0.5000000000000002, 1025)
+    >>> _frexp_gen(2**1024 + 5 * 2**971)  # round-ties-to-even check
+    (0.5000000000000002, 1025)
+    >>> _frexp_gen(2**1024 + 7 * 2**971)  # round-ties-to-even check
+    (0.5000000000000004, 1025)
+    >>> _frexp_gen(0)
+    (0.0, 0)
+    >>> _frexp_gen(-2**1024)
+    (-0.5, 1025)
+    """
+    try:
+        return math.frexp(n)
+    except OverflowError:
+        # Does correct rounding assuming IEEE 754 binary64 format, but
+        # shouldn't give badly wrong results for other formats. For abs(n) in
+        # the range [2**53, 2**1024 - 2**970), the code below should give
+        # identical results to math.frexp(n). For the purposes of _product,
+        # the abs(m) == 1.0 check can be omitted - it's safe to simply return
+        # (m, e).
+        e = n.bit_length()
+        m = ((n >> e - 54) - (-n >> e - 54)) / (1 << 55)
+        return (m / 2.0, e + 1) if abs(m) == 1.0 else (m, e)
+
+
 def _product(values):
     """Return product of values as (exponent, mantissa)."""
     errmsg = 'mixed Decimal and float is not supported'
@@ -499,11 +548,11 @@
     #
     #   x1*x2 = 2**p1*m1 * 2**p2*m2 = 2**(p1+p2)*(m1*m2)
     #
-    mant, scale = 1, 0  #math.frexp(prod)  # FIXME
+    mant, scale = _frexp_gen(prod)
     for y in chain([x], values):
         if isinstance(y, Decimal):
             raise TypeError(errmsg)
-        m1, e1 = math.frexp(y)
+        m1, e1 = _frexp_gen(y)
         m2, e2 = math.frexp(mant)
         scale += (e1 + e2)
         mant = m1*m2