在Python上实现Miller-Rabin算法

如何解决在Python上实现Miller-Rabin算法

我正在尝试在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。

解决方法

这是另一种效果很好的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,r   
       q,r = divmod(a,b) 
       x,y = y,x - q * y 
       if withstats == True: 
         print(f"a = {a},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,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,t,withstats) == False:
                return False
    return True