Perfect Numbers
In my quest to undertake various programming projects, for the sheer fun of it, I came across a definition
of perfect numbers:
A perfect number is a positive integer which is the sum of all its divisors excluding itself.
For example, 6 is a perfect number as 6 = 1 + 2 + 3, and 1, 2, and 3 are all the whole number divisors of
6, except 6 itself. (Divisor - a number which divides into another number a whole number of times).
Alternatively, we can define a perfect number as a number which is equal to half of the sum of all its
divisors (including itself) e.g. 6 = (1 + 2 + 3 + 6)/2.

This little study illustrates the limitations of computer power and ways to overcome it. I used Java
(NetBeans) for this project.

The Brute Force Approach

This seemed like a good way to test the power of a computer - by getting it to calculate perfect numbers.
The most obvious way to do this (superficially) is to take each positive integer in turn and find all of its
divisors and then check to see if the sum of these divisors is equal to that integer. Hence I wrote the
following code:

package perfectnumbers;

import java.util.*;
import java.math.BigDecimal;
import java.math.BigInteger;

// @BotRejectsInc.

public class Main
{
public static void main(String[] args)
{    
  findPerfect(10000);
}

//Finds all the perfect numbers <= n
private static ArrayList<Integer> findPerfect(int n)
{
   ArrayList<Integer> perfects = new ArrayList<Integer>();
   double remainder;
   int accumulator;
   ArrayList<Integer> perfect = new ArrayList<Integer>();
   int maxDivisor, breakCount = 0;
   
   //All perfect numbers are even
   for(int i = 2; i <= n; i += 2)
   {
       maxDivisor = (int)Math.floor(i/2);
       accumulator = 0;
       //Find divisors
       accumulator = 1 + 2; // we know that 1 and 2 are divisors
       for(int j = 3; j <= maxDivisor; j++)
       {
           remainder = i % j;
           if(remainder == 0)
           {
               //Divisor found
               accumulator += j;
               if(accumulator > i)
               {
                   //Sum of divisors is greater than the number so not a perfect number
                   break;
               }
           }
       }
       if(accumulator != i)
       {
           //Sum of divisors less than number so number not perfect
           continue;
       }
       if(accumulator == i)
       {
           System.out.println(i);
           perfect.add(i);
       }
   }
   return perfects;
}
}
Note the following short-cuts:

  • All perfect numbers are even (well it is not known if any are odd so we shall assume none are), so 1
    and 2 are divisors.
  • When testing the divisors there is no point testing a number more than half the number being tested
    for perfectness!

All very obvious so far, and the output is:

6
28
496
8128
BUILD SUCCESSFUL (total time: 4 seconds)

Compare this to a list of the first 8 known perfect numbers:

6
28
496
8128
33550336
8589869056
137438691328
2305843008139952128

Our first attempt has been quite successful, we have obtained the first 4 which were all those known to the
ancient Greeks. However we hit a problem - we cannot find number greater than about 1 billion (more
strictly 2^32 = 4 294 967 296, minus one if we include zero) since an integer is 32 bit and cannot represent
any number larger than about 4 thousand million. Besides, our method has become too slow, using:

    findPerfect(1000000000);

takes a long time to compute!

We need a short-cut and one is available in the form of Mersenne Primes.

Mersenne Primes - The Brute Force Approach Mk2

All (even) perfect numbers can be obtained from a mathematical series called the series of Mersenne
Prime numbers. Mersenne Primes are prime number satisfying:
So, find the first 100 Mersenne Primes, say, and we can then easily obtain the first 100 perfect numbers,
however, we will need to calculate lots of prime numbers (index p) as not all prime numbers, p, give a
Mersenne Prime, M. In fact very few of them do.

Finding the prime numbers

Recall that:

A prime number is a number with no other (whole number) divisors other than 1 and itself.

E.g. 2 is prime, as it can only be divided (without remainder) by 1 and 2;
11 is prime as it can only be divided by 1 and 11,
10 is not prime, as it can be divided by 2 and 5 as well as 1 and 10.

