Реализация алгоритма Миллера-Рабина на Python

Question

Я пытаюсь реализовать алгоритм Миллера-Рабина на Python.
Я закодировал так, как говорит псевдокод моего учебника, но по какой-то причине он не работает так, как я ожидал.
Чтобы быть особенным, функция test иногда возвращает true, когда переходит к проверке Ферма.

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

Я боролся из-за этого в течение нескольких часов и до сих пор не могу найти, какая часть кода является проблемой. Пожалуйста, дайте мне знать, если у вас есть какие-либо советы. PS exp (a,b,c) возвращает a^b mod c.

Answer 1

Вот другая реализация Миллера Рабина, которая хорошо работает. Функция MillerRabin может быть тем, что вас больше всего интересует:

      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