40 lines
874 B
Text
40 lines
874 B
Text
go2 =>
|
|
println("Using the formula: 2^(p-1)*(2^p-1) for prime p"),
|
|
foreach(P in primes(32))
|
|
PF=perfectf(P),
|
|
% Check that it is really a perfect number
|
|
if perfect2(PF) then
|
|
printf("%w (prime %w)\n",PF,P)
|
|
end
|
|
end,
|
|
nl.
|
|
|
|
% Formula for perfect number candidates:
|
|
% 2^(p-1)*(2^p-1) where p is a prime
|
|
%
|
|
perfectf(P) = (2**(P-1))*((2**P)-1).
|
|
|
|
% Faster check of a perfect number
|
|
perfect2(N) => sum_divisors(N) == N.
|
|
|
|
% Sum of divisors
|
|
table
|
|
sum_divisors(N) = Sum =>
|
|
sum_divisors(2,N,1,Sum).
|
|
|
|
sum_divisors(I,N,Sum0,Sum), I > floor(sqrt(N)) =>
|
|
Sum = Sum0.
|
|
|
|
% I is a divisor of N
|
|
sum_divisors(I,N,Sum0,Sum), N mod I == 0 =>
|
|
Sum1 = Sum0 + I,
|
|
(I != N div I ->
|
|
Sum2 = Sum1 + N div I
|
|
;
|
|
Sum2 = Sum1
|
|
),
|
|
sum_divisors(I+1,N,Sum2,Sum).
|
|
|
|
% I is not a divisor of N.
|
|
sum_divisors(I,N,Sum0,Sum) =>
|
|
sum_divisors(I+1,N,Sum0,Sum).
|