[StanFromIreland]

This builds (and trims) on #153224.

As I suggest in #153224 (comment), it drops math.integer.primes().

And, the strong probable prime test now runs in Montgomery form. A context (inverse of n mod 2**64, 2**64 mod n, and the odd part of n - 1) is computed once per n, and each modular multiplication costs two 64-bit high products instead of a 128-bit division. The reduction is then written subtractively so no intermediate exceeds 64 bits, which keeps n close to 2**64. FLINT made the same change to n_is_prime in flintlib/flint#2519 and reported 1.2x to 1.5x speedups for 33 to 64 bit inputs. Benchmarking our isprime on an optimised build (compared to the original PR), I get a nice speedup:

isprime[small]: Mean +- std dev: [before] 305 us +- 11 us -> [after] 303 us +- 8 us: 1.01x faster
isprime[below_2**32]: Mean +- std dev: [before] 699 us +- 13 us -> [after] 606 us +- 19 us: 1.15x faster
isprime[below_2**63]: Mean +- std dev: [before] 1.63 ms +- 0.04 ms -> [after] 1.33 ms +- 0.03 ms: 1.23x faster
isprime[below_2**64]: Mean +- std dev: [before] 1.68 ms +- 0.03 ms -> [after] 1.51 ms +- 0.03 ms: 1.11x faster

Geometric mean: 1.12x faster

Benchmark script

import pyperf
from math.integer import isprime

N = 10_000

# name: (start, number of primes in [start, start + N))
RANGES = {
    'small': (2, 1229),
    'below_2**32': (2**32 - N, 447),
    'below_2**63': (2**63 - N, 216),
    'below_2**64': (2**64 - N, 218),
}


def count_primes(start, expected):
    count = sum(map(isprime, range(start, start + N)))
    assert count == expected, (start, count, expected)
    return count


runner = pyperf.Runner()
for name, (start, expected) in RANGES.items():
    runner.bench_func(f'isprime[{name}]', count_primes, start, expected)

CC @mdickinson @tim-one

• Issue: Add math.integer.isprime() and math.integer.primes() #153222

[read-the-docs-community]

Documentation build overview

📚 cpython-previews | 🛠️ Build #35409304 | 📁 Comparing ded3453 against main (bb7e16a)

  🔍 Preview build  

5 files changed · ± 5 modified

± Modified

• glossary.html
• library/math.integer.html
• reference/grammar.html
• whatsnew/3.16.html
• whatsnew/changelog.html

[skirpichev]

I think this change is not required, no? isprime() returns an integer (boolean) and it's exact.

[skirpichev]

I don't think it's an acceptable API.

1. it isn't compatible with gmpy2 (is_prime() returns True for probable primes)
2. can't be used as an alternative to gmpy2 in applications like SymPy or Mpmath, currently using own isprime() helpers.

I suggest using Miller-Rabin test without constraints on the number range. You can document that test is certain for integers < 2^64, perhaps as an implementation detail. Or introduce instead API like gmpy2.is_probab_prime (ternary logic: definitely prime, probably prime, or definitely composite).

[picnixz]

For primality testing we should definitely use MR testing and make sure that we just do not limit ourselves to 64 bits. I think we should rather keep the other PR.

[StanFromIreland]

For primality testing we should definitely use MR testing and make sure that we just do not limit ourselves to 64 bits. I think we should rather keep the other PR.

Ok, I added a probabilistic M-R path for arbitrary-size integers, but kept the exact Miller–Rabin test with Montgomery arithmetic for inputs below 2**64. I considered Baillie-PSW, but that is quite complex, so I'm not sure its worth the maintenance cost here.

[picnixz]

I still wonder why you created a new PR. Why not creating a PR againt Seehiy's?

[picnixz]

And in general, I weas still unconvinced by this API. I would rather have a distinction between provable and probable primes.

[StanFromIreland]

I still wonder why you created a new PR. Why not creating a PR againt Seehiy's?

It's a little easier because I have permissions in this repository, I can run the CI, request reviews etc. It also allows Serhiy to cherry-pick the commits (which is why I split the work into various steps).

[StanFromIreland]

And in general, I weas still unconvinced by this API. I would rather have a distinction between provable and probable primes.

Hmm, what would your preferred API be? The current one is the same as SymPy's isprime(), and I personally find it reasonable, since we document the range to which primes are provable.

Edit, also see: #153222 (comment)