1

我编写了一个程序，通过使用 ((2^n)-1)(2^(n-1)) 计算出 1-1000 的所有梅森素数的偶数，其中 n 是梅森素数。

这是程序：

def PrimeFinder(PotPrime):
    PlaceNum=1
    for x in range (int(PotPrime**0.5)):
        PlaceNum=PlaceNum+1
        if int(PotPrime/PlaceNum) == (PotPrime/PlaceNum):
            return False
    return True


TrialNum = 1

for x in range (1000):

    if PrimeFinder(TrialNum) == True:

        if PrimeFinder((2**TrialNum)-1) == True:
            print(TrialNum,"is the Mersenne Prime for the perfect number:",(2**(TrialNum-1))*((2**TrialNum)-1))

    TrialNum = TrialNum+1

该程序运行良好，直到 32 < TrialNum < 60，因为它正确地识别出 31 是梅森素数，但它不适用于 61（以及所有大于它的数字）。

我想知道 Python 是否根本无法进行那么大的计算，或者我对梅森素数的理解或编程是否存在缺陷。

4

4 回答 4

2

Python 的整数对它们的大小没有限制，所以你应该安排你的计算，使它们都是基于整数的。以下是可以在您的程序中进行的一些更改，以使用整数而不是浮点数。

代替

for x in range (int(PotPrime**0.5)):

使用类似的东西

while x*x < PotPrime:

代替

if int(PotPrime/PlaceNum) == (PotPrime/PlaceNum):

使用更简单的

if PotPrime % PlaceNum == 0:
于 2015-03-26T10:36:36.610 回答
2

我想四舍五入的错误：如果你调试你注意到它认为 2 是 2^61 -1 的除数（这没有意义）。如果你用它替换if int(PotPrime/PlaceNum) == (PotPrime/PlaceNum):它 if PotPrime % PlaceNum == 0:是固定的。但是你的算法效率很低，而且 2^61 - 1 是一个非常大的数字，所以预计它需要几个小时。（可能更长）

于 2015-03-24T01:19:19.663 回答
2

使用Lucas-Lehmer 素数检验来检查哪些素数会产生梅森素数 (Mp)，这样您就可以避免检查梅森数本身的素数，因为它们非常大，而且因为只有素数指数会产生 Mp，通过这个测试你只需要那个 p 是素数并通过这个测试，那么你可以建立你的完美数为 ((2^p)-1)(2^(p-1))。

通过这个测试，我可以找到产生 Mp en 12seg in我的机器有点旧，奔腾双核 3GHz

这是代码（在 python 3.3 中），您还可以使用更强大的素性测试，例如Miller-Rabin 测试或 Baillie-PSW，这些测试与试验师伤口不同，不会永远检查大量数字。

def primality_Test_Trial_Division(n:int) -> bool:
    if n >= 0 :
        if n<2:
            return False
        elif n<4:
            return True
        elif not n&1 or n%3==0:
            return False
        else:
            mid = 1 + int( n**0.5 )
            i = 5
            while i<mid and n%i:
                i += 2
            return i>=mid
    else:
        raise ValueError

isPrime = primality_Test_Trial_Division

def primality_Test_LL(p:int) -> bool:
    """Lucas–Lehmer primality test. Check if Mp = 2^p − 1 is prime.
       en.wikipedia.org/wiki/Lucas%E2%80%93Lehmer_primality_test"""
    if isPrime(p):
        if p==2:
            return True
        mersenne = (2**p)-1 #Mp
        s = 4
        for x in range( p-2 ):
            s = pow(s,2,mersenne)-2
            #Performing the mod Mp at each iteration ensures
            #that all intermediate results are at most p bits
            #(otherwise the number of bits would double each iteration).
            #The same strategy is used in modular exponentiation.
        return s==0
    else:
        return False

import itertools
def perfect_numbers(n:int):
    """print the first n perfect numbers""" 
    perfect = 0
    for p in itertools.count(2):
        if primality_Test_LL(p):
            print(p,"is the Mersenne Prime for the perfect number:",(2**(p-1))*((2**p)-1))   
            perfect += 1
            if perfect >= n:
                break
于 2015-12-06T01:13:47.363 回答
1

完美数具有以下形式，它们可以被相同数量的 2 整除，奇数因子的 bit_length 是素数。像 120 和 2016 这样的半完美数具有相同的形式，只是奇数因子不是素数。例如，496 的因数是：[2, 2, 2, 2, 31]，2k-1 表示 2 * (2 * 2 * 2 * 2)-1 应该等于 31。确实如此，所以 496 只是需要在 31 上进行最终的 lucas lehmer 测试以测试最终测试。

很容易将偶数分解为用于质数测试的奇数和用于对完美数进行 bit_length 测试的偶数，这就是下面显示的第二个示例。这很容易做到。

例如，496 具有因子和以下属性：


import gmpy2