I wrote the following couple of functions to find the prime numbers, p, to use in obtaining the Mersenne
Primes:

//Finds and returns the first n prime numbers
private static ArrayList<Integer> findPrime(int n)
{
    ArrayList<Integer> primes = new ArrayList<Integer>();
    //Count the number of primes found
    int count;
    //2 is the only even prime
    primes.add(2);
    count = 1;
    //We will continue searching odd numbers from 3 onwards
    int number = 3;
    
    while(true)
    {
        //Consider odd numbers only
        for(int i = 3; i <= (int)Math.floor(number/2); i++)
        {
            if(number % i == 0)
            {
                //Found a divisor other than 1 and number
                //ergo number is not prime
                number += 2;
                continue;
            }
        }
        //No divisors found so number is prime
        primes.add(number);
        count++;
        if(count >= n)
        {
            break;
        }
        number += 2;
    }
    return primes;
}

//Determines whether or not a number of arbitrary size is prime using 'brute force'
private static boolean isPrime(BigInteger number)
{
    BigInteger ZERO = new BigInteger("0");
    BigInteger ONE = new BigInteger("1");
    BigInteger TWO = new BigInteger("2");
    BigInteger max = number.divide(TWO);
    BigInteger divisor = new BigInteger("3");
    
    if(number.compareTo(TWO) == 0)
    {
        return true;
    }
    
    if(number.mod(TWO).compareTo(ZERO) == 0)
    {
        //number is even so can not be prime
        //System.out.println("Number is even");
        return false;
    }
    while(true)
    {
        if(number.mod(divisor).compareTo(ZERO) == 0)
        {
            //divisor found, so not a prime
            System.out.println("Divisor found = " + divisor);
            return false;
        }
        divisor = divisor.add(BigInteger.ONE);
        if(divisor.compareTo(max) >= 0)
        {
            //end of loop
            return true;
        }
    }
}

//Returns the prime number p of the first n Mersenne primes 2^(p-1)
private static ArrayList<Integer> merseneIndex(int n)
{
    ArrayList<Integer> indices = new ArrayList<Integer>();
    BigInteger a;
    boolean isPrime;
    
    //Count the number of primes found
    int count;
    //2 is the only even prime and 2^(p-1) = 2 is also prime
    indices.add(2);
    count = 1;
    //We will continue searching odd numbers from 3 onwards
    int number = 3;
    
    infiniteLoop:
        while(true)
        {
            //Consider odd numbers only
            //System.out.println(number);
            for(int i = 3; i <= (int)Math.floor(number/2); i++)
            {
                if(number % i == 0)
                {
                    //Found a divisor other than 1 and number
                    //ergo number is not prime
                    //System.out.println(number + " is not prime");
                    number += 2;
                    continue;
                }
            }
            //No divisors found so number is prime
            //System.out.println(number + " is prime");
            //Is 2^(p-1) prime?
            a = new BigInteger("2").pow(number).add(new BigInteger("-1"));
            //is 'a' prime?
            isPrime = isPrime(a);
            System.out.println(number + ": " + a + " is prime? " + isPrime);
            if(isPrime)
            {
                indices.add(number);
                
                count++;
                //System.out.println(count);
                if(count >= n)
                {
                    break infiniteLoop;
                }
                number += 2;
            }
            else
            {
                number += 2;
                continue;
            }
            
        }
        
        return indices;
}    

