python - 在 Python 中实现 Miller-Rabin 复合性时遇到问题

Question

我不确定这是否是发布此问题的正确位置，所以如果不让我知道！我正在尝试在 python 中实现 Miller Rabin 测试。测试是找到第一个证明 N 的合数，一个奇数。我的代码适用于长度稍小的数字，但当我输入一个大数字时停止工作。(The "challenge" wants to find the witness of N := 14779897919793955962530084256322859998604150108176966387469447864639173396414229372284183833167 in which my code returns that it is prime when it isn't) The first part of the test is to convert N into the form 2^k + q, where q 是一个素数。python是否有一些限制，不允许为此使用大量数字？这是我的那部分测试的代码。

def convertN(n): #this turns n into 2^x * q
placeholder = False
list = []
#this will be x in the equation
count = 1
while placeholder == False:
    #x = result of division of 2^count
    x = (n / (2**count))
    #y tells if we can divide by 2 again or not
    y = x%2
    #if y != 0, it means that we cannot divide by 2, loop exits
    if y != 0:
        placeholder = True
        list.append(count) #x
        list.append(x)     #q
    else:
        count += 1
#makes list to return
#print(list)
return list

实际测试的代码：

def test(N):
#if even return false
if N == 2 | N%2 == 0:
    return "even"
#convert number to 2^k+q and put into said variables
n = N - 1
nArray = convertN(n)
k = nArray[0]
q = int(nArray[1])
#this is the upper limit a witness can be 
limit = int(math.floor(2 * (math.log(N))**2))

#Checks when 2^q*k = 1 mod N
for a in range(2,limit):
    modu = pow(a,q,N)
    for i in range(k):
        print(a,i,modu)
        if i==0:
            if modu == 1:
                break
        elif modu == -1:
            break
        elif i != 0:
            if modu == 1:
                #print(i)
                return a
        #instead of recalculating 2^q*k+1, can square old result and modN that.
        modu = pow(modu,2,N)

任何反馈表示赞赏！

标签： pythonmathcryptography

Answer 1

梅，我用 python 写了一个米勒拉宾测试，米勒拉宾部分是线程的，所以它非常快，比 sympy 更快，对于更大的数字：

import math

def strailing(N):
   return N>>lars_last_powers_of_two_trailing(N)

def lars_last_powers_of_two_trailing(N):
  """ This utilizes a bit trick to find the trailing zeros in a number
      Finding the trailing number of zeros is simply a lookup for most
      numbers and only in the case of 1 do you have to shift to find the
      number of zeros, so there is no need to bit shift in 7 of 8 cases.
      In those 7 cases, it's simply a lookup to find the amount of zeros.
  """
  p,y=1,2
  orign = N
  N = N&15
  if N == 1: 
     if ((orign -1) & (orign -2)) == 0: return orign.bit_length()-1
     while orign&y == 0:
       p+=1
       y<<=1
     return p
  if N in [3, 7, 11, 15]: return 1
  if N in [5, 13]: return 2
  if N == 9: return 3
  return 0
  
def primes_sieve2(limit):
    a = [True] * limit
    a[0] = a[1] = False

    for (i, isprime) in enumerate(a):
        if isprime:
            yield i
            for n in range(i*i, limit, i):
                a[n] = False
                
def llinear_diophantinex(a, b, divmodx=1, x=1, y=0, offset=0, withstats=False, pow_mod_p2=False): 
  """ For the case we use here, using a 
      llinear_diophantinex(num, 1<<num.bit_length()) returns the 
      same result as a 
      pow(num, 1<<num.bit_length()-1, 1<<num.bit_length()). This 
      is 100 to 1000x times faster so we use this instead of a pow. 
      The extra code is worth it for the time savings.
  """
  origa, origb = a, b 
  r=a  
  q = a//b 
  prevq=1  
  #k = powp2x(a)
  if a == 1:
    return 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
  if withstats == True: 
    print(f"x = {x}, y = {y}")  
  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*divmodx)%origb, y, modx, (origa)%origb
  else: 
    if x < 0: return (modx*divmodx)%origb 
    else: return (x*divmodx)%origb

def MillerRabin(arglist): 
  """ This is a standard MillerRabin Test, but refactored so it can be
      used with multi threading, so you can run a pool of MillerRabin 
      tests at the same time.
  """
  N = arglist[0]
  primetest = arglist[1]
  iterx = arglist[2]
  powx = arglist[3]
  withstats = arglist[4]
  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 

# For trial division, we setup this global variable to hold primes
# up to 1,000,000
SFACTORINT_PRIMES=list(primes_sieve2(100000))
  
  
# Uses MillerRabin in a unique algorithimically deterministic way and
# also uses multithreading so all MillerRabin Tests are performed at 
# the same time, speeding up the isprime test by a factor of 5 or more. 
# More k tests can be performed than 5, but in my testing i've found 
# that's all you need.
def sfactorint_isprime(N, kn=5, trialdivision=True, withstats=False):

    from multiprocessing import Pool

    if N == 2:
      return True
    if N % 2 == 0:
      return False
    if N < 2:
        return False
        
    # Trial Division Factoring
    if trialdivision == True:
      for xx in SFACTORINT_PRIMES:
        if N%xx == 0 and N != xx:
          return False
    
    iterx = lars_last_powers_of_two_trailing(N)
    """ This k test is a deterministic 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. I've tested this against all numbers in 
        https://oeis.org/A001262/b001262.txt and all fail, plus other
        exhaustive testing comparing to other isprimes to confirm it's 
        accuracy.
    """
    k = llinear_diophantinex(N, 1<<N.bit_length(), pow_mod_p2=True) - 1
    t = N >> iterx
    tests = []
    
    if kn % 2 == 0: offset = 0
    else: offset = 1
    
    for ktest in range(-(kn//2), (kn//2)+offset):  
       tests.append(k+ktest)
    
    for primetest in range(len(tests)):
      if tests[primetest] >= N:
         tests[primetest] %= N
  
    arglist = []
    for primetest in range(len(tests)):
      if tests[primetest] >= 2:
        arglist.append([N, tests[primetest], iterx, t, withstats])
     
    with Pool(kn) as p:
       s=p.map(MillerRabin, arglist)    
    
    if s.count(True) == len(arglist): return True
    else: return False

sinn=14779897919793955962530084256322859998604150108176966387469447864639173396414229372284183833167

print(sfactorint_isprime(sinn))