%Create perfect numbers

%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).

On Jan 24, 12:02=A0pm, Ricardo Nuno <9...@bigfoot.com> wrote:
> %Create perfect numbers
>
> %divisible(10,2).
> divisible(X,Y):- N is Y*Y,N =3D< X,X mod Y =3D:=3D 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,
> =A0 =A0 =A0 =A0 =A0 =A0 =A0 power(N,K1,R1),
> =A0 =A0 =A0 =A0 =A0 =A0 =A0 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,
> =A0 =A0 =A0 =A0 =A0 =A0 =A0power(2,K-1,R2),
> =A0 =A0 =A0 =A0 =A0 =A0 =A0R 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,
> =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0X is N+1,
> =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0N1 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).

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

On Jan 24, 12:42=A0pm, Chip Eastham <hardm...@gmail.com> wrote:
> On Jan 24, 12:02=A0pm, Ricardo Nuno <9...@bigfoot.com> wrote:
>
>
>
> > %Create perfect numbers
>
> > %divisible(10,2).
> > divisible(X,Y):- N is Y*Y,N =3D< X,X mod Y =3D:=3D 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,
> > =A0 =A0 =A0 =A0 =A0 =A0 =A0 power(N,K1,R1),
> > =A0 =A0 =A0 =A0 =A0 =A0 =A0 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,
> > =A0 =A0 =A0 =A0 =A0 =A0 =A0power(2,K-1,R2),
> > =A0 =A0 =A0 =A0 =A0 =A0 =A0R 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,
> > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0X is N+1,
> > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0N1 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).
>
> 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

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 =3D 2096128.  This corresponds
to p =3D 11, but 2^11 - 1 =3D 2047 is composite:

2047 =3D 23*89

You'll need to check not only that p is prime,
but that 2^p - 1 is also prime.

regards, chip

On Sun, 24 Jan 2010 09:02:59 -0800, Ricardo Nuno wrote:

> %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

On Jan 24, 1:33=A0pm, Chip Eastham <hardm...@gmail.com> wrote:
> On Jan 24, 12:42=A0pm, Chip Eastham <hardm...@gmail.com> wrote:
>
>
>
>
>
> > On Jan 24, 12:02=A0pm, Ricardo Nuno <9...@bigfoot.com> wrote:
>
> > > %Create perfect numbers
>
> > > %divisible(10,2).
> > > divisible(X,Y):- N is Y*Y,N =3D< X,X mod Y =3D:=3D 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,
> > > =A0 =A0 =A0 =A0 =A0 =A0 =A0 power(N,K1,R1),
> > > =A0 =A0 =A0 =A0 =A0 =A0 =A0 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,
> > > =A0 =A0 =A0 =A0 =A0 =A0 =A0power(2,K-1,R2),
> > > =A0 =A0 =A0 =A0 =A0 =A0 =A0R 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,
> > > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0X is N+1,
> > > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0N1 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).
>
> > 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
>
> 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). =A0[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 =3D 2096128. =A0This corresponds
> to p =3D 11, but 2^11 - 1 =3D 2047 is composite:
>
> 2047 =3D 23*89
>
> You'll need to check not only that p is prime,
> but that 2^p - 1 is also prime.
>
> regards, chip

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

Thank you
I  didnt knowed well ,the first part of this formula. I will use the
Lucas-Lehmer  test for Mersenne primes.

regards, Ricardo

correction of last program.
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).