//Finds the first n perfect numbers from prime numbers
//Problem - not all prime numbers give perfect numbers in this way
//2^(p-1) must also be prime
private static ArrayList<BigInteger> perfectSeries(ArrayList<Integer> primes, int n)
throws IllegalArgumentException
{
    if(n > primes.size() || primes.isEmpty())
    {
        //Illegal arguments
        throw new IllegalArgumentException("Not enough primes in supplied list to complete request");
    }
    
    ArrayList<BigInteger> perfectNumbers = new ArrayList<BigInteger>();
    int number;
    BigInteger perfect;
    BigDecimal a, b, c;
    
    for(int i = 0; i < n; i++)
    {
        //Use BigDecimals in case powers become too large
        number = primes.get(i);
        a = new BigDecimal(2).pow(number-1);
        b = new BigDecimal(2).pow(number).subtract(new BigDecimal(1));
        c = a.multiply(b);
        
        //Only first 6 (?) small enough to fit into Int32
        perfect = c.toBigInteger();
        perfectNumbers.add(perfect);
    }
    
    return perfectNumbers;
}
I am not suggesting that this is the best approach, it was merely experimental.

Note that unless BigIntegers are used we can only obtain the first 6 or so, any larger perfect
numbers will be in error due to an overflow error - an integer can only hold a number as large as
2^32, a long as large as 2^64. Since perfect numbers are few and far between, they rapidly become
enormous and so our function must be able to cope with arbitrarily large numbers.

The code can be tested in the Main function as follows:
Where we have used a brute-force approach to find the prime numbers by testing if a divisor of each
test number exists. We can then feed the numbers obtained in to our equation for a Mersenne Prime
and see which, if any, yield Mersenne Primes:

//Find the first 100 prime numbers       
   ArrayList<Integer> primes = findPrime(100);
   for(int prime : primes)
   {
       System.out.println(prime);
   }

//Test whether the number 2^p - 1 for large p, e.g. p = 31 is prime
   int number = 31;
   BigInteger a = new BigInteger("2").pow(number).add(new BigInteger("-1"));
   System.out.println(a);
    
   System.out.println(a + "is prime? " + isPrime(a));

//Find the first 10 primes p for the Mersene primes             
   ArrayList<Integer> merceneIndices = merseneIndex(10);
   for(int index : merceneIndices)
   {
       System.out.println(index);
   }

//Find the first 10 perfect numbers using these primes
   ArrayList<BigInteger> perfects = perfectSeries(merceneIndices, 10);
   for(BigInteger perfect : perfects)
   {
       System.out.println(perfect);
   }

}
The second test returns true after about 11 minutes on my PC! However, for p = 61, the Mersenne
Prime becomes so large that testing whether or not it is prime, by this brute force approach, will take
something like 100 000 times as long, that is about 20 years on my PC! (With the program thread
running in the background so that I can still use my PC for other things!).

The third test obtains the first 10 Mersenne Primes and the fourth test returns the first ten perfect
numbers corresponding to them. However, the code simply takes too long to finish - requiring, as the
second test shows, some 20 years to find the 9th Mersenne Prime. We can, however, quite readily find
the first 8 perfect numbers by this method, which puts us level with the great Leonhard Euler in the 19th
century (who used trial division without an electronic computer of course!).

The first test gives the first 100 prime numbers correctly as:
Some Mersenne Primes are known, for example the first few primes, p, that correspond to the first 9
Mersenne Primes, M, are:

p = 2: M = 3
3: 7
5: 31
7: 127
13: 8191
17: 131071
19: 524287
31: 2147483647
61: 2305843009213693951

So for example, for p = 3, 2^p - 1 = 2^3 - 1 = 8 - 1 = 7, and 7 is a prime number.
The prime numbers 11, give no Mersenne Primes.     
2
3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
61
67
71
73
79
83
89
95
97
101
103
107
109
113
121
127
131
137
139
149
151
157
163
167
173
179
181
191
193
197
199
209
211
221
223
227
229
233
239
241
251
257
263
269
271
277
281
283
293
305
307
311
313
317
323
329
331
337
347
349
353
359
367
373
379
383
389
395
397
401
409
419
421
431
433
439
443
449
455
457
461
463
467
BUILD SUCCESSFUL (total time: 0 seconds)
So far, so good for computing the primes. Using a brute force algorithm to find the primes is OK,
since we would have to get very far indeed into the perfect number series for the prime indices, p, to
become too large to handle. Before we look at how to get around the problem of calculating such
large Mersenne Primes, we side-track to another algorithm for obtaining prime numbers.

