Amicable pairs: Difference between revisions

From Rosetta Code
Content added Content deleted
m (→‎Direct approach: compiler switch for freepascal, put filehandling after summation , stop showing overflow)
Line 78: Line 78:
</pre>
</pre>


=={{header|Awk}}==
=={{header|AWK}}==
<lang awk>
<lang awk>
#!/bin/awk -f
#!/bin/awk -f
Line 111: Line 111:
17296 18416
17296 18416
</pre>
</pre>

=={{header|C}}==
=={{header|C}}==
Remark:
Remark:

Revision as of 23:33, 11 April 2015

Task
Amicable pairs
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.

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

AWK

<lang awk>

  1. !/bin/awk -f

function sumprop(num, i,sum) { sum=0 for ( i=1; i < num; i++) { 

   if (num % i == 0 )
   { 
   sum = sum + i
   }
   }

return sum } BEGIN{ limit=20000 for (n=1; n < limit+1; n++) 

   {
   m=sumprop(n)
   if (n == sumprop(m) && n < m) print n,m
   }

}</lang>

Output:
# ./amicable 
220 284
1184 1210
2620 2924
5020 5564
6232 6368
10744 10856
12285 14595
17296 18416

C

Remark: Look at Pascal Alternative [[1]].You are using the same principle, so here too both numbers of the pair must be < top.

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>

  1. 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
 // only sieve within lower half , the modification starts at 2*p
 for (uint p = 2; p+p <= top; p++) {
   if (divsum[p] > 1) {
     divsum[p] -= p;// subtract number itself from divisor sum ('proper')
     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];

   //counter where n is not a*p with a = ?*p, useful for most p.
   //think of p>31 seldom divisions or p>sqrt(top) than no division is needed
   //n = 2*p, so the prime itself is left unchanged => k=p-1
   uint k= p-1;
   for (uint n = p+p; n <= top; n += p) {
     uint s=1+pows[1];
     k--;
     // search the right power only if needed
     if ( k==0) {
       for (uint i = 2; i < x && !(n%pows[i]); s += pows[i++]);
       k = p; }
     divsum[n] *= s;
   }
 }

 //now correct the upper half
 for (uint p = (top >> 1)+1; p <= top; p++) {
   if (divsum[p] > 1){
     divsum[p] -= p;}
 }

 uint cnt = 0;
 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);
     cnt++;}
 }
 printf("\nTop %u count : %u\n",top,cnt);
 return 0;

}</lang>

Output:
% ./a.out 20000
220 284
1184 1210
2620 2924
5020 5564
6232 6368
10744 10856
12285 14595
17296 18416

Top 20000 count : 8


% ./a.out 524000000
..
475838415 514823985
491373104 511419856
509379344 523679536

Top 524000000 count : 442

real  0m16.285s
user  0m16.156s

D

Translation of: Python

<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>

Julia

Given factor, it is not necessary to calculate the individual divisors to compute their sum. See Abundant, deficient and perfect number classifications for the details.

Functions <lang Julia> function pcontrib(p::Int64, a::Int64) 

   n = one(p)
   pcon = one(p)
   for i in 1:a
       n *= p
       pcon += n
   end
   return pcon

end

function divisorsum(n::Int64) 

   dsum = one(n)
   for (p, a) in factor(n)
       dsum *= pcontrib(p, a)
   end
   dsum -= n

end </lang> pcontrib is a good candidate for Memoization should the performance of divisorsum become an issue.

Main

It is safe to exclude primes from consideration; their proper divisor sum is always 1. Also, this code uses a minor trick to ensure that none of the numbers identified are above the limit. All numbers in the range are checked for an amicable partner, but the pair is cataloged only when the greater member is reached. <lang Julia> const L = 2*10^4 acnt = 0

println("Amicable pairs not greater than ", L)

for i in 2:L 

   !isprime(i) || continue
   j = divisorsum(i)
   j < i && divisorsum(j) == i || continue
   acnt += 1
   println(@sprintf("%4d", acnt), " => ", j, ", ", i)

end </lang>

Output:
Amicable pairs not greater than 20000
   1 => 220, 284
   2 => 1184, 1210
   3 => 2620, 2924
   4 => 5020, 5564
   5 => 6232, 6368
   6 => 10744, 10856
   7 => 12285, 14595
   8 => 17296, 18416

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

Oforth

Using properDivs implementation tasks without optimization (calculating proper divisors sum without returning a list for instance) :

