Prime Generator

Back

Recently I wrote about the RSA algorithm and showed some python code to experiment with that. The starting point of RSA is having 2 large prime numbers, which I took from a website. I was interested in how one can generate these big numbers. I came across the Miller-Rabin primality test.

Miller-Rabin Test

It's a probabilistic test - it can have false positives but not false negatives. There is a parameter k such that an inner loop does k iterations. Complexity is O(k log^3 n) where n is the input number and error rate is at most 4^(-k).

Here is a simple python implementation. The idea is pick a random number with d bits using the python function random.getrandbits and then test it using the Miller-Rabin test. Choosing k > 10 should be good enough.

import random

def prime_test(n,k):
    r = 0
    n2 = n - 1
    while n2 & 1 == 0:
        r += 1
        n2 = n2 >> 1        
    d = n2 << 1
    
    for _ in range(k):
        a = random.randint(2, n - 2)

        x = pow(a, d, n)
        if x == 1 or x == n - 1:
            continue

        for _ in range(r-1):
            x = pow(x,2,n)
            if x == n - 1:
                break

        return "Composite"
    
    return "Probably Prime"

You run the test like this:

bit=500
while True:
    num = random.getrandbits(bit)
    if prime_test(num, 20) == 'Probably Prime':
        print(num)
        break    

On my laptop it takes about 5ms for 100bit number, 200ms for 500bit and a few seconds for 1000bit. I think that a C implementation could be at least 10x faster or more. Maybe I'll try implement it in Julia, for fun.

-------------------------------------

Addendum

As I thought, naive implementation in Julia runs about 10x faster, when comparing the same amount of "while True" iterations. I never coded in Julia before; it's really like Python, but without the ":" at the end of lines and with "end" at the end of blocks. pow is powermod, and I had to improvise when randomly choosing large integers. In Fact Python is very nice in that regard: integers are actually arbitrary-precision; in Julia they have BigInt but in Python it's transparent. In C, C++ I would have had to write it myself.