The Sieve of Eratosthenes

As an experiment I looked into another method: the Sieve of Eratosthenes, a very good explanation
of which can be found on Wikipedia
here (the animation is especially useful). Essentially a series or
grid (or 1D array) of positive integers of size equal to the largest prime is established and then it is
determined which numbers within that grid are prime, by beginning with the first number and
eliminating all multiples of it, then moving to the next number, which will be the next prime, and
eliminating all of its multiples, which are not prime, and so on until every number in the grid has been
checked. I wrote the following code to implement the Sieve of Eratosthenes:
package perfectnumbers;

import java.util.*;
import java.math.BigDecimal;
import java.math.BigInteger;

// @author BotRejectsInc.

public class Main
{
public static void main(String[] args)
{
 //Sieve of Eratosthenes       
   boolean[] myPrimes = sieveOfEratosthenes(100);
   for(int i = 0; i <= myPrimes.length - 1; i++)
   {
       System.out.println(i + " isPrime? " + myPrimes[i]);
   }
 }

//Uses the sieve of Eratosthenes to return an array of small prime numbers
private static boolean[] sieveOfEratosthenes(int prime)
{
   boolean[] primes = new boolean[prime+1];
   //We know that 0 and 1 are not prime        
   primes[0] = false;
   primes[1] = false;
   //Initialise rest of array
   for(int i = 2; i <= prime; i++)
   {
       primes[i] = true;
   }
   //Start sieve with 2 and multiples of 2
   int number = 2;
   for(int i = number + number; i <= prime; i = i + number)
   {
       primes[i] = false;
   }
   boolean nextFound = true;
   //Find next prime
   do
   {
       for(int i = number + 1; i <= prime + 1; i ++)
       {
           if(primes[i] == true)
           {
               number = i;
               break;
           }
           nextFound = false;
       }
       //Sieve out multiples of number
       for(int i = number + number; i <= prime; i = i + number)
       {
           primes[i] = false;
       }
   } while(nextFound);
   
   return primes;
}
}
Note we use a boolean array to store whether or not each element is a prime, as this is economical
on memory.

This algorithm is certainly fast, though for the size of primes we need it really makes no significant
time difference whether or not we use the brute force approach or the Sieve of Eratosthenes. The
latter has the advantage of returning all primes within the range specified by the one integer
parameter, which is the largest number we wish to test. However, it has limits - an array can only hold
up to 2^32 elements, since its elements are indexed by a 32 bit integer. We could use an ArrayList,
but in any case we would run out of heap memory before too long (in fact my PC ran out of heap
space even before the capacity of the boolean array could be filled). It might be possible to use the
individual bits of a BigInteger to represent each number (1 for prime, 0 otherwise) but otherwise
using brute force is very fast in obtaining the primes we need for our indices, p. There are other
algorithms, besides the Sieve of Eratosthenes, but again we seem likely to gain nothing significant
here for our purposes.

We still have the problem, however, that brute force fails in obtaining the large Mersenne Primes, M,
that we need. It simply takes too long to confirm which numbers obtained from our formula are
actually Mersenne Primes.

Incidentally, rather than determining whether or not a number has divisors, we could determine
whether or not it is a multiple of any smaller integer, like this:

