Simple implementation of Miller-Rabin primality test

main.py

from miller_rabin import *

## find the largest 1024-bit prime
## naive

n = 2 ** 1024 - 1
s = 40
for i in range(n, 2, -2):
    if miller_rabin(i, 40):
        print(f"{i} is the largest prime.")
        break

miller_rabin.py

import random
def miller_rabin(n, s):

    #####################################################
    ###  simple implementation of miller-rabin test
    ###  n: candidate prime
    ###  s: number of rounds
    #####################################################

    ### see http://stackoverflow.com/questions/6325576/how-many-iterations-of-rabin-miller-should-i-use-for-cryptographic-safe-primes
    ### for optimal rounds

    ### openssl prime $n
        
    if n == 2:
        return True 
    if n % 2 == 0:
        return False 

    ## n-1 = 2^u x r
    ## caculate u, r
    t = n-1
    u = 0
    while t % 2 == 0:
        u += 1
        t //= 2
    r = t

    ## start testing
    for i in range(s):
        a = random.randint(2, n-2)
        z = pow(a, r, n) # instead of a ** r % n
        if z != 1 and z != n-1:
            for j in range(1, u):
                z = z ** 2 % n 
                if z == n-1:
                    break 
                if z == 1:
                    return False 
            if z != n-1:
                return False 
          
    return True