From 81e40f0446eaab13387ba90b747683f8819fd906 Mon Sep 17 00:00:00 2001 From: Serhiy Storchaka Date: Mon, 6 Jul 2026 20:29:09 +0300 Subject: [PATCH 1/4] gh-153222: Add math.integer.isprime() and math.integer.primes() Draft implementation of prime-number functionality for the math.integer module (deferred out of the initial PEP 791 scope), supporting integers less than 2**64. isprime() uses the deterministic Miller-Rabin test: bases {2, 7, 61} below 4759123141, Jim Sinclair's seven bases up to 2**64 (verified against the Feitsma-Galway list of base-2 strong pseudoprimes). The result is exact for the whole supported range; larger arguments raise OverflowError. primes(start=2, stop=None) returns an iterator of the primes in [start, stop), unbounded if stop is None. Neither function is mirrored into the math namespace (PEP 791). Co-Authored-By: Claude Fable 5 --- Doc/library/math.integer.rst | 37 +- Doc/whatsnew/3.16.rst | 9 + .../pycore_global_objects_fini_generated.h | 1 + Include/internal/pycore_global_strings.h | 1 + .../internal/pycore_runtime_init_generated.h | 1 + .../internal/pycore_unicodeobject_generated.h | 4 + Lib/test/test_math_integer.py | 290 ++++++++++++++ ...-07-06-14-00-00.gh-issue-153222.pRiMe1.rst | 3 + Modules/clinic/mathintegermodule.c.h | 111 +++++- Modules/mathintegermodule.c | 372 +++++++++++++++++- 10 files changed, 821 insertions(+), 8 deletions(-) create mode 100644 Misc/NEWS.d/next/Library/2026-07-06-14-00-00.gh-issue-153222.pRiMe1.rst diff --git a/Doc/library/math.integer.rst b/Doc/library/math.integer.rst index c3f34cdfd85410c..d34e7f821c4df9a 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,20 @@ computed exactly and are integers. returns ``0``. +.. function:: isprime(n, /) + + Return ``True`` if *n* is a prime number, ``False`` otherwise. + + A prime number is a natural number greater than 1 + that is not a product of two smaller natural numbers. + Negative numbers, ``0`` and ``1`` are not prime. + + The argument must be less than 2\ :sup:`64`; + raises :exc:`OverflowError` otherwise. + + .. versionadded:: next + + .. function:: isqrt(n, /) Return the integer square root of the nonnegative integer *n*. This is the @@ -83,3 +98,21 @@ computed exactly and are integers. and the function returns ``n!``. Raises :exc:`ValueError` if either of the arguments are negative. + + +.. function:: primes(start=2, stop=None) + + Return an iterator of the prime numbers *p* with ``start <= p < stop``, + in increasing order. + If *stop* is ``None`` (the default), the iteration does not stop. + + Roughly equivalent to + ``(p for p in itertools.count(start) if isprime(p))`` + for the unbounded form. + + The bounds must be less than 2\ :sup:`64`; + raises :exc:`OverflowError` otherwise. + Unbounded iteration raises :exc:`OverflowError` + if the candidates reach that limit. + + .. versionadded:: next diff --git a/Doc/whatsnew/3.16.rst b/Doc/whatsnew/3.16.rst index e8c75d2f571061a..eb2537d15ce6e58 100644 --- a/Doc/whatsnew/3.16.rst +++ b/Doc/whatsnew/3.16.rst @@ -279,6 +279,15 @@ math (Contributed by Jeff Epler in :gh:`150534`.) +math.integer +------------ + +* Added :func:`math.integer.isprime` for primality testing of integers + less than 2**64 and + :func:`math.integer.primes` for iterating over prime numbers. + (Contributed by Serhiy Storchaka in :gh:`153222`.) + + os -- diff --git a/Include/internal/pycore_global_objects_fini_generated.h b/Include/internal/pycore_global_objects_fini_generated.h index 7d952a1e52561a6..be16ea2d4aa522c 100644 --- a/Include/internal/pycore_global_objects_fini_generated.h +++ b/Include/internal/pycore_global_objects_fini_generated.h @@ -2097,6 +2097,7 @@ _PyStaticObjects_CheckRefcnt(PyInterpreterState *interp) { _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(stdout)); _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(step)); _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(steps)); + _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(stop)); _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(store_name)); _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(strategy)); _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(strftime)); diff --git a/Include/internal/pycore_global_strings.h b/Include/internal/pycore_global_strings.h index 8a8bbc3b6d05bf7..91f2b1c979ee3c8 100644 --- a/Include/internal/pycore_global_strings.h +++ b/Include/internal/pycore_global_strings.h @@ -820,6 +820,7 @@ struct _Py_global_strings { STRUCT_FOR_ID(stdout) STRUCT_FOR_ID(step) STRUCT_FOR_ID(steps) + STRUCT_FOR_ID(stop) STRUCT_FOR_ID(store_name) STRUCT_FOR_ID(strategy) STRUCT_FOR_ID(strftime) diff --git a/Include/internal/pycore_runtime_init_generated.h b/Include/internal/pycore_runtime_init_generated.h index 366d2d300fb4780..273164bc882688c 100644 --- a/Include/internal/pycore_runtime_init_generated.h +++ b/Include/internal/pycore_runtime_init_generated.h @@ -2095,6 +2095,7 @@ extern "C" { INIT_ID(stdout), \ INIT_ID(step), \ INIT_ID(steps), \ + INIT_ID(stop), \ INIT_ID(store_name), \ INIT_ID(strategy), \ INIT_ID(strftime), \ diff --git a/Include/internal/pycore_unicodeobject_generated.h b/Include/internal/pycore_unicodeobject_generated.h index 00d6297432b5fca..179a47deb2dc899 100644 --- a/Include/internal/pycore_unicodeobject_generated.h +++ b/Include/internal/pycore_unicodeobject_generated.h @@ -3060,6 +3060,10 @@ _PyUnicode_InitStaticStrings(PyInterpreterState *interp) { _PyUnicode_InternStatic(interp, &string); assert(_PyUnicode_CheckConsistency(string, 1)); assert(PyUnicode_GET_LENGTH(string) != 1); + string = &_Py_ID(stop); + _PyUnicode_InternStatic(interp, &string); + assert(_PyUnicode_CheckConsistency(string, 1)); + assert(PyUnicode_GET_LENGTH(string) != 1); string = &_Py_ID(store_name); _PyUnicode_InternStatic(interp, &string); assert(_PyUnicode_CheckConsistency(string, 1)); diff --git a/Lib/test/test_math_integer.py b/Lib/test/test_math_integer.py index d1bc776dc42391a..ef0ed0b7791a4f5 100644 --- a/Lib/test/test_math_integer.py +++ b/Lib/test/test_math_integer.py @@ -1,5 +1,7 @@ from decimal import Decimal from fractions import Fraction +import itertools +import random import unittest from test import support @@ -78,6 +80,47 @@ 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. More bases are used for larger inputs, for +# which the test is probabilistic (but independent of the Baillie-PSW test +# used in the implementation). +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 +449,246 @@ class MathTests(IntMathTests): import math as module +# isprime() and primes() exist only in math.integer, not in math, so their +# 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 strong Lucas test. + isprime = self.module.isprime + for n in [2047, 3277, 4033, 4681, 8321, 15841, 29341, 42799, 49141, + 52633, 65281, 74665, 80581, 85489, 88357, 90751, + 3825123056546413051]: + with self.subTest(n=n): + self.assertIs(isprime(n), False) + + def test_isprime_strong_lucas_pseudoprimes(self): + # Composites that pass the strong Lucas test with Selfridge's + # parameters (A217255). + 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]: + with self.subTest(k=k): + self.assertIs(isprime(k*k), False) + + def test_isprime_word_boundaries(self): + # Exercise the boundaries between the base sets. + isprime = self.module.isprime + for boundary in [2**32, 4759123141]: + 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): + 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_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) + # Arguments not less than 2**64 are not supported. + for n in [2**64, 2**64 + 13, 2**89 - 1, 10**100]: + with self.subTest(n=n): + self.assertRaises(OverflowError, 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, 128]: + n = rng.getrandbits(bits) | (1 << (bits - 1)) + with self.subTest(n=n): + self.assertRaises(OverflowError, 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) + + 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) + + 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 PrimesIterTests(unittest.TestCase): + import math.integer as module + + def test_primes(self): + primes = self.module.primes + expected = primes_below(10**4) + self.assertEqual(list(primes(stop=10**4)), expected) + self.assertEqual(len(expected), 1229) + self.assertEqual(list(itertools.islice(primes(), 25)), expected[:25]) + + def test_primes_start(self): + primes = self.module.primes + for start in [-10**100, -100, -1, 0, 1, 2]: + with self.subTest(start=start): + self.assertEqual(list(primes(start, 10)), [2, 3, 5, 7]) + self.assertEqual(list(primes(3, 10)), [3, 5, 7]) + self.assertEqual(list(primes(4, 10)), [5, 7]) + self.assertEqual(list(primes(8, 12)), [11]) + self.assertEqual(list(primes(9, 12)), [11]) + self.assertEqual(list(primes(7, 8)), [7]) + + def test_primes_stop(self): + primes = self.module.primes + # The range is half-open. + self.assertEqual(list(primes(2, 2)), []) + self.assertEqual(list(primes(2, 3)), [2]) + self.assertEqual(list(primes(3, 3)), []) + self.assertEqual(list(primes(7, 7)), []) + self.assertEqual(list(primes(2, -10)), []) + self.assertEqual(list(primes(10, 5)), []) + self.assertEqual(list(primes(stop=0)), []) + + def test_primes_unbounded(self): + primes = self.module.primes + it = primes() + self.assertEqual([next(it) for _ in range(5)], [2, 3, 5, 7, 11]) + it = primes(10**6) + self.assertEqual(next(it), 1000003) + + def test_primes_huge(self): + primes = self.module.primes + boundary = 10**18 + expected = [n for n in range(boundary - 200, boundary + 200) + if py_isprime(n)] + self.assertEqual(list(primes(boundary - 200, boundary + 200)), + expected) + # The top of the supported range. + expected = [n for n in range(2**64 - 200, 2**64) if py_isprime(n)] + self.assertEqual(list(primes(2**64 - 200, 2**64 - 1)), expected) + + def test_primes_overflow(self): + primes = self.module.primes + # The bounds must be less than 2**64. + self.assertRaises(OverflowError, primes, 2**64) + self.assertRaises(OverflowError, primes, 2**64 + 100, 2**64 + 200) + self.assertRaises(OverflowError, primes, 0, 2**64) + self.assertRaises(OverflowError, primes, 10**100) + # An unbounded iterator raises when it runs out of the + # supported range. + it = primes(2**64 - 60) + self.assertEqual(next(it), 2**64 - 59) + self.assertRaises(OverflowError, next, it) + self.assertRaises(StopIteration, next, it) + + def test_primes_iterator_protocol(self): + primes = self.module.primes + it = primes(2, 10) + self.assertIs(iter(it), it) + self.assertEqual(list(it), [2, 3, 5, 7]) + # An exhausted iterator stays exhausted. + self.assertEqual(list(it), []) + self.assertRaises(StopIteration, next, it) + # The iterator type cannot be instantiated directly. + self.assertRaises(TypeError, type(primes())) + + def test_primes_keywords(self): + primes = self.module.primes + self.assertEqual(list(primes(start=10, stop=30)), [11, 13, 17, 19, 23, 29]) + self.assertEqual(list(primes(10, stop=30)), [11, 13, 17, 19, 23, 29]) + + def test_primes_integer_like(self): + primes = self.module.primes + self.assertEqual(list(primes(True, 10)), [2, 3, 5, 7]) + self.assertEqual(list(primes(IntSubclass(3), IntSubclass(10))), [3, 5, 7]) + self.assertEqual(list(primes(MyIndexable(3), MyIndexable(10))), [3, 5, 7]) + # The yielded values are exact ints. + for p in primes(IntSubclass(3), IntSubclass(10)): + self.assertIs(type(p), int) + + def test_primes_int_subclass_operators(self): + # Overridden operators of an int subclass must not affect the + # iteration. + primes = self.module.primes + self.assertEqual(list(primes(BadIntSubclass(3), BadIntSubclass(10))), + [3, 5, 7]) + big = 10**18 + self.assertEqual(list(primes(BadIntSubclass(big), big + 100)), + [n for n in range(big, big + 100) if py_isprime(n)]) + + def test_primes_non_integers(self): + primes = self.module.primes + self.assertRaises(TypeError, primes, 2.5) + self.assertRaises(TypeError, primes, 2.5, 10) + self.assertRaises(TypeError, primes, 2, 10.5) + self.assertRaises(TypeError, primes, '2') + self.assertRaises(TypeError, primes, 2, '10') + self.assertRaises(TypeError, primes, 2, 10, 3) + + class MiscTests(unittest.TestCase): def test_module_name(self): @@ -416,6 +699,13 @@ def test_module_name(self): obj = getattr(math.integer, name) self.assertEqual(obj.__module__, 'math.integer') + def test_math_namespace(self): + # New functions are added only to math.integer, not to math + # (PEP 791). + import math + self.assertFalse(hasattr(math, 'isprime')) + self.assertFalse(hasattr(math, 'primes')) + 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..85891ca4150003f --- /dev/null +++ b/Misc/NEWS.d/next/Library/2026-07-06-14-00-00.gh-issue-153222.pRiMe1.rst @@ -0,0 +1,3 @@ +Add :func:`math.integer.isprime` for primality testing of integers +less than 2**64 and +:func:`math.integer.primes` for iterating over prime numbers. diff --git a/Modules/clinic/mathintegermodule.c.h b/Modules/clinic/mathintegermodule.c.h index 29c2a0ac902ea39..e7796fc6c20baf8 100644 --- a/Modules/clinic/mathintegermodule.c.h +++ b/Modules/clinic/mathintegermodule.c.h @@ -2,7 +2,11 @@ preserve [clinic start generated code]*/ -#include "pycore_modsupport.h" // _PyArg_CheckPositional() +#if defined(Py_BUILD_CORE) && !defined(Py_BUILD_CORE_MODULE) +# include "pycore_gc.h" // PyGC_Head +# include "pycore_runtime.h" // _Py_ID() +#endif +#include "pycore_modsupport.h" // _PyArg_UnpackKeywords() PyDoc_STRVAR(math_integer_gcd__doc__, "gcd($module, /, *integers)\n" @@ -67,6 +71,109 @@ 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 a prime number, False otherwise.\n" +"\n" +"The argument must be less than 2**64."); + +#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_primes__doc__, +"primes($module, /, start=2, stop=None)\n" +"--\n" +"\n" +"Return an iterator of the prime numbers in the range [start, stop).\n" +"\n" +"If stop is None, the iteration does not stop.\n" +"The bounds must be less than 2**64."); + +#define MATH_INTEGER_PRIMES_METHODDEF \ + {"primes", _PyCFunction_CAST(math_integer_primes), METH_FASTCALL|METH_KEYWORDS, math_integer_primes__doc__}, + +static PyObject * +math_integer_primes_impl(PyObject *module, PyObject *start, PyObject *stop); + +static PyObject * +math_integer_primes(PyObject *module, PyObject *const *args, Py_ssize_t nargs, PyObject *kwnames) +{ + PyObject *return_value = NULL; + #if defined(Py_BUILD_CORE) && !defined(Py_BUILD_CORE_MODULE) + + #define NUM_KEYWORDS 2 + static struct { + PyGC_Head _this_is_not_used; + PyObject_VAR_HEAD + Py_hash_t ob_hash; + PyObject *ob_item[NUM_KEYWORDS]; + } _kwtuple = { + .ob_base = PyVarObject_HEAD_INIT(&PyTuple_Type, NUM_KEYWORDS) + .ob_hash = -1, + .ob_item = { &_Py_ID(start), &_Py_ID(stop), }, + }; + #undef NUM_KEYWORDS + #define KWTUPLE (&_kwtuple.ob_base.ob_base) + + #else // !Py_BUILD_CORE + # define KWTUPLE NULL + #endif // !Py_BUILD_CORE + + static const char * const _keywords[] = {"start", "stop", NULL}; + static _PyArg_Parser _parser = { + .keywords = _keywords, + .fname = "primes", + .kwtuple = KWTUPLE, + }; + #undef KWTUPLE + PyObject *argsbuf[2]; + Py_ssize_t noptargs = nargs + (kwnames ? PyTuple_GET_SIZE(kwnames) : 0) - 0; + PyObject *start = NULL; + PyObject *stop = Py_None; + + args = _PyArg_UnpackKeywords(args, nargs, NULL, kwnames, &_parser, + /*minpos*/ 0, /*maxpos*/ 2, /*minkw*/ 0, /*varpos*/ 0, argsbuf); + if (!args) { + goto exit; + } + if (!noptargs) { + goto skip_optional_pos; + } + if (args[0]) { + start = args[0]; + if (!--noptargs) { + goto skip_optional_pos; + } + } + stop = args[1]; +skip_optional_pos: + return_value = math_integer_primes_impl(module, start, stop); + +exit: + return return_value; +} + PyDoc_STRVAR(math_integer_factorial__doc__, "factorial($module, n, /)\n" "--\n" @@ -156,4 +263,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=172b302a40542e1c input=a9049054013a1b77]*/ diff --git a/Modules/mathintegermodule.c b/Modules/mathintegermodule.c index 0f660d461e349f8..7978e8b5e9d4944 100644 --- a/Modules/mathintegermodule.c +++ b/Modules/mathintegermodule.c @@ -5,12 +5,25 @@ #endif #include "Python.h" -#include "pycore_abstract.h" // _PyNumber_Index() -#include "pycore_bitutils.h" // _Py_bit_length() -#include "pycore_long.h" // _PyLong_GetZero() +#include "pycore_abstract.h" // _PyNumber_Index() +#include "pycore_bitutils.h" // _Py_bit_length() +#include "pycore_critical_section.h" // Py_BEGIN_CRITICAL_SECTION() +#include "pycore_long.h" // _PyLong_GetZero() #include "clinic/mathintegermodule.c.h" +typedef struct { + PyTypeObject *primes_iter_type; +} math_integer_state; + +static inline math_integer_state * +get_math_integer_state(PyObject *module) +{ + void *state = PyModule_GetState(module); + assert(state != NULL); + return (math_integer_state *)state; +} + /*[clinic input] module math module math.integer @@ -482,6 +495,323 @@ math_integer_isqrt(PyObject *module, PyObject *n) } +/* Primality testing 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). +*/ + +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 +# 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 + +static uint64_t +powmod_u64(uint64_t b, uint64_t e, uint64_t n) +{ + uint64_t result = 1; + b %= n; + while (e) { + if (e & 1) { + result = mulmod_u64(result, b, n); + } + b = mulmod_u64(b, b, n); + e >>= 1; + } + return result; +} + +/* Strong probable prime test to the given base. n must be odd and + coprime to the base (a base divisible by n gives no information and + is treated as passing). */ +static int +sprp_u64(uint64_t n, uint64_t base) +{ + base %= n; + if (base == 0) { + return 1; + } + uint64_t d = n - 1; + int s = 0; + while ((d & 1) == 0) { + d >>= 1; + s++; + } + uint64_t a = powmod_u64(base, d, n); + if (a == 1 || a == n - 1) { + return 1; + } + for (int r = 1; r < s; r++) { + a = mulmod_u64(a, a, n); + if (a == n - 1) { + 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; + } + 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(n, bases[i])) { + return 0; + } + } + return 1; +} + +/*[clinic input] +math.integer.isprime -> bool + + n: object + / + +Return True if n is a prime number, False otherwise. + +The argument must be less than 2**64. +[clinic start generated code]*/ + +static int +math_integer_isprime_impl(PyObject *module, PyObject *n) +/*[clinic end generated code: output=c808dc14d8d86875 input=bab96f73b765d9cf]*/ +{ + 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; + } + else { + result = isprime_u64(m); + } + Py_DECREF(n); + return result; +} + + +/* The iterator returned by primes(). */ + +typedef struct { + PyObject_HEAD + uint64_t candidate; /* next candidate to test; 2 or odd */ + uint64_t stop; /* exclusive upper bound; ignored if unbounded */ + char unbounded; + char done; + char overflowed; /* unbounded iteration ran out of the domain */ +} primesiterobject; + +static void +primes_iter_dealloc(PyObject *op) +{ + PyTypeObject *tp = Py_TYPE(op); + freefunc free_func = PyType_GetSlot(tp, Py_tp_free); + free_func(op); + Py_DECREF(tp); +} + +static PyObject * +primes_iter_next(PyObject *op) +{ + primesiterobject *it = (primesiterobject *)op; + PyObject *result = NULL; + Py_BEGIN_CRITICAL_SECTION(op); + for (;;) { + if (it->done) { + if (it->overflowed) { + /* The whole supported range has been iterated without + reaching a bound. Only reachable with stop=None. */ + it->overflowed = 0; + PyErr_SetString(PyExc_OverflowError, + "primes() argument must be less than 2**64"); + } + break; + } + if (PyErr_CheckSignals() < 0) { + break; + } + uint64_t cand = it->candidate; + if (!it->unbounded && cand >= it->stop) { + it->done = 1; + break; + } + /* Advance to the next candidate: 2 -> 3, odd -> odd + 2. */ + if (cand == 2) { + it->candidate = 3; + } + else if (cand == UINT64_MAX) { + it->done = 1; + it->overflowed = it->unbounded; + } + else { + it->candidate = cand + 2; + } + if (isprime_u64(cand)) { + result = PyLong_FromUnsignedLongLong(cand); + break; + } + } + Py_END_CRITICAL_SECTION(); + return result; +} + +static PyType_Slot primes_iter_slots[] = { + {Py_tp_dealloc, primes_iter_dealloc}, + {Py_tp_iter, PyObject_SelfIter}, + {Py_tp_iternext, primes_iter_next}, + {0, NULL}, +}; + +static PyType_Spec primes_iter_spec = { + .name = "math.integer.primes_iterator", + .basicsize = sizeof(primesiterobject), + .flags = (Py_TPFLAGS_DEFAULT | Py_TPFLAGS_IMMUTABLETYPE + | Py_TPFLAGS_DISALLOW_INSTANTIATION), + .slots = primes_iter_slots, +}; + +/*[clinic input] +math.integer.primes + + start: object(c_default="NULL") = 2 + stop: object = None + +Return an iterator of the prime numbers in the range [start, stop). + +If stop is None, the iteration does not stop. +The bounds must be less than 2**64. +[clinic start generated code]*/ + +static PyObject * +math_integer_primes_impl(PyObject *module, PyObject *start, PyObject *stop) +/*[clinic end generated code: output=b92270dab45a6b84 input=b5abd91fb58a8d90]*/ +{ + uint64_t candidate = 2; + uint64_t stopval = 0; + char unbounded = 0; + char done = 0; + + if (stop == Py_None) { + unbounded = 1; + } + else { + PyObject *s = _PyNumber_Index(stop); + if (s == NULL) { + return NULL; + } + if (_PyLong_IsNegative((PyLongObject *)s)) { + done = 1; + } + else if (PyLong_AsUInt64(s, &stopval) < 0) { + Py_DECREF(s); + return NULL; + } + Py_DECREF(s); + } + if (start != NULL) { + PyObject *s = _PyNumber_Index(start); + if (s == NULL) { + return NULL; + } + if (_PyLong_IsNegative((PyLongObject *)s)) { + /* The first prime is 2. */ + } + else if (PyLong_AsUInt64(s, &candidate) < 0) { + Py_DECREF(s); + return NULL; + } + Py_DECREF(s); + if (candidate < 2) { + candidate = 2; + } + else if (candidate > 2) { + /* Round an even start up to the next odd number. */ + candidate |= 1; + } + } + + math_integer_state *state = get_math_integer_state(module); + primesiterobject *it = PyObject_New(primesiterobject, + state->primes_iter_type); + if (it == NULL) { + return NULL; + } + it->candidate = candidate; + it->stop = stopval; + it->unbounded = unbounded; + it->done = done; + it->overflowed = 0; + return (PyObject *)it; +} + + static unsigned long count_set_bits(unsigned long n) { @@ -1234,15 +1564,24 @@ 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 + MATH_INTEGER_PRIMES_METHODDEF {NULL, NULL} /* sentinel */ }; static int math_integer_exec(PyObject *module) { + math_integer_state *state = get_math_integer_state(module); + state->primes_iter_type = (PyTypeObject *)PyType_FromModuleAndSpec( + module, &primes_iter_spec, NULL); + if (state->primes_iter_type == NULL) { + return -1; + } + /* Fix the __name__ attribute of the module and the __module__ attribute * of its functions. */ @@ -1271,6 +1610,28 @@ math_integer_exec(PyObject *module) return 0; } +static int +math_integer_traverse(PyObject *module, visitproc visit, void *arg) +{ + math_integer_state *state = get_math_integer_state(module); + Py_VISIT(state->primes_iter_type); + return 0; +} + +static int +math_integer_clear(PyObject *module) +{ + math_integer_state *state = get_math_integer_state(module); + Py_CLEAR(state->primes_iter_type); + return 0; +} + +static void +math_integer_free(void *module) +{ + (void)math_integer_clear((PyObject *)module); +} + static PyModuleDef_Slot math_integer_slots[] = { _Py_ABI_SLOT, {Py_mod_exec, math_integer_exec}, @@ -1286,9 +1647,12 @@ static struct PyModuleDef math_integer_module = { PyModuleDef_HEAD_INIT, .m_name = "math.integer", .m_doc = module_doc, - .m_size = 0, + .m_size = sizeof(math_integer_state), .m_methods = math_integer_methods, .m_slots = math_integer_slots, + .m_traverse = math_integer_traverse, + .m_clear = math_integer_clear, + .m_free = math_integer_free, }; PyMODINIT_FUNC From 651f262cfd9b5f0d44a8755277b171df023336be Mon Sep 17 00:00:00 2001 From: Stan Ulbrych Date: Sun, 25 Oct 2026 17:56:23 +0000 Subject: [PATCH 2/4] Drop `math.primes` --- Doc/library/math.integer.rst | 18 -- Doc/whatsnew/3.16.rst | 3 +- .../pycore_global_objects_fini_generated.h | 1 - Include/internal/pycore_global_strings.h | 1 - .../internal/pycore_runtime_init_generated.h | 1 - .../internal/pycore_unicodeobject_generated.h | 4 - Lib/test/test_math_integer.py | 117 +--------- ...-07-06-14-00-00.gh-issue-153222.pRiMe1.rst | 3 +- Modules/clinic/mathintegermodule.c.h | 81 +------ Modules/mathintegermodule.c | 209 +----------------- 10 files changed, 11 insertions(+), 427 deletions(-) diff --git a/Doc/library/math.integer.rst b/Doc/library/math.integer.rst index d34e7f821c4df9a..16959a8a9f84f2d 100644 --- a/Doc/library/math.integer.rst +++ b/Doc/library/math.integer.rst @@ -98,21 +98,3 @@ all return values are computed exactly and are integers. and the function returns ``n!``. Raises :exc:`ValueError` if either of the arguments are negative. - - -.. function:: primes(start=2, stop=None) - - Return an iterator of the prime numbers *p* with ``start <= p < stop``, - in increasing order. - If *stop* is ``None`` (the default), the iteration does not stop. - - Roughly equivalent to - ``(p for p in itertools.count(start) if isprime(p))`` - for the unbounded form. - - The bounds must be less than 2\ :sup:`64`; - raises :exc:`OverflowError` otherwise. - Unbounded iteration raises :exc:`OverflowError` - if the candidates reach that limit. - - .. versionadded:: next diff --git a/Doc/whatsnew/3.16.rst b/Doc/whatsnew/3.16.rst index eb2537d15ce6e58..9ba2d5f94b75bec 100644 --- a/Doc/whatsnew/3.16.rst +++ b/Doc/whatsnew/3.16.rst @@ -283,8 +283,7 @@ math.integer ------------ * Added :func:`math.integer.isprime` for primality testing of integers - less than 2**64 and - :func:`math.integer.primes` for iterating over prime numbers. + less than 2**64. (Contributed by Serhiy Storchaka in :gh:`153222`.) diff --git a/Include/internal/pycore_global_objects_fini_generated.h b/Include/internal/pycore_global_objects_fini_generated.h index be16ea2d4aa522c..7d952a1e52561a6 100644 --- a/Include/internal/pycore_global_objects_fini_generated.h +++ b/Include/internal/pycore_global_objects_fini_generated.h @@ -2097,7 +2097,6 @@ _PyStaticObjects_CheckRefcnt(PyInterpreterState *interp) { _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(stdout)); _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(step)); _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(steps)); - _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(stop)); _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(store_name)); _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(strategy)); _PyStaticObject_CheckRefcnt((PyObject *)&_Py_ID(strftime)); diff --git a/Include/internal/pycore_global_strings.h b/Include/internal/pycore_global_strings.h index 91f2b1c979ee3c8..8a8bbc3b6d05bf7 100644 --- a/Include/internal/pycore_global_strings.h +++ b/Include/internal/pycore_global_strings.h @@ -820,7 +820,6 @@ struct _Py_global_strings { STRUCT_FOR_ID(stdout) STRUCT_FOR_ID(step) STRUCT_FOR_ID(steps) - STRUCT_FOR_ID(stop) STRUCT_FOR_ID(store_name) STRUCT_FOR_ID(strategy) STRUCT_FOR_ID(strftime) diff --git a/Include/internal/pycore_runtime_init_generated.h b/Include/internal/pycore_runtime_init_generated.h index 273164bc882688c..366d2d300fb4780 100644 --- a/Include/internal/pycore_runtime_init_generated.h +++ b/Include/internal/pycore_runtime_init_generated.h @@ -2095,7 +2095,6 @@ extern "C" { INIT_ID(stdout), \ INIT_ID(step), \ INIT_ID(steps), \ - INIT_ID(stop), \ INIT_ID(store_name), \ INIT_ID(strategy), \ INIT_ID(strftime), \ diff --git a/Include/internal/pycore_unicodeobject_generated.h b/Include/internal/pycore_unicodeobject_generated.h index 179a47deb2dc899..00d6297432b5fca 100644 --- a/Include/internal/pycore_unicodeobject_generated.h +++ b/Include/internal/pycore_unicodeobject_generated.h @@ -3060,10 +3060,6 @@ _PyUnicode_InitStaticStrings(PyInterpreterState *interp) { _PyUnicode_InternStatic(interp, &string); assert(_PyUnicode_CheckConsistency(string, 1)); assert(PyUnicode_GET_LENGTH(string) != 1); - string = &_Py_ID(stop); - _PyUnicode_InternStatic(interp, &string); - assert(_PyUnicode_CheckConsistency(string, 1)); - assert(PyUnicode_GET_LENGTH(string) != 1); string = &_Py_ID(store_name); _PyUnicode_InternStatic(interp, &string); assert(_PyUnicode_CheckConsistency(string, 1)); diff --git a/Lib/test/test_math_integer.py b/Lib/test/test_math_integer.py index ef0ed0b7791a4f5..69d0b3bc631eb04 100644 --- a/Lib/test/test_math_integer.py +++ b/Lib/test/test_math_integer.py @@ -1,6 +1,5 @@ from decimal import Decimal from fractions import Fraction -import itertools import random import unittest from test import support @@ -449,8 +448,8 @@ class MathTests(IntMathTests): import math as module -# isprime() and primes() exist only in math.integer, not in math, so their -# tests are not in IntMathTests (which is re-run against math above). +# 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 @@ -580,115 +579,6 @@ def test_isprime_non_integers(self): self.assertRaises(TypeError, isprime, 7, 11) -class PrimesIterTests(unittest.TestCase): - import math.integer as module - - def test_primes(self): - primes = self.module.primes - expected = primes_below(10**4) - self.assertEqual(list(primes(stop=10**4)), expected) - self.assertEqual(len(expected), 1229) - self.assertEqual(list(itertools.islice(primes(), 25)), expected[:25]) - - def test_primes_start(self): - primes = self.module.primes - for start in [-10**100, -100, -1, 0, 1, 2]: - with self.subTest(start=start): - self.assertEqual(list(primes(start, 10)), [2, 3, 5, 7]) - self.assertEqual(list(primes(3, 10)), [3, 5, 7]) - self.assertEqual(list(primes(4, 10)), [5, 7]) - self.assertEqual(list(primes(8, 12)), [11]) - self.assertEqual(list(primes(9, 12)), [11]) - self.assertEqual(list(primes(7, 8)), [7]) - - def test_primes_stop(self): - primes = self.module.primes - # The range is half-open. - self.assertEqual(list(primes(2, 2)), []) - self.assertEqual(list(primes(2, 3)), [2]) - self.assertEqual(list(primes(3, 3)), []) - self.assertEqual(list(primes(7, 7)), []) - self.assertEqual(list(primes(2, -10)), []) - self.assertEqual(list(primes(10, 5)), []) - self.assertEqual(list(primes(stop=0)), []) - - def test_primes_unbounded(self): - primes = self.module.primes - it = primes() - self.assertEqual([next(it) for _ in range(5)], [2, 3, 5, 7, 11]) - it = primes(10**6) - self.assertEqual(next(it), 1000003) - - def test_primes_huge(self): - primes = self.module.primes - boundary = 10**18 - expected = [n for n in range(boundary - 200, boundary + 200) - if py_isprime(n)] - self.assertEqual(list(primes(boundary - 200, boundary + 200)), - expected) - # The top of the supported range. - expected = [n for n in range(2**64 - 200, 2**64) if py_isprime(n)] - self.assertEqual(list(primes(2**64 - 200, 2**64 - 1)), expected) - - def test_primes_overflow(self): - primes = self.module.primes - # The bounds must be less than 2**64. - self.assertRaises(OverflowError, primes, 2**64) - self.assertRaises(OverflowError, primes, 2**64 + 100, 2**64 + 200) - self.assertRaises(OverflowError, primes, 0, 2**64) - self.assertRaises(OverflowError, primes, 10**100) - # An unbounded iterator raises when it runs out of the - # supported range. - it = primes(2**64 - 60) - self.assertEqual(next(it), 2**64 - 59) - self.assertRaises(OverflowError, next, it) - self.assertRaises(StopIteration, next, it) - - def test_primes_iterator_protocol(self): - primes = self.module.primes - it = primes(2, 10) - self.assertIs(iter(it), it) - self.assertEqual(list(it), [2, 3, 5, 7]) - # An exhausted iterator stays exhausted. - self.assertEqual(list(it), []) - self.assertRaises(StopIteration, next, it) - # The iterator type cannot be instantiated directly. - self.assertRaises(TypeError, type(primes())) - - def test_primes_keywords(self): - primes = self.module.primes - self.assertEqual(list(primes(start=10, stop=30)), [11, 13, 17, 19, 23, 29]) - self.assertEqual(list(primes(10, stop=30)), [11, 13, 17, 19, 23, 29]) - - def test_primes_integer_like(self): - primes = self.module.primes - self.assertEqual(list(primes(True, 10)), [2, 3, 5, 7]) - self.assertEqual(list(primes(IntSubclass(3), IntSubclass(10))), [3, 5, 7]) - self.assertEqual(list(primes(MyIndexable(3), MyIndexable(10))), [3, 5, 7]) - # The yielded values are exact ints. - for p in primes(IntSubclass(3), IntSubclass(10)): - self.assertIs(type(p), int) - - def test_primes_int_subclass_operators(self): - # Overridden operators of an int subclass must not affect the - # iteration. - primes = self.module.primes - self.assertEqual(list(primes(BadIntSubclass(3), BadIntSubclass(10))), - [3, 5, 7]) - big = 10**18 - self.assertEqual(list(primes(BadIntSubclass(big), big + 100)), - [n for n in range(big, big + 100) if py_isprime(n)]) - - def test_primes_non_integers(self): - primes = self.module.primes - self.assertRaises(TypeError, primes, 2.5) - self.assertRaises(TypeError, primes, 2.5, 10) - self.assertRaises(TypeError, primes, 2, 10.5) - self.assertRaises(TypeError, primes, '2') - self.assertRaises(TypeError, primes, 2, '10') - self.assertRaises(TypeError, primes, 2, 10, 3) - - class MiscTests(unittest.TestCase): def test_module_name(self): @@ -700,11 +590,10 @@ def test_module_name(self): self.assertEqual(obj.__module__, 'math.integer') def test_math_namespace(self): - # New functions are added only to math.integer, not to math + # The new function is added only to math.integer, not to math # (PEP 791). import math self.assertFalse(hasattr(math, 'isprime')) - self.assertFalse(hasattr(math, 'primes')) if __name__ == '__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 index 85891ca4150003f..d2957788ee02448 100644 --- 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 @@ -1,3 +1,2 @@ Add :func:`math.integer.isprime` for primality testing of integers -less than 2**64 and -:func:`math.integer.primes` for iterating over prime numbers. +less than 2**64. diff --git a/Modules/clinic/mathintegermodule.c.h b/Modules/clinic/mathintegermodule.c.h index e7796fc6c20baf8..e9e6e53f4b0f785 100644 --- a/Modules/clinic/mathintegermodule.c.h +++ b/Modules/clinic/mathintegermodule.c.h @@ -2,11 +2,7 @@ preserve [clinic start generated code]*/ -#if defined(Py_BUILD_CORE) && !defined(Py_BUILD_CORE_MODULE) -# include "pycore_gc.h" // PyGC_Head -# include "pycore_runtime.h" // _Py_ID() -#endif -#include "pycore_modsupport.h" // _PyArg_UnpackKeywords() +#include "pycore_modsupport.h" // _PyArg_CheckPositional() PyDoc_STRVAR(math_integer_gcd__doc__, "gcd($module, /, *integers)\n" @@ -101,79 +97,6 @@ math_integer_isprime(PyObject *module, PyObject *n) return return_value; } -PyDoc_STRVAR(math_integer_primes__doc__, -"primes($module, /, start=2, stop=None)\n" -"--\n" -"\n" -"Return an iterator of the prime numbers in the range [start, stop).\n" -"\n" -"If stop is None, the iteration does not stop.\n" -"The bounds must be less than 2**64."); - -#define MATH_INTEGER_PRIMES_METHODDEF \ - {"primes", _PyCFunction_CAST(math_integer_primes), METH_FASTCALL|METH_KEYWORDS, math_integer_primes__doc__}, - -static PyObject * -math_integer_primes_impl(PyObject *module, PyObject *start, PyObject *stop); - -static PyObject * -math_integer_primes(PyObject *module, PyObject *const *args, Py_ssize_t nargs, PyObject *kwnames) -{ - PyObject *return_value = NULL; - #if defined(Py_BUILD_CORE) && !defined(Py_BUILD_CORE_MODULE) - - #define NUM_KEYWORDS 2 - static struct { - PyGC_Head _this_is_not_used; - PyObject_VAR_HEAD - Py_hash_t ob_hash; - PyObject *ob_item[NUM_KEYWORDS]; - } _kwtuple = { - .ob_base = PyVarObject_HEAD_INIT(&PyTuple_Type, NUM_KEYWORDS) - .ob_hash = -1, - .ob_item = { &_Py_ID(start), &_Py_ID(stop), }, - }; - #undef NUM_KEYWORDS - #define KWTUPLE (&_kwtuple.ob_base.ob_base) - - #else // !Py_BUILD_CORE - # define KWTUPLE NULL - #endif // !Py_BUILD_CORE - - static const char * const _keywords[] = {"start", "stop", NULL}; - static _PyArg_Parser _parser = { - .keywords = _keywords, - .fname = "primes", - .kwtuple = KWTUPLE, - }; - #undef KWTUPLE - PyObject *argsbuf[2]; - Py_ssize_t noptargs = nargs + (kwnames ? PyTuple_GET_SIZE(kwnames) : 0) - 0; - PyObject *start = NULL; - PyObject *stop = Py_None; - - args = _PyArg_UnpackKeywords(args, nargs, NULL, kwnames, &_parser, - /*minpos*/ 0, /*maxpos*/ 2, /*minkw*/ 0, /*varpos*/ 0, argsbuf); - if (!args) { - goto exit; - } - if (!noptargs) { - goto skip_optional_pos; - } - if (args[0]) { - start = args[0]; - if (!--noptargs) { - goto skip_optional_pos; - } - } - stop = args[1]; -skip_optional_pos: - return_value = math_integer_primes_impl(module, start, stop); - -exit: - return return_value; -} - PyDoc_STRVAR(math_integer_factorial__doc__, "factorial($module, n, /)\n" "--\n" @@ -263,4 +186,4 @@ math_integer_comb(PyObject *module, PyObject *const *args, Py_ssize_t nargs) exit: return return_value; } -/*[clinic end generated code: output=172b302a40542e1c input=a9049054013a1b77]*/ +/*[clinic end generated code: output=bb518a64a75a30f1 input=a9049054013a1b77]*/ diff --git a/Modules/mathintegermodule.c b/Modules/mathintegermodule.c index 7978e8b5e9d4944..5f682e870c1d418 100644 --- a/Modules/mathintegermodule.c +++ b/Modules/mathintegermodule.c @@ -5,25 +5,12 @@ #endif #include "Python.h" -#include "pycore_abstract.h" // _PyNumber_Index() -#include "pycore_bitutils.h" // _Py_bit_length() -#include "pycore_critical_section.h" // Py_BEGIN_CRITICAL_SECTION() -#include "pycore_long.h" // _PyLong_GetZero() +#include "pycore_abstract.h" // _PyNumber_Index() +#include "pycore_bitutils.h" // _Py_bit_length() +#include "pycore_long.h" // _PyLong_GetZero() #include "clinic/mathintegermodule.c.h" -typedef struct { - PyTypeObject *primes_iter_type; -} math_integer_state; - -static inline math_integer_state * -get_math_integer_state(PyObject *module) -{ - void *state = PyModule_GetState(module); - assert(state != NULL); - return (math_integer_state *)state; -} - /*[clinic input] module math module math.integer @@ -657,161 +644,6 @@ math_integer_isprime_impl(PyObject *module, PyObject *n) } -/* The iterator returned by primes(). */ - -typedef struct { - PyObject_HEAD - uint64_t candidate; /* next candidate to test; 2 or odd */ - uint64_t stop; /* exclusive upper bound; ignored if unbounded */ - char unbounded; - char done; - char overflowed; /* unbounded iteration ran out of the domain */ -} primesiterobject; - -static void -primes_iter_dealloc(PyObject *op) -{ - PyTypeObject *tp = Py_TYPE(op); - freefunc free_func = PyType_GetSlot(tp, Py_tp_free); - free_func(op); - Py_DECREF(tp); -} - -static PyObject * -primes_iter_next(PyObject *op) -{ - primesiterobject *it = (primesiterobject *)op; - PyObject *result = NULL; - Py_BEGIN_CRITICAL_SECTION(op); - for (;;) { - if (it->done) { - if (it->overflowed) { - /* The whole supported range has been iterated without - reaching a bound. Only reachable with stop=None. */ - it->overflowed = 0; - PyErr_SetString(PyExc_OverflowError, - "primes() argument must be less than 2**64"); - } - break; - } - if (PyErr_CheckSignals() < 0) { - break; - } - uint64_t cand = it->candidate; - if (!it->unbounded && cand >= it->stop) { - it->done = 1; - break; - } - /* Advance to the next candidate: 2 -> 3, odd -> odd + 2. */ - if (cand == 2) { - it->candidate = 3; - } - else if (cand == UINT64_MAX) { - it->done = 1; - it->overflowed = it->unbounded; - } - else { - it->candidate = cand + 2; - } - if (isprime_u64(cand)) { - result = PyLong_FromUnsignedLongLong(cand); - break; - } - } - Py_END_CRITICAL_SECTION(); - return result; -} - -static PyType_Slot primes_iter_slots[] = { - {Py_tp_dealloc, primes_iter_dealloc}, - {Py_tp_iter, PyObject_SelfIter}, - {Py_tp_iternext, primes_iter_next}, - {0, NULL}, -}; - -static PyType_Spec primes_iter_spec = { - .name = "math.integer.primes_iterator", - .basicsize = sizeof(primesiterobject), - .flags = (Py_TPFLAGS_DEFAULT | Py_TPFLAGS_IMMUTABLETYPE - | Py_TPFLAGS_DISALLOW_INSTANTIATION), - .slots = primes_iter_slots, -}; - -/*[clinic input] -math.integer.primes - - start: object(c_default="NULL") = 2 - stop: object = None - -Return an iterator of the prime numbers in the range [start, stop). - -If stop is None, the iteration does not stop. -The bounds must be less than 2**64. -[clinic start generated code]*/ - -static PyObject * -math_integer_primes_impl(PyObject *module, PyObject *start, PyObject *stop) -/*[clinic end generated code: output=b92270dab45a6b84 input=b5abd91fb58a8d90]*/ -{ - uint64_t candidate = 2; - uint64_t stopval = 0; - char unbounded = 0; - char done = 0; - - if (stop == Py_None) { - unbounded = 1; - } - else { - PyObject *s = _PyNumber_Index(stop); - if (s == NULL) { - return NULL; - } - if (_PyLong_IsNegative((PyLongObject *)s)) { - done = 1; - } - else if (PyLong_AsUInt64(s, &stopval) < 0) { - Py_DECREF(s); - return NULL; - } - Py_DECREF(s); - } - if (start != NULL) { - PyObject *s = _PyNumber_Index(start); - if (s == NULL) { - return NULL; - } - if (_PyLong_IsNegative((PyLongObject *)s)) { - /* The first prime is 2. */ - } - else if (PyLong_AsUInt64(s, &candidate) < 0) { - Py_DECREF(s); - return NULL; - } - Py_DECREF(s); - if (candidate < 2) { - candidate = 2; - } - else if (candidate > 2) { - /* Round an even start up to the next odd number. */ - candidate |= 1; - } - } - - math_integer_state *state = get_math_integer_state(module); - primesiterobject *it = PyObject_New(primesiterobject, - state->primes_iter_type); - if (it == NULL) { - return NULL; - } - it->candidate = candidate; - it->stop = stopval; - it->unbounded = unbounded; - it->done = done; - it->overflowed = 0; - return (PyObject *)it; -} - - static unsigned long count_set_bits(unsigned long n) { @@ -1568,20 +1400,12 @@ static PyMethodDef math_integer_methods[] = { MATH_INTEGER_ISQRT_METHODDEF MATH_INTEGER_LCM_METHODDEF MATH_INTEGER_PERM_METHODDEF - MATH_INTEGER_PRIMES_METHODDEF {NULL, NULL} /* sentinel */ }; static int math_integer_exec(PyObject *module) { - math_integer_state *state = get_math_integer_state(module); - state->primes_iter_type = (PyTypeObject *)PyType_FromModuleAndSpec( - module, &primes_iter_spec, NULL); - if (state->primes_iter_type == NULL) { - return -1; - } - /* Fix the __name__ attribute of the module and the __module__ attribute * of its functions. */ @@ -1610,28 +1434,6 @@ math_integer_exec(PyObject *module) return 0; } -static int -math_integer_traverse(PyObject *module, visitproc visit, void *arg) -{ - math_integer_state *state = get_math_integer_state(module); - Py_VISIT(state->primes_iter_type); - return 0; -} - -static int -math_integer_clear(PyObject *module) -{ - math_integer_state *state = get_math_integer_state(module); - Py_CLEAR(state->primes_iter_type); - return 0; -} - -static void -math_integer_free(void *module) -{ - (void)math_integer_clear((PyObject *)module); -} - static PyModuleDef_Slot math_integer_slots[] = { _Py_ABI_SLOT, {Py_mod_exec, math_integer_exec}, @@ -1647,12 +1449,9 @@ static struct PyModuleDef math_integer_module = { PyModuleDef_HEAD_INIT, .m_name = "math.integer", .m_doc = module_doc, - .m_size = sizeof(math_integer_state), + .m_size = 0, .m_methods = math_integer_methods, .m_slots = math_integer_slots, - .m_traverse = math_integer_traverse, - .m_clear = math_integer_clear, - .m_free = math_integer_free, }; PyMODINIT_FUNC From 0392f1391b7ca9425171ef729e526773706ec6ce Mon Sep 17 00:00:00 2001 From: Stan Ulbrych Date: Sun, 25 Oct 2026 18:57:35 +0000 Subject: [PATCH 3/4] Use Montgomery multiplication in `isprime()` --- Doc/whatsnew/3.16.rst | 2 +- Lib/test/test_math_integer.py | 21 +++++-- Modules/mathintegermodule.c | 114 +++++++++++++++++++++++++--------- 3 files changed, 103 insertions(+), 34 deletions(-) diff --git a/Doc/whatsnew/3.16.rst b/Doc/whatsnew/3.16.rst index 9ba2d5f94b75bec..b28109a12f8dccc 100644 --- a/Doc/whatsnew/3.16.rst +++ b/Doc/whatsnew/3.16.rst @@ -284,7 +284,7 @@ math.integer * Added :func:`math.integer.isprime` for primality testing of integers less than 2**64. - (Contributed by Serhiy Storchaka in :gh:`153222`.) + (Contributed by Serhiy Storchaka and Stan Ulbrych in :gh:`153222`.) os diff --git a/Lib/test/test_math_integer.py b/Lib/test/test_math_integer.py index 69d0b3bc631eb04..d0c8ada9da6c86f 100644 --- a/Lib/test/test_math_integer.py +++ b/Lib/test/test_math_integer.py @@ -478,7 +478,7 @@ def test_isprime_carmichael(self): 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 strong Lucas test. + # (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, @@ -487,8 +487,8 @@ def test_isprime_strong_pseudoprimes_base_2(self): self.assertIs(isprime(n), False) def test_isprime_strong_lucas_pseudoprimes(self): - # Composites that pass the strong Lucas test with Selfridge's - # parameters (A217255). + # 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]: @@ -519,9 +519,9 @@ def test_isprime_perfect_squares(self): self.assertIs(isprime(k*k), False) def test_isprime_word_boundaries(self): - # Exercise the boundaries between the base sets. + # Exercise word boundaries and the boundary between the base sets. isprime = self.module.isprime - for boundary in [2**32, 4759123141]: + 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)) @@ -530,6 +530,17 @@ def test_isprime_word_boundaries(self): 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 diff --git a/Modules/mathintegermodule.c b/Modules/mathintegermodule.c index 5f682e870c1d418..7a5e63dbe09dd2c 100644 --- a/Modules/mathintegermodule.c +++ b/Modules/mathintegermodule.c @@ -491,6 +491,11 @@ math_integer_isqrt(PyObject *module, PyObject *n) 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[] = { @@ -540,44 +545,96 @@ mulmod_u64(uint64_t a, uint64_t b, uint64_t n) } #endif -static uint64_t -powmod_u64(uint64_t b, uint64_t e, uint64_t n) +/* High 64 bits of the full product a * b. */ +static inline uint64_t +mulhi_u64(uint64_t a, uint64_t b) { - uint64_t result = 1; - b %= n; - while (e) { - if (e & 1) { - result = mulmod_u64(result, b, n); - } - b = mulmod_u64(b, b, n); - e >>= 1; - } - return result; +#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 } -/* Strong probable prime test to the given base. n must be odd and - coprime to the base (a base divisible by n gives no information and - is treated as passing). */ +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(uint64_t n, uint64_t base) +sprp_u64(const montgomery_context *ctx, uint64_t base) { - base %= n; + base %= ctx->n; if (base == 0) { return 1; } - uint64_t d = n - 1; - int s = 0; - while ((d & 1) == 0) { - d >>= 1; - s++; + // 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); } - uint64_t a = powmod_u64(base, d, n); - if (a == 1 || a == n - 1) { + if (a == ctx->one || a == ctx->n - ctx->one) { return 1; } - for (int r = 1; r < s; r++) { - a = mulmod_u64(a, a, n); - if (a == n - 1) { + for (int r = 1; r < ctx->s; r++) { + a = montgomery_mul(a, a, ctx); + if (a == ctx->n - ctx->one) { return 1; } } @@ -595,6 +652,7 @@ isprime_u64(uint64_t n) 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) { @@ -602,7 +660,7 @@ isprime_u64(uint64_t n) num_bases = Py_ARRAY_LENGTH(mr_bases_small); } for (size_t i = 0; i < num_bases; i++) { - if (!sprp_u64(n, bases[i])) { + if (!sprp_u64(&ctx, bases[i])) { return 0; } } From ded345313ba96b60e6aaafe2be555fb2c69b947d Mon Sep 17 00:00:00 2001 From: Stan Ulbrych Date: Mon, 26 Oct 2026 09:54:51 +0000 Subject: [PATCH 4/4] Support arbitrary-size integers --- Doc/library/math.integer.rst | 15 +- Doc/whatsnew/3.16.rst | 4 +- Lib/test/test_math_integer.py | 51 +++-- ...-07-06-14-00-00.gh-issue-153222.pRiMe1.rst | 4 +- Modules/clinic/mathintegermodule.c.h | 7 +- Modules/mathintegermodule.c | 180 +++++++++++++++++- 6 files changed, 230 insertions(+), 31 deletions(-) diff --git a/Doc/library/math.integer.rst b/Doc/library/math.integer.rst index 16959a8a9f84f2d..752790dde8e11f1 100644 --- a/Doc/library/math.integer.rst +++ b/Doc/library/math.integer.rst @@ -49,14 +49,17 @@ all return values are computed exactly and are integers. .. function:: isprime(n, /) - Return ``True`` if *n* is a prime number, ``False`` otherwise. + Return ``True`` if *n* is prime or probably prime, ``False`` otherwise. - A prime number is a natural number greater than 1 - that is not a product of two smaller natural numbers. - Negative numbers, ``0`` and ``1`` are not prime. + 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. - The argument must be less than 2\ :sup:`64`; - raises :exc:`OverflowError` otherwise. + 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 diff --git a/Doc/whatsnew/3.16.rst b/Doc/whatsnew/3.16.rst index 44a90323a1e20ac..c2bf230b2416e87 100644 --- a/Doc/whatsnew/3.16.rst +++ b/Doc/whatsnew/3.16.rst @@ -568,8 +568,8 @@ math math.integer ------------ -* Added :func:`math.integer.isprime` for primality testing of integers - less than 2**64. +* 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`.) diff --git a/Lib/test/test_math_integer.py b/Lib/test/test_math_integer.py index d0c8ada9da6c86f..7cf53bc185a0a31 100644 --- a/Lib/test/test_math_integer.py +++ b/Lib/test/test_math_integer.py @@ -93,9 +93,8 @@ def primes_below(n): 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. More bases are used for larger inputs, for -# which the test is probabilistic (but independent of the Baillie-PSW test -# used in the implementation). +# 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): @@ -482,7 +481,9 @@ def test_isprime_strong_pseudoprimes_base_2(self): isprime = self.module.isprime for n in [2047, 3277, 4033, 4681, 8321, 15841, 29341, 42799, 49141, 52633, 65281, 74665, 80581, 85489, 88357, 90751, - 3825123056546413051]: + 3825123056546413051, + 318665857834031151167461, + 3317044064679887385961981]: with self.subTest(n=n): self.assertIs(isprime(n), False) @@ -514,7 +515,7 @@ def test_isprime_perfect_squares(self): 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]: + 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) @@ -525,7 +526,7 @@ def test_isprime_word_boundaries(self): 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): + 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 @@ -546,10 +547,28 @@ def test_isprime_large_values(self): self.assertIs(isprime(2**61 - 1), True) # Mersenne prime self.assertIs(isprime(2**62 - 1), False) self.assertIs(isprime(10**19), False) - # Arguments not less than 2**64 are not supported. - for n in [2**64, 2**64 + 13, 2**89 - 1, 10**100]: + # 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.assertRaises(OverflowError, isprime, 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 @@ -559,10 +578,11 @@ def test_isprime_random(self): n = rng.getrandbits(bits) with self.subTest(n=n): self.assertIs(isprime(n), py_isprime(n)) - for bits in [65, 80, 128]: - n = rng.getrandbits(bits) | (1 << (bits - 1)) - with self.subTest(n=n): - self.assertRaises(OverflowError, 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 @@ -572,6 +592,9 @@ def test_isprime_integer_like(self): 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 @@ -580,6 +603,8 @@ def test_isprime_int_subclass_operators(self): 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 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 index d2957788ee02448..27016e8dd7ccfac 100644 --- 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 @@ -1,2 +1,2 @@ -Add :func:`math.integer.isprime` for primality testing of integers -less than 2**64. +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 e9e6e53f4b0f785..9230faaaf3c3cee 100644 --- a/Modules/clinic/mathintegermodule.c.h +++ b/Modules/clinic/mathintegermodule.c.h @@ -71,9 +71,10 @@ PyDoc_STRVAR(math_integer_isprime__doc__, "isprime($module, n, /)\n" "--\n" "\n" -"Return True if n is a prime number, False otherwise.\n" +"Return True if n is prime or probably prime, False otherwise.\n" "\n" -"The argument must be less than 2**64."); +"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__}, @@ -186,4 +187,4 @@ math_integer_comb(PyObject *module, PyObject *const *args, Py_ssize_t nargs) exit: return return_value; } -/*[clinic end generated code: output=bb518a64a75a30f1 input=a9049054013a1b77]*/ +/*[clinic end generated code: output=4e47e334e7fe1156 input=a9049054013a1b77]*/ diff --git a/Modules/mathintegermodule.c b/Modules/mathintegermodule.c index 054d3e9745073c1..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,7 +482,7 @@ math_integer_isqrt(PyObject *module, PyObject *n) } -/* Primality testing with the deterministic Miller-Rabin test. +/* 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 @@ -666,22 +667,187 @@ isprime_u64(uint64_t n) 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 a prime number, False otherwise. +Return True if n is prime or probably prime, False otherwise. -The argument must be less than 2**64. +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=bab96f73b765d9cf]*/ +/*[clinic end generated code: output=c808dc14d8d86875 input=e771fc2bdb5e29fb]*/ { - n = _PyNumber_Index(n); + /* The arbitrary-precision path must not call int subclass operators. */ + n = PyNumber_Index(n); if (n == NULL) { return -1; } @@ -692,6 +858,10 @@ math_integer_isprime_impl(PyObject *module, PyObject *n) } 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);