Deterministic Miller-Rabin primality test

miller_rabin.py

""" Deterministic Miller-Rabin primality test

    See https://en.m.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
    and https://miller-rabin.appspot.com/

    Written by PM 2Ring 2015.04.29
    Updated 2022.04.20
"""

small_primes = [
  2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,
]

# Smallest strong pseudoprime to the first k prime bases
# From https://oeis.org/A014233
psi = [
  2047, # 2
  1373653, # 3
  25326001, # 5
  3215031751, # 7
  2152302898747, # 11
  3474749660383, # 13
  341550071728321, # 17
  341550071728321, # 19
  3825123056546413051, # 23
  3825123056546413051, # 29
  3825123056546413051, # 31
  318665857834031151167461, # 37
  3317044064679887385961981, # 41
]

# Miller-Rabin primality test.
def is_prime_mr(n):
    if n < 2:
        return False

    # Test small prime factors. This also ensures
    # that n is coprime to every witness
    for p in small_primes:
        if n % p == 0:
            return n == p
    if n < 43 * 43:
       return True

    # Find d & s: m = d * 2**s, d odd
    m = n - 1
    # m & -m is the highest power of 2 dividing m
    s = (m & -m).bit_length() - 1
    d = m >> s
    srange = range(s - 1)

    # Loop over witnesses & perform strong pseudoprime tests
    for a, spp in zip(small_primes, psi):
        x = pow(a, d, n)
        if 1 < x < m:
            for _ in srange:
                x = x * x % n
                if x == 1:
                    # Previous (x+1)(x-1) == 0 mod n
                    return False
                if x == m:
                    # Looks good
                    break
            else:
                return False
        # else x == 1 or x == m
        if n < spp:
            # n must be prime
            break
    return True

Live demo on SageMathCell