diff --git a/Doc/library/math.integer.rst b/Doc/library/math.integer.rst
index c3f34cdfd85410c..752790dde8e11f1 100644
--- a/Doc/library/math.integer.rst
+++ b/Doc/library/math.integer.rst
@@ -13,8 +13,9 @@ These functions accept integers and objects that implement the
 :meth:`~object.__index__` method which is used to convert the object to an integer
 number.
 
-The following functions are provided by this module.  All return values are
-computed exactly and are integers.
+The following functions are provided by this module.
+Unless stated otherwise below,
+all return values are computed exactly and are integers.
 
 
 .. function:: comb(n, k, /)
@@ -46,6 +47,23 @@ computed exactly and are integers.
    returns ``0``.
 
 
+.. function:: isprime(n, /)
+
+   Return ``True`` if *n* is prime or probably prime, ``False`` otherwise.
+
+   For ``n < 2**64``, the result is exact. For larger integers, ``True``
+   indicates that *n* is probably prime. However, a composite number may be
+   accepted. A ``False`` result always means that *n* is not prime.
+
+   For large inputs, the implementation uses 64 Miller--Rabin rounds with
+   independently and uniformly chosen random bases. For any fixed composite
+   input, the probability of returning ``True`` is at most ``2**-128``.
+   This bound is over the random choices made by the test, not the probability
+   that an input is composite given a ``True`` result.
+
+   .. versionadded:: next
+
+
 .. function:: isqrt(n, /)
 
    Return the integer square root of the nonnegative integer *n*. This is the
diff --git a/Doc/whatsnew/3.16.rst b/Doc/whatsnew/3.16.rst
index e342aa71c7924ea..c2bf230b2416e87 100644
--- a/Doc/whatsnew/3.16.rst
+++ b/Doc/whatsnew/3.16.rst
@@ -565,6 +565,14 @@ math
   (Contributed by Jeff Epler in :gh:`150534`.)
 
 
+math.integer
+------------
+
+* Added :func:`math.integer.isprime` for primality testing of integers.
+  Results are exact below ``2**64``, and larger inputs use a probable-prime test.
+  (Contributed by Serhiy Storchaka and Stan Ulbrych in :gh:`153222`.)
+
+
 multiprocessing
 ---------------
 
diff --git a/Lib/test/test_math_integer.py b/Lib/test/test_math_integer.py
index d1bc776dc42391a..7cf53bc185a0a31 100644
--- a/Lib/test/test_math_integer.py
+++ b/Lib/test/test_math_integer.py
@@ -1,5 +1,6 @@
 from decimal import Decimal
 from fractions import Fraction
+import random
 import unittest
 from test import support
 
@@ -78,6 +79,46 @@ def py_factorial(n):
         outer *= inner
     return outer << (n - count_set_bits(n))
 
+# Reference implementations for primality testing.
+
+def primes_below(n):
+    """List of the primes below n, by the sieve of Eratosthenes."""
+    sieve = bytearray([1]) * n
+    sieve[:2] = bytes(2)
+    i = 2
+    while i * i < n:
+        if sieve[i]:
+            sieve[i*i::i] = bytes(len(range(i*i, n, i)))
+        i += 1
+    return [i for i in range(n) if sieve[i]]
+
+# The Miller-Rabin test with the first 13 primes as bases is known to be
+# exact for n < 3.3 * 10**24. These fixed bases provide an independent
+# cross-check within that range; results beyond it are not guaranteed.
+MILLER_RABIN_BASES = primes_below(100)
+
+def py_isprime(n):
+    """Miller-Rabin primality test, for cross-checking."""
+    if n < 2:
+        return False
+    for p in MILLER_RABIN_BASES:
+        if n % p == 0:
+            return n == p
+    d = n - 1
+    s = (d & -d).bit_length() - 1
+    d >>= s
+    for a in MILLER_RABIN_BASES:
+        x = pow(a, d, n)
+        if x == 1 or x == n - 1:
+            continue
+        for _ in range(s - 1):
+            x = x * x % n
+            if x == n - 1:
+                break
+        else:
+            return False
+    return True
+
 
 class IntMathTests(unittest.TestCase):
     import math.integer as module
@@ -406,6 +447,174 @@ class MathTests(IntMathTests):
     import math as module
 
 
