🚀 go-pugleaf

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.

>   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