//Determines whether or not a number is a multiple (has factors) using 'brute force'
private static boolean isMultiple(BigInteger multiple)
throws IllegalArgumentException
{
    if(multiple.compareTo(BigInteger.ZERO) < 0)
    {
        //Negative argument not allowed
        throw new IllegalArgumentException("Negative argument not allowed!");
    }
    
    boolean isMultiple = false;
    
    BigInteger number = new BigInteger("2");
    BigInteger TWO = new BigInteger("2");
    if(multiple.compareTo(number) < 0)
    {
        //Number is less than 2
        return false;
    }
    while(number.compareTo(multiple.divide(TWO)) <= 0)
    {
        if(multiple.mod(number).compareTo(BigInteger.ZERO) == 0)
        {
            //multiple found so number not prime
            System.out.println(multiple + " Is a multiple of " + number);
            return true;
        }
        else
        {
            number = number.add(BigInteger.ONE);
        }
    }
    
    return isMultiple;
}
However, this alternative brute force approach gives no speed advantage over testing divisors. (It
takes about 20 minutes on my PC to determine whether or not M = 2147483647 (for p = 31) is a
prime number.

What we need is a clever mathematical short-cut. At least one such shortcut exists - the LLT
(Lucas-Lehmer primality test) which is explained quite well on wikipedia,
here and efficiently
determines whether or not a number is a Mercenne Prime. My Java implementation of this algorithm
is shown below (complete program):

package perfectnumbers;

import java.util.*;
import java.math.BigDecimal;
import java.math.BigInteger;


// @author BotRejectsInc.

public class Main
{

/** Creates a new instance of Main */
public Main()
{
}

/**
 * @param args the command line arguments
 */
public static void main(String[] args)
{
    int n = 11;
    ArrayList<BigInteger> perfectSeries = findPerfects(n);
    for(BigInteger perfect : perfectSeries)
    {
        System.out.println(perfect);
    }
 }

//Finds the first n perfects using LLT to obtain Mersenne primes
private static ArrayList<BigInteger> findPerfects(int n)
{
    ArrayList<BigInteger> perfects = new ArrayList<BigInteger>();
    perfects.add(new BigInteger("6"));
    
    int count = 1;
    int p = 1;
    while(count < n+1)
    {
        if(isMersennePrime(p))
        {
            count ++;
            BigInteger mersenne = new BigInteger("2").pow(p).subtract(BigInteger.ONE);
            BigInteger factor = new BigInteger("2").pow(p-1);
            BigInteger perfect = mersenne.multiply(factor);
            perfects.add(perfect);
        }
        p++;
    }
    
    return perfects;
}  

//Determines whether a number is a Mersenne prime using LLT
static boolean isMersennePrime(int p)
{        
    BigDecimal seriesSum = new BigDecimal(4);
    BigDecimal TWO = new BigDecimal(2);
    BigDecimal ZERO = new BigDecimal(0);
    
    BigDecimal mersenne = TWO.pow(p).subtract(new BigDecimal(1));
    //System.out.println("Testing " + mersenne + " with p = " + p + ": ");
    
    //System.out.println("0: " + seriesSum);
    
    for(int i = 1; i <= p - 2; i++)
    {
        seriesSum = seriesSum.multiply(seriesSum).subtract(TWO);            
        seriesSum = seriesSum.remainder(mersenne);            
        //System.out.println(i + ": " + seriesSum);     
    }
     
    //System.out.println(seriesSum);
    if(seriesSum.compareTo(ZERO) == 0)
    {
        return true;
    }
    
    return false;
}
}
This program calculates the first eleven perfect numbers quite rapidly. Here is the output:

6
28
496
8128
33550336
8589869056
137438691328
2305843008139952128
2658455991569831744654692615953842176
191561942608236107294793378084303638130997321548169216
13164036458569648337239753460458722910223472318386943117783728128
14474011154664524427946373126085988481573677491474835889066354349131199152128
BUILD SUCCESSFUL (total time: 7 seconds)

Going further becomes increasingly difficult, computations slow as the numbers become larger and the
perfect numbers further apart. The 11th perfect number was first discovered in 1914, so we are doing
quite well now. We could do much better with a supercomputer.

GIMPS

The search for perfect numbers goes on. It is not known whether or not they are infinite in number - no
proof is currently known. The biggest supercomputer is the Internet and the
GIMPS (Great Internet
Mersenne Prime Search) project enlists computers via the WWW to find more perfect numbers.

To read more about GIMPS or partake in their project using your own computer, go to
www.mersenne.org.