在 Python 上实现 Miller-Rabin 算法

Question

我正在尝试在 Python 上实现 Miller-Rabin 算法。
我已经按照教科书的伪代码进行了编码，但由于某种原因，它没有按我的预期工作。
具体来说，函数“test”在进行 Fermat 测试时有时会返回“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

我为此苦苦挣扎了几个小时，但仍然找不到代码的哪一部分是问题所在。如果您有任何建议，请告诉我。 附言exp (a,b,c) 返回 a^b mod c。

Answer 1

这是另一个运行良好的 Miller Rabin 实现。 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