primes.py

import math

class PrimesConstants:
  SMALL_PRIMES = [
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 
    53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 
    109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 
    173, 179, 181, 191, 193, 197, 199
  ]

  PRIME_GAPS = [
    10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 
    4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 
    4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2
  ]

def is_prime(n):
  if n < 2: 
    return False
  
  for a in PrimesConstants.SMALL_PRIMES:
    b = n // a
    if b < a: return True
    if a * b == n: return False

  a = 211
  while True:
    for gap in PrimesConstants.PRIME_GAPS:
      b = n // a
      if b < a: return True
      if a * b == n: return False
      a += gap

def get_factors(n):
  for a in PrimesConstants.SMALL_PRIMES:
    while True:
      b = n // a
      if b < a: yield n; return
      if b * a != n: break
      yield a; n = b

  a = 211
  while True:
    for gap in PrimesConstants.PRIME_GAPS:
      while True:
        b = n // a
        if b < a: yield n; return
        if b * a != n: break
        yield a; n = b
      a += gap

def list_primes():
  for a in PrimesConstants.SMALL_PRIMES:
    yield a

  a = 211
  while True:
    for gap in PrimesConstants.PRIME_GAPS:
      if is_prime(a):
        yield a
      a += gap

def miller_rabin_test_bases(n, bases):
  def exp_mod(base, power, mod):
    base %= mod
    result = 1
  
    while power > 0:
      if power % 2 == 1:
        result = (result * base) % mod
      base = (base * base) % mod
      power >>= 1
  
    return result
  
  def test_base(n, base):
    if n == base:
      return is_prime(n)
    
    d = n - 1
    s = 0
  
    while d % 2 == 0:
      s += 1
      d //= 2
  
    x = exp_mod(base, d, n)
    if x == 1 or x == n - 1:
      return True
  
    for _ in range(1, s):
      x *= x
      x %= n
      if x == n - 1:
        return True
  
    return False

  for base in bases:
    if not test_base(n, base):
      return False
  return True

def miller_rabin_test(n):
  for m, l in [
    (2047, [2]),
    (1373653, [2, 3]),
    (9080191, [31, 73]),
    (25326001, [2, 3, 5]),
    (3215031751, [2, 3, 5, 7]),
    (4759123141, [2, 7, 61]),
    (1122004669633, [2, 13, 23, 1662803]),
    (2152302898747, [2, 3, 5, 7, 11]),
    (3474749660383, [2, 3, 5, 7, 11, 13]),
    (341550071728321, [2, 3, 5, 7, 11, 13, 17]),
    (3825123056546413051, [2, 3, 5, 7, 11, 13, 17, 19, 23]),
    (318665857834031151167461, [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]),
    (3317044064679887385961981, [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41])
  ]:
    if n < m:
      return miller_rabin_test_bases(n, l)
  return miller_rabin_test_bases(n, range(2, int(2 * math.log(n) ** 2)))