diff --git a/Lib/test/ieee754.txt b/Lib/test/ieee754.txt
new file mode 100644
index 0000000000..3e986cdb10
--- /dev/null
+++ b/Lib/test/ieee754.txt
@@ -0,0 +1,183 @@
+======================================
+Python IEEE 754 floating point support
+======================================
+
+>>> from sys import float_info as FI
+>>> from math import *
+>>> PI = pi
+>>> E = e
+
+You must never compare two floats with == because you are not going to get
+what you expect. We treat two floats as equal if the difference between them
+is small than epsilon.
+>>> EPS = 1E-15
+>>> def equal(x, y):
+...     """Almost equal helper for floats"""
+...     return abs(x - y) < EPS
+
+
+NaNs and INFs
+=============
+
+In Python 2.6 and newer NaNs (not a number) and infinity can be constructed
+from the strings 'inf' and 'nan'.
+
+>>> INF = float('inf')
+>>> NINF = float('-inf')
+>>> NAN = float('nan')
+
+>>> INF
+inf
+>>> NINF
+-inf
+>>> NAN
+nan
+
+The math module's ``isnan`` and ``isinf`` functions can be used to detect INF
+and NAN:
+>>> isinf(INF), isinf(NINF), isnan(NAN)
+(True, True, True)
+>>> INF == -NINF
+True
+
+Infinity
+--------
+
+Ambiguous operations like ``0 * inf`` or ``inf - inf`` result in NaN.
+>>> INF * 0
+nan
+>>> INF - INF
+nan
+>>> INF / INF
+nan
+
+However unambiguous operations with inf return inf:
+>>> INF * INF
+inf
+>>> 1.5 * INF
+inf
+>>> 0.5 * INF
+inf
+>>> INF / 1000
+inf
+
+Not a Number
+------------
+
+NaNs are never equal to another number, even itself
+>>> NAN == NAN
+False
+>>> NAN < 0
+False
+>>> NAN >= 0
+False
+
+All operations involving a NaN return a NaN except for nan**0 and 1**nan.
+>>> 1 + NAN
+nan
+>>> 1 * NAN
+nan
+>>> 0 * NAN
+nan
+>>> 1 ** NAN
+1.0
+>>> NAN ** 0
+1.0
+>>> 0 ** NAN
+nan
+>>> (1.0 + FI.epsilon) * NAN
+nan
+
+Misc Functions
+==============
+
+The power of 1 raised to x is always 1.0, even for special values like 0,
+infinity and NaN.
+
+>>> pow(1, 0)
+1.0
+>>> pow(1, INF)
+1.0
+>>> pow(1, -INF)
+1.0
+>>> pow(1, NAN)
+1.0
+
+The power of 0 raised to x is defined as 0, if x is positive. Negative
+finite values are a domain error or zero division error and NaN result in a
+silent NaN.
+
+>>> pow(0, 0)
+1.0
+>>> pow(0, INF)
+0.0
+>>> pow(0, -INF)
+inf
+>>> 0 ** -1
+Traceback (most recent call last):
+...
+ZeroDivisionError: 0.0 cannot be raised to a negative power
+>>> pow(0, NAN)
+nan
+
+
+Trigonometric Functions
+=======================
+
+>>> sin(INF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> sin(NINF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> sin(NAN)
+nan
+>>> cos(INF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> cos(NINF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> cos(NAN)
+nan
+>>> tan(INF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> tan(NINF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> tan(NAN)
+nan
+
+Neither pi nor tan are exact, but you can assume that tan(pi/2) is a large value
+and tan(pi) is a very small value:
+>>> tan(PI/2) > 1E10
+True
+>>> -tan(-PI/2) > 1E10
+True
+>>> tan(PI) < 1E-15
+True
+
+>>> asin(NAN), acos(NAN), atan(NAN)
+(nan, nan, nan)
+>>> asin(INF), asin(NINF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> acos(INF), acos(NINF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> equal(atan(INF), PI/2), equal(atan(NINF), -PI/2)
+(True, True)
+
+
+Hyberbolic Functions
+====================
+
diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py
index 2f64180652..fa79456ed4 100644
--- a/Lib/test/test_math.py
+++ b/Lib/test/test_math.py
@@ -4,6 +4,7 @@
 from test.support import verbose, requires_IEEE_754
 from test import support
 import unittest
+import fractions
 import itertools
 import decimal
 import math
@@ -186,6 +187,9 @@ def result_check(expected, got, ulp_tol=5, abs_tol=0.0):
 
     # Check exactly equal (applies also to strings representing exceptions)
     if got == expected:
+        if not got and not expected:
+            if math.copysign(1, got) != math.copysign(1, expected):
+                return f"expected {expected}, got {got} (zero has wrong sign)"
         return None
 
     failure = "not equal"
@@ -234,6 +238,10 @@ def __init__(self, value):
     def __index__(self):
         return self.value
 
+class BadDescr:
+    def __get__(self, obj, objtype=None):
+        raise ValueError
+
 class MathTests(unittest.TestCase):
 
     def ftest(self, name, got, expected, ulp_tol=5, abs_tol=0.0):
@@ -323,6 +331,7 @@ def testAtan2(self):
         self.ftest('atan2(0, 1)', math.atan2(0, 1), 0)
         self.ftest('atan2(1, 1)', math.atan2(1, 1), math.pi/4)
         self.ftest('atan2(1, 0)', math.atan2(1, 0), math.pi/2)
+        self.ftest('atan2(1, -1)', math.atan2(1, -1), 3*math.pi/4)
 
         # math.atan2(0, x)
         self.ftest('atan2(0., -inf)', math.atan2(0., NINF), math.pi)
@@ -416,16 +425,22 @@ def __ceil__(self):
                 return 42
         class TestNoCeil:
             pass
+        class TestBadCeil:
+            __ceil__ = BadDescr()
         self.assertEqual(math.ceil(TestCeil()), 42)
         self.assertEqual(math.ceil(FloatCeil()), 42)
         self.assertEqual(math.ceil(FloatLike(42.5)), 43)
         self.assertRaises(TypeError, math.ceil, TestNoCeil())
+        self.assertRaises(ValueError, math.ceil, TestBadCeil())
 
         t = TestNoCeil()
         t.__ceil__ = lambda *args: args
         self.assertRaises(TypeError, math.ceil, t)
         self.assertRaises(TypeError, math.ceil, t, 0)
 
+        self.assertEqual(math.ceil(FloatLike(+1.0)), +1.0)
+        self.assertEqual(math.ceil(FloatLike(-1.0)), -1.0)
+
     @requires_IEEE_754
     def testCopysign(self):
         self.assertEqual(math.copysign(1, 42), 1.0)
@@ -566,16 +581,22 @@ def __floor__(self):
                 return 42
         class TestNoFloor:
             pass
+        class TestBadFloor:
+            __floor__ = BadDescr()
         self.assertEqual(math.floor(TestFloor()), 42)
         self.assertEqual(math.floor(FloatFloor()), 42)
         self.assertEqual(math.floor(FloatLike(41.9)), 41)
         self.assertRaises(TypeError, math.floor, TestNoFloor())
+        self.assertRaises(ValueError, math.floor, TestBadFloor())
 
         t = TestNoFloor()
         t.__floor__ = lambda *args: args
         self.assertRaises(TypeError, math.floor, t)
         self.assertRaises(TypeError, math.floor, t, 0)
 
+        self.assertEqual(math.floor(FloatLike(+1.0)), +1.0)
+        self.assertEqual(math.floor(FloatLike(-1.0)), -1.0)
+
     def testFmod(self):
         self.assertRaises(TypeError, math.fmod)
         self.ftest('fmod(10, 1)', math.fmod(10, 1), 0.0)
@@ -597,6 +618,7 @@ def testFmod(self):
         self.assertEqual(math.fmod(-3.0, NINF), -3.0)
         self.assertEqual(math.fmod(0.0, 3.0), 0.0)
         self.assertEqual(math.fmod(0.0, NINF), 0.0)
+        self.assertRaises(ValueError, math.fmod, INF, INF)
 
     def testFrexp(self):
         self.assertRaises(TypeError, math.frexp)
@@ -638,7 +660,7 @@ def testFsum(self):
         def msum(iterable):
             """Full precision summation.  Compute sum(iterable) without any
             intermediate accumulation of error.  Based on the 'lsum' function
-            at http://code.activestate.com/recipes/393090/
+            at https://code.activestate.com/recipes/393090-binary-floating-point-summation-accurate-to-full-p/
 
             """
             tmant, texp = 0, 0
@@ -666,6 +688,7 @@ def msum(iterable):
             ([], 0.0),
             ([0.0], 0.0),
             ([1e100, 1.0, -1e100, 1e-100, 1e50, -1.0, -1e50], 1e-100),
+            ([1e100, 1.0, -1e100, 1e-100, 1e50, -1, -1e50], 1e-100),
             ([2.0**53, -0.5, -2.0**-54], 2.0**53-1.0),
             ([2.0**53, 1.0, 2.0**-100], 2.0**53+2.0),
             ([2.0**53+10.0, 1.0, 2.0**-100], 2.0**53+12.0),
@@ -713,6 +736,22 @@ def msum(iterable):
             s = msum(vals)
             self.assertEqual(msum(vals), math.fsum(vals))
 
+        self.assertEqual(math.fsum([1.0, math.inf]), math.inf)
+        self.assertTrue(math.isnan(math.fsum([math.nan, 1.0])))
+        self.assertEqual(math.fsum([1e100, FloatLike(1.0), -1e100, 1e-100,
+                                    1e50, FloatLike(-1.0), -1e50]), 1e-100)
+        self.assertRaises(OverflowError, math.fsum, [1e+308, 1e+308])
+        self.assertRaises(ValueError, math.fsum, [math.inf, -math.inf])
+        self.assertRaises(TypeError, math.fsum, ['spam'])
+        self.assertRaises(TypeError, math.fsum, 1)
+        self.assertRaises(OverflowError, math.fsum, [10**1000])
+
+        def bad_iter():
+            yield 1.0
+            raise ZeroDivisionError
+
+        self.assertRaises(ZeroDivisionError, math.fsum, bad_iter())
+
     def testGcd(self):
         gcd = math.gcd
         self.assertEqual(gcd(0, 0), 0)
@@ -773,9 +812,13 @@ def testHypot(self):
         # Test allowable types (those with __float__)
         self.assertEqual(hypot(12.0, 5.0), 13.0)
         self.assertEqual(hypot(12, 5), 13)
+        self.assertEqual(hypot(0.75, -1), 1.25)
+        self.assertEqual(hypot(-1, 0.75), 1.25)
+        self.assertEqual(hypot(0.75, FloatLike(-1.)), 1.25)
+        self.assertEqual(hypot(FloatLike(-1.), 0.75), 1.25)
         self.assertEqual(hypot(Decimal(12), Decimal(5)), 13)
         self.assertEqual(hypot(Fraction(12, 32), Fraction(5, 32)), Fraction(13, 32))
-        self.assertEqual(hypot(bool(1), bool(0), bool(1), bool(1)), math.sqrt(3))
+        self.assertEqual(hypot(True, False, True, True, True), 2.0)
 
         # Test corner cases
         self.assertEqual(hypot(0.0, 0.0), 0.0)     # Max input is zero
@@ -830,6 +873,8 @@ def testHypot(self):
             scale = FLOAT_MIN / 2.0 ** exp
             self.assertEqual(math.hypot(4*scale, 3*scale), 5*scale)
 
+        self.assertRaises(TypeError, math.hypot, *([1.0]*18), 'spam')
+
     @requires_IEEE_754
     @unittest.skipIf(HAVE_DOUBLE_ROUNDING,
                      "hypot() loses accuracy on machines with double rounding")
@@ -922,12 +967,16 @@ def testDist(self):
         # Test allowable types (those with __float__)
         self.assertEqual(dist((14.0, 1.0), (2.0, -4.0)), 13.0)
         self.assertEqual(dist((14, 1), (2, -4)), 13)
+        self.assertEqual(dist((FloatLike(14.), 1), (2, -4)), 13)
+        self.assertEqual(dist((11, 1), (FloatLike(-1.), -4)), 13)
+        self.assertEqual(dist((14, FloatLike(-1.)), (2, -6)), 13)
+        self.assertEqual(dist((14, -1), (2, -6)), 13)
         self.assertEqual(dist((D(14), D(1)), (D(2), D(-4))), D(13))
         self.assertEqual(dist((F(14, 32), F(1, 32)), (F(2, 32), F(-4, 32))),
                          F(13, 32))
-        self.assertEqual(dist((True, True, False, True, False),
-                              (True, False, True, True, False)),
-                         sqrt(2.0))
+        self.assertEqual(dist((True, True, False, False, True, True),
+                              (True, False, True, False, False, False)),
+                         2.0)
 
         # Test corner cases
         self.assertEqual(dist((13.25, 12.5, -3.25),
@@ -965,6 +1014,8 @@ class T(tuple):
             dist((1, 2, 3, 4), (5, 6, 7))
         with self.assertRaises(ValueError):        # Check dimension agree
             dist((1, 2, 3), (4, 5, 6, 7))
+        with self.assertRaises(TypeError):
+            dist((1,)*17 + ("spam",), (1,)*18)
         with self.assertRaises(TypeError):         # Rejects invalid types
             dist("abc", "xyz")
         int_too_big_for_float = 10 ** (sys.float_info.max_10_exp + 5)
@@ -972,6 +1023,16 @@ class T(tuple):
             dist((1, int_too_big_for_float), (2, 3))
         with self.assertRaises((ValueError, OverflowError)):
             dist((2, 3), (1, int_too_big_for_float))
+        with self.assertRaises(TypeError):
+            dist((1,), 2)
+        with self.assertRaises(TypeError):
+            dist([1], 2)
+
+        class BadFloat:
+            __float__ = BadDescr()
+
+        with self.assertRaises(ValueError):
+            dist([1], [BadFloat()])
 
         # Verify that the one dimensional case is equivalent to abs()
         for i in range(20):
@@ -1110,6 +1171,7 @@ def test_lcm(self):
 
     def testLdexp(self):
         self.assertRaises(TypeError, math.ldexp)
+        self.assertRaises(TypeError, math.ldexp, 2.0, 1.1)
         self.ftest('ldexp(0,1)', math.ldexp(0,1), 0)
         self.ftest('ldexp(1,1)', math.ldexp(1,1), 2)
         self.ftest('ldexp(1,-1)', math.ldexp(1,-1), 0.5)
@@ -1142,6 +1204,7 @@ def testLdexp(self):
 
     def testLog(self):
         self.assertRaises(TypeError, math.log)
+        self.assertRaises(TypeError, math.log, 1, 2, 3)
         self.ftest('log(1/e)', math.log(1/math.e), -1)
         self.ftest('log(1)', math.log(1), 0)
         self.ftest('log(e)', math.log(math.e), 1)
@@ -1152,6 +1215,7 @@ def testLog(self):
                    2302.5850929940457)
         self.assertRaises(ValueError, math.log, -1.5)
         self.assertRaises(ValueError, math.log, -10**1000)
+        self.assertRaises(ValueError, math.log, 10, -10)
         self.assertRaises(ValueError, math.log, NINF)
         self.assertEqual(math.log(INF), INF)
         self.assertTrue(math.isnan(math.log(NAN)))
@@ -1202,6 +1266,277 @@ def testLog10(self):
         self.assertEqual(math.log(INF), INF)
         self.assertTrue(math.isnan(math.log10(NAN)))
 
+    def testSumProd(self):
+        sumprod = math.sumprod
+        Decimal = decimal.Decimal
+        Fraction = fractions.Fraction
+
+        # Core functionality
+        self.assertEqual(sumprod(iter([10, 20, 30]), (1, 2, 3)), 140)
+        self.assertEqual(sumprod([1.5, 2.5], [3.5, 4.5]), 16.5)
+        self.assertEqual(sumprod([], []), 0)
+        self.assertEqual(sumprod([-1], [1.]), -1)
+        self.assertEqual(sumprod([1.], [-1]), -1)
+
+        # Type preservation and coercion
+        for v in [
+            (10, 20, 30),
+            (1.5, -2.5),
+            (Fraction(3, 5), Fraction(4, 5)),
+            (Decimal(3.5), Decimal(4.5)),
+            (2.5, 10),             # float/int
+            (2.5, Fraction(3, 5)), # float/fraction
+            (25, Fraction(3, 5)),  # int/fraction
+            (25, Decimal(4.5)),    # int/decimal
+        ]:
+            for p, q in [(v, v), (v, v[::-1])]:
+                with self.subTest(p=p, q=q):
+                    expected = sum(p_i * q_i for p_i, q_i in zip(p, q, strict=True))
+                    actual = sumprod(p, q)
+                    self.assertEqual(expected, actual)
+                    self.assertEqual(type(expected), type(actual))
+
+        # Bad arguments
+        self.assertRaises(TypeError, sumprod)               # No args
+        self.assertRaises(TypeError, sumprod, [])           # One arg
+        self.assertRaises(TypeError, sumprod, [], [], [])   # Three args
+        self.assertRaises(TypeError, sumprod, None, [10])   # Non-iterable
+        self.assertRaises(TypeError, sumprod, [10], None)   # Non-iterable
+        self.assertRaises(TypeError, sumprod, ['x'], [1.0])
+
+        # Uneven lengths
+        self.assertRaises(ValueError, sumprod, [10, 20], [30])
+        self.assertRaises(ValueError, sumprod, [10], [20, 30])
+
+        # Overflows
+        self.assertEqual(sumprod([10**20], [1]), 10**20)
+        self.assertEqual(sumprod([1], [10**20]), 10**20)
+        self.assertEqual(sumprod([10**10], [10**10]), 10**20)
+        self.assertEqual(sumprod([10**7]*10**5, [10**7]*10**5), 10**19)
+        self.assertRaises(OverflowError, sumprod, [10**1000], [1.0])
+        self.assertRaises(OverflowError, sumprod, [1.0], [10**1000])
+
+        # Error in iterator
+        def raise_after(n):
+            for i in range(n):
+                yield i
+            raise RuntimeError
+        with self.assertRaises(RuntimeError):
+            sumprod(range(10), raise_after(5))
+        with self.assertRaises(RuntimeError):
+            sumprod(raise_after(5), range(10))
+
+        from test.test_iter import BasicIterClass
+
+        self.assertEqual(sumprod(BasicIterClass(1), [1]), 0)
+        self.assertEqual(sumprod([1], BasicIterClass(1)), 0)
+
+        # Error in multiplication
+        class BadMultiply:
+            def __mul__(self, other):
+                raise RuntimeError
+            def __rmul__(self, other):
+                raise RuntimeError
+        with self.assertRaises(RuntimeError):
+            sumprod([10, BadMultiply(), 30], [1, 2, 3])
+        with self.assertRaises(RuntimeError):
+            sumprod([1, 2, 3], [10, BadMultiply(), 30])
+
+        # Error in addition
+        with self.assertRaises(TypeError):
+            sumprod(['abc', 3], [5, 10])
+        with self.assertRaises(TypeError):
+            sumprod([5, 10], ['abc', 3])
+
+        # Special values should give the same as the pure python recipe
+        self.assertEqual(sumprod([10.1, math.inf], [20.2, 30.3]), math.inf)
+        self.assertEqual(sumprod([10.1, math.inf], [math.inf, 30.3]), math.inf)
+        self.assertEqual(sumprod([10.1, math.inf], [math.inf, math.inf]), math.inf)
+        self.assertEqual(sumprod([10.1, -math.inf], [20.2, 30.3]), -math.inf)
+        self.assertTrue(math.isnan(sumprod([10.1, math.inf], [-math.inf, math.inf])))
+        self.assertTrue(math.isnan(sumprod([10.1, math.nan], [20.2, 30.3])))
+        self.assertTrue(math.isnan(sumprod([10.1, math.inf], [math.nan, 30.3])))
+        self.assertTrue(math.isnan(sumprod([10.1, math.inf], [20.3, math.nan])))
+
+        # Error cases that arose during development
+        args = ((-5, -5, 10), (1.5, 4611686018427387904, 2305843009213693952))
+        self.assertEqual(sumprod(*args), 0.0)
+
+
+    @requires_IEEE_754
+    @unittest.skipIf(HAVE_DOUBLE_ROUNDING,
+                         "sumprod() accuracy not guaranteed on machines with double rounding")
+    @support.cpython_only    # Other implementations may choose a different algorithm
+    def test_sumprod_accuracy(self):
+        sumprod = math.sumprod
+        self.assertEqual(sumprod([0.1] * 10, [1]*10), 1.0)
+        self.assertEqual(sumprod([0.1] * 20, [True, False] * 10), 1.0)
+        self.assertEqual(sumprod([True, False] * 10, [0.1] * 20), 1.0)
+        self.assertEqual(sumprod([1.0, 10E100, 1.0, -10E100], [1.0]*4), 2.0)
+
+    @support.requires_resource('cpu')
+    def test_sumprod_stress(self):
+        sumprod = math.sumprod
+        product = itertools.product
+        Decimal = decimal.Decimal
+        Fraction = fractions.Fraction
+
+        class Int(int):
+            def __add__(self, other):
+                return Int(int(self) + int(other))
+            def __mul__(self, other):
+                return Int(int(self) * int(other))
+            __radd__ = __add__
+            __rmul__ = __mul__
+            def __repr__(self):
+                return f'Int({int(self)})'
+
+        class Flt(float):
+            def __add__(self, other):
+                return Int(int(self) + int(other))
+            def __mul__(self, other):
+                return Int(int(self) * int(other))
+            __radd__ = __add__
+            __rmul__ = __mul__
+            def __repr__(self):
+                return f'Flt({int(self)})'
+
+        def baseline_sumprod(p, q):
+            """This defines the target behavior including exceptions and special values.
+            However, it is subject to rounding errors, so float inputs should be exactly
+            representable with only a few bits.
+            """
+            total = 0
+            for p_i, q_i in zip(p, q, strict=True):
+                total += p_i * q_i
+            return total
+
+        def run(func, *args):
+            "Make comparing functions easier. Returns error status, type, and result."
+            try:
+                result = func(*args)
+            except (AssertionError, NameError):
+                raise
+            except Exception as e:
+                return type(e), None, 'None'
+            return None, type(result), repr(result)
+
+        pools = [
+            (-5, 10, -2**20, 2**31, 2**40, 2**61, 2**62, 2**80, 1.5, Int(7)),
+            (5.25, -3.5, 4.75, 11.25, 400.5, 0.046875, 0.25, -1.0, -0.078125),
+            (-19.0*2**500, 11*2**1000, -3*2**1500, 17*2*333,
+               5.25, -3.25, -3.0*2**(-333),  3, 2**513),
+            (3.75, 2.5, -1.5, float('inf'), -float('inf'), float('NaN'), 14,
+                9, 3+4j, Flt(13), 0.0),
+            (13.25, -4.25, Decimal('10.5'), Decimal('-2.25'), Fraction(13, 8),
+                 Fraction(-11, 16), 4.75 + 0.125j, 97, -41, Int(3)),
+            (Decimal('6.125'), Decimal('12.375'), Decimal('-2.75'), Decimal(0),
+                 Decimal('Inf'), -Decimal('Inf'), Decimal('NaN'), 12, 13.5),
+            (-2.0 ** -1000, 11*2**1000, 3, 7, -37*2**32, -2*2**-537, -2*2**-538,
+                 2*2**-513),
+            (-7 * 2.0 ** -510, 5 * 2.0 ** -520, 17, -19.0, -6.25),
+            (11.25, -3.75, -0.625, 23.375, True, False, 7, Int(5)),
+        ]
+
+        for pool in pools:
+            for size in range(4):
+                for args1 in product(pool, repeat=size):
+                    for args2 in product(pool, repeat=size):
+                        args = (args1, args2)
+                        self.assertEqual(
+                            run(baseline_sumprod, *args),
+                            run(sumprod, *args),
+                            args,
+                        )
+
+    @requires_IEEE_754
+    @unittest.skipIf(HAVE_DOUBLE_ROUNDING,
+                         "sumprod() accuracy not guaranteed on machines with double rounding")
+    @support.cpython_only    # Other implementations may choose a different algorithm
+    @support.requires_resource('cpu')
+    def test_sumprod_extended_precision_accuracy(self):
+        import operator
+        from fractions import Fraction
+        from itertools import starmap
+        from collections import namedtuple
+        from math import log2, exp2, fabs
+        from random import choices, uniform, shuffle
+        from statistics import median
+
+        DotExample = namedtuple('DotExample', ('x', 'y', 'target_sumprod', 'condition'))
+
+        def DotExact(x, y):
+            vec1 = map(Fraction, x)
+            vec2 = map(Fraction, y)
+            return sum(starmap(operator.mul, zip(vec1, vec2, strict=True)))
+
+        def Condition(x, y):
+            return 2.0 * DotExact(map(abs, x), map(abs, y)) / abs(DotExact(x, y))
+
+        def linspace(lo, hi, n):
+            width = (hi - lo) / (n - 1)
+            return [lo + width * i for i in range(n)]
+
+        def GenDot(n, c):
+            """ Algorithm 6.1 (GenDot) works as follows. The condition number (5.7) of
+            the dot product xT y is proportional to the degree of cancellation. In
+            order to achieve a prescribed cancellation, we generate the first half of
+            the vectors x and y randomly within a large exponent range. This range is
+            chosen according to the anticipated condition number. The second half of x
+            and y is then constructed choosing xi randomly with decreasing exponent,
+            and calculating yi such that some cancellation occurs. Finally, we permute
+            the vectors x, y randomly and calculate the achieved condition number.
+            """
+
+            assert n >= 6
+            n2 = n // 2
+            x = [0.0] * n
+            y = [0.0] * n
+            b = log2(c)
+
+            # First half with exponents from 0 to |_b/2_| and random ints in between
+            e = choices(range(int(b/2)), k=n2)
+            e[0] = int(b / 2) + 1
+            e[-1] = 0.0
+
+            x[:n2] = [uniform(-1.0, 1.0) * exp2(p) for p in e]
+            y[:n2] = [uniform(-1.0, 1.0) * exp2(p) for p in e]
+
+            # Second half
+            e = list(map(round, linspace(b/2, 0.0 , n-n2)))
+            for i in range(n2, n):
+                x[i] = uniform(-1.0, 1.0) * exp2(e[i - n2])
+                y[i] = (uniform(-1.0, 1.0) * exp2(e[i - n2]) - DotExact(x, y)) / x[i]
+
+            # Shuffle
+            pairs = list(zip(x, y))
+            shuffle(pairs)
+            x, y = zip(*pairs)
+
+            return DotExample(x, y, DotExact(x, y), Condition(x, y))
+
+        def RelativeError(res, ex):
+            x, y, target_sumprod, condition = ex
+            n = DotExact(list(x) + [-res], list(y) + [1])
+            return fabs(n / target_sumprod)
+
+        def Trial(dotfunc, c, n):
+            ex = GenDot(10, c)
+            res = dotfunc(ex.x, ex.y)
+            return RelativeError(res, ex)
+
+        times = 1000          # Number of trials
+        n = 20                # Length of vectors
+        c = 1e30              # Target condition number
+
+        # If the following test fails, it means that the C math library
+        # implementation of fma() is not compliant with the C99 standard
+        # and is inaccurate.  To solve this problem, make a new build
+        # with the symbol UNRELIABLE_FMA defined.  That will enable a
+        # slower but accurate code path that avoids the fma() call.
+        relative_err = median(Trial(math.sumprod, c, n) for i in range(times))
+        self.assertLess(relative_err, 1e-16)
+
     def testModf(self):
         self.assertRaises(TypeError, math.modf)
 
@@ -1235,6 +1570,7 @@ def testPow(self):
         self.assertTrue(math.isnan(math.pow(2, NAN)))
         self.assertTrue(math.isnan(math.pow(0, NAN)))
         self.assertEqual(math.pow(1, NAN), 1)
+        self.assertRaises(OverflowError, math.pow, 1e+100, 1e+100)
 
         # pow(0., x)
         self.assertEqual(math.pow(0., INF), 0.)
@@ -1550,7 +1886,7 @@ def testTan(self):
         try:
             self.assertTrue(math.isnan(math.tan(INF)))
             self.assertTrue(math.isnan(math.tan(NINF)))
-        except:
+        except ValueError:
             self.assertRaises(ValueError, math.tan, INF)
             self.assertRaises(ValueError, math.tan, NINF)
         self.assertTrue(math.isnan(math.tan(NAN)))
@@ -1591,6 +1927,8 @@ def __trunc__(self):
                 return 23
         class TestNoTrunc:
             pass
+        class TestBadTrunc:
+            __trunc__ = BadDescr()
 
         self.assertEqual(math.trunc(TestTrunc()), 23)
         self.assertEqual(math.trunc(FloatTrunc()), 23)
@@ -1599,6 +1937,7 @@ class TestNoTrunc:
         self.assertRaises(TypeError, math.trunc, 1, 2)
         self.assertRaises(TypeError, math.trunc, FloatLike(23.5))
         self.assertRaises(TypeError, math.trunc, TestNoTrunc())
+        self.assertRaises(ValueError, math.trunc, TestBadTrunc())
 
     def testIsfinite(self):
         self.assertTrue(math.isfinite(0.0))
@@ -1626,11 +1965,11 @@ def testIsinf(self):
         self.assertFalse(math.isinf(0.))
         self.assertFalse(math.isinf(1.))
 
-    @requires_IEEE_754
     def test_nan_constant(self):
+        # `math.nan` must be a quiet NaN with positive sign bit
         self.assertTrue(math.isnan(math.nan))
+        self.assertEqual(math.copysign(1., math.nan), 1.)
 
-    @requires_IEEE_754
     def test_inf_constant(self):
         self.assertTrue(math.isinf(math.inf))
         self.assertGreater(math.inf, 0.0)
@@ -1719,6 +2058,13 @@ def test_testfile(self):
             except OverflowError:
                 result = 'OverflowError'
 
+            # C99+ says for math.h's sqrt: If the argument is +∞ or ±0, it is
+            # returned, unmodified.  On another hand, for csqrt: If z is ±0+0i,
+            # the result is +0+0i.  Lets correct zero sign of er to follow
+            # first convention.
+            if id in ['sqrt0002', 'sqrt0003', 'sqrt1001', 'sqrt1023']:
+                er = math.copysign(er, ar)
+
             # Default tolerances
             ulp_tol, abs_tol = 5, 0.0
 
@@ -1802,6 +2148,8 @@ def test_mtestfile(self):
                       '\n  '.join(failures))
 
     def test_prod(self):
+        from fractions import Fraction as F
+
         prod = math.prod
         self.assertEqual(prod([]), 1)
         self.assertEqual(prod([], start=5), 5)
@@ -1813,6 +2161,14 @@ def test_prod(self):
         self.assertEqual(prod([1.0, 2.0, 3.0, 4.0, 5.0]), 120.0)
         self.assertEqual(prod([1, 2, 3, 4.0, 5.0]), 120.0)
         self.assertEqual(prod([1.0, 2.0, 3.0, 4, 5]), 120.0)
+        self.assertEqual(prod([1., F(3, 2)]), 1.5)
+
+        # Error in multiplication
+        class BadMultiply:
+            def __rmul__(self, other):
+                raise RuntimeError
+        with self.assertRaises(RuntimeError):
+            prod([10., BadMultiply()])
 
         # Test overflow in fast-path for integers
         self.assertEqual(prod([1, 1, 2**32, 1, 1]), 2**32)
@@ -2044,11 +2400,20 @@ def test_nextafter(self):
                          float.fromhex('0x1.fffffffffffffp-1'))
         self.assertEqual(math.nextafter(1.0, INF),
                          float.fromhex('0x1.0000000000001p+0'))
+        self.assertEqual(math.nextafter(1.0, -INF, steps=1),
+                         float.fromhex('0x1.fffffffffffffp-1'))
+        self.assertEqual(math.nextafter(1.0, INF, steps=1),
+                         float.fromhex('0x1.0000000000001p+0'))
+        self.assertEqual(math.nextafter(1.0, -INF, steps=3),
+                         float.fromhex('0x1.ffffffffffffdp-1'))
+        self.assertEqual(math.nextafter(1.0, INF, steps=3),
+                         float.fromhex('0x1.0000000000003p+0'))
 
         # x == y: y is returned
-        self.assertEqual(math.nextafter(2.0, 2.0), 2.0)
-        self.assertEqualSign(math.nextafter(-0.0, +0.0), +0.0)
-        self.assertEqualSign(math.nextafter(+0.0, -0.0), -0.0)
+        for steps in range(1, 5):
+            self.assertEqual(math.nextafter(2.0, 2.0, steps=steps), 2.0)
+            self.assertEqualSign(math.nextafter(-0.0, +0.0, steps=steps), +0.0)
+            self.assertEqualSign(math.nextafter(+0.0, -0.0, steps=steps), -0.0)
 
         # around 0.0
         smallest_subnormal = sys.float_info.min * sys.float_info.epsilon
@@ -2073,6 +2438,11 @@ def test_nextafter(self):
         self.assertIsNaN(math.nextafter(1.0, NAN))
         self.assertIsNaN(math.nextafter(NAN, NAN))
 
+        self.assertEqual(1.0, math.nextafter(1.0, INF, steps=0))
+        with self.assertRaises(ValueError):
+            math.nextafter(1.0, INF, steps=-1)
+
+
     @requires_IEEE_754
     def test_ulp(self):
         self.assertEqual(math.ulp(1.0), sys.float_info.epsilon)
@@ -2112,6 +2482,14 @@ def __float__(self):
             # argument to a float.
             self.assertFalse(getattr(y, "converted", False))
 
+    def test_input_exceptions(self):
+        self.assertRaises(TypeError, math.exp, "spam")
+        self.assertRaises(TypeError, math.erf, "spam")
+        self.assertRaises(TypeError, math.atan2, "spam", 1.0)
+        self.assertRaises(TypeError, math.atan2, 1.0, "spam")
+        self.assertRaises(TypeError, math.atan2, 1.0)
+        self.assertRaises(TypeError, math.atan2, 1.0, 2.0, 3.0)
+
     # Custom assertions.
 
     def assertIsNaN(self, value):
@@ -2252,7 +2630,7 @@ def test_fractions(self):
 
 def load_tests(loader, tests, pattern):
     from doctest import DocFileSuite
-    # tests.addTest(DocFileSuite("ieee754.txt"))
+    tests.addTest(DocFileSuite("ieee754.txt"))
     return tests
 
 if __name__ == '__main__':
diff --git a/common/src/float_ops.rs b/common/src/float_ops.rs
index 69ae8833a2..8e16e71c14 100644
--- a/common/src/float_ops.rs
+++ b/common/src/float_ops.rs
@@ -133,6 +133,69 @@ pub fn nextafter(x: f64, y: f64) -> f64 {
     }
 }
 
+#[allow(clippy::float_cmp)]
+pub fn nextafter_with_steps(x: f64, y: f64, steps: u64) -> f64 {
+    if x == y {
+        y
+    } else if x.is_nan() || y.is_nan() {
+        f64::NAN
+    } else if x >= f64::INFINITY {
+        f64::MAX
+    } else if x <= f64::NEG_INFINITY {
+        f64::MIN
+    } else if x == 0.0 {
+        f64::from_bits(1).copysign(y)
+    } else {
+        if steps == 0 {
+            return x;
+        }
+
+        if x.is_nan() {
+            return x;
+        }
+
+        if y.is_nan() {
+            return y;
+        }
+
+        let sign_bit: u64 = 1 << 63;
+
+        let mut ux = x.to_bits();
+        let uy = y.to_bits();
+
+        let ax = ux & !sign_bit;
+        let ay = uy & !sign_bit;
+
+        // If signs are different
+        if ((ux ^ uy) & sign_bit) != 0 {
+            return if ax + ay <= steps {
+                f64::from_bits(uy)
+            } else if ax < steps {
+                let result = (uy & sign_bit) | (steps - ax);
+                f64::from_bits(result)
+            } else {
+                ux -= steps;
+                f64::from_bits(ux)
+            };
+        }
+
+        // If signs are the same
+        if ax > ay {
+            if ax - ay >= steps {
+                ux -= steps;
+                f64::from_bits(ux)
+            } else {
+                f64::from_bits(uy)
+            }
+        } else if ay - ax >= steps {
+            ux += steps;
+            f64::from_bits(ux)
+        } else {
+            f64::from_bits(uy)
+        }
+    }
+}
+
 pub fn ulp(x: f64) -> f64 {
     if x.is_nan() {
         return x;
diff --git a/stdlib/src/math.rs b/stdlib/src/math.rs
index 30b088f714..fc4e740e40 100644
--- a/stdlib/src/math.rs
+++ b/stdlib/src/math.rs
@@ -629,7 +629,7 @@ mod math {
 
     #[pyfunction]
     fn fsum(seq: ArgIterable<ArgIntoFloat>, vm: &VirtualMachine) -> PyResult<f64> {
-        let mut partials = vec![];
+        let mut partials = Vec::with_capacity(32);
         let mut special_sum = 0.0;
         let mut inf_sum = 0.0;
 
@@ -637,11 +637,11 @@ mod math {
             let mut x = *obj?;
 
             let xsave = x;
-            let mut j = 0;
+            let mut i = 0;
             // This inner loop applies `hi`/`lo` summation to each
             // partial so that the list of partial sums remains exact.
-            for i in 0..partials.len() {
-                let mut y: f64 = partials[i];
+            for j in 0..partials.len() {
+                let mut y: f64 = partials[j];
                 if x.abs() < y.abs() {
                     std::mem::swap(&mut x, &mut y);
                 }
@@ -650,33 +650,33 @@ mod math {
                 let hi = x + y;
                 let lo = y - (hi - x);
                 if lo != 0.0 {
-                    partials[j] = lo;
-                    j += 1;
+                    partials[i] = lo;
+                    i += 1;
                 }
                 x = hi;
             }
 
-            if !x.is_finite() {
-                // a nonfinite x could arise either as
-                // a result of intermediate overflow, or
-                // as a result of a nan or inf in the
-                // summands
-                if xsave.is_finite() {
-                    return Err(vm.new_overflow_error("intermediate overflow in fsum".to_owned()));
-                }
-                if xsave.is_infinite() {
-                    inf_sum += xsave;
+            partials.truncate(i);
+            if x != 0.0 {
+                if !x.is_finite() {
+                    // a nonfinite x could arise either as
+                    // a result of intermediate overflow, or
+                    // as a result of a nan or inf in the
+                    // summands
+                    if xsave.is_finite() {
+                        return Err(
+                            vm.new_overflow_error("intermediate overflow in fsum".to_owned())
+                        );
+                    }
+                    if xsave.is_infinite() {
+                        inf_sum += xsave;
+                    }
+                    special_sum += xsave;
+                    // reset partials
+                    partials.clear();
+                } else {
+                    partials.push(x);
                 }
-                special_sum += xsave;
-                // reset partials
-                partials.clear();
-            }
-
-            if j >= partials.len() {
-                partials.push(x);
-            } else {
-                partials[j] = x;
-                partials.truncate(j + 1);
             }
         }
         if special_sum != 0.0 {
@@ -831,9 +831,38 @@ mod math {
         (x.fract(), x.trunc())
     }
 
-    #[pyfunction]
-    fn nextafter(x: ArgIntoFloat, y: ArgIntoFloat) -> f64 {
-        float_ops::nextafter(*x, *y)
+    #[derive(FromArgs)]
+    struct NextAfterArgs {
+        #[pyarg(positional)]
+        x: ArgIntoFloat,
+        #[pyarg(positional)]
+        y: ArgIntoFloat,
+        #[pyarg(named, optional)]
+        steps: OptionalArg<ArgIndex>,
+    }
+
+    #[pyfunction]
+    fn nextafter(arg: NextAfterArgs, vm: &VirtualMachine) -> PyResult<f64> {
+        let steps: Option<i64> = arg
+            .steps
+            .map(|v| v.try_to_primitive(vm))
+            .transpose()?
+            .into_option();
+        match steps {
+            Some(steps) => {
+                if steps < 0 {
+                    return Err(
+                        vm.new_value_error("steps must be a non-negative integer".to_string())
+                    );
+                }
+                Ok(float_ops::nextafter_with_steps(
+                    *arg.x,
+                    *arg.y,
+                    steps as u64,
+                ))
+            }
+            None => Ok(float_ops::nextafter(*arg.x, *arg.y)),
+        }
     }
 
     #[pyfunction]
@@ -917,10 +946,38 @@ mod math {
         // refer: https://github.com/python/cpython/blob/main/Modules/mathmodule.c#L3093-L3193
         for obj in iter.iter(vm)? {
             let obj = obj?;
+            result = vm._mul(&result, &obj)?;
+        }
+
+        Ok(result)
+    }
 
-            result = vm
-                ._mul(&result, &obj)
-                .map_err(|_| vm.new_type_error("math type error".to_owned()))?;
+    #[pyfunction]
+    fn sumprod(
+        p: ArgIterable<PyObjectRef>,
+        q: ArgIterable<PyObjectRef>,
+        vm: &VirtualMachine,
+    ) -> PyResult<PyObjectRef> {
+        let mut p_iter = p.iter(vm)?;
+        let mut q_iter = q.iter(vm)?;
+        // We cannot just create a float because the iterator may contain
+        // anything as long as it supports __add__ and __mul__.
+        let mut result = vm.new_pyobj(0);
+        loop {
+            let m_p = p_iter.next();
+            let m_q = q_iter.next();
+            match (m_p, m_q) {
+                (Some(r_p), Some(r_q)) => {
+                    let p = r_p?;
+                    let q = r_q?;
+                    let tmp = vm._mul(&p, &q)?;
+                    result = vm._add(&result, &tmp)?;
+                }
+                (None, None) => break,
+                _ => {
+                    return Err(vm.new_value_error("Inputs are not the same length".to_string()));
+                }
+            }
         }
 
         Ok(result)
diff --git a/vm/src/builtins/float.rs b/vm/src/builtins/float.rs
index 1cd041b7b9..27fdff12ee 100644
--- a/vm/src/builtins/float.rs
+++ b/vm/src/builtins/float.rs
@@ -124,8 +124,8 @@ fn inner_divmod(v1: f64, v2: f64, vm: &VirtualMachine) -> PyResult<(f64, f64)> {
 
 pub fn float_pow(v1: f64, v2: f64, vm: &VirtualMachine) -> PyResult {
     if v1.is_zero() && v2.is_sign_negative() {
-        let msg = format!("{v1} cannot be raised to a negative power");
-        Err(vm.new_zero_division_error(msg))
+        let msg = "0.0 cannot be raised to a negative power";
+        Err(vm.new_zero_division_error(msg.to_owned()))
     } else if v1.is_sign_negative() && (v2.floor() - v2).abs() > f64::EPSILON {
         let v1 = Complex64::new(v1, 0.);
         let v2 = Complex64::new(v2, 0.);