在Python中获取Miller-Rabin素性测试有困难

Question

我一直在尝试用Python实现Miller-Rabin素性测试。不幸的是，在一些质数上似乎存在问题。任何帮助都是非常感谢的。

代码：

def _isPrime(n):

if n % 2 == 0:
    return False

a = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
s, d = _Apart(n) # turns it into 2^s x d
print(s, d)

basePrime = False
isPrime = True

for base in a:
    if base >= n-2:
        break 

    if isPrime == False and basePrime == False:
        return False
    
    witness = pow(base, d)

    if (witness % n) == 1 or (witness % n) == (n-1):
        basePrime = True
        continue
    else:
        basePrime = False

    for r in range(1, s):
        witness = pow(base, pow(2, r) * d)

        if (witness % n) != ( n - 1 ):
            isPrime = False

return True

测试：

isPrime(17)

期望值：

True

结果：

False

Answer 1

我写了一个米勒·拉宾测试，它是确定性的，不需要随机数。这个实现是针对python 3.7的。在Python3.8中，llinear_diophantinex可以替换为pow(x，-1，y)。我还使用了gmpy2，因为它非常快，但是如果您不能使用它，您可以用普通的能力替换gmpy2语句，并删除gmpy2.mpz()包装器，因为这些包装器只是用来重载操作符。

import gmpy2

sinn = 2110229697309202254897383305762150945330987087513434511395506048950594976569434432057019507105035289374307720719984431280856161609820548842778454256113246763860786119268583367543952735347969627478873317341364209555365064365565504232770227619462128918701942169785585423104678142850200975026619010035331023744330713985615650556129731348659986462960062760308034462660525448390420668021248422741300646552941285862310410598374242189448623917196191138254637812716211329113836605859918549332304189053950819346551095911511755911832183789503704294770046935064469435830299623205136625543859303686699678929069468518950480476841246805908501510754550017255944080874819287974625925494008373883250410775902993163965873632474224574883242826458163446781002284368017611606202344050570737818087202137703099075773680753707346415849787963446390136517016131227807076254668461445862154978026041507116570585784569893773262639243954090283224759975513502582494002154146757110676408972377044584495342170277522887809749465855954126593100747444378301829661568735873345178089061677917127496915956539418931430313218084338374827152407795095072639044306222222695685778907958272820576498682506540189586657786292950574081739269257159839589987847266550007783514316481286222515710538845836151864127815058116482680058626451349913138908040817800742009650450811565324184631847563730941344941348929727603343965091116543702880556850922077216848669966268219928808236163268726995495688157209747596437162960244538054993785127947211290438554095851924381172697827312534174244295581184309147813790451951453564726742200569263225639113681905176376701339808868274637448606821696026703034737428319530072483125495383057919894902076566679023694181381398377144302767983385253577700652358431959604517728821603076762965129019244904679015099154368058005173028200266632883632953133017055122970338782493475762548347258351148037427739052271661340801912188203749647918379812483260399614599813650518046331670764766419886619324840045611486524123102046413946014624119568013100078163986683199814025915420877588778260860713148420321896163326473203441644820182490479899368048072263481024886708136521847014624735722333931331098969321911443978386868675912141648200500219168920887757573018380579532261231821382787339600631297820996466930957801607217549420247654458172818940238337170577825003408756362106088558651381993611741503374243481167926898332728164900189941804942580426055589622673679047058619682175301326905577843405270203660160407401675700528981573327582844330828861745574031416926871562443652858767649050943181353635950301154441954046214987718582670685455252774874198771086552440702483933126644594300464549471422237478151976561680719370424626162642534252062987911763456822609569209140676822858933588602318066530038691463577379331113471591913447226829868760176810195567325921301390329055242213842898142597360121925124635965685365925901913816717677946911762631634793638450106377437599347740569467683272089859392249351406815344105961234868327316964137925419770514177021722214309784062017826024217906664090209434553785436385765927274067126192143337589109608949427467825999057058702263715338956534536892852849984934736685814891286495169007648767081688963426768409476169071460997622740467533572971356017575900999100928776382541052696124463195981888715845688808970103527288822088031150716134784735332326775370417950625124642515148342694377095213470544739900830244879573205335578256682901821773047071352497997708791157012233232529777513203024818391621220967964874173106990772425289900446640237659116713251437567138729645677868024033209183367071421651937808005637679844370347367922676824239404492688418047080583797577102267329067247758368597488680401670673861120323439239792549053895366970423259196919428554146265587250617656401028722578111927104663315250291888502226235291264834968061065817079511872899991276288365723969841290984981957389126603952133124328219936785870274843554107325931034103072894378438818494802517594594270034007832922248742746517915210656205746338575621725899098414488628833412591266637224507533934158213117522503993423240638893845121918647788013


def llinear_diophantinex(a, b, divmodx=1, x=1, y=0, offset=0, withstats=False, pow_mod_p2=False):  
  origa, origb = a, gmpy2.mpz(b)  
  r=a   
  q = a//b  
  prevq=1   
  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 = gmpy2.mpz(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 ffs(x): 
    """Returns the index, counting from 0, of the 
    least significant set bit in `x`. 
    """ 
    return (x&-x).bit_length()-1 
 
def MillerRabin(arglist):  
  N = arglist[0] 
  primetest = arglist[1] 
  iterx = arglist[2] 
  powx = arglist[3] 
  withstats = arglist[4] 
  primetest = gmpy2.powmod(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):  
       primetest = gmpy2.powmod(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): 
 
    N = gmpy2.mpz(N) 
    from multiprocessing import Pool 
 
    if N <= 1: return False 
    if N == 2: 
      return True 
    if N % 2 == 0: 
      return False 
    if N < 2: 
        return False 
         
    # Add Trial Factoring here to speed up smaller factored number testing 
 
     
    iterx = ffs(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.  
    """ 
    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 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(5) as p: 
       s=p.map(MillerRabin, arglist)     
     
    if s.count(True) == len(arglist): return True
    else: return False 
     
    return s   

Answer 2

最近，我还尝试以确定性的形式实现此算法，如Miller Test中所示，但遇到了同样的问题。我不明白为什么它在一些数字上失败了，所以我决定实现“完整”的Miller-Rabin Test，并且它起作用了。

我应该警告一下，随着n的增加，它会变得很慢。建议使用像Numpy这样的数值库来优化计算。

通过对代码进行一些调整，可以实现以下目标：

import random

def _isPrime(n):
    if n % 2 == 0:
        return False

    rounds = 40  # Determines the accuracy of the test
    s, d = _Apart(n)  # Turns it into 2^s x d + 1

    for _ in range(rounds):
        base = random.randrange(2, n-1)  # Randomly selects a base
        witness = pow(base, d)

        if (witness % n) == 1 or (witness % n) == (n-1):
            continue

        for _ in range(1, s):
            witness = pow(witness, 2) % n

            if witness == (n-1):
                break
        else:  # if the for-loop is not 'broken', n is not prime
            return False

    return True

Answer 3

def is_prime(n):
    """returns True if n is prime else False"""
    if n < 5 or n & 1 == 0 or n % 3 == 0:
        return 2 <= n <= 3
    s = ((n - 1) & (1 - n)).bit_length() - 1
    d = n >> s
    for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
        p = pow(a, d, n)
        if p == 1 or p == n - 1 or a % n == 0:
            continue
        for _ in range(s):
            p = (p * p) % n
            if p == n - 1:
                break
        else:
            return False
    return True