Perfect Numbers |

In my quest to undertake various programming projects, for the sheer fun of it, I came across a definition

of perfect numbers:

of perfect numbers:

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:

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

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

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:

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.

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:

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:

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:

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.

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

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

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

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)

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:

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.

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:

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):

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.

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.

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.