Mersenne数的素性检验

Question

请协助解决此问题；

对于给定的指数为p的Mersenne数，如果Lucas-Lehmer级数在p- 2的位置为0，则该数是素数。编写一个函数来检验指数为p的Mersenne数是否为素数。检验素数p在3和65之间的Mersenne数是否为素数。你的最终答案应该是一个元组列表，它由(Mersenne指数，0)或(1)组成，对于你测试的每个Mersenne数，其中0和1分别是'False‘和'True’的替代。下面显示的是我尝试的解决方案

#function to define a mersenne number
def mersenne_number(p):
    return 2**p - 1

#function to generate the Lucas-Lehmer sequence
def lucas_lehmer(p):
    ll_seq = [4]
    if p > 2:
        for i in range(1, (p - 2) + 1):
            n_i = (ll_seq[i-1] ** 2 - 2)%(2**p - 1)
            ll_seq.append(n_i)
    return ll_seq

#To generate a list of Tuples consisting of (mersenne exponents, 1 or 0), 1 and 0 represent True and False respectively.
mersenne_primes = []

def ll_prime(num1, num2):
    for p in range(3, 65):
        if lucas_lehmer(p)[-1] == 0:
            mersenne_primes.append((p, 1))
        else:
            mersenne_primes.append((p, 0))

    print(mersenne_primes)

在这个阶段，我得到了一个错误的输出。函数返回素数mersenne指数，而不是素数，意思是0。

Answer 1

我从维基百科页面上的伪代码中实现了Lucas Lehmer测试：

def LucasLehmer(p):
   if p == 2:
     return True
   s = 4
   M = pow(2, p) - 1
   for x in range (1, (p-2)+1):
      s = ((s * s) - 2) % M
   if s == 0: return True
   else: return False

for x in range(2,10000): 
   if LucasLehmer(x): 
      print(f"2**{x}-1 is Prime") 
                                                                                                                                                           
2**2-1 is Prime
2**3-1 is Prime
2**5-1 is Prime
2**7-1 is Prime
2**13-1 is Prime
2**17-1 is Prime
2**19-1 is Prime
2**31-1 is Prime
2**61-1 is Prime
2**89-1 is Prime
2**107-1 is Prime
2**127-1 is Prime
2**521-1 is Prime
2**607-1 is Prime
2**1279-1 is Prime
2**2203-1 is Prime
2**2281-1 is Prime

我还改变了算法，使其成为因子分解引擎，如果你有兴趣在这里使用它进行因子分解：

import math

def PrimeFinderLucasLehmer4(N):
   p = 1<<N.bit_length()-1
   if N == 2:
     return 2
   if N == 3:
     return 3
   s = 4
   M = pow(p, 2) - 1
   for x in range (1, 100000):
     s = (((s * N ) - 2 )) % M
     xx = [math.gcd(s, N)] + [math.gcd(s*p+x,N) for x in range(7)] + [math.gcd(s*p-x,N) for x in range(1,7)] 
     try:
        prime = min(list(filter(lambda x: x not in set([1]),xx)))
     except:
        prime = 1
     if prime == 1:
        continue
     else:
        break
   #print (s, x, prime, xx)
   return prime


In [70]: PrimeFinderLucasLehmer4(1009732533765211)                                                                                                                                        
Out[70]: 11344301

And from https://stackoverflow.com/questions/4078902/cracking-short-rsa-keys

In [71]: PrimeFinderLucasLehmer4(10142789312725007)                                                                                                                                       
Out[71]: 100711423