Yet another perfect number program
Author Message Yet another perfect number program

When the Cat is away the mice play. Now that RO'K is ...

Following is a program for generating perfect numbers. It is different
from others as it uses the theorem cited by Thomas Sj|land.

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

NOTE: The isprime/1 algorithm is substantially different and faster than
what I posted earlier. It is *NOT* derived from any previous version.

The previous isprime algorithm checked for the existence of a divisor
for N between 1 .. N//2+1. This algorithm checks for divisors between
I .. N//I+1 starting from I = 1.

The algorithm is substantially faster and gives the 6, 28, 496, 8128,
and 33550336 perfect numbers in a wink. Numbers beyond that exceed the
integer bounds of SICStus and Quintus Prologs.

/* checking for primeness */

divides(Dividend,Divisor,Quotient) :-
Quotient is Dividend//Divisor,
Dividend is Quotient * Divisor.

isprime(2).
isprime(Number) :-
\+ divides(Number, 2),
isprime(Number, 1, Number).

isprime(_M, Lower, Upper) :-
Lower >= Upper,!.
isprime(Number, Lower, _Upper) :-
%       Lower < _Upper,
divides(Number, Lower, Quotient),!,
Lower = 1,
Lower1 is Lower +1,
isprime(Number, Lower1, Quotient).
isprime(Number, Lower, _Upper) :-
%       Lower < _Upper,
Lower1 is Lower +1,
Upper1 is Number//Lower+1,
isprime(Number, Lower1, Upper1).

/*     Perfect Numbers  */

perfect_numbers :-
write('P. No'), tab(4), write('Merseme Prime = 2^p-1'),nl,
perfect_number(PNo, Merseme/N),
write(PNo), tab(3), write((Merseme = 2^N-1)), nl,
fail.

perfect_number(PNo, Merseme/N) :-
candidate(PNo, Merseme/N),
isprime(Merseme).      % Is merseme prime

candidate(PNo, M/N) :-
int(N), N > 1,
pow2(N, P),         % P = 2^N
M is P-1, N1 is N-1,
pow2(N1, P1),       % 2^(N-1)
PNo is M * P1.      % 2^N-1 * 2^(N-1)

int(1).
int(I) :- int(I1), I is I1+1.

pow2(N, P) :- P is 1 << N.

=

Arun

Wed, 19 May 1993 05:54:00 GMT

 Page 1 of 1 [ 1 post ]

Relevant Pages