arun@dsrgsun.ces.cwru.edu (Arun Lakhotia) (12/22/88)
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
bimandre@kulcs.uucp (Andre Marien) (01/10/89)
The following program is about 4 times faster than all previous posted versions, using the normal search method. It would run a lot faster if you new perfect numbers are all even. I could not prove that in short time. Anyone ? Timings were done with BIM_Prolog. perfect(_X) :- gen(_X,1), sumd(_X,1,2) . sumd(_X,_PS,_Y) :- _X mod _Y =:= 0, !, _NS is _PS + _X/_Y + _Y, check(_X,_NS,_Y) . sumd(_X,_PS,_Y) :- _NY is _Y + 1, _NY*_NY =< _X, sumd(_X,_PS,_NY) . check(_X,_S,_Y) :- _S < _X, !, _NY is _Y + 1, _NY*_NY =< _X, sumd(_X,_S,_NY) . check(_X,_S,_Y) :- _S > _X, !, _Y*_Y =:= _X, _X =:= _S - _Y . check(_X,_S,_Y) :- _Y*_Y < _X, _NY is _Y + 1, nodiv(_X,_NY) . nodiv(_X,_Y) :- _Y*_Y =< _X, !, _X mod _Y =\= 0, _NY is _Y + 1, nodiv(_X,_NY) . nodiv(_X,_Y) . gen(_N,_N) . gen(_N,_S) :- _NS is _S + 1, gen(_N,_NS) . z :- perfect(_N), write(_N) . test :- perfect(_N), _N >= 8128, ! . timetest :- time(test) . timing results (SUN 3/50) 42.66 42.5 Andre' Marien B.I.M. mcvax!prlb2!kulcs!bimandre