Working with large primes in Python
Question: 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?


Answer: 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](http://en.wikipedia.org/wiki/Miller-Rabin_primality_test) algorithm to
be a good choice. It is a little slower the
[Fermat](http://en.wikipedia.org/wiki/Fermat_primality_test) 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 += 1
        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.

See: [fermat.py](https://gist.github.com/bnlucas/5857437),
[miller_rabin.py](https://gist.github.com/bnlucas/5857478), and
[solovay_strassen.py](https://gist.github.com/bnlucas/5857525) on
_gist.github.com_

