[Back to MATH SWAG index]  [Back to Main SWAG index]  [Original]

{
>Samiel (samiel@fastlane.net) wrote:
>: Here's my fast and elegant code... snarf it and add it to SWAG...
>: 
>: {****************************************************************
>: Perfect Numbers...
>: 
>: A Perfect number is a number whose divisors not including the original
>: number add up to the original number.  An example is 6, 6=3*2*1=3+2+1
>: and 28, 28=14*7*4*2*1=14+7+4+2+1.  The definition of a Perfect number
>: can also be defined as,
>: 2^(p-1)*(2^p-1) where p is prime and 2^p-1 is a Mersenne prime...
>: 
>: A Mersenne prime is defined as,
>: 2^p-1 where p is prime and 2^p-1 is prime, so 2^2-1=3 is a Mersenne
>: prime where p is 2, and 2^3-1=5 is a Mersenne prime where p is 3.
>: 
>: Mersenne primes under 2^32 (32-bit) are:
>:   [2] -> 3
>:   [3] -> 7
>:   [5] -> 31
>:   [7] -> 127
>: [13] -> 8191
>: [17] -> 131071
>: [19] -> 524287
>: [31] -> 2147483647
>: 
>: (Values of p are in braces [])
>: 
>: The Perfect numbers under 2^32 are:
>: Perfect numbers...
>:   [2] 6
>:   [3] 28
>:   [5] 496
>:   [7] 8128
>: [13] 33550336
>: 
>: (Values of p are in braces [])
>: 
>: Here is my code...
>:
>: ****************************************************************}

 PROGRAM Perfect;
 {Computes Perfect Numbers...}

 VAR
   tmp,num:longint; {Long Integer signed 32-bit (31-bit on each side)}
   j,k:byte;

 { Slow? Fast? way to find primes, dividing by odd numbers }
 Function IsPrime(num:longint):boolean;
 Var
   tmp:boolean;
   j:longint;
 Begin
   tmp:=true;
   if num mod 2=0 then
     tmp:=false;
   for j:=3 to round(sqrt(num)) do
     if odd(j) then
       if num mod j=0 then
         tmp:=false;
   if num=2 then
      tmp:=true;
   IsPrime:=tmp;
 End;

 BEGIN
   tmp:=2; {2^1}
   writeln('Perfect numbers...');
   for j:=2 to 31 do
     begin
       tmp:=tmp*2; {Raise 2 to another power}
       if IsPrime(j) then
         begin
           num:=tmp-1;
           if IsPrime(num) then
             begin
               num:=num*(tmp div 2);
               writeln(num:1,' is perfect'); {Ignore negatives}
             end;
         end;
     end;
 END.


>: - Samiel
>: samiel@fastlane.net
>: http://www.fastlane.net/~samiel
>: 

>Considering that 2^859433-1 is prime(and is the largest known prime), there
>are obviously better methods.  Use the Lucas-Lehmer test.

>B(n+1) = B(n)^2 - 2 mod (2^p-1) where B(0)=4

>if B(n-1) = 0 then 2^p-1 is prime.  The division is a special case and is
>done easily because of the all ONE's when 2^p-1 is written in binary.  All
>that is necessary is a fast implementation of multiplication and this can
>be done with FFT's.

>See http://www.utm.edu/research/primes/mersenne.shtml

>or for software
>http://ourworld.compuserve.com/homepages/justforfun/prime.htm

Well, considering we are looking at numbers under 32 bits, your code
would probably not get done as fast as mine, though in the long run,
(say over 64 bits or so) it would.  All the FFT's and rather long
multiplication would slow it down considerably.  Even if you could get
a really fast implementation of FFT's and multiplications, we're
talking about a net savings of about .625 seconds or so with numbers
under 32 bits on a computer with a Math Processor... so mine's good
enough... although a few of multiplications could be taken out...

- Samiel
samiel@fastlane.net
http://www.fastlane.net/~samiel

[Back to MATH SWAG index]  [Back to Main SWAG index]  [Original]