# If you don't have p2ecm, you can use any 
# factorization method, but p2ecm is not needed for the final program below 
# the example, we use it here just to show the factors.

def isoneminuspowersoftwo(N):
   N=N+1
   return (N & (N-1) == 0) and N != 0

In [2]: p2ecm(496)    # from https://github.com/oppressionslayer/primalitytest                                                           
Out[2]: [2, 2, 2, 2, 31]

In [3]: 2*2*2*2                                                                 
Out[3]: 16

In [5]: isoneminuspowersoftwo(31)                                               
Out[5]: True

In [7]: (2*2*2*2).bit_length() == (31).bit_length()                             
Out[7]: True

#and

In [8]: 2*(2*2*2*2) -1     # 2k-1                                                      
Out[8]: 31

# Perform Lucas Lehmer test on 31 above:
In [106]: s=4 
     ...: for x in range((31).bit_length()-1): 
     ...:    s = (s * s - 2) % 31 
     ...:    if s in [0,2]: print(True) 
     ...:                                                                       
True

# All three tests above passed, so the number 496 is a perfect number.

第二个例子

获得没有因子的因子的一种简单方法是使用简单的 python 程序来提取数字，您可以使用以下方法进行测试：

def ffs(x):
    return (x&-x).bit_length()-1

def extractoddfactor(N):
   return N//(2**ffs(N))
   
In [13]: extractoddfactor(496)                                                      
Out[13]: 31

In [14]: 496//31                                                                
Out[14]: 16

下面是 (2 ** 5)*(2 ** 6-1) 的半完美示例。尽管它通过了前两个测试，但它在素数测试中失败了，因此它是半完美的而不是完美的数字。

In [11]: (2**5)*(2**6-1)                                                        
Out[11]: 2016

In [21]: extractoddfactor(2016)                                                 
Out[21]: 63

In [23]: 2016//63                                                               
Out[23]: 32

In [24]: (32).bit_length() == (63).bit_length()                                 
Out[24]: True

In [25]: 2*(32) -1 == 63  # 2k-1 test                                                      
Out[25]: True

In [107]: s=4 
     ...: for x in range((63).bit_length()-1): 
     ...:    s = (s * s - 2) % 63 
     ...:    if s in [0,2]: print(True) 
     ...:  
# No output, so False

我编写了一个程序，通过上述测试和上述 LucasLehmer 素数测试找到完美数。如果您没有 gmpy2，只需删除 gmpy2.mpz 包装器并将 gmpy2 gmpy2.bit_length(x) 语句更改为 python 等价物，如下所示：x.bit_length()。我使用了 gmpy2，因为它比没有它快十倍，但它不是必需的，并且可以很容易地修改为不使用它。

这个程序执行上面的测试来测试一个数字是否完美。它测试上面列出的完美 num 的特征，然后以 lucas lehmer 测试结束。

如果删除 gmpy2 mpz 包装器并将 bit_length() 语句更改为 python 格式，则不需要外部库

import gmpy2

def ffs(x):
    x=gmpy2.mpz(x)
    return gmpy2.bit_length(x&-x)-1

def levelx(N, withstats=False):
  if N <= 1: return False
  N = gmpy2.mpz(N)
  zero, one, two = gmpy2.mpz(0), gmpy2.mpz(1), gmpy2.mpz(2)
  temp = gmpy2.mpz(N)
  newlevel = gmpy2.bit_length(temp)
  primetest = gmpy2.mpz(temp//(two**(ffs(temp))))
  offset = gmpy2.mpz(temp//primetest)
  s = gmpy2.mpz(4)
  
  nextlevel = newlevel // two
  
  check = temp // (two**nextlevel)

  prevtemp = two**nextlevel

  if withstats == True:
    print (newlevel, newlevel, temp)
    print (newlevel, nextlevel+one, prevtemp)
 
        
  if (prevtemp & (prevtemp-one) == zero) and prevtemp != zero:
       if gmpy2.bit_length(offset) == gmpy2.bit_length(primetest):
          if ((primetest+one) & ((primetest+one)-one) == zero):
             for x in range(gmpy2.bit_length(primetest)-1):
               s = (s * s - 2) % primetest
               if withstats == True:
                  print(s)
               if s in [0,2]: return True
             return False
          else: return False
       else: return False
  else: return False

     
In [69]: levelx((2**4)*(2**5-1))                                                
Out[69]: True

In [101]: for x in range(1, 10000):  
...:     if levelx((2**x)*(2**(x+1)-1)) == True:  
...:         print(x+1, (2**x)*(2**(x+1)-1))  
...:                                                                       
2 6
3 28
5 496
7 8128
13 33550336
17 8589869056
19 137438691328
31 2305843008139952128
61 2658455991569831744654692615953842176
89 191561942608236107294793378084303638130997321548169216
107 13164036458569648337239753460458722910223472318386943117783728128
127 14474011154664524427946373126085988481573677491474835889066354349131199152128
...
于 2021-01-18T00:55:25.400 回答