Perfect Numbers and Mersenne Primes - Python

Question

I've written a program which works out the even perfect numbers for all Mersenne Primes from 1-1000, by using ((2^n)-1)(2^(n-1)) where n is a Mersenne Prime number.

This is the program:

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

This program works fine, up until somewhere where 32 < TrialNum < 60, as it correctly identifies that 31 is a Mersenne Prime, however it does not for 61 (and all numbers greater than it).

I was wondering if Python simply cannot do calculations that large, or whether there is a flaw in either my understanding Mersenne Primes or programming.

Answer 1

Rounding errors I suppose: if you debug you notice it thinks 2 is a divisor of 2^61 -1 (which makes no sense). If you replace if int(PotPrime/PlaceNum) == (PotPrime/PlaceNum): with if PotPrime % PlaceNum == 0: it's fixed. But your algorithm is quite inefficient and 2^61 - 1 is a very large number so expect it to take hours. (probably even longer)

Answer 2

Python's integers do not have limits on their size, so you should arrange your calculations so they're all integer-based. Here are a few changes that could be made in your program to use integers instead of floating-point.

Instead of

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

use something like

while x*x < PotPrime:

Instead of

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

use the simpler

if PotPrime % PlaceNum == 0:

Answer 3

Use the Lucas-Lehmer primality test to check what primes produce Mersenne primes (Mp), that way you avoid check the primality of the Mersenne number itself because they are pretty big, and because only prime exponent produce Mp, with this test you only need that that p to be prime and pass this test, then you can build your perfect number as ((2^p)-1)(2^(p-1)).

With this test I can find the first 16 p [2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203] that produce a Mp en 12seg in my machine that is kind of old, a Pentium Dual-Core of 3GHz

Here is the code (in python 3.3) you can also use a more powerful primality Test like the Miller-Rabin test or the Baillie-PSW those unlike trial division wound't take forever to check large numbers.

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

Answer 4

Perfect Numbers have the following form, they are divisible by the same number of 2's as the bit_length of the odd factor, which is prime. semi perfect numbers like 120 and 2016, have this same form, except that the odd factor is not prime. For example, the factors of 496 are: [2, 2, 2, 2, 31], and 2k-1 says that 2 * (2 * 2 * 2 * 2)-1 should equal 31. It does, so 496 just needs a final lucas lehmer test on 31 to test for a final test.

It is very easy to break up an even number to an odd number that you use for prime testing and an even number that you use for bit_length testing on perfect numbers, and that is the second example shown below. it's very easy to do.

For example, 496 has the factors and following properties:

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.

2nd Example

A Simple way to get the factors without factors is using a simple python program to extract the numbers, which you can than test with:

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

A semiperfect example is below for (2 ** 5)*(2 ** 6-1). It fails on the prime test so is semiperfect and not a perfect number, althought it passes the first two tests.

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

I wrote a program that finds perfect numbers with the above tests and the above LucasLehmer prime test. If you don't have gmpy2, just remove the gmpy2.mpz wrappers and change the gmpy2 gmpy2.bit_length(x) statements to pythons equivelant like so: x.bit_length(). I used gmpy2 as it's ten time faster than without it, but it's not needed, and can be easily modified to not use it.

This program performs the tests above to test if a number is perfect or not. It tests for the features of a perfect num, listed above, and then finishes with a lucas lehmer test.

No external libraries are needed if you remove the gmpy2 mpz wrappers and change the bit_length() statements to python format

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
...