diff --git a/.cspell.dict/python-more.txt b/.cspell.dict/python-more.txt
index 7ea660a5d..f2383f15f 100644
--- a/.cspell.dict/python-more.txt
+++ b/.cspell.dict/python-more.txt
@@ -127,6 +127,7 @@ isabstractmethod
 isbytes
 iscased
 isfinal
+issubnormal
 istext
 itemiterator
 itemsize
diff --git a/Cargo.lock b/Cargo.lock
index 5f43ee577..c31bf98ce 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -3169,8 +3169,7 @@ dependencies = [
 [[package]]
 name = "pymath"
 version = "0.2.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "bc10e50b7a1f2cc3887e983721cb51fc7574be0066c84bff3ef9e5c096e8d6d5"
+source = "git+https://github.com/RustPython/pymath.git?rev=dd1e9a4cf26d68caf3b549c09476f52e356e16a3#dd1e9a4cf26d68caf3b549c09476f52e356e16a3"
 dependencies = [
  "libc",
  "libm",
diff --git a/Cargo.toml b/Cargo.toml
index 4b84f3093..72f06ffd7 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -266,7 +266,7 @@ pem-rfc7468 = "1.0"
 pkcs8 = "0.11"
 proc-macro2 = "1.0.105"
 psm = "0.1"
-pymath = { version = "0.2.0", features = ["mul_add", "malachite-bigint", "complex"] }
+pymath = { git = "https://github.com/RustPython/pymath.git", rev = "dd1e9a4cf26d68caf3b549c09476f52e356e16a3", features = ["mul_add", "malachite-bigint", "complex"] }
 pyo3 = "0.29"
 quote = "1.0.45"
 radium = "1.1.1"
diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py
index 943448d69..c16c25970 100644
--- a/Lib/test/test_math.py
+++ b/Lib/test/test_math.py
@@ -17,6 +17,7 @@ import sys
 
 eps = 1E-05
 NAN = float('nan')
+NNAN = float('-nan')
 INF = float('inf')
 NINF = float('-inf')
 FLOAT_MAX = sys.float_info.max
@@ -54,56 +55,6 @@ def to_ulps(x):
     return n
 
 
-# Here's a pure Python version of the math.factorial algorithm, for
-# documentation and comparison purposes.
-#
-# Formula:
-#
-#   factorial(n) = factorial_odd_part(n) << (n - count_set_bits(n))
-#
-# where
-#
-#   factorial_odd_part(n) = product_{i >= 0} product_{0 < j <= n >> i; j odd} j
-#
-# The outer product above is an infinite product, but once i >= n.bit_length,
-# (n >> i) < 1 and the corresponding term of the product is empty.  So only the
-# finitely many terms for 0 <= i < n.bit_length() contribute anything.
-#
-# We iterate downwards from i == n.bit_length() - 1 to i == 0.  The inner
-# product in the formula above starts at 1 for i == n.bit_length(); for each i
-# < n.bit_length() we get the inner product for i from that for i + 1 by
-# multiplying by all j in {n >> i+1 < j <= n >> i; j odd}.  In Python terms,
-# this set is range((n >> i+1) + 1 | 1, (n >> i) + 1 | 1, 2).
-
-def count_set_bits(n):
-    """Number of '1' bits in binary expansion of a nonnnegative integer."""
-    return 1 + count_set_bits(n & n - 1) if n else 0
-
-def partial_product(start, stop):
-    """Product of integers in range(start, stop, 2), computed recursively.
-    start and stop should both be odd, with start <= stop.
-
-    """
-    numfactors = (stop - start) >> 1
-    if not numfactors:
-        return 1
-    elif numfactors == 1:
-        return start
-    else:
-        mid = (start + numfactors) | 1
-        return partial_product(start, mid) * partial_product(mid, stop)
-
-def py_factorial(n):
-    """Factorial of nonnegative integer n, via "Binary Split Factorial Formula"
-    described at http://www.luschny.de/math/factorial/binarysplitfact.html
-
-    """
-    inner = outer = 1
-    for i in reversed(range(n.bit_length())):
-        inner *= partial_product((n >> i + 1) + 1 | 1, (n >> i) + 1 | 1)
-        outer *= inner
-    return outer << (n - count_set_bits(n))
-
 def ulp_abs_check(expected, got, ulp_tol, abs_tol):
     """Given finite floats `expected` and `got`, check that they're
     approximately equal to within the given number of ulps or the
@@ -238,18 +189,21 @@ class MyIndexable(object):
     def __index__(self):
         return self.value
 
-# int subclass with broken arithmetic operators; implementations must
-# convert their arguments to exact ints instead of using these.
-class BadIntSubclass(int):
-    def _binop(self, other='ignored', mod=None):
-        return 42
-    __add__ = __radd__ = __sub__ = __rsub__ = _binop
-    __mul__ = __rmul__ = __mod__ = __rmod__ = _binop
-    __divmod__ = __rdivmod__ = __pow__ = __rpow__ = _binop
-    __floordiv__ = __rfloordiv__ = _binop
-    __lshift__ = __rlshift__ = __rshift__ = __rrshift__ = _binop
-    __and__ = __rand__ = __or__ = __ror__ = __xor__ = __rxor__ = _binop
-    __lt__ = __le__ = __gt__ = __ge__ = _binop
+class IndexableFloatLike:
+    def __init__(self, float_value, index_value):
+        self.float_value = float_value
+        self.index_value = index_value
+
+    def __float__(self):
+        if isinstance(self.float_value, BaseException):
+            raise self.float_value
+        return self.float_value
+
+    def __index__(self):
+        if isinstance(self.index_value, BaseException):
+            raise self.index_value
+        return self.index_value
+
 
 class BadDescr:
     def __get__(self, obj, objtype=None):
@@ -490,6 +444,19 @@ class MathTests(unittest.TestCase):
         # similarly, copysign(2., NAN) could be 2. or -2.
         self.assertEqual(abs(math.copysign(2., NAN)), 2.)
 
+    def test_signbit(self):
+        self.assertRaises(TypeError, math.signbit)
+        self.assertRaises(TypeError, math.signbit, '1.0')
+
+        # C11, §7.12.3.6 requires signbit() to return a nonzero value
+        # if and only if the sign of its argument value is negative,
+        # but in practice, we are only interested in a boolean value.
+        self.assertIsInstance(math.signbit(1.0), bool)
+
+        for arg in [0., 1., INF, NAN]:
+            self.assertFalse(math.signbit(arg))
+            self.assertTrue(math.signbit(-arg))
+
     def testCos(self):
         self.assertRaises(TypeError, math.cos)
         self.ftest('cos(-pi/2)', math.cos(-math.pi/2), 0, abs_tol=math.ulp(1))
@@ -548,33 +515,6 @@ class MathTests(unittest.TestCase):
         self.ftest('fabs(0)', math.fabs(0), 0)
         self.ftest('fabs(1)', math.fabs(1), 1)
 
-    def testFactorial(self):
-        self.assertEqual(math.factorial(0), 1)
-        total = 1
-        for i in range(1, 1000):
-            total *= i
-            self.assertEqual(math.factorial(i), total)
-            self.assertEqual(math.factorial(i), py_factorial(i))
-        self.assertRaises(ValueError, math.factorial, -1)
-        self.assertRaises(ValueError, math.factorial, -10**100)
-
-    def testFactorialNonIntegers(self):
-        self.assertRaises(TypeError, math.factorial, 5.0)
-        self.assertRaises(TypeError, math.factorial, 5.2)
-        self.assertRaises(TypeError, math.factorial, -1.0)
-        self.assertRaises(TypeError, math.factorial, -1e100)
-        self.assertRaises(TypeError, math.factorial, decimal.Decimal('5'))
-        self.assertRaises(TypeError, math.factorial, decimal.Decimal('5.2'))
-        self.assertRaises(TypeError, math.factorial, "5")
-
-    # Other implementations may place different upper bounds.
-    @support.cpython_only
-    def testFactorialHugeInputs(self):
-        # Currently raises OverflowError for inputs that are too large
-        # to fit into a C long.
-        self.assertRaises(OverflowError, math.factorial, 10**100)
-        self.assertRaises(TypeError, math.factorial, 1e100)
-
     def testFloor(self):
         self.assertRaises(TypeError, math.floor)
         self.assertEqual(int, type(math.floor(0.5)))
@@ -638,6 +578,92 @@ class MathTests(unittest.TestCase):
         self.assertEqual(math.fmod(0.0, NINF), 0.0)
         self.assertRaises(ValueError, math.fmod, INF, INF)
 
+    def test_fmax(self):
+        self.assertRaises(TypeError, math.fmax)
+        self.assertRaises(TypeError, math.fmax, 'x', 'y')
+
+        self.assertEqual(math.fmax(0., 0.), 0.)
+        self.assertEqual(math.fmax(1., 2.), 2.)
+        self.assertEqual(math.fmax(2., 1.), 2.)
+
+        self.assertEqual(math.fmax(+1., +0.), 1.)
+        self.assertEqual(math.fmax(+0., +1.), 1.)
+        self.assertEqual(math.fmax(+1., -0.), 1.)
+        self.assertEqual(math.fmax(-0., +1.), 1.)
+
+        self.assertEqual(math.fmax(-1., +0.), 0.)
+        self.assertEqual(math.fmax(+0., -1.), 0.)
+        self.assertEqual(math.fmax(-1., -0.), 0.)
+        self.assertEqual(math.fmax(-0., -1.), 0.)
+
+        for x in [NINF, -1., -0., 0., 1., INF]:
+            self.assertFalse(math.isnan(x))
+
+            with self.subTest(x=x, is_negative=math.copysign(1, x) < 0):
+                self.assertEqual(math.fmax(INF, x), INF)
+                self.assertEqual(math.fmax(x, INF), INF)
+                self.assertEqual(math.fmax(NINF, x), x)
+                self.assertEqual(math.fmax(x, NINF), x)
+
+    @requires_IEEE_754
+    def test_fmax_nans(self):
+        # When exactly one operand is NaN, the other is returned.
+        for x in [NINF, -1., -0., 0., 1., INF]:
+            with self.subTest(x=x, is_negative=math.copysign(1, x) < 0):
+                self.assertFalse(math.isnan(math.fmax(NAN, x)))
+                self.assertFalse(math.isnan(math.fmax(x, NAN)))
+                self.assertFalse(math.isnan(math.fmax(NNAN, x)))
+                self.assertFalse(math.isnan(math.fmax(x, NNAN)))
+        # When both operands are NaNs, fmax() returns NaN (see C11, F.10.9.2)
+        # whose sign is implementation-defined (see C11, F.10.0.3).
+        self.assertTrue(math.isnan(math.fmax(NAN, NAN)))
+        self.assertTrue(math.isnan(math.fmax(NNAN, NNAN)))
+        self.assertTrue(math.isnan(math.fmax(NAN, NNAN)))
+        self.assertTrue(math.isnan(math.fmax(NNAN, NAN)))
+
+    def test_fmin(self):
+        self.assertRaises(TypeError, math.fmin)
+        self.assertRaises(TypeError, math.fmin, 'x', 'y')
+
+        self.assertEqual(math.fmin(0., 0.), 0.)
+        self.assertEqual(math.fmin(1., 2.), 1.)
+        self.assertEqual(math.fmin(2., 1.), 1.)
+
+        self.assertEqual(math.fmin(+1., +0.), 0.)
+        self.assertEqual(math.fmin(+0., +1.), 0.)
+        self.assertEqual(math.fmin(+1., -0.), 0.)
+        self.assertEqual(math.fmin(-0., +1.), 0.)
+
+        self.assertEqual(math.fmin(-1., +0.), -1.)
+        self.assertEqual(math.fmin(+0., -1.), -1.)
+        self.assertEqual(math.fmin(-1., -0.), -1.)
+        self.assertEqual(math.fmin(-0., -1.), -1.)
+
+        for x in [NINF, -1., -0., 0., 1., INF]:
+            self.assertFalse(math.isnan(x))
+
+            with self.subTest(x=x, is_negative=math.copysign(1, x) < 0):
+                self.assertEqual(math.fmin(INF, x), x)
+                self.assertEqual(math.fmin(x, INF), x)
+                self.assertEqual(math.fmin(NINF, x), NINF)
+                self.assertEqual(math.fmin(x, NINF), NINF)
+
+    @requires_IEEE_754
+    def test_fmin_nans(self):
+        # When exactly one operand is NaN, the other is returned.
+        for x in [NINF, -1., -0., 0., 1., INF]:
+            with self.subTest(x=x, is_negative=math.copysign(1, x) < 0):
+                self.assertFalse(math.isnan(math.fmin(NAN, x)))
+                self.assertFalse(math.isnan(math.fmin(x, NAN)))
+                self.assertFalse(math.isnan(math.fmin(NNAN, x)))
+                self.assertFalse(math.isnan(math.fmin(x, NNAN)))
+        # When both operands are NaNs, fmin() returns NaN (see C11, F.10.9.3)
+        # whose sign is implementation-defined (see C11, F.10.0.3).
+        self.assertTrue(math.isnan(math.fmin(NAN, NAN)))
+        self.assertTrue(math.isnan(math.fmin(NNAN, NNAN)))
+        self.assertTrue(math.isnan(math.fmin(NAN, NNAN)))
+        self.assertTrue(math.isnan(math.fmin(NNAN, NAN)))
+
     def testFrexp(self):
         self.assertRaises(TypeError, math.frexp)
 
@@ -1090,78 +1116,6 @@ class MathTests(unittest.TestCase):
         with self.assertRaises(ValueError):
             math.dist([1, 2], [3, 4, 5])
 
-    def testIsqrt(self):
-        # Test a variety of inputs, large and small.
-        test_values = (
-            list(range(1000))
-            + list(range(10**6 - 1000, 10**6 + 1000))
-            + [2**e + i for e in range(60, 200) for i in range(-40, 40)]
-            + [3**9999, 10**5001]
-        )
-
-        for value in test_values:
-            with self.subTest(value=value):
-                s = math.isqrt(value)
-                self.assertIs(type(s), int)
-                self.assertLessEqual(s*s, value)
-                self.assertLess(value, (s+1)*(s+1))
-
-        # Negative values
-        with self.assertRaises(ValueError):
-            math.isqrt(-1)
-
-        # Integer-like things
-        s = math.isqrt(True)
-        self.assertIs(type(s), int)
-        self.assertEqual(s, 1)
-
-        s = math.isqrt(False)
-        self.assertIs(type(s), int)
-        self.assertEqual(s, 0)
-
-        class IntegerLike(object):
-            def __init__(self, value):
-                self.value = value
-
-            def __index__(self):
-                return self.value
-
-        s = math.isqrt(IntegerLike(1729))
-        self.assertIs(type(s), int)
-        self.assertEqual(s, 41)
-
-        # Overridden operators of an int subclass must not affect the
-        # result.
-        s = math.isqrt(BadIntSubclass(10**20))
-        self.assertIs(type(s), int)
-        self.assertEqual(s, 10**10)
-
-        s = math.isqrt(BadIntSubclass(10**20 - 1))
-        self.assertIs(type(s), int)
-        self.assertEqual(s, 10**10 - 1)
-
-        with self.assertRaises(ValueError):
-            math.isqrt(IntegerLike(-3))
-
-        # Non-integer-like things
-        bad_values = [
-            3.5, "a string", decimal.Decimal("3.5"), 3.5j,
-            100.0, -4.0,
-        ]
-        for value in bad_values:
-            with self.subTest(value=value):
-                with self.assertRaises(TypeError):
-                    math.isqrt(value)
-
-    @support.bigmemtest(2**32, memuse=0.85)
-    def test_isqrt_huge(self, size):
-        if size & 1:
-            size += 1
-        v = 1 << size
-        w = math.isqrt(v)
-        self.assertEqual(w.bit_length(), size // 2 + 1)
-        self.assertEqual(w.bit_count(), 1)
-
     def test_lcm(self):
         lcm = math.lcm
         self.assertEqual(lcm(0, 0), 0)
@@ -1256,13 +1210,32 @@ class MathTests(unittest.TestCase):
         self.ftest('log(10**40, 10**20)', math.log(10**40, 10**20), 2)
         self.ftest('log(10**1000)', math.log(10**1000),
                    2302.5850929940457)
+        self.ftest('log(10**2000, 10**1000)', math.log(10**2000, 10**1000), 2)
+        self.ftest('log(MyIndexable(32), MyIndexable(2))',
+                   math.log(MyIndexable(32), MyIndexable(2)), 5)
+        self.ftest('log(MyIndexable(10**1000))',
+                   math.log(MyIndexable(10**1000)),
+                   2302.5850929940457)
+        self.ftest('log(MyIndexable(10**2000), MyIndexable(10**1000))',
+                   math.log(MyIndexable(10**2000), MyIndexable(10**1000)),
+                   2)
+        self.assertRaises(ValueError, math.log, 0.0)
+        self.assertRaises(ValueError, math.log, 0)
+        self.assertRaises(ValueError, math.log, MyIndexable(0))
         self.assertRaises(ValueError, math.log, -1.5)
+        self.assertRaises(ValueError, math.log, -1)
+        self.assertRaises(ValueError, math.log, MyIndexable(-1))
         self.assertRaises(ValueError, math.log, -10**1000)
+        self.assertRaises(ValueError, math.log, MyIndexable(-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)))
 
+        self.assertEqual(math.log(IndexableFloatLike(math.e, 10**1000)), 1.0)
+        self.assertAlmostEqual(math.log(IndexableFloatLike(OverflowError(), 10**1000)),
+                               2302.5850929940457)
+
     def testLog1p(self):
         self.assertRaises(TypeError, math.log1p)
         for n in [2, 2**90, 2**300]:
@@ -1278,16 +1251,28 @@ class MathTests(unittest.TestCase):
         self.assertEqual(math.log2(1), 0.0)
         self.assertEqual(math.log2(2), 1.0)
         self.assertEqual(math.log2(4), 2.0)
+        self.assertEqual(math.log2(MyIndexable(4)), 2.0)
 
         # Large integer values
         self.assertEqual(math.log2(2**1023), 1023.0)
         self.assertEqual(math.log2(2**1024), 1024.0)
         self.assertEqual(math.log2(2**2000), 2000.0)
+        self.assertEqual(math.log2(MyIndexable(2**2000)), 2000.0)
 
+        self.assertRaises(ValueError, math.log2, 0.0)
+        self.assertRaises(ValueError, math.log2, 0)
+        self.assertRaises(ValueError, math.log2, MyIndexable(0))
         self.assertRaises(ValueError, math.log2, -1.5)
+        self.assertRaises(ValueError, math.log2, -1)
+        self.assertRaises(ValueError, math.log2, MyIndexable(-1))
+        self.assertRaises(ValueError, math.log2, -2**2000)
+        self.assertRaises(ValueError, math.log2, MyIndexable(-2**2000))
         self.assertRaises(ValueError, math.log2, NINF)
         self.assertTrue(math.isnan(math.log2(NAN)))
 
+        self.assertEqual(math.log2(IndexableFloatLike(8.0, 2**2000)), 3.0)
+        self.assertEqual(math.log2(IndexableFloatLike(OverflowError(), 2**2000)), 2000.0)
+
     @requires_IEEE_754
     # log2() is not accurate enough on Mac OS X Tiger (10.4)
     @support.requires_mac_ver(10, 5)
@@ -1303,12 +1288,24 @@ class MathTests(unittest.TestCase):
         self.ftest('log10(1)', math.log10(1), 0)
         self.ftest('log10(10)', math.log10(10), 1)
         self.ftest('log10(10**1000)', math.log10(10**1000), 1000.0)
+        self.ftest('log10(MyIndexable(10))', math.log10(MyIndexable(10)), 1)
+        self.ftest('log10(MyIndexable(10**1000))',
+                   math.log10(MyIndexable(10**1000)), 1000.0)
+        self.assertRaises(ValueError, math.log10, 0.0)
+        self.assertRaises(ValueError, math.log10, 0)
+        self.assertRaises(ValueError, math.log10, MyIndexable(0))
         self.assertRaises(ValueError, math.log10, -1.5)
+        self.assertRaises(ValueError, math.log10, -1)
+        self.assertRaises(ValueError, math.log10, MyIndexable(-1))
         self.assertRaises(ValueError, math.log10, -10**1000)
+        self.assertRaises(ValueError, math.log10, MyIndexable(-10**1000))
         self.assertRaises(ValueError, math.log10, NINF)
         self.assertEqual(math.log(INF), INF)
         self.assertTrue(math.isnan(math.log10(NAN)))
 
+        self.assertEqual(math.log10(IndexableFloatLike(100.0, 10**1000)), 2.0)
+        self.assertEqual(math.log10(IndexableFloatLike(OverflowError(), 10**1000)), 1000.0)
+
     @support.bigmemtest(2**32, memuse=0.2)
     def test_log_huge_integer(self, size):
         v = 1 << size
@@ -1412,7 +1409,6 @@ class MathTests(unittest.TestCase):
         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")
@@ -1999,6 +1995,28 @@ class MathTests(unittest.TestCase):
         self.assertFalse(math.isfinite(float("inf")))
         self.assertFalse(math.isfinite(float("-inf")))
 
+    def testIsnormal(self):
+        self.assertTrue(math.isnormal(1.25))
+        self.assertTrue(math.isnormal(-1.0))
+        self.assertFalse(math.isnormal(0.0))
+        self.assertFalse(math.isnormal(-0.0))
+        self.assertFalse(math.isnormal(INF))
+        self.assertFalse(math.isnormal(NINF))
+        self.assertFalse(math.isnormal(NAN))
+        self.assertFalse(math.isnormal(FLOAT_MIN/2))
+        self.assertFalse(math.isnormal(-FLOAT_MIN/2))
+
+    def testIssubnormal(self):
+        self.assertFalse(math.issubnormal(1.25))
+        self.assertFalse(math.issubnormal(-1.0))
+        self.assertFalse(math.issubnormal(0.0))
+        self.assertFalse(math.issubnormal(-0.0))
+        self.assertFalse(math.issubnormal(INF))
+        self.assertFalse(math.issubnormal(NINF))
+        self.assertFalse(math.issubnormal(NAN))
+        self.assertTrue(math.issubnormal(FLOAT_MIN/2))
+        self.assertTrue(math.issubnormal(-FLOAT_MIN/2))
+
     def testIsnan(self):
         self.assertTrue(math.isnan(float("nan")))
         self.assertTrue(math.isnan(float("-nan")))
@@ -2297,140 +2315,6 @@ class MathTests(unittest.TestCase):
         self.assertEqual(type(prod([1, decimal.Decimal(2.0), 3, 4, 5, 6])),
                          decimal.Decimal)
 
-    def testPerm(self):
-        perm = math.perm
-        factorial = math.factorial
-        # Test if factorial definition is satisfied
-        for n in range(500):
-            for k in (range(n + 1) if n < 100 else range(30) if n < 200 else range(10)):
-                self.assertEqual(perm(n, k),
-                                 factorial(n) // factorial(n - k))
-
-        # Test for Pascal's identity
-        for n in range(1, 100):
-            for k in range(1, n):
-                self.assertEqual(perm(n, k), perm(n - 1, k - 1) * k + perm(n - 1, k))
-
-        # Test corner cases
-        for n in range(1, 100):
-            self.assertEqual(perm(n, 0), 1)
-            self.assertEqual(perm(n, 1), n)
-            self.assertEqual(perm(n, n), factorial(n))
-
-        # Test one argument form
-        for n in range(20):
-            self.assertEqual(perm(n), factorial(n))
-            self.assertEqual(perm(n, None), factorial(n))
-
-        # Raises TypeError if any argument is non-integer or argument count is
-        # not 1 or 2
-        self.assertRaises(TypeError, perm, 10, 1.0)
-        self.assertRaises(TypeError, perm, 10, decimal.Decimal(1.0))
-        self.assertRaises(TypeError, perm, 10, "1")
-        self.assertRaises(TypeError, perm, 10.0, 1)
-        self.assertRaises(TypeError, perm, decimal.Decimal(10.0), 1)
-        self.assertRaises(TypeError, perm, "10", 1)
-
-        self.assertRaises(TypeError, perm)
-        self.assertRaises(TypeError, perm, 10, 1, 3)
-        self.assertRaises(TypeError, perm)
-
-        # Raises Value error if not k or n are negative numbers
-        self.assertRaises(ValueError, perm, -1, 1)
-        self.assertRaises(ValueError, perm, -2**1000, 1)
-        self.assertRaises(ValueError, perm, 1, -1)
-        self.assertRaises(ValueError, perm, 1, -2**1000)
-
-        # Returns zero if k is greater than n
-        self.assertEqual(perm(1, 2), 0)
-        self.assertEqual(perm(1, 2**1000), 0)
-
-        n = 2**1000
-        self.assertEqual(perm(n, 0), 1)
-        self.assertEqual(perm(n, 1), n)
-        self.assertEqual(perm(n, 2), n * (n-1))
-        if support.check_impl_detail(cpython=True):
-            self.assertRaises(OverflowError, perm, n, n)
-
-        for n, k in (True, True), (True, False), (False, False):
-            self.assertEqual(perm(n, k), 1)
-            self.assertIs(type(perm(n, k)), int)
-        self.assertEqual(perm(IntSubclass(5), IntSubclass(2)), 20)
-        self.assertEqual(perm(MyIndexable(5), MyIndexable(2)), 20)
-        for k in range(3):
-            self.assertIs(type(perm(IntSubclass(5), IntSubclass(k))), int)
-            self.assertIs(type(perm(MyIndexable(5), MyIndexable(k))), int)
-
-    def testComb(self):
-        comb = math.comb
-        factorial = math.factorial
-        # Test if factorial definition is satisfied
-        for n in range(500):
-            for k in (range(n + 1) if n < 100 else range(30) if n < 200 else range(10)):
-                self.assertEqual(comb(n, k), factorial(n)
-                    // (factorial(k) * factorial(n - k)))
-
-        # Test for Pascal's identity
-        for n in range(1, 100):
-            for k in range(1, n):
-                self.assertEqual(comb(n, k), comb(n - 1, k - 1) + comb(n - 1, k))
-
-        # Test corner cases
-        for n in range(100):
-            self.assertEqual(comb(n, 0), 1)
-            self.assertEqual(comb(n, n), 1)
-
-        for n in range(1, 100):
-            self.assertEqual(comb(n, 1), n)
-            self.assertEqual(comb(n, n - 1), n)
-
-        # Test Symmetry
-        for n in range(100):
-            for k in range(n // 2):
-                self.assertEqual(comb(n, k), comb(n, n - k))
-
-        # Raises TypeError if any argument is non-integer or argument count is
-        # not 2
-        self.assertRaises(TypeError, comb, 10, 1.0)
-        self.assertRaises(TypeError, comb, 10, decimal.Decimal(1.0))
-        self.assertRaises(TypeError, comb, 10, "1")
-        self.assertRaises(TypeError, comb, 10.0, 1)
-        self.assertRaises(TypeError, comb, decimal.Decimal(10.0), 1)
-        self.assertRaises(TypeError, comb, "10", 1)
-
-        self.assertRaises(TypeError, comb, 10)
-        self.assertRaises(TypeError, comb, 10, 1, 3)
-        self.assertRaises(TypeError, comb)
-
-        # Raises Value error if not k or n are negative numbers
-        self.assertRaises(ValueError, comb, -1, 1)
-        self.assertRaises(ValueError, comb, -2**1000, 1)
-        self.assertRaises(ValueError, comb, 1, -1)
-        self.assertRaises(ValueError, comb, 1, -2**1000)
-
-        # Returns zero if k is greater than n
-        self.assertEqual(comb(1, 2), 0)
-        self.assertEqual(comb(1, 2**1000), 0)
-
-        n = 2**1000
-        self.assertEqual(comb(n, 0), 1)
-        self.assertEqual(comb(n, 1), n)
-        self.assertEqual(comb(n, 2), n * (n-1) // 2)
-        self.assertEqual(comb(n, n), 1)
-        self.assertEqual(comb(n, n-1), n)
-        self.assertEqual(comb(n, n-2), n * (n-1) // 2)
-        if support.check_impl_detail(cpython=True):
-            self.assertRaises(OverflowError, comb, n, n//2)
-
-        for n, k in (True, True), (True, False), (False, False):
-            self.assertEqual(comb(n, k), 1)
-            self.assertIs(type(comb(n, k)), int)
-        self.assertEqual(comb(IntSubclass(5), IntSubclass(2)), 10)
-        self.assertEqual(comb(MyIndexable(5), MyIndexable(2)), 10)
-        for k in range(3):
-            self.assertIs(type(comb(IntSubclass(5), IntSubclass(k))), int)
-            self.assertIs(type(comb(MyIndexable(5), MyIndexable(k))), int)
-
     @requires_IEEE_754
     def test_nextafter(self):
         # around 2^52 and 2^63
@@ -2490,7 +2374,6 @@ class MathTests(unittest.TestCase):
         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)
@@ -2805,11 +2688,8 @@ class FMATests(unittest.TestCase):
     @unittest.skipIf(
         sys.platform.startswith(("freebsd", "wasi", "netbsd", "emscripten"))
         or (sys.platform == "android" and platform.machine() == "x86_64")
-        or support.linked_to_musl(),  # gh-131032
+        or (support.linked_to_musl() and support.linked_to_musl() < (1, 2, 6)),  # gh-131032
         f"this platform doesn't implement IEE 754-2008 properly")
-    # gh-131032: musl is fixed but the fix is not yet released; when the fixed
-    # version is known change this to:
-    #   or support.linked_to_musl() < (1, <m>, <p>)
     def test_fma_zero_result(self):
         nonnegative_finites = [0.0, 1e-300, 2.3, 1e300]
 
diff --git a/Lib/test/test_math_integer.py b/Lib/test/test_math_integer.py
new file mode 100644
index 000000000..d1bc776dc
--- /dev/null
+++ b/Lib/test/test_math_integer.py
@@ -0,0 +1,421 @@
+from decimal import Decimal
+from fractions import Fraction
+import unittest
+from test import support
+
+
+class IntSubclass(int):
+    pass
+
+# Class providing an __index__ method.
+class MyIndexable(object):
+    def __init__(self, value):
+        self.value = value
+
+    def __index__(self):
+        return self.value
+
+# int subclass with broken arithmetic operators; implementations must
+# convert their arguments to exact ints instead of using these.
+class BadIntSubclass(int):
+    def _binop(self, other='ignored', mod=None):
+        return 42
+    __add__ = __radd__ = __sub__ = __rsub__ = _binop
+    __mul__ = __rmul__ = __mod__ = __rmod__ = _binop
+    __divmod__ = __rdivmod__ = __pow__ = __rpow__ = _binop
+    __floordiv__ = __rfloordiv__ = _binop
+    __lshift__ = __rlshift__ = __rshift__ = __rrshift__ = _binop
+    __and__ = __rand__ = __or__ = __ror__ = __xor__ = __rxor__ = _binop
+    __lt__ = __le__ = __gt__ = __ge__ = _binop
+
+# Here's a pure Python version of the math.integer.factorial algorithm, for
+# documentation and comparison purposes.
+#
+# Formula:
+#
+#   factorial(n) = factorial_odd_part(n) << (n - count_set_bits(n))
+#
+# where
+#
+#   factorial_odd_part(n) = product_{i >= 0} product_{0 < j <= n >> i; j odd} j
+#
+# The outer product above is an infinite product, but once i >= n.bit_length,
+# (n >> i) < 1 and the corresponding term of the product is empty.  So only the
+# finitely many terms for 0 <= i < n.bit_length() contribute anything.
+#
+# We iterate downwards from i == n.bit_length() - 1 to i == 0.  The inner
+# product in the formula above starts at 1 for i == n.bit_length(); for each i
+# < n.bit_length() we get the inner product for i from that for i + 1 by
+# multiplying by all j in {n >> i+1 < j <= n >> i; j odd}.  In Python terms,
+# this set is range((n >> i+1) + 1 | 1, (n >> i) + 1 | 1, 2).
+
+def count_set_bits(n):
+    """Number of '1' bits in binary expansion of a nonnnegative integer."""
+    return 1 + count_set_bits(n & n - 1) if n else 0
+
+def partial_product(start, stop):
+    """Product of integers in range(start, stop, 2), computed recursively.
+    start and stop should both be odd, with start <= stop.
+
+    """
+    numfactors = (stop - start) >> 1
+    if not numfactors:
+        return 1
+    elif numfactors == 1:
+        return start
+    else:
+        mid = (start + numfactors) | 1
+        return partial_product(start, mid) * partial_product(mid, stop)
+
+def py_factorial(n):
+    """Factorial of nonnegative integer n, via "Binary Split Factorial Formula"
+    described at http://www.luschny.de/math/factorial/binarysplitfact.html
+
+    """
+    inner = outer = 1
+    for i in reversed(range(n.bit_length())):
+        inner *= partial_product((n >> i + 1) + 1 | 1, (n >> i) + 1 | 1)
+        outer *= inner
+    return outer << (n - count_set_bits(n))
+
+
+class IntMathTests(unittest.TestCase):
+    import math.integer as module
+
+    def assertIntEqual(self, actual, expected):
+        self.assertEqual(actual, expected)
+        self.assertIs(type(actual), int)
+
+    def test_factorial(self):
+        factorial = self.module.factorial
+        self.assertEqual(factorial(0), 1)
+        total = 1
+        for i in range(1, 1000):
+            total *= i
+            self.assertEqual(factorial(i), total)
+            self.assertEqual(factorial(i), py_factorial(i))
+
+        self.assertIntEqual(factorial(False), 1)
+        self.assertIntEqual(factorial(True), 1)
+        for i in range(3):
+            expected = factorial(i)
+            self.assertIntEqual(factorial(IntSubclass(i)), expected)
+            self.assertIntEqual(factorial(MyIndexable(i)), expected)
+
+        self.assertRaises(ValueError, factorial, -1)
+        self.assertRaises(ValueError, factorial, -10**1000)
+
+    def test_factorial_non_integers(self):
+        factorial = self.module.factorial
+        self.assertRaises(TypeError, factorial, 5.0)
+        self.assertRaises(TypeError, factorial, 5.2)
+        self.assertRaises(TypeError, factorial, -1.0)
+        self.assertRaises(TypeError, factorial, -1e100)
+        self.assertRaises(TypeError, factorial, Decimal('5'))
+        self.assertRaises(TypeError, factorial, Decimal('5.2'))
+        self.assertRaises(TypeError, factorial, Fraction(5, 1))
+        self.assertRaises(TypeError, factorial, "5")
+
+    # Other implementations may place different upper bounds.
+    @support.cpython_only
+    def test_factorial_huge_inputs(self):
+        factorial = self.module.factorial
+        # Currently raises OverflowError for inputs that are too large
+        # to fit into a C long.
+        self.assertRaises(OverflowError, factorial, 10**100)
+        self.assertRaises(TypeError, factorial, 1e100)
+
+    def test_gcd(self):
+        gcd = self.module.gcd
+        self.assertEqual(gcd(0, 0), 0)
+        self.assertEqual(gcd(1, 0), 1)
+        self.assertEqual(gcd(-1, 0), 1)
+        self.assertEqual(gcd(0, 1), 1)
+        self.assertEqual(gcd(0, -1), 1)
+        self.assertEqual(gcd(7, 1), 1)
+        self.assertEqual(gcd(7, -1), 1)
+        self.assertEqual(gcd(-23, 15), 1)
+        self.assertEqual(gcd(120, 84), 12)
+        self.assertEqual(gcd(84, -120), 12)
+        self.assertEqual(gcd(1216342683557601535506311712,
+                             436522681849110124616458784), 32)
+        c = 652560
+        x = 434610456570399902378880679233098819019853229470286994367836600566
+        y = 1064502245825115327754847244914921553977
+        a = x * c
+        b = y * c
+        self.assertEqual(gcd(a, b), c)
+        self.assertEqual(gcd(b, a), c)
+        self.assertEqual(gcd(-a, b), c)
+        self.assertEqual(gcd(b, -a), c)
+        self.assertEqual(gcd(a, -b), c)
+        self.assertEqual(gcd(-b, a), c)
+        self.assertEqual(gcd(-a, -b), c)
+        self.assertEqual(gcd(-b, -a), c)
+        c = 576559230871654959816130551884856912003141446781646602790216406874
+        a = x * c
+        b = y * c
+        self.assertEqual(gcd(a, b), c)
+        self.assertEqual(gcd(b, a), c)
+        self.assertEqual(gcd(-a, b), c)
+        self.assertEqual(gcd(b, -a), c)
+        self.assertEqual(gcd(a, -b), c)
+        self.assertEqual(gcd(-b, a), c)
+        self.assertEqual(gcd(-a, -b), c)
+        self.assertEqual(gcd(-b, -a), c)
+
+        self.assertRaises(TypeError, gcd, 120.0, 84)
+        self.assertRaises(TypeError, gcd, 120, 84.0)
+        self.assertIntEqual(gcd(IntSubclass(120), IntSubclass(84)), 12)
+        self.assertIntEqual(gcd(MyIndexable(120), MyIndexable(84)), 12)
+
+    def test_lcm(self):
+        lcm = self.module.lcm
+        self.assertEqual(lcm(0, 0), 0)
+        self.assertEqual(lcm(1, 0), 0)
+        self.assertEqual(lcm(-1, 0), 0)
+        self.assertEqual(lcm(0, 1), 0)
+        self.assertEqual(lcm(0, -1), 0)
+        self.assertEqual(lcm(7, 1), 7)
+        self.assertEqual(lcm(7, -1), 7)
+        self.assertEqual(lcm(-23, 15), 345)
+        self.assertEqual(lcm(120, 84), 840)
+        self.assertEqual(lcm(84, -120), 840)
+        self.assertEqual(lcm(1216342683557601535506311712,
+                             436522681849110124616458784),
+                             16592536571065866494401400422922201534178938447014944)
+
+        x = 43461045657039990237
+        y = 10645022458251153277
+        for c in (652560,
+                  57655923087165495981):
+            a = x * c
+            b = y * c
+            d = x * y * c
+            self.assertEqual(lcm(a, b), d)
+            self.assertEqual(lcm(b, a), d)
+            self.assertEqual(lcm(-a, b), d)
+            self.assertEqual(lcm(b, -a), d)
+            self.assertEqual(lcm(a, -b), d)
+            self.assertEqual(lcm(-b, a), d)
+            self.assertEqual(lcm(-a, -b), d)
+            self.assertEqual(lcm(-b, -a), d)
+
+        self.assertEqual(lcm(), 1)
+        self.assertEqual(lcm(120), 120)
+        self.assertEqual(lcm(-120), 120)
+        self.assertEqual(lcm(120, 84, 102), 14280)
+        self.assertEqual(lcm(120, 0, 84), 0)
+
+        self.assertRaises(TypeError, lcm, 120.0)
+        self.assertRaises(TypeError, lcm, 120.0, 84)
+        self.assertRaises(TypeError, lcm, 120, 84.0)
+        self.assertRaises(TypeError, lcm, 120, 0, 84.0)
+        self.assertEqual(lcm(MyIndexable(120), MyIndexable(84)), 840)
+
+    def test_isqrt(self):
+        isqrt = self.module.isqrt
+        # Test a variety of inputs, large and small.
+        test_values = (
+            list(range(1000))
+            + list(range(10**6 - 1000, 10**6 + 1000))
+            + [2**e + i for e in range(60, 200) for i in range(-40, 40)]
+            + [3**9999, 10**5001]
+        )
+
+        for value in test_values:
+            with self.subTest(value=value):
+                s = isqrt(value)
+                self.assertIs(type(s), int)
+                self.assertLessEqual(s*s, value)
+                self.assertLess(value, (s+1)*(s+1))
+
+        # Negative values
+        with self.assertRaises(ValueError):
+            isqrt(-1)
+
+        # Integer-like things
+        self.assertIntEqual(isqrt(True), 1)
+        self.assertIntEqual(isqrt(False), 0)
+        self.assertIntEqual(isqrt(MyIndexable(1729)), 41)
+
+        # Overridden operators of an int subclass must not affect the
+        # result.
+        self.assertIntEqual(isqrt(BadIntSubclass(10**20)), 10**10)
+        self.assertIntEqual(isqrt(BadIntSubclass(10**20 - 1)), 10**10 - 1)
+
+        with self.assertRaises(ValueError):
+            isqrt(MyIndexable(-3))
+
+        # Non-integer-like things
+        bad_values = [
+            3.5, "a string", Decimal("3.5"), 3.5j,
+            100.0, -4.0,
+        ]
+        for value in bad_values:
+            with self.subTest(value=value):
+                with self.assertRaises(TypeError):
+                    isqrt(value)
+
+    @support.bigmemtest(2**32, memuse=0.85)
+    def test_isqrt_huge(self, size):
+        isqrt = self.module.isqrt
+        if size & 1:
+            size += 1
+        v = 1 << size
+        w = isqrt(v)
+        self.assertEqual(w.bit_length(), size // 2 + 1)
+        self.assertEqual(w.bit_count(), 1)
+
+    def test_perm(self):
+        perm = self.module.perm
+        factorial = self.module.factorial
+        # Test if factorial definition is satisfied
+        for n in range(500):
+            for k in (range(n + 1) if n < 100 else range(30) if n < 200 else range(10)):
+                self.assertEqual(perm(n, k),
+                                 factorial(n) // factorial(n - k))
+
+        # Test for Pascal's identity
+        for n in range(1, 100):
+            for k in range(1, n):
+                self.assertEqual(perm(n, k), perm(n - 1, k - 1) * k + perm(n - 1, k))
+
+        # Test corner cases
+        for n in range(1, 100):
+            self.assertEqual(perm(n, 0), 1)
+            self.assertEqual(perm(n, 1), n)
+            self.assertEqual(perm(n, n), factorial(n))
+
+        # Test one argument form
+        for n in range(20):
+            self.assertEqual(perm(n), factorial(n))
+            self.assertEqual(perm(n, None), factorial(n))
+
+        # Raises TypeError if any argument is non-integer or argument count is
+        # not 1 or 2
+        self.assertRaises(TypeError, perm, 10, 1.0)
+        self.assertRaises(TypeError, perm, 10, Decimal(1.0))
+        self.assertRaises(TypeError, perm, 10, Fraction(1, 1))
+        self.assertRaises(TypeError, perm, 10, "1")
+        self.assertRaises(TypeError, perm, 10.0, 1)
+        self.assertRaises(TypeError, perm, Decimal(10.0), 1)
+        self.assertRaises(TypeError, perm, Fraction(10, 1), 1)
+        self.assertRaises(TypeError, perm, "10", 1)
+
+        self.assertRaises(TypeError, perm)
+        self.assertRaises(TypeError, perm, 10, 1, 3)
+        self.assertRaises(TypeError, perm)
+
+        # Raises Value error if not k or n are negative numbers
+        self.assertRaises(ValueError, perm, -1, 1)
+        self.assertRaises(ValueError, perm, -2**1000, 1)
+        self.assertRaises(ValueError, perm, 1, -1)
+        self.assertRaises(ValueError, perm, 1, -2**1000)
+
+        # Returns zero if k is greater than n
+        self.assertEqual(perm(1, 2), 0)
+        self.assertEqual(perm(1, 2**1000), 0)
+
+        n = 2**1000
+        self.assertEqual(perm(n, 0), 1)
+        self.assertEqual(perm(n, 1), n)
+        self.assertEqual(perm(n, 2), n * (n-1))
+        if support.check_impl_detail(cpython=True):
+            self.assertRaises(OverflowError, perm, n, n)
+
+        for n, k in (True, True), (True, False), (False, False):
+            self.assertIntEqual(perm(n, k), 1)
+        self.assertEqual(perm(IntSubclass(5), IntSubclass(2)), 20)
+        self.assertEqual(perm(MyIndexable(5), MyIndexable(2)), 20)
+        for k in range(3):
+            self.assertIs(type(perm(IntSubclass(5), IntSubclass(k))), int)
+            self.assertIs(type(perm(MyIndexable(5), MyIndexable(k))), int)
+
+    def test_comb(self):
+        comb = self.module.comb
+        factorial = self.module.factorial
+        # Test if factorial definition is satisfied
+        for n in range(500):
+            for k in (range(n + 1) if n < 100 else range(30) if n < 200 else range(10)):
+                self.assertEqual(comb(n, k), factorial(n)
+                    // (factorial(k) * factorial(n - k)))
+
+        # Test for Pascal's identity
+        for n in range(1, 100):
+            for k in range(1, n):
+                self.assertEqual(comb(n, k), comb(n - 1, k - 1) + comb(n - 1, k))
+
+        # Test corner cases
+        for n in range(100):
+            self.assertEqual(comb(n, 0), 1)
+            self.assertEqual(comb(n, n), 1)
+
+        for n in range(1, 100):
+            self.assertEqual(comb(n, 1), n)
+            self.assertEqual(comb(n, n - 1), n)
+
+        # Test Symmetry
+        for n in range(100):
+            for k in range(n // 2):
+                self.assertEqual(comb(n, k), comb(n, n - k))
+
+        # Raises TypeError if any argument is non-integer or argument count is
+        # not 2
+        self.assertRaises(TypeError, comb, 10, 1.0)
+        self.assertRaises(TypeError, comb, 10, Decimal(1.0))
+        self.assertRaises(TypeError, comb, 10, "1")
+        self.assertRaises(TypeError, comb, 10.0, 1)
+        self.assertRaises(TypeError, comb, Decimal(10.0), 1)
+        self.assertRaises(TypeError, comb, "10", 1)
+
+        self.assertRaises(TypeError, comb, 10)
+        self.assertRaises(TypeError, comb, 10, 1, 3)
+        self.assertRaises(TypeError, comb)
+
+        # Raises Value error if not k or n are negative numbers
+        self.assertRaises(ValueError, comb, -1, 1)
+        self.assertRaises(ValueError, comb, -2**1000, 1)
+        self.assertRaises(ValueError, comb, 1, -1)
+        self.assertRaises(ValueError, comb, 1, -2**1000)
+
+        # Returns zero if k is greater than n
+        self.assertEqual(comb(1, 2), 0)
+        self.assertEqual(comb(1, 2**1000), 0)
+
+        n = 2**1000
+        self.assertEqual(comb(n, 0), 1)
+        self.assertEqual(comb(n, 1), n)
+        self.assertEqual(comb(n, 2), n * (n-1) // 2)
+        self.assertEqual(comb(n, n), 1)
+        self.assertEqual(comb(n, n-1), n)
+        self.assertEqual(comb(n, n-2), n * (n-1) // 2)
+        if support.check_impl_detail(cpython=True):
+            self.assertRaises(OverflowError, comb, n, n//2)
+
+        for n, k in (True, True), (True, False), (False, False):
+            self.assertIntEqual(comb(n, k), 1)
+        self.assertEqual(comb(IntSubclass(5), IntSubclass(2)), 10)
+        self.assertEqual(comb(MyIndexable(5), MyIndexable(2)), 10)
+        for k in range(3):
+            self.assertIs(type(comb(IntSubclass(5), IntSubclass(k))), int)
+            self.assertIs(type(comb(MyIndexable(5), MyIndexable(k))), int)
+
+
+class MathTests(IntMathTests):
+    import math as module
+
+
+class MiscTests(unittest.TestCase):
+
+    def test_module_name(self):
+        import math.integer
+        self.assertEqual(math.integer.__name__, 'math.integer')
+        for name in dir(math.integer):
+            if not name.startswith('_'):
+                obj = getattr(math.integer, name)
+                self.assertEqual(obj.__module__, 'math.integer')
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/crates/stdlib/src/lib.rs b/crates/stdlib/src/lib.rs
index bc3350231..5eb8062ce 100644
--- a/crates/stdlib/src/lib.rs
+++ b/crates/stdlib/src/lib.rs
@@ -232,6 +232,7 @@ pub fn stdlib_module_defs(ctx: &Context) -> Vec<&'static builtins::PyModuleDef>
         locale::module_def(ctx),
         lzma::module_def(ctx),
         math::module_def(ctx),
+        math::integer_module_def(ctx),
         md5::module_def(ctx),
         #[cfg(all(feature = "host_env", any(unix, windows)))]
         mmap::module_def(ctx),
diff --git a/crates/stdlib/src/math.rs b/crates/stdlib/src/math.rs
index 97e34efcc..9ea0be522 100644
--- a/crates/stdlib/src/math.rs
+++ b/crates/stdlib/src/math.rs
@@ -1,3 +1,4 @@
+pub(crate) use _math_integer::module_def as integer_module_def;
 pub(crate) use math::module_def;
 
 use crate::vm::{VirtualMachine, builtins::PyBaseExceptionRef};
@@ -5,19 +6,31 @@ use crate::vm::{VirtualMachine, builtins::PyBaseExceptionRef};
 #[pymodule]
 mod math {
     use crate::vm::{
-        AsObject, PyObject, PyObjectRef, PyRef, PyResult, VirtualMachine,
-        builtins::{PyFloat, PyInt, PyIntRef, PyStrInterned, try_bigint_to_f64, try_f64_to_bigint},
+        AsObject, Py, PyObject, PyObjectRef, PyRef, PyResult, VirtualMachine,
+        builtins::{
+            PyFloat, PyInt, PyIntRef, PyModule, PyStrInterned, try_bigint_to_f64, try_f64_to_bigint,
+        },
         function::{
-            ArgIndex, ArgIntoFloat, ArgIterable, Either, NameCoordinates, NameIntegers,
-            OptionalArg, PosArgs,
+            ArgIndex, ArgIntoFloat, ArgIterable, Either, NameCoordinates, OptionalArg, PosArgs,
         },
         identifier,
     };
     use malachite_bigint::BigInt;
-    use num_traits::{Signed, ToPrimitive};
 
     use super::{float_repr, pymath_exception};
 
+    pub(crate) fn module_exec(vm: &VirtualMachine, module: &Py<PyModule>) -> PyResult<()> {
+        __module_exec(vm, module);
+        let integer = vm.import("_math_integer", 0)?;
+        for name in super::INTEGER_FUNCTIONS {
+            module.set_attr(name, integer.get_attr(name, vm)?, vm)?;
+        }
+        vm.sys_module
+            .get_attr("modules", vm)?
+            .set_item("math.integer", integer.clone(), vm)?;
+        module.set_attr("integer", integer, vm)
+    }
+
     // Constants
     #[pyattr]
     use core::f64::consts::{E as e, PI as pi, TAU as tau};
@@ -48,6 +61,36 @@ mod math {
         pymath::math::isnan(x.into_float())
     }
 
+    /// Return True if x is normal, and False otherwise.
+    #[pyfunction]
+    fn isnormal(x: ArgIntoFloat) -> bool {
+        pymath::math::isnormal(x.into_float())
+    }
+
+    /// Return True if x is subnormal, and False otherwise.
+    #[pyfunction]
+    fn issubnormal(x: ArgIntoFloat) -> bool {
+        pymath::math::issubnormal(x.into_float())
+    }
+
+    /// Return True if the sign of x is negative and False otherwise.
+    #[pyfunction]
+    fn signbit(x: ArgIntoFloat) -> bool {
+        pymath::math::signbit(x.into_float())
+    }
+
+    /// Return the larger of two floating-point arguments.
+    #[pyfunction]
+    fn fmax(x: ArgIntoFloat, y: ArgIntoFloat) -> f64 {
+        pymath::math::fmax(x.into_float(), y.into_float())
+    }
+
+    /// Return the smaller of two floating-point arguments.
+    #[pyfunction]
+    fn fmin(x: ArgIntoFloat, y: ArgIntoFloat) -> f64 {
+        pymath::math::fmin(x.into_float(), y.into_float())
+    }
+
     #[derive(FromArgs)]
     struct IsCloseArgs {
         #[pyarg(any)]
@@ -773,7 +816,18 @@ mod math {
             pymath::Error::ERANGE => vm.new_overflow_error("overflow in fma"),
         })
     }
+}
 
+const INTEGER_FUNCTIONS: [&str; 6] = ["comb", "factorial", "gcd", "isqrt", "lcm", "perm"];
+
+#[pymodule(sub, name = "math")]
+mod integer_functions {
+    use crate::vm::{
+        PyResult, VirtualMachine,
+        function::{ArgIndex, NameIntegers, PosArgs},
+    };
+    use malachite_bigint::BigInt;
+    use num_traits::{Signed, ToPrimitive};
     // Integer functions:
 
     #[pyfunction]
@@ -960,3 +1014,21 @@ fn float_repr(value: f64) -> String {
         }
     }
 }
+
+/// This module provides access to integer related mathematical functions.
+#[pymodule(with(integer_functions))]
+mod _math_integer {
+    use crate::vm::{Py, PyResult, VirtualMachine, builtins::PyModule};
+
+    pub(crate) fn module_exec(vm: &VirtualMachine, module: &Py<PyModule>) -> PyResult<()> {
+        __module_exec(vm, module);
+        let name = vm.ctx.new_str("math.integer");
+        module.set_attr("__name__", name.clone(), vm)?;
+        for function_name in super::INTEGER_FUNCTIONS {
+            module
+                .get_attr(function_name, vm)?
+                .set_attr("__module__", name.clone(), vm)?;
+        }
+        Ok(())
+    }
+}
diff --git a/extra_tests/snippets/stdlib_math.py b/extra_tests/snippets/stdlib_math.py
index 6e4e74b04..b06bf8446 100644
--- a/extra_tests/snippets/stdlib_math.py
+++ b/extra_tests/snippets/stdlib_math.py
@@ -7,6 +7,81 @@ NAN = float("nan")
 INF = float("inf")
 NINF = float("-inf")
 
+
+def test_math_315_apis():
+    import math.integer
+    import sys
+
+    import _math_integer
+
+    assert math.integer is _math_integer
+    assert sys.modules["math.integer"] is math.integer
+    assert math.integer.__name__ == "math.integer"
+    for name, args, expected in (
+        ("comb", (5, 2), 10),
+        ("factorial", (5,), 120),
+        ("gcd", (12, 18), 6),
+        ("isqrt", (17,), 4),
+        ("lcm", (4, 6), 12),
+        ("perm", (5, 2), 20),
+    ):
+        function = getattr(math, name)
+        assert function is getattr(math.integer, name)
+        assert function.__module__ == "math.integer"
+        assert function(*args) == expected
+
+    for value in (1.0, -1.0, sys.float_info.min, sys.float_info.max):
+        assert math.isnormal(value) is True
+        assert math.issubnormal(value) is False
+    for value in (sys.float_info.min / 2, -sys.float_info.min / 2):
+        assert math.isnormal(value) is False
+        assert math.issubnormal(value) is True
+    for value in (0.0, -0.0, INF, NINF, NAN):
+        assert math.isnormal(value) is False
+        assert math.issubnormal(value) is False
+    for value in (0.0, 1.0, INF, NAN):
+        assert math.signbit(math.copysign(value, 1.0)) is False
+        assert math.signbit(math.copysign(value, -1.0)) is True
+
+    class FloatArgument:
+        def __float__(self):
+            return -2.0
+
+    class IndexArgument:
+        def __index__(self):
+            return 2
+
+    for function in (math.isnormal, math.issubnormal, math.signbit):
+        assert function(FloatArgument()) == function(-2.0)
+        assert function(IndexArgument()) == function(2.0)
+        assert_raises(TypeError, function, object())
+        assert_raises(TypeError, function, x=1.0)
+        assert_raises(TypeError, function)
+        assert_raises(TypeError, function, 1.0, 2.0)
+
+    for function, expected in ((math.fmax, 2.0), (math.fmin, -2.0)):
+        assert function(-2.0, 2.0) == expected
+        assert function(2.0, -2.0) == expected
+        assert type(function(1, 2)) is float
+        assert function(FloatArgument(), IndexArgument()) == expected
+        assert math.isnan(function(NAN, NAN))
+        for value in (0.0, -0.0, 1.0, -1.0, INF, NINF):
+            for result in (function(NAN, value), function(value, NAN)):
+                assert result == value
+                assert math.signbit(result) == math.signbit(value)
+        assert function(0.0, -0.0) == 0.0
+        assert function(-0.0, 0.0) == 0.0
+        assert_raises(TypeError, function, object(), 1.0)
+        assert_raises(TypeError, function, 1.0, object())
+        assert_raises(TypeError, function, x=1.0, y=2.0)
+        assert_raises(TypeError, function, 1.0)
+        assert_raises(TypeError, function, 1.0, 2.0, 3.0)
+    assert math.fmax(NINF, INF) == INF
+    assert math.fmin(NINF, INF) == NINF
+
+
+test_math_315_apis()
+
 # assert(math.exp(2) == math.exp(2.0))
 # assert(math.exp(True) == math.exp(1.0))
 #
diff --git a/scripts/update_lib/deps.py b/scripts/update_lib/deps.py
index da49490c7..8c1a28ef8 100644
--- a/scripts/update_lib/deps.py
+++ b/scripts/update_lib/deps.py
@@ -651,6 +651,7 @@ DEPENDENCIES = {
         "lib": [],
         "test": [
             "test_math.py",
+            "test_math_integer.py",
             "test_math_property.py",
         ],
     },