Amicable pairs: Difference between revisions
m (Corrected level of header) |
m (→version 2: corrected a comment in the REXX (version 2) section header.) |
||
Line 1,107: | Line 1,107: | ||
===version 2=== |
===version 2=== |
||
This REXX version allows the specification of the upper limit (for the searching of amicable pairs), some optimization was |
This REXX version allows the specification of the upper limit (for the searching of amicable pairs), some optimization was |
||
<br>incorporated into the ''proper |
<br>incorporated into the ''proper divisors'' subroutine as well as the search itself, the subroutine is a modified version of the |
||
<br>subroutine taken from the REXX language entry for Rosetta code task ''integer factors''. |
<br>subroutine taken from the REXX language entry for Rosetta code task ''integer factors''. |
||
Revision as of 21:54, 27 January 2015
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
C
The program will overflow and error in all sorts of ways when given a commandline argument >= UINT_MAX/2 (generally 2^31) <lang c>#include <stdio.h>
- include <stdlib.h>
typedef unsigned int uint;
int main(int argc, char **argv) { uint top = atoi(argv[1]); uint *divsum = malloc((top + 1) * sizeof(*divsum)); uint pows[32] = {1, 0};
for (uint i = 0; i <= top; i++) divsum[i] = 1;
// sieve for (uint p = 2; p <= top; p++) { if (divsum[p] > 1) continue; // p not prime
uint x; // highest power of p we need
// checking x <= top/y instead of x*y <= top to avoid overflow for (x = 1; pows[x - 1] <= top/p; x++) pows[x] = p*pows[x - 1];
for (uint n = p; n <= top; n += p) { uint s; for (uint i = s = 1; i < x && !(n%pows[i]); s += pows[i++]); divsum[n] *= s; } }
// subtract number itself from divisor sum ('proper') for (uint i = 0; i <= top; i++) divsum[i] -= i;
for (uint a = 1; a <= top; a++) { uint b = divsum[a]; if (b > a && b <= top && divsum[b] == a) printf("%u %u\n", a, b); }
return 0; }</lang>
- Output:
% ./a.out 20000 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]
Haskell
<lang Haskell>divisors :: (Integral a) => a -> [a] divisors n = filter ((0 ==) . (n `mod`)) [1 .. (n `div` 2)]
main :: IO () main = do
let range = [1 .. 20000 :: Int] divs = zip range $ map (sum . divisors) range pairs = [(n, m) | (n, nd) <- divs, (m, md) <- divs, n < m, nd == m, md == n] print pairs</lang>
- Output:
[(220,284),(1184,1210),(2620,2924),(5020,5564),(6232,6368),(10744,10856),(12285,14595),(17296,18416)]
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># unordered def proper_divisors:
. as $n | if $n > 1 then 1, (sqrt|floor as $s | range(2; $s+1) as $i | if ($n % $i) == 0 then $i, (if $i * $i == $n then empty else ($n / $i) end)
else empty end)
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
PARI/GP
<lang parigp>for(x=1,20000,my(y=sigma(x)-x); if(y>x && x == sigma(y)-y,print(x" "y)))</lang>
- Output:
220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416
Pascal
a "normal" Version. Nearly fast as perl using nTheory. <lang pascal>program AmicablePairs; {$IFDEF FPC}
{$MODE DELPHI} {$H+}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF} uses
sysutils;
const
MAX = 20000;
//MAX = 20*1000*1000; type
tValue = LongWord; tpValue = ^tValue; tPower = array[0..31] of tValue; tIndex = record idxI, idxS : Uint64; end;
var
Indices : array[0..511] of tIndex; //primes up to 65536 enough until 2^32 primes : array[0..6542] of tValue;
procedure InitPrimes; // sieve of erathosthenes without multiples of 2 type
tSieve = array[0..(65536-1) div 2] of char;
var
ESieve : ^tSieve; idx,i,j,p : LongINt;
Begin
new(ESieve); fillchar(ESieve^[0],SizeOF(tSieve),#1); primes[0] := 2; idx := 1;
//sieving j := 1; p := 2*j+1; repeat if Esieve^[j] = #1 then begin i := (2*j+2)*j;// i := (sqr(p) -1) div 2; if i > High(tSieve) then BREAK; repeat ESIeve^[i] := #0; inc(i,p); until i > High(tSieve); end; inc(j); inc(p,2); until j >High(tSieve);
//collecting For i := 1 to High(tSieve) do IF Esieve^[i] = #1 then Begin primes[idx] := 2*i+1; inc(idx); IF idx>High(primes) then BREAK; end; dispose(Esieve);
end;
procedure Su_append(n,factor:tValue;var su:string); var
q,p : tValue;
begin
p := 0; repeat q := n div factor; IF q*factor<>n then Break; inc(p); n := q; until false; IF p > 0 then IF p= 1 then su:= su+IntToStr(factor)+'*' else su:= su+IntToStr(factor)+'^'+IntToStr(p)+'*';
end;
procedure ProperDivs(n: Uint64); //output of prime factorization var
su : string; primNo : tValue; p:tValue;
begin
str(n:8,su); su:= su +' ['; primNo := 0; p := primes[0]; repeat Su_Append(n,p,su); inc(primNo); p := primes[primNo]; until (p=0) OR (p*p >= n); p := n; Su_Append(n,p,su); su[length(su)] := ']'; writeln(su);
end;
procedure AmPairOutput(cnt:tValue); var
i : tValue; r_max,r_min,r : double;
begin
r_max := 1.0; r_min := 16.0; For i := 0 to cnt-1 do with Indices[i] do begin r := IdxS/IDxI; writeln(i+1:4,IdxI:16,IDxS:16,' ratio ',r:10:7); IF r < 1 then begin writeln(i); readln; halt; end; if r_max < r then r_max := r else if r_min > r then r_min := r; IF cnt < 20 then begin ProperDivs(IdxI); ProperDivs(IdxS); end; end; writeln(' min ratio ',r_min:12:10); writeln(' max ratio ',r_max:12:10);
end;
procedure SumOFProperDiv(n: tValue;var SumOfProperDivs:tValue); // calculated by prime factorization var
i,q, primNo, Prime,pot : tValue; SumOfDivs: tValue;
begin
i := N; SumOfDivs := 1; primNo := 0; Prime := Primes[0]; q := i DIV Prime; repeat if q*Prime = i then Begin pot := 1; repeat i := q; q := i div Prime; Pot := Pot * Prime+1; until q*Prime <> i; SumOfDivs := SumOfDivs * pot; end; Inc(primNo); Prime := Primes[primNo]; q := i DIV Prime;
{check if i already prime} if Prime > q then begin prime := i; q := 1; end; until i = 1; SumOfProperDivs := SumOfDivs - N;
end;
function Check:tValue; const
//going backwards DIV23 : array[0..5] of byte = //== 5,4,3,2,1,0 (1,0,0,0,1,0);
var
i,s,k,n : tValue; idx : nativeInt;
begin
n := 0; idx := 3; For i := 2 to MAX do begin //must be divisble by 2 or 3 ( n < High(tValue) < 1e14 ) IF DIV23[idx] = 0 then begin SumOFProperDiv(i,s); //only 24.7...% IF s>i then Begin SumOFProperDiv(s,k); IF k = i then begin With indices[n] do begin idxI := i; idxS := s; end; inc(n); end; end; end; dec(idx); IF idx < 0 then idx := high(DIV23); end; result := n;
end;
var
T2,T1: TDatetime; APcnt: tValue;
begin
InitPrimes; T1:= time; APCnt:= Check; T2:= time; AmPairOutput(APCnt); writeln('Time to find amicable pairs ',FormatDateTime('HH:NN:SS.ZZZ' ,T2-T1)); {$IFNDEF UNIX} readln;{$ENDIF}
end.</lang> Output
1 220 284 ratio 1.2909091 220 [2^2*5*11*220] 284 [2^2*284] 2 1184 1210 ratio 1.0219595 1184 [2^5*1184] 1210 [2*5*11^2*1210] 3 2620 2924 ratio 1.1160305 2620 [2^2*5*2620] 2924 [2^2*17*43*2924] 4 5020 5564 ratio 1.1083665 5020 [2^2*5*5020] 5564 [2^2*13*5564] 5 6232 6368 ratio 1.0218228 6232 [2^3*19*41*6232] 6368 [2^5*6368] 6 10744 10856 ratio 1.0104244 10744 [2^3*17*79*10744] 10856 [2^3*23*59*10856] 7 12285 14595 ratio 1.1880342 12285 [3^3*5*7*13*12285] 14595 [3*5*7*14595] 8 17296 18416 ratio 1.0647549 17296 [2^4*23*47*17296] 18416 [2^4*18416]
Alternative
about 25-times faster. This will not output the amicable number unless both! numbers are under the given limit.
So there will be differences to "Table of n, a(n) for n=1..39374" https://oeis.org/A002025/b002025.txt Up to 524'000'000 the pairs found are only correct by number up to no. 437 460122410 and only 442 out of 455 are found, because some pairs exceed the limit. The limits of the ratio between the numbers of the amicable pair up to 1E14 are, based on b002025.txt:
No. lower upper 31447 52326552030976 52326637800704 ratio 1.0000016 52326552030976 [2^8*563*6079*59723] 52326637800704 [2^8*797*1439*178223] 38336 92371445691525 154378742017851 ratio 1.6712821 92371445691525 [3^2*5^2*7^2*11*13^2*23*29^2*233] 154378742017851 [3^2*13^2*53*337*5682671]
The distance check is being corrected, the lower number is now not limited.
The used method is not useful for very high limits.
n = p[1]^a[1]*p[2]^a[2]*...p[l]^a[l]
sum of divisors(n) = 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 beware that >both< numbers must be smaller than MAX there are 455 amicable pairs up to 524*1000*1000 correct up to
- 437 460122410
} //optimized for freepascal 2.6.4 32-Bit {$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,peephole,cse,asmcse,regvar} {$CODEALIGN loop=1,proc=8}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF}
uses
sysutils;
const //MAX = 20000; {$IFDEF UNIX} MAX = 524*1000*1000;{$ELSE}MAX = 499*1000*1000;{$ENDIF} type
tValue = LongWord; tpValue = ^tValue; tPower = array[0..31] of tValue; tIndex = record idxI, idxS : tValue; end; tdpa = array[0..2] of LongWord;
var
power : tPower; PowerFac : tPower; DivSumField : array[0..MAX] of tValue; Indices : array[0..511] of tIndex; DpaCnt : tdpa;
procedure Init; var
i : LongInt;
begin
DivSumField[0]:= 0; For i := 1 to MAX do DivSumField[i]:= 1;
end;
procedure ProperDivs(n: tValue); //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 AmPairOutput(cnt:tValue); var
i : tValue; r : double;
begin
r := 1.0; For i := 0 to cnt-1 do with Indices[i] do begin writeln(i+1:4,IdxI:12,IDxS:12,' ratio ',IdxS/IDxI:10:7); if r < IdxS/IDxI then r := IdxS/IDxI; IF cnt < 20 then begin ProperDivs(IdxI); ProperDivs(IdxS); end; end; writeln(' max ratio ',r:10:4);
end;
function Check:tValue; var
i,s,n : tValue;
begin
fillchar(DpaCnt,SizeOf(dpaCnt),#0); n := 0; For i := 1 to MAX do begin //s = sum of proper divs (I) == sum of divs (I) - I s := DivSumField[i]-i; IF (s <=MAX) AND (s>i) then begin IF DivSumField[s]-s = i then begin With indices[n] do begin idxI := i; idxS := s; end; inc(n); end; end; inc(DpaCnt[Ord(s>=i)-Ord(s<=i)+1]); end; result := n;
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;
procedure InitPW(prim:tValue); begin
fillchar(power,SizeOf(power),#0); CalcPotfactor(prim);
end;
function NextPotCnt(p: tValue):tValue;inline; //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;
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
T2,T1,T0: TDatetime; APcnt: tValue;
begin
T0:= time; Init; Sieve(2); T1:= time; APCnt := Check; T2:= time; AmPairOutput(APCnt); writeln(DpaCnt[0]:10,' deficient'); writeln(DpaCnt[1]:10,' perfect'); writeln(DpaCnt[2]:10,' abundant'); writeln(DpaCnt[2]/DpaCnt[0]:14:10,' ratio abundant/deficient '); writeln('Time to calc sum of divs ',FormatDateTime('HH:NN:SS.ZZZ' ,T1-T0)); writeln('Time to find amicable pairs ',FormatDateTime('HH:NN:SS.ZZZ' ,T2-T1)); {$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 {Win7 nearly nothing else running. MAX = 499*1000*1000 435 460122410 484817110 ratio 1.0536698 436 466389344 472453792 ratio 1.0130030 max ratio 1.3537 375440837 deficient 5 perfect 123559158 abundant 0.3291042045 ratio abundant/deficient Time to calc sum of divs 00:00:17.818 Time to find amicable pairs 00:00:04.493
Perl
Not particularly clever, but instant for this example, and does up to 20 million in 11 seconds.
<lang perl>use ntheory qw/divisor_sum/; for my $x (1..20000) {
my $y = divisor_sum($x)-$x; say "$x $y" if $y > $x && $x == divisor_sum($y)-$y;
}</lang>
- Output:
220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416
Perl 6
<lang perl6>sub propdivsum (\x) {
[+] x > 1, gather for 2 .. x.sqrt.floor -> \d { my \y = x div d; if y * d == x { take d; take y unless y == d } }
}
for 1..20000 -> $i {
my $j = propdivsum($i); say "$i $j" if $j > $i and $i == propdivsum($j);
}</lang>
- Output:
220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416
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
version 1
<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
version 2
This REXX version allows the specification of the upper limit (for the searching of amicable pairs), some optimization was
incorporated into the proper divisors subroutine as well as the search itself, the subroutine is a modified version of the
subroutine taken from the REXX language entry for Rosetta code task integer factors.
The generation/summation is about forty times faster, searching is about four times faster. <lang rexx>/*REXX program finds/displays all amicable pairs up to a given number.*/ parse arg H .; if H== then H=20000 /*get optional arg (high limit).*/ w=length(H) ; H.=H || . /*for columnar aligned output. */ @.=0
do k=1 for H; _=Pdivs(k); #=words(_) /*gen proper divs.*/ do i=1 for #; @.k=@.k + word(_,i) /*gen Pdivs sums. */ end /*i*/ /* [↑] sum the proper divisors.*/ end /*k*/ /* [↑] process a range of ints.*/
- =0 /*number of amicable pairs found.*/
do m=220 for H-220+1 /*start search at lowest number. */ do n=m+1 for H-m if m==@.n then if n==@.m then do; #=#+1 /*bump the counter.*/ say right(m,w) ' and ' right(n,w) " are amicable pairs." end end /*p*/ end /*n*/ /*DO loop FORs: faster than TOs.*/
say say # 'amicable pairs found up to' H. /*display count of amicable pairs*/ exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────PDIVS subroutine────────────────────*/ Pdivs: procedure; parse arg x,b; odd=x//2 /* [↑] modified for amicable*/ a=1 /* [↓] use only EVEN|ODD integers*/
do j=2+odd by 1+odd while j*j<x /*divide by all integers up to √x*/ if x//j==0 then do; a=a j; b=x%j b; end /*add divs to α&ß lists if ÷*/ end /*j*/ /* [↑] % is REXX integer divide*/ /* [↓] adjust for square. _ */
if j*j==x then return a j b /*Was X a square? If so, add √x.*/
return a b /*return divisors (both lists). */</lang>
output when the default input is used:
220 and 284 are amicable pairs. 1184 and 1210 are amicable pairs. 2620 and 2924 are amicable pairs. 5020 and 5564 are amicable pairs. 6232 and 6368 are amicable pairs. 10744 and 10856 are amicable pairs. 12285 and 14595 are amicable pairs. 17296 and 18416 are amicable pairs. 8 amicable pairs found up to 20000.
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
Swift
<lang Swift>import func Darwin.sqrt
func sqrt(x:Int) -> Int { return Int(sqrt(Double(x))) }
func properDivs(n: Int) -> [Int] {
if n == 1 { return [] } var result = [Int]() for div in filter (1...sqrt(n), { n % $0 == 0 }) { result.append(div)
if n/div != div && n/div != n { result.append(n/div) } } return sorted(result)
}
func sumDivs(n:Int) -> Int {
struct Cache { static var sum = [Int:Int]() } if let sum = Cache.sum[n] { return sum } let sum = properDivs(n).reduce(0) { $0 + $1 } Cache.sum[n] = sum return sum
}
func amicable(n:Int, m:Int) -> Bool {
if n == m { return false } if sumDivs(n) != m || sumDivs(m) != n { return false } return true
}
var pairs = [(Int, Int)]()
for n in 1 ..< 20_000 {
for m in n+1 ... 20_000 { if amicable(n, m) { pairs.append(n, m) println("\(n, m)") } }
}</lang>
Alternative
about 800 times faster.<lang Swift>import func Darwin.sqrt
func sqrt(x:Int) -> Int { return Int(sqrt(Double(x))) }
func sigma(n: Int) -> Int {
if n == 1 { return 0 } // definition of aliquot sum var result = 1 let root = sqrt(n) for var div = 2; div <= root; ++div { if n % div == 0 { result += div + n/div } } if root*root == n { result -= root } return (result)
}
func amicables (upTo: Int) -> () {
var aliquot = Array(count: upTo+1, repeatedValue: 0)
for i in 1 ... upTo { // fill lookup array aliquot[i] = sigma(i) } for i in 1 ... upTo { let a = aliquot[i] if a > upTo {continue} //first candidate out-of-bounds
if a == i {continue} //skip perfect numbers let b = aliquot[a] if b > upTo {continue} //second candidate out-of-bounds if i == b { println("\(i, a)") aliquot[a] = upTo+1 //prevent second display of pair } }
}
amicables(20_000)</lang>
- Output:
(220, 284) (1184, 1210) (2620, 2924) (5020, 5564) (6232, 6368) (10744, 10856) (12285, 14595) (17296, 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))