diff -r 5673cf852daa Lib/statistics.py
--- a/Lib/statistics.py	Sat Oct 01 01:06:51 2016 -0500
+++ b/Lib/statistics.py	Sat Oct 01 14:22:11 2016 +0100
@@ -87,6 +87,7 @@
 import decimal
 import math
 import numbers
+import sys
 
 from fractions import Fraction
 from decimal import Decimal
@@ -478,6 +479,53 @@
 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:
+        # Assuming IEEE 754 binary64 format, the code below gives
+        # correctly-rounded results (using round-ties-to-even, as usual),
+        # and should agree with math.frexp(n) for abs(n) in the range
+        # [2**53, 2**1024 - 2**970).
+        shift = n.bit_length() - sys.float_info.mant_dig - 1
+        m, e = math.frexp((n >> shift) - (-n >> shift))
+        return m, e + shift - 1
+
+
 def _product(values):
     """Return product of values as (exponent, mantissa)."""
     errmsg = 'mixed Decimal and float is not supported'
@@ -499,11 +547,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