Amicable pairs: Difference between revisions
| Line 189: | Line 189: | ||
=={{header|Pascal}}== |
=={{header|Pascal}}== |
||
extremly fast. |
extremly fast. |
||
up to |
up to 524'000'000 in 11.3 secs (linux 32-Bit fpc 2.6.4. i4330 3.5 Ghz ) |
||
n = p[1]^a[1]*p[2]^a[2]*...p[l]^a[l] |
|||
sum of divisors= s(n) = (p[1]^(a[1]+1) -1) / (p[1] -1) * ... * (p[l]^(a[l]+1) -1) / (p[l] -1) with |
|||
p[k]^(a[k]+1) -1) / (p[k] -1) = sum (i= [1..a[k]])(p[k]^i) |
|||
Using "Sieve of Erathosthenes"-style |
Using "Sieve of Erathosthenes"-style |
||
| Line 302: | Line 305: | ||
function NextNPotCnt(pot,prim: tValue):tValue; |
function NextNPotCnt(pot,prim: tValue):tValue; |
||
//return the first power <> 0 |
//return the first power[pot+?] <> 0 |
||
//power == n to base prim |
//power == n to base prim |
||
var |
var |
||
| Line 494: | Line 497: | ||
T0:= time; |
T0:= time; |
||
Init; |
Init; |
||
IF MAX < 1000*1000 then |
|||
Sieve(2) |
|||
| ⚫ | |||
else |
|||
| ⚫ | |||
T1:= time; |
T1:= time; |
||
Check; |
Check; |
||
| Line 541: | Line 546: | ||
00:00:00.000 |
00:00:00.000 |
||
..... |
..... |
||
{.. |
|||
{.. linux 32-Bit |
|||
509379344 523679536 |
509379344 523679536 |
||
[1,2,4,8,16,23,46,92,184,347,368,694,1388,2776,3989,5552,7978,7981,15956,15962,31912,31924,63824,63848,91747,127696,183494,366988,733976,1384183,1467952,2768366,5536732,11073464,22146928,31836209,63672418,127344836,254689672] |
[1,2,4,8,16,23,46,92,184,347,368,694,1388,2776,3989,5552,7978,7981,15956,15962,31912,31924,63824,63848,91747,127696,183494,366988,733976,1384183,1467952,2768366,5536732,11073464,22146928,31836209,63672418,127344836,254689672] |
||
| Line 548: | Line 553: | ||
442 amicable numbers upto 524000000 |
442 amicable numbers upto 524000000 |
||
00:00:11.213 |
00:00:11.213 |
||
12 |
|||
real 0m15.222s |
real 0m15.222s |
||
user 0m14.977s |
user 0m14.977s |
||
sys 0m0.223s |
sys 0m0.223s |
||
} |
|||
} |
} |
||
Revision as of 15:02, 31 December 2014
You are encouraged to solve this task according to the task description, using any language you may know.
Two integers and are said to be amicable pairs if and the sum of the proper divisors of () as well as .
For example 1184 and 1210 are an amicable pair (with proper divisors 1, 2, 4, 8, 16, 32, 37, 74, 148, 296, 592 and 1, 2, 5, 10, 11, 22, 55, 110, 121, 242, 605 respectively).
- Task
Calculate and show here the Amicable pairs below 20,000; (there are eight).
- Cf.
- Proper divisors
- Abundant, deficient and perfect number classifications
- Aliquot sequence classifications and its amicable classification.
AutoHotkey
<lang d>SetBatchLines -1 Loop, 20000 { m := A_index
; Getting factors loop % floor(sqrt(m)) { if ( mod(m, A_index) = 0 ) { if ( A_index ** 2 == m ) { sum += A_index continue } else if ( A_index != 1 ) { sum += A_index + m//A_index } else if ( A_index = 1 ) { sum += A_index } } } ; Factors obtained
; Checking factors of sum if ( sum > 1 ) { loop % floor(sqrt(sum)) { if ( mod(sum, A_index) = 0 ) { if ( A_index ** 2 == sum ) { sum2 += A_index continue } else if ( A_index != 1 ) { sum2 += A_index + sum//A_index } else if ( A_index = 1 ) { sum2 += A_index } } } if ( m = sum2 ) && ( m != sum ) && ( m < sum ) final .= m . ":" . sum . "`n" } ; Checked
sum := 0 sum2 := 0 } MsgBox % final ExitApp</lang>
- Output:
220:284 1184:1210 2620:2924 5020:5564 6232:6368 10744:10856 12285:14595 17296:18416
D
<lang d>void main() /*@safe @nogc*/ {
import std.stdio, std.algorithm, std.range, std.typecons, std.array;
immutable properDivs = (in uint n) pure nothrow @safe /*@nogc*/ =>
iota(1, (n + 1) / 2 + 1).filter!(x => n % x == 0);
enum rangeMax = 20_000; auto n2d = iota(1, rangeMax + 1).map!(n => properDivs(n).sum);
foreach (immutable n, immutable divSum; n2d.enumerate(1))
if (n < divSum && divSum <= rangeMax && n2d[divSum - 1] == n)
writefln("Amicable pair: %d and %d with proper divisors:\n %s\n %s",
n, divSum, properDivs(n), properDivs(divSum));
}</lang>
- Output:
Amicable pair: 220 and 284 with proper divisors:
[1, 2, 4, 5, 10, 11, 20, 22, 44, 55, 110]
[1, 2, 4, 71, 142]
Amicable pair: 1184 and 1210 with proper divisors:
[1, 2, 4, 8, 16, 32, 37, 74, 148, 296, 592]
[1, 2, 5, 10, 11, 22, 55, 110, 121, 242, 605]
Amicable pair: 2620 and 2924 with proper divisors:
[1, 2, 4, 5, 10, 20, 131, 262, 524, 655, 1310]
[1, 2, 4, 17, 34, 43, 68, 86, 172, 731, 1462]
Amicable pair: 5020 and 5564 with proper divisors:
[1, 2, 4, 5, 10, 20, 251, 502, 1004, 1255, 2510]
[1, 2, 4, 13, 26, 52, 107, 214, 428, 1391, 2782]
Amicable pair: 6232 and 6368 with proper divisors:
[1, 2, 4, 8, 19, 38, 41, 76, 82, 152, 164, 328, 779, 1558, 3116]
[1, 2, 4, 8, 16, 32, 199, 398, 796, 1592, 3184]
Amicable pair: 10744 and 10856 with proper divisors:
[1, 2, 4, 8, 17, 34, 68, 79, 136, 158, 316, 632, 1343, 2686, 5372]
[1, 2, 4, 8, 23, 46, 59, 92, 118, 184, 236, 472, 1357, 2714, 5428]
Amicable pair: 12285 and 14595 with proper divisors:
[1, 3, 5, 7, 9, 13, 15, 21, 27, 35, 39, 45, 63, 65, 91, 105, 117, 135, 189, 195, 273, 315, 351, 455, 585, 819, 945, 1365, 1755, 2457, 4095]
[1, 3, 5, 7, 15, 21, 35, 105, 139, 417, 695, 973, 2085, 2919, 4865]
Amicable pair: 17296 and 18416 with proper divisors:
[1, 2, 4, 8, 16, 23, 46, 47, 92, 94, 184, 188, 368, 376, 752, 1081, 2162, 4324, 8648]
[1, 2, 4, 8, 16, 1151, 2302, 4604, 9208]
J
Proper Divisor implementation:
<lang J>factors=: [: /:~@, */&>@{@((^ i.@>:)&.>/)@q:~&__ properDivisors=: factors -. -.&1</lang>
Amicable pairs:
<lang J> 1+0 20000 #:I.,(</~@i.@#*(*|:))(=/ +/@properDivisors@>) 1+i.20000
220 284 1184 1210 2620 2924 5020 5564 6232 6368
10744 10856 12285 14595 17296 18416</lang>
jq
<lang jq>def proper_divisors:
. as $n | range(1; 1 + ($n / 2 | floor)) as $i | if ($n % $i) == 0 then $i else empty end;
def addup(stream): reduce stream as $i (0; . + $i);
def task(n):
(reduce range(0; n+1) as $n ( []; . + [$n | addup(proper_divisors)] )) as $listing | range(1;n+1) as $j | range(1;$j) as $k | if $listing[$j] == $k and $listing[$k] == $j then "\($k) and \($j) are amicable" else empty end ;
task(20000)</lang>
- Output:
<lang sh>$ jq -c -n -f amicable_pairs.jq 220 and 284 are amicable 1184 and 1210 are amicable 2620 and 2924 are amicable 5020 and 5564 are amicable 6232 and 6368 are amicable 10744 and 10856 are amicable 12285 and 14595 are amicable 17296 and 18416 are amicable</lang>
Mathematica / Wolfram Language
<lang Mathematica>amicableQ[n_] :=
Module[{sum = Total[Most@Divisors@n]},
sum != n && n == Total[Most@Divisors@sum]]
Grid@Partition[Cases[Range[4, 20000], _?(amicableQ@# &)], 2]</lang>
- Output:
220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416
Pascal
extremly fast. up to 524'000'000 in 11.3 secs (linux 32-Bit fpc 2.6.4. i4330 3.5 Ghz )
n = p[1]^a[1]*p[2]^a[2]*...p[l]^a[l]
sum of divisors= s(n) = (p[1]^(a[1]+1) -1) / (p[1] -1) * ... * (p[l]^(a[l]+1) -1) / (p[l] -1) with
p[k]^(a[k]+1) -1) / (p[k] -1) = sum (i= [1..a[k]])(p[k]^i)
Using "Sieve of Erathosthenes"-style
<lang pascal>program AmicablePairs; //find amicable pairs in a limited region 2..MAX //optimized for freepascal 2.6.4 {$IFDEF FPC}
{$MODE DELPHI}
{$OPTIMIZATION ON}
{$OPTIMIZATION PeepHole}
{$OPTIMIZATION CSE}
{$OPTIMIZATION ASMCSE}
{$CODEALIGN loop=1,proc=8}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF}
uses
sysutils;
const // MAX = 200*1000*1000; {$IFDEF UNIX} MAX = 524*1000*1000;{$ELSE}MAX = 100*1000*1000;{$ENDIF}
type
tValue = LongWord; tpValue = ^tValue; tPower = array[0..31] of tValue;
var
power : tPower; PowerFac : tPower; SmallPF : array[0..17*17*17-1]of tValue; DivSumField : array[0..MAX] of tValue;
procedure Init; var
i : LongInt;
begin
DivSumField[0]:= 0; DivSumField[1]:= 0; For i := 2 to MAX do DivSumField[i]:= 1;
end;
procedure ProperDivs(n: LongWord); //Only for output, normally a factorication would do var
su,so : string; i,q : tValue;
begin
su:= '1';
so:= ;
i := 2;
while i*i <= n do
begin
q := n div i;
IF q*i -n = 0 then
begin
su:= su+','+IntToStr(i);
IF q <> i then
so:= ','+IntToStr(q)+so;
end;
inc(i);
end;
writeln(' [',su+so,']');
end;
procedure Check; var
i,s,n : LongInt;
begin
n := 0;
For i := 2 to MAX do
begin
s := DivSumField[i]-i;
//s never above MAX
IF (i > s) AND (s > 0) then
begin
IF DivSumField[s]= DivSumField[i] then
begin
inc(n);
writeln(s:10,i:10);
ProperDivs(s);ProperDivs(i);writeln;
end;
end;
end;
writeln(n,' amicable numbers upto ',MAX);
end;
Procedure CalcPotfactor(prim:tValue); //PowerFac[k] = (prim^(k+1)-1)/(prim-1) == Sum (i=1..k) prim^i var
k: tValue; Pot, //== prim^k PFac : Int64;
begin
Pot := prim;
PFac := 1;
For k := 0 to High(PowerFac) do
begin
PFac := PFac+Pot;
IF (POT > MAX) then
BREAK;
PowerFac[k] := PFac;
Pot := Pot*prim;
end;
end;
function NextNPotCnt(pot,prim: tValue):tValue; //return the first power[pot+?] <> 0 //power == n to base prim var
i : tValue;
begin
result := pot;
repeat
i := power[result];
Inc(i);
IF i < prim then
BREAK
else
begin
i := 0;
power[result] := 0;
inc(result);
end;
until false;
power[result] := i;
end;
function NextPotCnt(p: tValue):tValue; //return the first power <> 0 //power == n to base prim var
i : tValue;
begin
result := 0;
repeat
i := power[result];
Inc(i);
IF i < p then
BREAK
else
begin
i := 0;
power[result] := 0;
inc(result);
end;
until false;
power[result] := i;
end;
procedure InitPW(prim:tValue); begin
fillchar(power,SizeOf(power),#0); CalcPotfactor(prim);
end;
function FillSmallPF(prim:tValue):tValue; var
i : tValue; pSPF : tpValue;
begin
InitPW(prim); i :=1; result := 0; Repeat i := i*prim; inc(result); until i > High(SmallPF)+1; i := i DIV prim; dec(result); pSPF := @SmallPF[0]; for i := i-1 downto 0 do begin pSPF^ := PowerFac[NextPotCnt(prim)]; inc(pSPF); end; fillchar(power,SizeOf(power),#0);
// writeln((PtrUint(pSPF) - PtrUint(@SmallPF[0])) DIV SizeOF(tValue)); end;
function SieveReallySmall(prim: tValue):tValue; // precalculate values up to prim^n var
actNumber,n,i,pot,k : tValue;
begin
while prim*prim <= MAX do
begin
actNumber := prim;
n := FillSmallPF(prim);
IF n < 1 then
BREAK;
pot := 1;
For i := 1 to n do
Pot := Pot*Prim;
k := MAX DIV (prim*Pot);
dec(Pot);
while k >=1 do
begin
For i := 0 to Pot-1 do
begin
DivSumField[actNumber] := DivSumField[actNumber] *SmallPF[i];
inc(actNumber,prim);
end;
DivSumField[actNumber] := DivSumField[actNumber]*PowerFac[NextNPotCnt(n,prim)];
inc(actNumber,prim);
dec(k);
end;
For i := 1 to Pot do
begin
IF actNumber>= MAX then
BREAK;
DivSumField[actNumber] := DivSumField[actNumber] *SmallPF[i];
inc(actNumber,prim);
end;
//next prime
repeat
inc(prim);
until (DivSumField[prim] = 1);
end;
result := prim;
end;
function SieveSmall(prim: tValue):tValue; var
actNumber : tValue;
begin
//small primes
while prim*prim <= MAX do
begin
InitPW(prim);
//actNumber = actual number = n*prim
//power == n to base prim
actNumber := prim;
while actNumber < MAX do
begin
DivSumField[actNumber] := DivSumField[actNumber] *PowerFac[NextPotCnt(prim)];
inc(actNumber,prim);
end;
//next prime
repeat
inc(prim);
until (DivSumField[prim] = 1);
end;
result := prim;
end;
function SieveBig(prim: tValue):tValue; var
actNumber,p1 : tValue;
Begin
while prim <= MAX do
begin
actNumber := prim;
p1 := prim+1;
while actNumber <= MAX do
begin
DivSumField[actNumber] := DivSumField[actNumber]*p1;
inc(actNumber,prim);
end;
//next prime
repeat
inc(prim);
until (prim > MAX) OR (DivSumField[prim] = 1);
end;
result := prim;
end;
function Sieve(prim: tValue):tValue; //simple version var
actNumber : tValue;
begin
while prim <= MAX do
begin
InitPW(prim);
//actNumber = actual number = n*prim
//power == n to base prim
actNumber := prim;
while actNumber < MAX do
begin
DivSumField[actNumber] := DivSumField[actNumber] *PowerFac[NextPotCnt(prim)];
inc(actNumber,prim);
end;
//next prime
repeat
inc(prim);
until (DivSumField[prim] = 1);
end;
result := prim;
end;
var
T1,T0: TDatetime;
begin
T0:= time; Init; IF MAX < 1000*1000 then Sieve(2) else SieveBig(SieveSmall(SieveReallySmall(2))); T1:= time; Check;
writeln(FormatDateTime('HH:NN:SS.ZZZ' ,T1-T0));
{$IFNDEF UNIX}
readln;
{$ENDIF}
end.</lang> output
220 284
[1,2,4,5,10,11,20,22,44,55,110]
[1,2,4,71,142]
1184 1210
[1,2,4,8,16,32,37,74,148,296,592]
[1,2,5,10,11,22,55,110,121,242,605]
2620 2924
[1,2,4,5,10,20,131,262,524,655,1310]
[1,2,4,17,34,43,68,86,172,731,1462]
5020 5564
[1,2,4,5,10,20,251,502,1004,1255,2510]
[1,2,4,13,26,52,107,214,428,1391,2782]
6232 6368
[1,2,4,8,19,38,41,76,82,152,164,328,779,1558,3116]
[1,2,4,8,16,32,199,398,796,1592,3184]
10744 10856
[1,2,4,8,17,34,68,79,136,158,316,632,1343,2686,5372]
[1,2,4,8,23,46,59,92,118,184,236,472,1357,2714,5428]
12285 14595
[1,3,5,7,9,13,15,21,27,35,39,45,63,65,91,105,117,135,189,195,273,315,351,455,585,819,945,1365,1755,2457,4095]
[1,3,5,7,15,21,35,105,139,417,695,973,2085,2919,4865]
17296 18416
[1,2,4,8,16,23,46,47,92,94,184,188,368,376,752,1081,2162,4324,8648]
[1,2,4,8,16,1151,2302,4604,9208]
8 amicable numbers up to 20000
00:00:00.000
.....
{..
509379344 523679536
[1,2,4,8,16,23,46,92,184,347,368,694,1388,2776,3989,5552,7978,7981,15956,15962,31912,31924,63824,63848,91747,127696,183494,366988,733976,1384183,1467952,2768366,5536732,11073464,22146928,31836209,63672418,127344836,254689672]
[1,2,4,8,16,83,166,179,332,358,664,716,1328,1432,2203,2864,4406,8812,14857,17624,29714,35248,59428,118856,182849,237712,365698,394337,731396,788674,1462792,1577348,2925584,3154696,6309392,32729971,65459942,130919884,261839768]
442 amicable numbers upto 524000000
00:00:11.213
real 0m15.222s
user 0m14.977s
sys 0m0.223s
}
}
PL/I
<lang pli>*process source xref;
ami: Proc Options(main); p9a=time(); Dcl (p9a,p9b,p9c) Pic'(9)9'; Dcl sumpd(20000) Bin Fixed(31); Dcl pd(300) Bin Fixed(31); Dcl npd Bin Fixed(31); Dcl (x,y) Bin Fixed(31);
Do x=1 To 20000;
Call proper_divisors(x,pd,npd);
sumpd(x)=sum(pd,npd);
End;
p9b=time();
Put Edit('sum(pd) computed in',(p9b-p9a)/1000,' seconds elapsed')
(Skip,col(7),a,f(6,3),a);
Do x=1 To 20000;
Do y=x+1 To 20000;
If y=sumpd(x) &
x=sumpd(y) Then
Put Edit(x,y,' found after ',elapsed(),' seconds')
(Skip,2(f(6)),a,f(6,3),a);
End;
End;
Put Edit(elapsed(),' seconds total search time')(Skip,f(6,3),a);
proper_divisors: Proc(n,pd,npd);
Dcl (n,pd(300),npd) Bin Fixed(31);
Dcl (d,delta) Bin Fixed(31);
npd=0;
If n>1 Then Do;
If mod(n,2)=1 Then /* odd number */
delta=2;
Else /* even number */
delta=1;
Do d=1 To n/2 By delta;
If mod(n,d)=0 Then Do;
npd+=1;
pd(npd)=d;
End;
End;
End;
End;
sum: Proc(pd,npd) Returns(Bin Fixed(31)); Dcl (pd(300),npd) Bin Fixed(31); Dcl sum Bin Fixed(31) Init(0); Dcl i Bin Fixed(31); Do i=1 To npd; sum+=pd(i); End; Return(sum); End;
elapsed: Proc Returns(Dec Fixed(6,3)); p9c=time(); Return((p9c-p9b)/1000); End; End;</lang>
- Output:
sum(pd) computed in 0.510 seconds elapsed 220 284 found after 0.010 seconds 1184 1210 found after 0.060 seconds 2620 2924 found after 0.110 seconds 5020 5564 found after 0.210 seconds 6232 6368 found after 0.260 seconds 10744 10856 found after 2.110 seconds 12285 14595 found after 2.150 seconds 17296 18416 found after 2.240 seconds 2.250 seconds total search time
Python
Importing Proper divisors from prime factors: <lang python>from proper_divisors import proper_divs
def amicable(rangemax=20000):
n2divsum = {n: sum(proper_divs(n)) for n in range(1, rangemax + 1)}
for num, divsum in n2divsum.items():
if num < divsum and divsum <= rangemax and n2divsum[divsum] == num:
yield num, divsum
if __name__ == '__main__':
for num, divsum in amicable():
print('Amicable pair: %i and %i With proper divisors:\n %r\n %r'
% (num, divsum, sorted(proper_divs(num)), sorted(proper_divs(divsum))))</lang>
- Output:
Amicable pair: 220 and 284 With proper divisors:
[1, 2, 4, 5, 10, 11, 20, 22, 44, 55, 110]
[1, 2, 4, 71, 142]
Amicable pair: 1184 and 1210 With proper divisors:
[1, 2, 4, 8, 16, 32, 37, 74, 148, 296, 592]
[1, 2, 5, 10, 11, 22, 55, 110, 121, 242, 605]
Amicable pair: 2620 and 2924 With proper divisors:
[1, 2, 4, 5, 10, 20, 131, 262, 524, 655, 1310]
[1, 2, 4, 17, 34, 43, 68, 86, 172, 731, 1462]
Amicable pair: 5020 and 5564 With proper divisors:
[1, 2, 4, 5, 10, 20, 251, 502, 1004, 1255, 2510]
[1, 2, 4, 13, 26, 52, 107, 214, 428, 1391, 2782]
Amicable pair: 6232 and 6368 With proper divisors:
[1, 2, 4, 8, 19, 38, 41, 76, 82, 152, 164, 328, 779, 1558, 3116]
[1, 2, 4, 8, 16, 32, 199, 398, 796, 1592, 3184]
Amicable pair: 10744 and 10856 With proper divisors:
[1, 2, 4, 8, 17, 34, 68, 79, 136, 158, 316, 632, 1343, 2686, 5372]
[1, 2, 4, 8, 23, 46, 59, 92, 118, 184, 236, 472, 1357, 2714, 5428]
Amicable pair: 12285 and 14595 With proper divisors:
[1, 3, 5, 7, 9, 13, 15, 21, 27, 35, 39, 45, 63, 65, 91, 105, 117, 135, 189, 195, 273, 315, 351, 455, 585, 819, 945, 1365, 1755, 2457, 4095]
[1, 3, 5, 7, 15, 21, 35, 105, 139, 417, 695, 973, 2085, 2919, 4865]
Amicable pair: 17296 and 18416 With proper divisors:
[1, 2, 4, 8, 16, 23, 46, 47, 92, 94, 184, 188, 368, 376, 752, 1081, 2162, 4324, 8648]
[1, 2, 4, 8, 16, 1151, 2302, 4604, 9208]
Racket
With Proper_divisors#Racket in place: <lang racket>#lang racket (require "proper-divisors.rkt") (define SCOPE 20000)
(define P
(let ((P-v (vector)))
(λ (n)
(set! P-v (fold-divisors P-v n 0 +))
(vector-ref P-v n))))
- returns #f if not an amicable number, amicable pairing otherwise
(define (amicable? n)
(define m (P n))
(define m-sod (P m))
(and (= m-sod n)
(< m n) ; each pair exactly once, also eliminates perfect numbers
m))
(void (amicable? SCOPE)) ; prime the memoisation
(for* ((n (in-range 1 (add1 SCOPE)))
(m (in-value (amicable? n)))
#:when m)
(printf #<<EOS
amicable pair: ~a, ~a
~a: divisors: ~a ~a: divisors: ~a
EOS
n m n (proper-divisors n) m (proper-divisors m)))
</lang>
- Output:
amicable pair: 284, 220 284: divisors: (1 2 4 71 142) 220: divisors: (1 2 4 5 10 11 20 22 44 55 110) amicable pair: 1210, 1184 1210: divisors: (1 2 5 10 11 22 55 110 121 242 605) 1184: divisors: (1 2 4 8 16 32 37 74 148 296 592) amicable pair: 2924, 2620 2924: divisors: (1 2 4 17 34 43 68 86 172 731 1462) 2620: divisors: (1 2 4 5 10 20 131 262 524 655 1310) amicable pair: 5564, 5020 5564: divisors: (1 2 4 13 26 52 107 214 428 1391 2782) 5020: divisors: (1 2 4 5 10 20 251 502 1004 1255 2510) amicable pair: 6368, 6232 6368: divisors: (1 2 4 8 16 32 199 398 796 1592 3184) 6232: divisors: (1 2 4 8 19 38 41 76 82 152 164 328 779 1558 3116) amicable pair: 10856, 10744 10856: divisors: (1 2 4 8 23 46 59 92 118 184 236 472 1357 2714 5428) 10744: divisors: (1 2 4 8 17 34 68 79 136 158 316 632 1343 2686 5372) amicable pair: 14595, 12285 14595: divisors: (1 3 5 7 15 21 35 105 139 417 695 973 2085 2919 4865) 12285: divisors: (1 3 5 7 9 13 15 21 27 35 39 45 63 65 91 105 117 135 189 195 273 315 351 455 585 819 945 1365 1755 2457 4095) amicable pair: 18416, 17296 18416: divisors: (1 2 4 8 16 1151 2302 4604 9208) 17296: divisors: (1 2 4 8 16 23 46 47 92 94 184 188 368 376 752 1081 2162 4324 8648)
REXX
<lang rexx>Call time 'R' Do x=1 To 20000
pd=proper_divisors(x) sumpd.x=sum(pd) End
Say 'sum(pd) computed in' time('E') 'seconds' Call time 'R' Do x=1 To 20000
/* If x//1000=0 Then Say x time() */
Do y=x+1 To 20000
If y=sumpd.x &,
x=sumpd.y Then
Say x y 'found after' time('E') 'seconds'
End
End
Say time('E') 'seconds total search time' Exit
proper_divisors: Procedure Parse Arg n Pd= If n=1 Then Return If n//2=1 Then /* odd number */
delta=2
Else /* even number */
delta=1
Do d=1 To n%2 By delta
If n//d=0 Then pd=pd d End
Return space(pd)
sum: Procedure Parse Arg list sum=0 Do i=1 To words(list)
sum=sum+word(list,i) End
Return sum</lang>
- Output:
sum(pd) computed in 48.502000 seconds 220 284 found after 3.775000 seconds 1184 1210 found after 21.611000 seconds 2620 2924 found after 46.817000 seconds 5020 5564 found after 84.296000 seconds 6232 6368 found after 100.918000 seconds 10744 10856 found after 150.126000 seconds 12285 14595 found after 162.124000 seconds 17296 18416 found after 185.600000 seconds 188.836000 seconds total search time
Ruby
With proper_divisors#Ruby in place: <lang ruby>h = {} (1..20_000).each{|n| h[n] = n.proper_divisors.inject(:+)} h.select{|k,v| h[v] == k && k < v}.each do |key,val| # k<v filters out doubles and perfects
puts "#{key} and #{val}"
end </lang>
- Output:
220 and 284 1184 and 1210 2620 and 2924 5020 and 5564 6232 and 6368 10744 and 10856 12285 and 14595 17296 and 18416
Tcl
<lang tcl>proc properDivisors {n} {
if {$n == 1} return
set divs 1
set sum 1
for {set i 2} {$i*$i <= $n} {incr i} {
if {!($n % $i)} { lappend divs $i incr sum $i if {$i*$i < $n} { lappend divs [set d [expr {$n / $i}]] incr sum $d } }
} return [list $sum $divs]
}
proc amicablePairs {limit} {
set result {}
set sums [set divs {{}}]
for {set n 1} {$n < $limit} {incr n} {
lassign [properDivisors $n] sum d lappend sums $sum lappend divs [lsort -integer $d]
}
for {set n 1} {$n < $limit} {incr n} {
set nsum [lindex $sums $n] for {set m 1} {$m < $n} {incr m} { if {$n==[lindex $sums $m] && $m==$nsum} { lappend result $m $n [lindex $divs $m] [lindex $divs $n] } }
} return $result
}
foreach {m n md nd} [amicablePairs 20000] {
puts "$m and $n are an amicable pair with these proper divisors" puts "\t$m : $md" puts "\t$n : $nd"
}</lang>
- Output:
220 and 284 are an amicable pair with these proper divisors 220 : 1 2 4 5 10 11 20 22 44 55 110 284 : 1 2 4 71 142 1184 and 1210 are an amicable pair with these proper divisors 1184 : 1 2 4 8 16 32 37 74 148 296 592 1210 : 1 2 5 10 11 22 55 110 121 242 605 2620 and 2924 are an amicable pair with these proper divisors 2620 : 1 2 4 5 10 20 131 262 524 655 1310 2924 : 1 2 4 17 34 43 68 86 172 731 1462 5020 and 5564 are an amicable pair with these proper divisors 5020 : 1 2 4 5 10 20 251 502 1004 1255 2510 5564 : 1 2 4 13 26 52 107 214 428 1391 2782 6232 and 6368 are an amicable pair with these proper divisors 6232 : 1 2 4 8 19 38 41 76 82 152 164 328 779 1558 3116 6368 : 1 2 4 8 16 32 199 398 796 1592 3184 10744 and 10856 are an amicable pair with these proper divisors 10744 : 1 2 4 8 17 34 68 79 136 158 316 632 1343 2686 5372 10856 : 1 2 4 8 23 46 59 92 118 184 236 472 1357 2714 5428 12285 and 14595 are an amicable pair with these proper divisors 12285 : 1 3 5 7 9 13 15 21 27 35 39 45 63 65 91 105 117 135 189 195 273 315 351 455 585 819 945 1365 1755 2457 4095 14595 : 1 3 5 7 15 21 35 105 139 417 695 973 2085 2919 4865 17296 and 18416 are an amicable pair with these proper divisors 17296 : 1 2 4 8 16 23 46 47 92 94 184 188 368 376 752 1081 2162 4324 8648 18416 : 1 2 4 8 16 1151 2302 4604 9208
VbScript
Not at all optimal. :-( <lang VbScript>start = Now Set nlookup = CreateObject("Scripting.Dictionary") Set uniquepair = CreateObject("Scripting.Dictionary")
For i = 1 To 20000 sum = 0 For n = 1 To 20000 If n < i Then If i Mod n = 0 Then sum = sum + n End If End If Next nlookup.Add i,sum Next
For j = 1 To 20000 sum = 0 For m = 1 To 20000 If m < j Then If j Mod m = 0 Then sum = sum + m End If End If Next If nlookup.Exists(sum) And nlookup.Item(sum) = j And j <> sum _ And uniquepair.Exists(sum) = False Then uniquepair.Add j,sum End If Next
For Each key In uniquepair.Keys WScript.Echo key & ":" & uniquepair.Item(key) Next
WScript.Echo "Execution Time: " & DateDiff("s",Start,Now) & " seconds"</lang>
- Output:
220:284 1184:1210 2620:2924 5020:5564 6232:6368 10744:10856 12285:14595 17296:18416 Execution Time: 162 seconds
zkl
Slooooow <lang zkl>fcn properDivs(n){ [1.. (n + 1)/2 + 1].filter('wrap(x){ n%x==0 and n!=x }) } const N=20000; sums:=[1..N].pump(T(-1),fcn(n){ properDivs(n).sum(0) }); [0..].zip(sums).filter('wrap([(n,s)]){ (n<s<=N) and sums[s]==n }).println();</lang>
- Output:
L(L(220,284),L(1184,1210),L(2620,2924),L(5020,5564),L(6232,6368),L(10744,10856),L(12285,14595),L(17296,18416))