What is an efficient way for working with large prime numbers with Python? You search on here or on google, and you find many different methods for doing so... sieves, primality test algorithms... Which ways work for larger primes?

解决方案

For determining if a number is a prime, there a sieves and primality tests.

# for large numbers, xrange will throw an error.

# OverflowError: Python int too large to convert to C long

# to get over this:

def mrange(start, stop, step):

while start < stop:

yield start

start += step

# benchmarked on an old single-core system with 2GB RAM.

from math import sqrt

def is_prime(num):

if num == 2:

return True

if (num < 2) or (num % 2 == 0):

return False

return all(num % i for i in mrange(3, int(sqrt(num)) + 1, 2))

# benchmark is_prime(100**10-1) using mrange

# 10000 calls, 53191 per second.

# 60006 function calls in 0.190 seconds.

This seems to be the fastest. There is another version using not any that you see,

def is_prime(num)

# ...

return not any(num % i == 0 for i in mrange(3, int(sqrt(num)) + 1, 2))

However, in the benchmarks I got 70006 function calls in 0.272 seconds. over the use of all 60006 function calls in 0.190 seconds. while testing if 100**10-1 was prime.

If you're needing to find the next highest prime, this method will not work for you. You need to go with a primality test, I have found the Miller-Rabin algorithm to be a good choice. It is a little slower the Fermat method, but more accurate against pseudoprimes. Using the above mentioned method takes +5 minutes on this system.

Miller-Rabin algorithm:

from random import randrange

def is_prime(n, k=10):

if n == 2:

return True

if not n & 1:

return False

def check(a, s, d, n):

x = pow(a, d, n)

if x == 1:

return True

for i in xrange(s - 1):

if x == n - 1:

return True

x = pow(x, 2, n)

return x == n - 1

s = 0

d = n - 1

while d % 2 == 0:

d >>= 1

s += 1

for i in xrange(k):

a = randrange(2, n - 1)

if not check(a, s, d, n):

return False

return True

Fermat algoithm:

def is_prime(num):

if num == 2:

return True

if not num & 1:

return False

return pow(2, num-1, num) == 1

To get the next highest prime:

def next_prime(num):

if (not num & 1) and (num != 2):

num += 1

if is_prime(num):

num += 2

while True:

if is_prime(num):

break

num += 2

return num

print next_prime(100**10-1) # returns `100000000000000000039`

# benchmark next_prime(100**10-1) using Miller-Rabin algorithm.

1000 calls, 337 per second.

258669 function calls in 2.971 seconds

Using the Fermat test, we got a benchmark of 45006 function calls in 0.885 seconds., but you run a higher chance of pseudoprimes.

So, if just needing to check if a number is prime or not, the first method for is_prime works just fine. It is the fastest, if you use the mrange method with it.

Ideally, you would want to store the primes generated by next_prime and just read from that.

For example, using next_prime with the Miller-Rabin algorithm:

print next_prime(10^301)

# prints in 2.9s on the old single-core system, opposed to fermat's 2.8s

1000000000000000000000000000000000000000000000000000000000000000000000000000

0000000000000000000000000000000000000000000000000000000000000000000000000000

0000000000000000000000000000000000000000000000000000000000000000000000000000

00000000000000000000000000000000000000000000000000000000000000000000000531

You wouldn't be able to do this with return all(num % i for i in mrange(3, int(sqrt(num)) + 1, 2)) in a timely fashion. I can't even do it on this old system.

And to be sure that next_prime(10^301) and Miller-Rabin yields a correct value, this was also tested using the Fermat and the Solovay-Strassen algorithms.

Edit: Fixed a bug in next_prime