generate N perfect numbers.
Ricardo Nuno
%divisible(10,2).
divisible(X,Y):- N is Y*Y,N =< X,X mod Y =:= 0.
divisible(X,Y):- Y < X, Y1 is Y + 1, divisible(X,Y1).
%isprime([3],Z).
isprime([X|_],X):-Y is 2, X >1, \+divisible(X,Y).
isprime([_|T],Z):-isprime(T,Z).
%Calculate the power of one number
%ie power(2,10,R)
power(_,0,1):-!.
power(N,K,R):-K1 is K-1,
power(N,K1,R1),
R is R1*N.
%formula of perfect numbers 2^(p-1)*(2^p-1)
%ie calc(2,10,R)
calc(2,K,R):-power(2,K,X),R1 is X-1,
power(2,K-1,R2),
R is R1 * R2.
%using lists
%ie calc([2,3,4],R).
listperf([K|_],R):-calc(2,K,R).
listperf([_|T],Z):-listperf(T,Z).
%generate one list of N numbers.
%genList(10,L).
generateList(0,[]).
generateList(N,[X|Xs]):- N > 0,
X is N+1,
N1 is N-1,generateList(N1,Xs).
%list of N perfect numbers
%perfect(100,C)
perfect(N,C):-generateList(N,R),findall(L,isprime(R,L),P),listperf
(P,C).
Chip Eastham
Strictly speaking the formula you have:
(2^(p-1))*((2^p)-1)
where the last factor is a (Mersenne) prime
is for _even_ perfect numbers.
It is not known whether an odd perfect number
exists, though none do below 10^300 (Brent,
Cohen, and te Riele ~ 1989).
I've not checked running your code, but it
appears you intend to prove primality of X
from trial division by 2:
> %isprime([3],Z).
> isprime([X|_],X):-Y is 2, X >1, \+divisible(X,Y).
> isprime([_|T],Z):-isprime(T,Z).
What will happen with the goal isprime([9],Z)?
regards, chip
Chip Eastham
With respect to what happens with isprime([9],Z),
I tried and it correctly reports "No", because
divisible(X,Y) checks whether X is divisible
by any integer between Y and sqrt(X). [The
iteration actually runs through all integers up
to X, but only checks divisibility by those
below sqrt(X), so there's a bit of room for
improvement.]
However perfect(10,C) returns (with backtracking)
five values for C, the first of which is _not_ a
perfect number, C = 2096128. This corresponds
to p = 11, but 2^11 - 1 = 2047 is composite:
2047 = 23*89
You'll need to check not only that p is prime,
but that 2^p - 1 is also prime.
regards, chip
bart demoen
> %Create perfect numbers
[...]
Thank you for another benchmark :-) [I am not mocking you, but I do
at this point not care about the correctness of your program versus
your intentions]
I tested/timed the query ?- perfect(5000,_), fail.
on
SWI -O: 2489 msecs
ECLiPSe: 1860
SICStus: 1360
hProlog: 1090
XSB, B-Prolog: wrong results because they have no bigints
Yap: bug
:-)
Cheers
Bart Demoen
Chip Eastham
Since 2^p - 1 grows much faster than p, it is not
as practical to test the primality of 2^p - 1 by
trial division as it is to test exponents p in that
way. Fortunately we have the Lucas-Lehmer primality
test for Mersenne primes:
[Lucas-Lehmer primality test -- Wikipedia]
http://en.wikipedia.org/wiki/Lucas%E2%80%93Lehmer_primality_test
which determines if a candidate 2^p - 1 is a Mersenne
prime in p-2 steps of arithmetic mod (2^p - 1). This
is the efficient test used by GIMPS:
[GIMPS - The Great Internet Mersenne Prime Search]
http://www.mersenne.org/
regards, chip
Ricardo Nuno
Thank you
I didnt knowed well ,the first part of this formula. I will use the
Lucas-Lehmer test for Mersenne primes.
regards, Ricardo
Ricardo Nuno
This program need to improve performance with a prime numbers test
faster...
%----------------
%Find perfect numbers
hasFactor(_,0):-!,fail.
hasFactor(N,L):-L<N,N mod L=:=0.
hasFactor(N,L):-L<N,L2 is L+1,hasFactor(N,L2).
%test of prime number
%isprime([3],Z).
isprime([X|_],X):-Y is 2, X >1, \+hasFactor(X,Y).
isprime([_|T],Z):-isprime(T,Z).
%multiply with sums
%product of natural numbers with successives sums:
% |0 se y=0
% mul(x,y)|
% |x+mul(x,y-1) se y>0
mul(_,0,0):-!.
mul(X,Y,Z):-X>0,Y>0,
Y1 is Y-1,
mul(X,Y1,Z1),Z is X+Z1.
%Calculate the power of one number
%ie power(2,10,R)
power(_,0,1):-!.
power(N,K,R):-K1 is K-1,
power(N,K1,R1),
mul(R1,N,R).
%generate one list of N numbers.
%genList(10,L).
generateList(0,[]).
generateList(N,[X|Xs]):- N > 0,
X is N+1,
N1 is N-1,generateList(N1,Xs).
%formula of perfect numbers (2^p-1)*2^(p-1)
%ie calc(2,10,R)
calc(2,K,R):-power(2,K,X),R1 is X-1,primo(R1),
power(2,K-1,R2),
mul(R1,R2,R).
primo(X):-Y is 2, X >1, \+hasFactor(X,Y).
%output the list
noLst([X|_],X).
noLst([_|T],Z):-noLst(T,Z).
%calculate perfect numbers
%ie perfect(12,N).
perfect(N,R):-generateList(N,Xs),noLst(Xs,K),calc(2,K,R).