Implementing Miller-Rabin Algorithm on Python

Question

I'm trying to implement the Miller-Rabin algorithm on Python.
I have coded just as my textbook's pseudocode says, but for some reason, it it not working as I have expected.
To be sepcific, the function 'test' returns sometimes return 'true' when it goes to Fermat test.

def miller_rabin(n, s):
    if n == 2:
        return Prime
    elif n % 2 == 0:
        return Composite

    for _ in range(s):
        a = random.randint(1, n-1)
        if test(a, n) == True:
            return Composite

    return Prime

def test(a, n):
    t, u = 0, n-1
    while (u % 2 == 0):
        t += 1                  #t >= 1, u is odd, n=1 = 2^t * u
        u //= 2                 #initialization
    x = exp(a, u, n)           #initializing x0 = a^u mod n

    for _ in range(t-1):        #for i = 1 to t
        x_prev = x              #xi-1
        x = exp(x_prev, 2, n)   #xi = (xi-1)^2 mod n
        if x == 1 and x_prev != 1 and x_prev != (n-1):      #NSR test
            return True

    if x != 1:                  #Fermat test
        return True
    return False

I have been struggling because of this for few hours, and still cannot find which part of the code is the problem. Please let me know if you have any advice. P.S. exp (a,b,c) returns a^b mod c.

Answer 1

Here is a different Miller Rabin Implementation that works well. The MillerRabin function may be what you're most interested in:

def llinear_diophantinex(a, b, divmodx=1, x=1, y=0, withstats=False, pow_mod_p2=False): 
  origa, origb = a, b 
  r=a  
  q = a//b 
  prevq=1  
  if withstats == True: 
    print(f"a = {a}, b = {b}, q = {q}, r = {r}")   
  while r != 0:  
       prevr = r  
       a,r,b = b, b, r   
       q,r = divmod(a,b) 
       x, y = y, x - q * y 
       if withstats == True: 
         print(f"a = {a}, b = {b}, q = {q}, r = {r}, x = {x}, y = {y}")  
  y = 1 - origb*x // origa - 1 
  x,y=y,x 
  modx = (-abs(x)*divmodx)%origb 
  if withstats == True: 
    print(f"x = {x}, y = {y}, modx = {modx}") 
  if pow_mod_p2==False:   
    return x, y, modx 
  else: 
    if x < 0: return modx 
    else: return x 

def ltrailing(N):
    return len(str(bin(N))) - len(str(bin(N)).rstrip('0'))

def MillerRabin(N, primetest, iterx, powx, withstats=False): 
  primetest = pow(primetest, powx, N) 
  if withstats == True:
     print("first: ",primetest) 
  if primetest == 1 or primetest == N - 1: 
    return True 
  else: 
    for x in range(0, iterx-1): 
       primetest = pow(primetest, 2, N) 
       if withstats == True:
          print("else: ", primetest) 
       if primetest == N - 1: return True 
       if primetest == 1: return False 
  return False 

def sfactorint_isprime(N, withstats=False):
    if N == 2:
      return True
    if N % 2 == 0:
      return False
    if N < 2:
        return False
    iterx = ltrailing(N - 1)
    """ This k test is an algorithmic test builder instead of using
        random numbers. The offset of k, from -2 to +2 produces pow tests
        that fail or pass instead of having to use random numbers and more
        iterations. All you need are those 5 numbers from k to get a 
        primality answer. This is the same as doing:
        pow(N, (1<<N.bit_length()) - 1, 1<<N.bit_length()) but much faster
        using a linear diophantine formula which gives the same result for 
        powers of 2
    """
    k = llinear_diophantinex(N, 1<<N.bit_length(), pow_mod_p2=True) - 1
    t = N >> iterx
    tests = [k-2, k-1, k, k+1, k+2]
    for primetest in tests:
        if primetest >= N:
            primetest %= N
        if primetest >= 2:
            if MillerRabin(N, primetest, iterx, t, withstats) == False:
                return False
    return True