<lang Oforth>Integer method: properDivs { seq(self 2 / ) filter(#[ self swap rem 0 == ]) }

func: amicables { | i j | 

  ListBuffer new
  20000 loop: i [
     i properDivs sum dup ->j i <= ifTrue: [ continue ]
     j properDivs sum i <> ifTrue: [ continue ]
     [ i, j ] over add
     ]

}</lang>

Output:
amicables println
[[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

Direct approach

Works with: Turbo Pascal
Works with: Free Pascal

This version mutates the Sieve of Eratoshenes from striking out factors into summing factors. The Pascal source compiles with Turbo Pascal (7, patched to avoid the zero divide problem for cpu speeds better than ~150MHz) except that the array limit is too large: 15,000 works but does not reach 20,000. The Free Pascal compiler however can handle an array of 20,000 elements. Because the sum of factors of N can exceed N an ad-hoc SumF procedure is provided, thus the search could continue past the table limit, but at a cost in calculation time.

Output is 

Chasing Chains of Sums of Factors of Numbers.
Perfect!! 6,
Perfect!! 28,
Amicable! 220,284,
Perfect!! 496,
Amicable! 1184,1210,
Amicable! 2620,2924,
Amicable! 5020,5564,
Amicable! 6232,6368,
Perfect!! 8128,
Amicable! 10744,10856,
Amicable! 12285,14595,
Sociable: 12496,14288,15472,14536,14264,
Sociable: 14316,19116,31704,47616,83328,177792,295488,629072,589786,294896,358336,418904,366556,274924,275444,243760,376736,381028,285778,152990,122410,97946,48976,45946,22976,22744,19916,17716,
Amicable! 17296,18416,

Source file:<lang pascal> 

Program SumOfFactors; uses crt; {Perpetrated by R.N.McLean, December MCMXCV}

//{$DEFINE ShowOverflow} {$IFDEF FPC} 

 {$MODE DELPHI}//tested with lots = 524*1000*1000 takes 75 secs generating KnownSum

{$ENDIF} 

 var outf: text;
 const Limit = 2147483647;
 const lots = 20000;       {This should be much bigger, but problems apply.}
 var KnownSum: array[1..lots] of longint;
 Function SumF(N: Longint): Longint;
  var f,f2,s,ulp: longint;
  Begin
   if n <= lots then SumF:=KnownSum[N] {Hurrah!}
    else
     begin      {This is really crude...}
      s:=1;     {1 is always a factor, but N is not.}
      f:=2;
      f2:=f*f;
      while f2 < N do
       begin
        if N mod f = 0 then
         begin  {We have a divisor, and its friend.}
          ulp:=f + (N div f);
          if s > Limit - ulp then begin SumF:=-666; exit; end;
          s:=s + ulp;
         end;
        f:=f + 1;
        f2:=f*f;
       end;
       if f2 = N then {A perfect square gets its factor in once only.}
        if s <= Limit - f then s:=s + f
         else begin SumF:=-667; exit; end;
      SumF:=s;
     end;
  End;
 var i,j,l,sf,fs: LongInt;
 const enuff = 666; {Only so much sociability.}
 var trail: array[0..enuff] of longint;
 BEGIN
  ClrScr;
  WriteLn('Chasing Chains of Sums of Factors of Numbers.');
  for i:=1 to lots do KnownSum[i]:=1; {Sigh. KnownSum:=1;}

{start summing every divisor } 

  for i:=2 to lots do
   begin
    j:=i + i;
    While j <= lots do    {Sigh. For j:=i + i:Lots:i do KnownSum[j]:=KnownSum[j] + i;}
    begin
      KnownSum[j]:=KnownSum[j] + i;
      j:=j + i;
     end;
   end; 

{Enough preparation.}
  Assign(outf,'Factors.txt'); ReWrite(Outf);
  WriteLn(Outf,'Chasing Chains of Sums of Factors of Numbers.');

  for i:=2 to lots do    {Search.}
   begin
    l:=0;
    sf:=SumF(i);
    while (sf > i) and (l < enuff) do
     begin
      l:=l + 1;
      trail[l]:=sf;
      sf:=SumF(sf);
     end;
    if l >= enuff then writeln('Rope ran out! ',i);

{$IFDEF ShowOverflow} 

    if sf < 0 then writeln('Overflow with ',i);

{$ENDIF} 

    if i = sf then      {A loop?}
     begin              {Yes. Reveal its members.}
      trail[0]:=i;      {The first.}
      if l = 0 then write('Perfect!! ')
       else if l = 1 then write('Amicable! ')
        else write('Sociable: ');
      for j:=0 to l do Write(Trail[j],',');
      WriteLn;
      if l = 0 then write(outf,'Perfect!! ')
       else if l = 1 then write(outf,'Amicable! ')
        else write(outf,'Sociable: ');
      for j:=0 to l do write(outf,Trail[j],',');
      WriteLn(outf);
     end;
   end;
  Close (outf);
 END.</lang>

More expansive.

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

  1. 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.

Library: ntheory

<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

Translation of: REXX

<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 by using a   sigma   function   (which was a re-coded   proper divisors   (Pdivs)   function.

The   Pdivs   function was taken from the REXX language entry for Rosetta code task   integer factors.

Other optimizations were incorporated which took advantage of several well-known generalizations about amicable pairs.

The generation/summation is about forty times faster;   searching is about two orders of magnitude 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)  ; low=220 /*W: used for columnar alignment,*/ @.=. /* [↑] LOW is lowest amicable #.*/ 

    do k=low  for H-low;  _=sigma(k)  /*generate sigma sums for a range*/
    if _>=low  then @.k=_             /*only keep pertinent sigma sums.*/
    end   /*k*/                       /* [↑]   process a range of ints.*/
  1. =0 /*number of amicable pairs found.*/
    do   m=low  to  H;    n=@.m       /*start search at lowest number. */
    if n==.  then iterate             /*if not pertinent,  then ignore.*/
    if n==m  then iterate             /*Is sigma(M)=M?  Skip perfect #s*/
    if m==@.n  then do;   #=#+1       /*bump the amicable pair counter.*/
                    say right(m,w) ' and ' right(n,w) " are an amicable pair."
                    m=n               /*start  M (do loop index) from N*/
                    end
    end     /*m*/                     /*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.*/ /*──────────────────────────────────SIGMA subroutine────────────────────*/ sigma: procedure; parse arg x; if x<2 then return 0; odd=x//2 s=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  s=s+j+ x%j       /*add the two divisors to the sum*/
   end   /*j*/                        /* [↑]  %  is REXX integer divide*/
                                      /* [↓]  adjust for square.     _ */

if j*j==x then s=s+j /*Was X a square? If so, add √x.*/ return s /*return the sum of the divisors.*/</lang> output   when the default input is used: 

  220  and    284  are an amicable pair.
 1184  and   1210  are an amicable pair.
 2620  and   2924  are an amicable pair.
 5020  and   5564  are an amicable pair.
 6232  and   6368  are an amicable pair.
10744  and  10856  are an amicable pair.
12285  and  14595  are an amicable pair.
17296  and  18416  are an amicable pair.

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


Rust

<lang rust> fn sum_of_divisors(val:i64) -> i64 { 

       let mut sum = 0i64;

       for i in 1..val {
               if val % i == 0 {
                       sum += i;
               }
       }
               
       sum

}

fn main() { 

       for i in 1..20000i64 {
               let sum1 = sum_of_divisors(i);

               if sum1 <= i {
                       /*
                        * 1) Amicable pairs can not be the same
                        * number.
                        * 2) If sum1 is less than i then we would
                        * have already examined it for amicable pair.
                        */
                       continue;
               }

               let sum2 = sum_of_divisors(sum1);

               if sum2 == i {
                       println!("{} {}", i, sum1);
               }
       }

} </lang>

Output:
220 284
1184 1210
2620 2924
5020 5564
6232 6368
10744 10856
12285 14595
17296 18416


Scala

<lang Scala>def properDivisors(n: Int) = (1 to n/2).filter(i => n % i == 0) val divisorsSum = (1 to 20000).map(i => i -> properDivisors(i).sum).toMap val result = divisorsSum.filter(v => v._1 < v._2 && divisorsSum.get(v._2) == Some(v._1))

println( result mkString ", " )</lang>

Output:
5020 -> 5564, 220 -> 284, 6232 -> 6368, 17296 -> 18416, 2620 -> 2924, 10744 -> 10856, 12285 -> 14595, 1184 -> 1210

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}      //second part of pair out-of-bounds

       if a == i {continue}        //skip perfect numbers
       
       if i == aliquot[a] {
           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))
Retrieved from "https://rosettacode.org/wiki/Amicable_pairs?oldid=252714"