+# isprime() exists only in math.integer, not in math, so its tests
+# are not in IntMathTests (which is re-run against math above).
+
+class IsPrimeTests(unittest.TestCase):
+    import math.integer as module
+
+    def test_isprime_small(self):
+        isprime = self.module.isprime
+        sieve = set(primes_below(10**4))
+        for n in range(-10, 10**4):
+            with self.subTest(n=n):
+                self.assertIs(isprime(n), n in sieve)
+
+    def test_isprime_negative(self):
+        isprime = self.module.isprime
+        self.assertIs(isprime(-1), False)
+        self.assertIs(isprime(-2), False)
+        self.assertIs(isprime(-3), False)
+        self.assertIs(isprime(-10**100), False)
+
+    def test_isprime_carmichael(self):
+        # Carmichael numbers are composite (A002997).
+        isprime = self.module.isprime
+        for n in [561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841,
+                  29341, 41041, 46657, 52633, 62745, 63973, 75361]:
+            with self.subTest(n=n):
+                self.assertIs(isprime(n), False)
+
+    def test_isprime_strong_pseudoprimes_base_2(self):
+        # Composites that pass the strong probable prime test to base 2
+        # (A001262); they must be caught by the other bases.
+        isprime = self.module.isprime
+        for n in [2047, 3277, 4033, 4681, 8321, 15841, 29341, 42799, 49141,
+                  52633, 65281, 74665, 80581, 85489, 88357, 90751,
+                  3825123056546413051,
+                  318665857834031151167461,
+                  3317044064679887385961981]:
+            with self.subTest(n=n):
+                self.assertIs(isprime(n), False)
+
+    def test_isprime_strong_lucas_pseudoprimes(self):
+        # Composites that are strong Lucas pseudoprimes (A217255); they
+        # are unrelated to the Miller-Rabin bases and must still fail.
+        isprime = self.module.isprime
+        for n in [5459, 5777, 10877, 16109, 18971, 22499, 24569, 25199,
+                  40309, 58519, 75077, 97439]:
+            with self.subTest(n=n):
+                self.assertIs(isprime(n), False)
+
+    def test_isprime_base_divisors(self):
+        # Divisors of the Miller-Rabin bases exercise the case where a
+        # base is divisible by the tested number.
+        isprime = self.module.isprime
+        for base in [2, 7, 61, 325, 9375, 28178, 450775, 9780504,
+                     1795265022]:
+            d = 1
+            while d * d <= base:
+                if base % d == 0:
+                    for n in [d, base // d]:
+                        with self.subTest(base=base, n=n):
+                            self.assertIs(isprime(n), py_isprime(n))
+                d += 1
+
+    def test_isprime_perfect_squares(self):
+        isprime = self.module.isprime
+        for k in range(1000):
+            with self.subTest(k=k):
+                self.assertIs(isprime(k*k), False)
+        for k in [2**31 - 1, 2**32 - 5, 2**61 - 1, 2**89 - 1]:
+            with self.subTest(k=k):
+                self.assertIs(isprime(k*k), False)
+
+    def test_isprime_word_boundaries(self):
+        # Exercise word boundaries and the boundary between the base sets.
+        isprime = self.module.isprime
+        for boundary in [2**32, 4759123141, 2**63]:
+            for n in range(boundary - 200, boundary + 200):
+                with self.subTest(n=n):
+                    self.assertIs(isprime(n), py_isprime(n))
+        for n in range(2**64 - 200, 2**64 + 200):
+            with self.subTest(n=n):
+                self.assertIs(isprime(n), py_isprime(n))
+        self.assertIs(isprime(2**64 - 59), True)   # largest prime < 2**64
+
+    def test_isprime_many_squarings(self):
+        # Exercise repeated squaring when n - 1 has many factors of two,
+        # including large moduli such as the prime 2**64 - 2**32 + 1.
+        isprime = self.module.isprime
+        for s in range(1, 64):
+            for d in [1, 3, 5, 2**(64 - s) - 1]:
+                n = (d << s) + 1
+                if n < 2**64:
+                    with self.subTest(n=n):
+                        self.assertIs(isprime(n), py_isprime(n))
+
+    def test_isprime_large_values(self):
+        isprime = self.module.isprime
+        self.assertIs(isprime(2**61 - 1), True)    # Mersenne prime
+        self.assertIs(isprime(2**62 - 1), False)
+        self.assertIs(isprime(10**19), False)
+        # Known primes above the exact 64-bit range.
+        for n in [2**64 + 13, 2**89 - 1, 2**107 - 1,
+                  2**127 - 1, 2**521 - 1]:
+            with self.subTest(n=n):
+                self.assertIs(isprime(n), True)
+        for n in [2**64, 2**64 + 1, 2**128 + 1, 10**100,
+                  (2**61 - 1) * (2**89 - 1)]:
+            with self.subTest(n=n):
+                self.assertIs(isprime(n), False)
+
+    def test_isprime_large_small_factors(self):
+        for p in primes_below(67):
+            with self.subTest(p=p):
+                self.assertIs(self.module.isprime(p * (2**127 - 1)), False)
+
+    def test_isprime_large_many_squarings(self):
+        # Stay within the reference test's proven range.
+        for s in [64, 65, 70, 75]:
+            for d in [1, 3, 5, 7]:
+                n = (d << s) + 1
+                with self.subTest(n=n):
+                    self.assertIs(self.module.isprime(n), py_isprime(n))
+
+    def test_isprime_random(self):
+        isprime = self.module.isprime
+        rng = random.Random(1729)
+        for bits in [32, 34, 63, 64]:
+            for _ in range(300):
+                n = rng.getrandbits(bits)
+                with self.subTest(n=n):
+                    self.assertIs(isprime(n), py_isprime(n))
+        for bits in [65, 80]:
+            for _ in range(30):
+                n = rng.getrandbits(bits) | (1 << (bits - 1)) | 1
+                with self.subTest(n=n):
+                    self.assertIs(isprime(n), py_isprime(n))
+
+    def test_isprime_integer_like(self):
+        isprime = self.module.isprime
+        self.assertIs(isprime(False), False)
+        self.assertIs(isprime(True), False)
+        self.assertIs(isprime(IntSubclass(7)), True)
+        self.assertIs(isprime(IntSubclass(8)), False)
+        self.assertIs(isprime(MyIndexable(97)), True)
+        self.assertIs(isprime(MyIndexable(-97)), False)
+        self.assertIs(isprime(IntSubclass(2**89 - 1)), True)
+        self.assertIs(isprime(MyIndexable(2**89 - 1)), True)
+        self.assertIs(isprime(MyIndexable(2**89)), False)
+
+    def test_isprime_int_subclass_operators(self):
+        # Overridden operators of an int subclass must not affect the
+        # result.
+        isprime = self.module.isprime
+        self.assertIs(isprime(BadIntSubclass(97)), True)
+        self.assertIs(isprime(BadIntSubclass(2**61 - 1)), True)
+        self.assertIs(isprime(BadIntSubclass(2**62 - 1)), False)
+        self.assertIs(isprime(BadIntSubclass(2**89 - 1)), True)
+        self.assertIs(isprime(BadIntSubclass(2**89)), False)
+
+    def test_isprime_non_integers(self):
+        isprime = self.module.isprime
+        for value in [7.0, 7.5, Decimal('7'), Fraction(7, 1), '7', 7.5j]:
+            with self.subTest(value=value):
+                self.assertRaises(TypeError, isprime, value)
+        self.assertRaises(TypeError, isprime)
+        self.assertRaises(TypeError, isprime, 7, 11)
+
+
 class MiscTests(unittest.TestCase):
 
     def test_module_name(self):
@@ -416,6 +625,12 @@ def test_module_name(self):
                 obj = getattr(math.integer, name)
                 self.assertEqual(obj.__module__, 'math.integer')
 
+    def test_math_namespace(self):
+        # The new function is added only to math.integer, not to math
+        # (PEP 791).
+        import math
+        self.assertFalse(hasattr(math, 'isprime'))
+
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/Misc/NEWS.d/next/Library/2026-07-06-14-00-00.gh-issue-153222.pRiMe1.rst b/Misc/NEWS.d/next/Library/2026-07-06-14-00-00.gh-issue-153222.pRiMe1.rst
new file mode 100644
index 000000000000000..27016e8dd7ccfac
--- /dev/null
+++ b/Misc/NEWS.d/next/Library/2026-07-06-14-00-00.gh-issue-153222.pRiMe1.rst
@@ -0,0 +1,2 @@
+Add :func:`math.integer.isprime` for primality testing of integers.
+Results are exact below ``2**64``; larger inputs use a probable-prime test.
diff --git a/Modules/clinic/mathintegermodule.c.h b/Modules/clinic/mathintegermodule.c.h
index 29c2a0ac902ea39..9230faaaf3c3cee 100644
--- a/Modules/clinic/mathintegermodule.c.h
+++ b/Modules/clinic/mathintegermodule.c.h
@@ -67,6 +67,37 @@ PyDoc_STRVAR(math_integer_isqrt__doc__,
 #define MATH_INTEGER_ISQRT_METHODDEF    \
     {"isqrt", (PyCFunction)math_integer_isqrt, METH_O, math_integer_isqrt__doc__},
 
+PyDoc_STRVAR(math_integer_isprime__doc__,
+"isprime($module, n, /)\n"
+"--\n"
+"\n"
+"Return True if n is prime or probably prime, False otherwise.\n"
+"\n"
+"The result is exact for n < 2**64. For larger integers, True means\n"
+"probably prime. False always means that n is not prime.");
+
+#define MATH_INTEGER_ISPRIME_METHODDEF    \
+    {"isprime", (PyCFunction)math_integer_isprime, METH_O, math_integer_isprime__doc__},
+
+static int
+math_integer_isprime_impl(PyObject *module, PyObject *n);
+
+static PyObject *
+math_integer_isprime(PyObject *module, PyObject *n)
+{
+    PyObject *return_value = NULL;
+    int _return_value;
+
+    _return_value = math_integer_isprime_impl(module, n);
+    if ((_return_value == -1) && PyErr_Occurred()) {
+        goto exit;
+    }
+    return_value = PyBool_FromLong((long)_return_value);
+
+exit:
+    return return_value;
+}
+
 PyDoc_STRVAR(math_integer_factorial__doc__,
 "factorial($module, n, /)\n"
 "--\n"
@@ -156,4 +187,4 @@ math_integer_comb(PyObject *module, PyObject *const *args, Py_ssize_t nargs)
 exit:
     return return_value;
 }
-/*[clinic end generated code: output=34697570c923a3af input=a9049054013a1b77]*/
+/*[clinic end generated code: output=4e47e334e7fe1156 input=a9049054013a1b77]*/
diff --git a/Modules/mathintegermodule.c b/Modules/mathintegermodule.c
index 6d35c825349e451..ef2469250dd0003 100644
--- a/Modules/mathintegermodule.c
+++ b/Modules/mathintegermodule.c
@@ -8,6 +8,7 @@
 #include "pycore_abstract.h"      // _PyNumber_Index()
 #include "pycore_bitutils.h"      // _Py_bit_length()
 #include "pycore_long.h"          // _PyLong_GetZero()
+#include "pycore_pylifecycle.h"   // _PyOS_URandom()
 
 #include "clinic/mathintegermodule.c.h"
 
@@ -481,6 +482,395 @@ math_integer_isqrt(PyObject *module, PyObject *n)
 }
 
 
+/* Primality testing below 2**64 with the deterministic Miller-Rabin test.
+
+   A strong probable prime test is exact for n < 4759123141 with the
+   bases {2, 7, 61} (G. Jaeschke, "On strong pseudoprimes to several
+   bases", Mathematics of Computation 61 (1993), 915-926), and for all
+   n < 2**64 with the seven bases below, found by Jim Sinclair and
+   verified against the Feitsma-Galway exhaustive list of base-2 strong
+   pseudoprimes below 2**64 (see http://miller-rabin.appspot.com/ and
+   http://ntheory.org/pseudoprimes.html).
+
+   Modular exponentiation and squaring use Montgomery arithmetic with
+   R = 2**64 to avoid division in their inner loops (P. L. Montgomery,
+   "Modular multiplication without trial division", Mathematics of
+   Computation 44 (1985), 519-521).
+*/
+
+static const uint8_t small_primes[] = {
+    2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
+};
+
+static const uint32_t mr_bases_small[] = {2, 7, 61};
+static const uint32_t mr_bases[] = {
+    2, 325, 9375, 28178, 450775, 9780504, 1795265022,
+};
+
+#if defined(HAVE_GCC_UINT128_T)
+static inline uint64_t
+mulmod_u64(uint64_t a, uint64_t b, uint64_t n)
+{
+    return (uint64_t)((__uint128_t)a * b % n);
+}
+#elif defined(_MSC_VER) && defined(_M_X64)
+#  include <intrin.h>
+#  pragma intrinsic(_umul128, _udiv128)
+static inline uint64_t
+mulmod_u64(uint64_t a, uint64_t b, uint64_t n)
+{
+    /* Since a, b < n, the high 64 bits of a*b are less than n and
+       _udiv128() cannot fault. */
+    uint64_t hi, rem;
+    uint64_t lo = _umul128(a, b, &hi);
+    (void)_udiv128(hi, lo, n, &rem);
+    return rem;
+}
+#else
+/* No fast 128-bit multiplication: fall back to shift-and-add, with the
+   modular additions written to avoid overflow for n close to 2**64.
+   Assumes a < n. */
+static inline uint64_t
+mulmod_u64(uint64_t a, uint64_t b, uint64_t n)
+{
+    uint64_t result = 0;
+    while (b) {
+        if (b & 1) {
+            result = (result >= n - a) ? result - (n - a) : result + a;
+        }
+        a = (a >= n - a) ? a - (n - a) : a + a;
+        b >>= 1;
+    }
+    return result;
+}
+#endif
+
+/* High 64 bits of the full product a * b. */
+static inline uint64_t
+mulhi_u64(uint64_t a, uint64_t b)
+{
+#if defined(HAVE_GCC_UINT128_T)
+    return (uint64_t)(((__uint128_t)a * b) >> 64);
+#elif defined(_MSC_VER) && defined(_M_X64)
+    uint64_t hi;
+    (void)_umul128(a, b, &hi);
+    return hi;
+#else
+    /* Split the operands into 32-bit limbs. Each sum fits in uint64_t. */
+    uint64_t a0 = (uint32_t)a, a1 = a >> 32;
+    uint64_t b0 = (uint32_t)b, b1 = b >> 32;
+    uint64_t t = a1 * b0 + ((a0 * b0) >> 32);
+    uint64_t carry = t >> 32;
+    t = a0 * b1 + (uint32_t)t;
+    return a1 * b1 + carry + (t >> 32);
+#endif
+}
+
+typedef struct {
+    uint64_t n;
+    uint64_t inverse;  // n**-1 modulo 2**64
+    uint64_t one;      // 2**64 modulo n
+    uint64_t d;        // Odd part of n - 1
+    int s;             // n - 1 = d * 2**s
+} montgomery_context;
+
+/* Precompute once for all bases. n must be odd and greater than 1. */
+static montgomery_context
+montgomery_init(uint64_t n)
+{
+    montgomery_context ctx = {
+        .n = n,
+        .inverse = 1,
+        .one = (0 - n) % n,
+        .d = n - 1,
+        .s = 0,
+    };
+    /* Newton iteration doubles the number of correct inverse bits,
+       starting with one correct bit because n is odd. */
+    for (int i = 0; i < 6; i++) {
+        ctx.inverse *= 2 - n * ctx.inverse;
+    }
+    while ((ctx.d & 1) == 0) {
+        ctx.d >>= 1;
+        ctx.s++;
+    }
+    return ctx;
+}
+
+/* Return a * b * R**-1 modulo n, where R = 2**64 and a, b < n.
+   With t = a*b and m = (t mod R) * n**-1 mod R, t - m*n is divisible
+   by R, and (t - m*n) / R lies strictly between -n and n. Subtracting
+   avoids the possible 129-bit carry of the additive Montgomery form.
+   All uint64_t arithmetic intentionally wraps modulo R. */
+static inline uint64_t
+montgomery_mul(uint64_t a, uint64_t b, const montgomery_context *ctx)
+{
+    uint64_t lo = a * b, hi = mulhi_u64(a, b);
+    uint64_t mhi = mulhi_u64(lo * ctx->inverse, ctx->n);
+    uint64_t result = hi - mhi;
+    return hi < mhi ? result + ctx->n : result;
+}
+
+/* Strong probable prime test to the given base, in Montgomery form.
+   A base divisible by n gives no information and is treated as passing. */
+static int
+sprp_u64(const montgomery_context *ctx, uint64_t base)
+{
+    base %= ctx->n;
+    if (base == 0) {
+        return 1;
+    }
+    // Represent x by x*R modulo n throughout exponentiation.
+    uint64_t b = mulmod_u64(base, ctx->one, ctx->n);
+    uint64_t a = ctx->one;
+    for (uint64_t e = ctx->d; e; e >>= 1) {
+        if (e & 1) {
+            a = montgomery_mul(a, b, ctx);
+        }
+        b = montgomery_mul(b, b, ctx);
+    }
+    if (a == ctx->one || a == ctx->n - ctx->one) {
+        return 1;
+    }
+    for (int r = 1; r < ctx->s; r++) {
+        a = montgomery_mul(a, a, ctx);
+        if (a == ctx->n - ctx->one) {
+            return 1;
+        }
+    }
+    return 0;
+}
+
+static int
+isprime_u64(uint64_t n)
+{
+    for (size_t i = 0; i < Py_ARRAY_LENGTH(small_primes); i++) {
+        if (n % small_primes[i] == 0) {
+            return n == small_primes[i];
+        }
+    }
+    if (n < 67 * 67) {
+        return n > 1;
+    }
+    montgomery_context ctx = montgomery_init(n);
+    const uint32_t *bases = mr_bases;
+    size_t num_bases = Py_ARRAY_LENGTH(mr_bases);
+    if (n < 4759123141) {
+        bases = mr_bases_small;
+        num_bases = Py_ARRAY_LENGTH(mr_bases_small);
+    }
+    for (size_t i = 0; i < num_bases; i++) {
+        if (!sprp_u64(&ctx, bases[i])) {
+            return 0;
+        }
+    }
+    return 1;
+}
+
+/* Strong probable prime test using arbitrary-precision arithmetic.
+   n - 1 = d * 2**s, where d is odd. Return -1 on error. */
+static int
+sprp_long(PyObject *n, PyObject *base, PyObject *d, int64_t s,
+          PyObject *n_minus_one)
+{
+    PyObject *x = PyNumber_Power(base, d, n);
+    if (x == NULL) {
+        return -1;
+    }
+    int result = PyObject_RichCompareBool(x, _PyLong_GetOne(), Py_EQ);
+    if (result != 0) {
+        goto done;
+    }
+    for (int64_t r = 0; r < s; r++) {
+        result = PyObject_RichCompareBool(x, n_minus_one, Py_EQ);
+        if (result != 0 || r == s - 1) {
+            goto done;
+        }
+        if (PyErr_CheckSignals() < 0) {
+            result = -1;
+            goto done;
+        }
+        Py_SETREF(x, PyNumber_Multiply(x, x));
+        if (x == NULL) {
+            return -1;
+        }
+        Py_SETREF(x, PyNumber_Remainder(x, n));
+        if (x == NULL) {
+            return -1;
+        }
+        result = PyObject_RichCompareBool(x, _PyLong_GetOne(), Py_EQ);
+        if (result != 0) {
+            /* Reaching 1 before -1 proves that n is composite. */
+            if (result > 0) {
+                result = 0;
+            }
+            goto done;
+        }
+    }
+
+  done:
+    Py_DECREF(x);
+    return result;
+}
+
+/* Probable primality for n >= 2**64. A base-2 test rejects most composites
+   before requesting randomness. Then 64 independent, uniformly sampled
+   bases give a false-positive probability at most 4**-64 = 2**-128 for
+   any fixed composite n. */
+static int
+isprime_long(PyObject *n)
+{
+    assert(PyLong_CheckExact(n));
+    assert(_PyLong_NumBits(n) > 64);
+
+    for (size_t i = 0; i < Py_ARRAY_LENGTH(small_primes); i++) {
+        PyObject *p = PyLong_FromLong(small_primes[i]);
+        if (p == NULL) {
+            return -1;
+        }
+        PyObject *rem = PyNumber_Remainder(n, p);
+        Py_DECREF(p);
+        if (rem == NULL) {
+            return -1;
+        }
+        int divisible = _PyLong_IsZero((PyLongObject *)rem);
+        Py_DECREF(rem);
+        if (divisible) {
+            return 0;
+        }
+    }
+
+    int result = -1;
+    PyObject *n_minus_one = NULL, *d = NULL, *base = NULL, *lowbit = NULL;
+    unsigned char *buffer = NULL;
+    n_minus_one = PyNumber_Subtract(n, _PyLong_GetOne());
+    if (n_minus_one == NULL) {
+        goto done;
+    }
+    /* Extract the largest power of two dividing n - 1. */
+    lowbit = PyNumber_Negative(n_minus_one);
+    if (lowbit == NULL) {
+        goto done;
+    }
+    Py_SETREF(lowbit, PyNumber_And(n_minus_one, lowbit));
+    if (lowbit == NULL) {
+        goto done;
+    }
+    int64_t s = _PyLong_NumBits(lowbit) - 1;
+    Py_CLEAR(lowbit);
+    d = _PyLong_Rshift(n_minus_one, s);
+    if (d == NULL) {
+        goto done;
+    }
+    base = PyLong_FromLong(2);
+    if (base == NULL) {
+        goto done;
+    }
+    int passed = sprp_long(n, base, d, s, n_minus_one);
+    if (passed <= 0) {
+        result = passed;
+        goto done;
+    }
+
+    int64_t bits = _PyLong_NumBits(n);
+    int64_t bytes = (bits - 1) / 8 + 1;
+    if (bytes > PY_SSIZE_T_MAX) {
+        PyErr_NoMemory();
+        goto done;
+    }
+    Py_ssize_t size = (Py_ssize_t)bytes;
+    buffer = PyMem_Malloc(size);
+    if (buffer == NULL) {
+        PyErr_NoMemory();
+        goto done;
+    }
+    unsigned char mask = (unsigned char)(0xff >> (7 - (bits - 1) % 8));
+    for (int i = 0; i < 64; i++) {
+        /* Rejection sampling avoids the bias of reduction modulo n. */
+        for (;;) {
+            if (PyErr_CheckSignals() < 0 || _PyOS_URandom(buffer, size) < 0) {
+                goto done;
+            }
+            buffer[size - 1] &= mask;
+            Py_SETREF(base, _PyLong_FromByteArray(buffer, size, 1, 0));
+            if (base == NULL) {
+                goto done;
+            }
+            int above_one = PyObject_RichCompareBool(
+                base, _PyLong_GetOne(), Py_GT);
+            if (above_one < 0) {
+                goto done;
+            }
+            if (!above_one) {
+                continue;
+            }
+            int below_limit = PyObject_RichCompareBool(
+                base, n_minus_one, Py_LT);
+            if (below_limit < 0) {
+                goto done;
+            }
+            if (below_limit) {
+                break;
+            }
+        }
+        passed = sprp_long(n, base, d, s, n_minus_one);
+        if (passed <= 0) {
+            result = passed;
+            goto done;
+        }
+    }
+    result = 1;
+
+  done:
+    Py_XDECREF(n_minus_one);
+    Py_XDECREF(d);
+    Py_XDECREF(base);
+    Py_XDECREF(lowbit);
+    PyMem_Free(buffer);
+    return result;
+}
+
+/*[clinic input]
+math.integer.isprime -> bool
+
+    n: object
+    /
+
+Return True if n is prime or probably prime, False otherwise.
+
+The result is exact for n < 2**64. For larger integers, True means
+probably prime. False always means that n is not prime.
+[clinic start generated code]*/
+
+static int
+math_integer_isprime_impl(PyObject *module, PyObject *n)
+/*[clinic end generated code: output=c808dc14d8d86875 input=e771fc2bdb5e29fb]*/
+{
+    /* The arbitrary-precision path must not call int subclass operators. */
+    n = PyNumber_Index(n);
+    if (n == NULL) {
+        return -1;
+    }
+    int result;
+    uint64_t m;
+    if (_PyLong_IsNegative((PyLongObject *)n)) {
+        result = 0;
+    }
+    else if (PyLong_AsUInt64(n, &m) < 0) {
+        result = -1;
+        if (PyErr_ExceptionMatches(PyExc_OverflowError)) {
+            PyErr_Clear();
+            result = isprime_long(n);
+        }
+    }
+    else {
+        result = isprime_u64(m);
+    }
+    Py_DECREF(n);
+    return result;
+}
+
+
 static unsigned long
 count_set_bits(unsigned long n)
 {
@@ -1233,6 +1623,7 @@ static PyMethodDef math_integer_methods[] = {
     MATH_INTEGER_COMB_METHODDEF
     MATH_INTEGER_FACTORIAL_METHODDEF
     MATH_INTEGER_GCD_METHODDEF
+    MATH_INTEGER_ISPRIME_METHODDEF
     MATH_INTEGER_ISQRT_METHODDEF
     MATH_INTEGER_LCM_METHODDEF
     MATH_INTEGER_PERM_METHODDEF