>> Regardless, how can we best move forward with this idea?

> Provide a pull request.

Hi, I looked into what scientific programs where doing. Most of them uses a form of the Baillie–PSW primality test which is a probabilistic primality test that's never wrong up to 2**64 and for which their is currently no known pseudoprime.

This first version I wrote uses a deterministic variant of the Miller-Rabin test for n < 3317044064679887385961981. For larger n, a probabilistic Miller-Rabin test is used with s=25 random bases. The probabilistic Miller-Rabin test is never wrong when it concludes that n is composite and has a probability of error of 4^(-s) when it concludes that n is prime.

The implementations of next_prime() and previous_prime() are straightforward and factorise() uses the Phollard's rho heuristic which gives satisfactory results for numbers with small factors. It's a generator as it may hang when n has very large factors e.g. 2**(2**8)+1.

I implemented all functions Steven D'Aprano suggested but did not bother with the sieve as Tim Peters already provided one in the python-ideas thread. The code is in imath for now but I can move it.

Hopefully this is enough to bikeshed the API, if this proposal is accepted I will write more tests and fix any bug found.