🚀 go-pugleaf

/*

There seems to be some suggestion that the question of perfect numbers has
been worked to death.  Not so!  Previous programs have used the theorem:

>   if (2^p-1) is a (mersenne) prime 
>   then (2^p-1) * 2^(p-1) is a perfect number.

The program below makes use of the fact that there are not many values of p
for which 2^p-1 is a prime, and the first 24 of them are listed in Knuth's
"The Art of Computer Programming", Volume 2.  (Look in the index under
"Mersenne primes".)  These values of p are tabulated below.

This program uses the Quintus library package "long" to handle the large
integers which are involved.  I'm sorry to say that it runs up against a
current limitation in "long" that exponents over 9999 are not allowed.  All
the same, it does manage to generate 22 perfect numbers, the last of which
is a 5985 digit number (if I counted correctly).  That number (p=9941) took
199.5 seconds to generate on a Sun-3/60.  The first 12 numbers, with the
times taken to generate them on a Sun-3/60, are:

6
[P=2, 0.017 sec]

28
[P=3, 0.033 sec]

496
[P=5, 0.033 sec]

8128
[P=7, 0.017 sec]

33550336
[P=13, 0.016 sec]

8589869056
[P=17, 0.033 sec]

137438691328
[P=19, 0.017 sec]

2305843008139952128
[P=31, 0.033 sec]

2658455991569831744654692615953842176
[P=61, 0.033 sec]

191561942608236107294793378084303638130997321548169216
[P=89, 0.033 sec]

13164036458569648337239753460458722910223472318386943117783728128
[P=107, 0.050 sec]

14474011154664524427946373126085988481573677491474835889066354349131199152128
[P=127, 0.100 sec]

*/

:- ensure_loaded(library(long)).    % for eval/[1,2], portray_number/1

portray(N) :- portray_number(N).

perfect_numbers :-
	mersenne_prime_exponent(P),
	statistics(runtime, [T0|_]),
	eval((2^P)-1, M),	    % M = 2^P-1 is a Mersenne prime
	eval(M*(2^(P-1)), N),	    % N = M*2^(P-1)
	statistics(runtime, [T1|_]),
	Time is T1 - T0,
	format('~p~n[P=~p, ~3D sec]~n~n', [N, P, Time]),
	fail.
perfect_numbers.

mersenne_prime_exponent(2).
mersenne_prime_exponent(3).
mersenne_prime_exponent(5).
mersenne_prime_exponent(7).
mersenne_prime_exponent(13).
mersenne_prime_exponent(17).
mersenne_prime_exponent(19).
mersenne_prime_exponent(31).
mersenne_prime_exponent(61).
mersenne_prime_exponent(89).
mersenne_prime_exponent(107).
mersenne_prime_exponent(127).
mersenne_prime_exponent(521).
mersenne_prime_exponent(607).
mersenne_prime_exponent(1279).
mersenne_prime_exponent(2203).
mersenne_prime_exponent(2281).
mersenne_prime_exponent(3217).
mersenne_prime_exponent(4253).
mersenne_prime_exponent(4423).
mersenne_prime_exponent(9689).
mersenne_prime_exponent(9941).
mersenne_prime_exponent(11213).
mersenne_prime_exponent(19937).