**[**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**]**