Two integers and are said to beamicable pairs if and the sum of theproper divisors of () as well as.
1184 and1210 are an amicable pair, with proper divisors:
Calculate and show here the Amicable pairs below 20,000; (there are eight).
F sum_proper_divisors(n) R I n < 2 {0} E sum((1 .. n I/ 2).filter(it -> (@n % it) == 0))L(n) 1..20000 V m = sum_proper_divisors(n) I m > n & sum_proper_divisors(m) == n print(n"\t"m)org100h;;;Calculate proper divisors of 2..20000lxih,pdiv + 4; 2 bytes per entrylxid,19999; [2 .. 20000] means 19999 entrieslxib,1; Initialize each entry to 1init:movm,cinxhmovm,binxhdcxdmova,doraejnzinitlxib,1; BC = outer loop variableiouter:inxblxih,-10001; Are we there yet?dadbjcidone; If so, we've calculated all of themmovh,bmovl,cdadhxchg; DE = inner loop variableiinner:pushd; save DExchgdad h; calculate *pdiv[DE]lxid,pdivdaddmove,m; DE = pdiv[DE]inxhmovd,mxchg; pdiv[DE] += BCdad bxchg; store it backmovm,ddcxhmovm,epoph; restore DE (into HL)dadb; add BClxid,-20001; are we there yet?daddjciouter; then continue with outer looplxid,20001; otherwise continue with inner loopdaddxchgjmpiinneridone:lxib,1; BC = outer loop variabletouter:inxblxih,-20001; Are we there yet?dadbrc; If so, stopmovd,b; DE = outer loop variablemove,ctinner:inxdlxih,-20001; Are we there yet?daddjctouter; If so continue with outer looppushd; Store the variablespushbmovh,b; find *pdiv[BC]movl,cdad b lxib,pdivdadbmova,m; Compare low byte (to E) cmpejnztnext1; Not equal = not amicableinxhmova,mcmpd; Compare high byte (to B) jnztnext1; Not equal = not amicablepopb; Restore BCxchg; find *pdiv[DE]dadhlxid,pdivdaddmova,m; Compare low byte (to C)cmpcjnztnext2; Not equal = not amicableinxhmova,m; Compare high byte (to B) cmp bjnztnext2; Not equal = not amicablepopd; Restore DEpushd; Save them both on the stack againpushbpush dmovh,b; Print the first numbermovl,ccallprhlpop h; And the second numbercallprhllxid,nl; And a newlinemvic,9call 5tnext1:popb; Restore Btnext2:popd; Restore Djmptinner; Continue;;;Print the number in HLprhl:lxid,nbuf; Store buffer pointer on stackpush dlxib,-10; Divisorpdgt:lxid,-1; Quotientpdivlp:inxddadbjcpdivlpmvia,'0'+10; Make ASCII digitaddlpoph; Store in output buffer dcxhmovm,apushhxchg; Keep going with rest of numbermova,h; if not zerooraljnzpdgtmvic,9; CP/M call to print stringpopd; Get buffer pointerjmp5db'*****'nbuf:db' $'nl:db13,10,'$'pdiv:equ$; base
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
LIMIT:equ20000; Maximum valuecpu8086org100hsection.textmovax,final; Set DS and ES to point just beyond themovcl,4; program. We're just going to assume MS-DOSshrax,cl; gave us enough memory. (Generally the case,incax; a .COM gets a 64K segment and we need ~40K.)movcx,csaddax,cxmovds,axmoves,axcalc:movax,1; Calculate proper divisors for 2..20000movdi,4; Initially, set each entry to 1.movcx,LIMIT-1; 2 to 20000 inclusive = 19999 entriesrepstoswmovax,2; AX = outer loop countermovcl,2movdx,LIMIT*2; Keep inner loop limit ready in DXmovbp,LIMIT/2; And outer loop limit in BP.outer:movbx,ax; BX = inner loop counter (multiplied by two)shlbx,cl; Each entry is 2 bytes wide.inner:add[bx],ax; divsum[BX/2] += AXaddbx,ax; Advance to next entryaddbx,ax; Twice, because each entry is 2 bytes widecmpbx,dx; Are we there yet?jbe.inner; If not, keep goingincaxcmpax,bp; Is the outer loop done yet?jbe.outer; If not, keep goingshow:movdx,LIMIT; Keep limit ready in DXmovax,2; AX = outer loop countermovsi,4; SI = address for outer loop.outer:movcx,ax; CX = inner loop counterinccxmovdi,cx; DI = address for inner loopshldi,1movbx,[si]; Preload divsum[AX].inner:cmpcx,bx; CX == divsum[AX]?jne.next; If not, the pair is not amicablecmpax,[di]; AX == divsum[CX]?jne.next; If not, the pair is not amicablepushax; Keep the registerspushbxpushcxpushdxpushcx; And CX twice because we need to print itcallprax; Print the first numberpopaxcallprax; And the second numbermovdx,nl; And a newlinecallpstrpopdx; Restore the registerspopcxpopbxpopax.next:incdi; Increment inner loop variable and addressincdi; Address twice because each entry has 2 bytesinccxcmpcx,dx; Are we done yet?jbe.inner; If not, keep goingincsi; Increment outer loop variable and addressincsi; Address twice because each entry has 2 bytesincaxcmpax,dx; Are we done yet?jbe.outer; If not, keep going.ret;;;Print the number in AX. Destroys AX, BX, CX, DX.prax:movcx,10; Divisormovbx,nbuf; Buffer pointer.digit:xordx,dxdivcx; Divide by 10 and extract digitadddl,'0'; Add ASCII 0 to digitdecbxmov[cs:bx],dl; Store in stringtestax,ax; Any more?jnz.digit; If so, keep goingmovdx,bx; If not, print the result;;;Print string from CS.pstr:pushds; Save DSmovax,cs; Set DS to CSmovds,axmovah,9; Print string using MS-DOSint21hpopds; Restore DSretdb'*****'nbuf:db'$'nl:db13,10,'$'final:equ$
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
/* ARM assembly AARCH64 Raspberry PI 3B *//* program amicable64.s */ /*******************************************//* Constantes file *//*******************************************//* for this file see task include a file in language AArch64 assembly*/.include "../includeConstantesARM64.inc" .equ NMAXI, 20000.equ TABMAXI, 100 /*********************************//* Initialized data *//*********************************/.datasMessResult: .asciz " @ : @\n"szCarriageReturn: .asciz "\n"szMessErr1: .asciz "Array too small !!"/*********************************//* UnInitialized data *//*********************************/.bss sZoneConv: .skip 24tResult: .skip 8 * TABMAXI/*********************************//* code section *//*********************************/.text.global main main: // entry of program ldr x3,qNMaxi // load limit mov x4,#2 // number begin1: mov x0,x4 // number bl decFactor // compute sum factors cmp x0,x4 // equal ? beq 2f mov x2,x0 // factor sum 1 bl decFactor cmp x0,x4 // equal number ? bne 2f mov x0,x4 // yes -> search in array mov x1,x2 // and store sum bl searchRes cmp x0,#0 // find ? bne 2f // yes mov x0,x4 // no -> display number ans sum mov x1,x2 bl displayResult2: add x4,x4,#1 // increment number cmp x4,x3 // end ? ble 1b 100: // standard end of the program mov x0, #0 // return code mov x8, #EXIT // request to exit program svc #0 // perform the system callqAdrszCarriageReturn: .quad szCarriageReturnqNMaxi: .quad NMAXI/***************************************************//* display message number *//***************************************************//* x0 contains number 1 *//* x1 contains number 2 */displayResult: stp x1,lr,[sp,-16]! // save registers stp x2,x3,[sp,-16]! // save registers mov x2,x1 ldr x1,qAdrsZoneConv bl conversion10 // call décimal conversion ldr x0,qAdrsMessResult ldr x1,qAdrsZoneConv // insert conversion in message bl strInsertAtCharInc mov x3,x0 mov x0,x2 ldr x1,qAdrsZoneConv bl conversion10 // call décimal conversion mov x0,x3 ldr x1,qAdrsZoneConv // insert conversion in message bl strInsertAtCharInc bl affichageMess // display message ldp x2,x3,[sp],16 // restaur 2 registers ldp x1,lr,[sp],16 // restaur 2 registers ret // return to address lr x30qAdrsMessResult: .quad sMessResultqAdrsZoneConv: .quad sZoneConv/***************************************************//* compute factors sum *//***************************************************//* x0 contains the number */decFactor: stp x1,lr,[sp,-16]! // save registers stp x2,x3,[sp,-16]! // save registers stp x4,x5,[sp,-16]! // save registers mov x4,#1 // init sum mov x1,#2 // start factor -> divisor1: udiv x2,x0,x1 msub x3,x2,x1,x0 // remainder cmp x1,x2 // divisor > quotient ? bgt 3f cmp x3,#0 // remainder = 0 ? bne 2f add x4,x4,x1 // add divisor to sum cmp x1,x2 // divisor = quotient ? beq 3f // yes -> end add x4,x4,x2 // no -> add quotient to sum2: add x1,x1,#1 // increment factor b 1b // and loop3: mov x0,x4 // return sum ldp x4,x5,[sp],16 // restaur 2 registers ldp x2,x3,[sp],16 // restaur 2 registers ldp x1,lr,[sp],16 // restaur 2 registers ret // return to address lr x30/***************************************************//* search and store result in array *//***************************************************//* x0 contains the number *//* x1 contains factors sum *//* x0 return 1 if find 0 else -1 if error */searchRes: stp x1,lr,[sp,-16]! // save registers stp x2,x3,[sp,-16]! // save registers stp x4,x5,[sp,-16]! // save registers ldr x4,qAdrtResult // array address mov x2,#0 // indice begin1: ldr x3,[x4,x2,lsl #3] // load one result array cmp x3,#0 // if 0 store new result beq 2f cmp x3,x0 // equal ? beq 3f // find -> return 1 add x2,x2,#1 // increment indice cmp x2,#TABMAXI // maxi array ? blt 1b ldr x0,qAdrszMessErr1 // error bl affichageMess mov x0,#-1 b 100f2: str x1,[x4,x2,lsl #3] mov x0,#0 // not find -> store and retun 0 b 100f3: mov x0,#1100: ldp x4,x5,[sp],16 // restaur 2 registers ldp x2,x3,[sp],16 // restaur 2 registers ldp x1,lr,[sp],16 // restaur 2 registers ret // return to address lr x30qAdrtResult: .quad tResultqAdrszMessErr1: .quad szMessErr1/********************************************************//* File Include fonctions *//********************************************************//* for this file see task include a file in language AArch64 assembly */.include "../includeARM64.inc"
220 : 284 1184 : 1210 2620 : 2924 5020 : 5564 6232 : 6368 10744 : 10856 12285 : 14595 17296 : 18416
HOW TO RETURN proper.divisor.sum.table n: PUT {} IN propdivs FOR i IN {1..n}: PUT 1 IN propdivs[i] FOR i IN {2..floor (n/2)}: PUT i+i IN j WHILE j<=n: PUT propdivs[j] + i IN propdivs[j] PUT i + j IN j RETURN propdivsPUT 20000 IN maximumPUT proper.divisor.sum.table maximum IN propdivsFOR cand IN {1..maximum}: PUT propdivs[cand] IN other IF cand<other<maximum AND propdivs[other]=cand: WRITE cand, other/220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
Calculations on a real Atari 8-bit computer take quite long time. It is recommended to use an emulator capable with increasing speed of Atari CPU.
INCLUDE "H6:SIEVE.ACT"CARD FUNC SumDivisors(CARD x) CARD i,max,sum sum=1 i=2 max=x WHILE i<max DO IF x MOD i=0 THEN max=x/i IF i<max THEN sum==+i+max ELSEIF i=max THEN sum==+i FI FI i==+1 ODRETURN (sum)PROC Main() DEFINE MAXNUM="20000" BYTE ARRAY primes(MAXNUM+1) CARD m,n Put(125) PutE() ;clear the screen Sieve(primes,MAXNUM+1) FOR m=1 TO MAXNUM-1 DO IF primes(m)=0 THEN n=SumDivisors(m) IF n<MAXNUM AND primes(n)=0 AND n>m THEN IF m=SumDivisors(n) THEN PrintF("%U %U%E",m,n) FI FI FI ODRETURNScreenshot from Atari 8-bit computer
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
This solution uses the packageGeneric_Divisors from the Proper Divisors task [[1]].
withAda.Text_IO,Generic_Divisors;useAda.Text_IO;procedureAmicable_PairsisfunctionSame(P:Positive)returnPositiveis(P);packageDivisor_Sumis newGeneric_Divisors(Result_Type=> Natural,None=> 0,One=> Same,Add=> "+");Num2:Integer;beginforNum1in4..20_000loopNum2:=Divisor_Sum.Process(Num1);ifNum1<Num2thenifNum1=Divisor_Sum.Process(Num2)thenPut_Line(Integer'Image(Num1)&","&Integer'Image(Num2));endif;endif;endloop;endAmicable_Pairs;
220, 284 1184, 1210 2620, 2924 5020, 5564 6232, 6368 10744, 10856 12285, 14595 17296, 18416
begincomment - return n mod m;integer procedure mod(n, m); value n, m; integer n, m;begin mod := n - m * entier(n / m);end;comment - return sum of the proper divisors of n;integer procedure sumf(n); value n; integer n;begin integer sum, f1, f2; sum := 1; f1 := 2; for f1 := f1 while (f1 * f1) <= n do begin if mod(n, f1) = 0 then begin sum := sum + f1; f2 := n / f1; if f2 > f1 then sum := sum + f2; end; f1 := f1 + 1; end; sumf := sum;end;comment - main program begins here;integer a, b, found;outstring(1,"Searching up to 20000 for amicable pairs\n");found := 0;for a := 2 step 1 until 20000 do begin b := sumf(a); if b > a then begin if a = sumf(b) then begin found := found + 1; outinteger(1,a); outinteger(1,b); outstring(1,"\n"); end; end; end;outinteger(1,found);outstring(1,"pairs were found");end
Searching up to 20000 for amicable pairs 220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416 8 pairs were found
BEGIN # find amicable pairs p1, p2 where each is equal to the other's proper divisor sum # [ 1 : 20 000 ]INT pd sum; # table of proper divisors # FOR n TO UPB pd sum DO pd sum[ n ] := 1 OD; FOR i FROM 2 TO UPB pd sum DO FOR j FROM i + i BY i TO UPB pd sum DO pd sum[ j ] +:= i OD OD; # find the amicable pairs up to 20 000 # FOR p1 TO UPB pd sum - 1 DO INT pd sum p1 = pd sum[ p1 ]; IF pd sum p1 > p1 AND pd sum p1 <= UPB pd sum THEN IF pd sum[ pd sum p1 ] = p1 THEN print( ( whole( p1, -6 ), " and ", whole( pd sum p1, -6 ), " are an amicable pair", newline ) ) FI FI ODEND
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
begin % find amicable pairs p1, p2 where each is equal to the other's % % proper divisor sum % integer MAX_NUMBER; MAX_NUMBER := 20000; begin integer array pdSum( 1 :: MAX_NUMBER ); % table of proper divisors % for i := 1 until MAX_NUMBER do pdSum( i ) := 1; for i := 2 until MAX_NUMBER do begin for j := i + i step i until MAX_NUMBER do pdSum( j ) := pdSum( j ) + i end for_i ; % find the amicable pairs up to 20 000 % for p1 := 1 until MAX_NUMBER - 1 do begin integer pdSumP1; pdSumP1 := pdSum( p1 ); if pdSumP1 > p1 and pdSumP1 <= MAX_NUMBER and pdSum( pdSumP1 ) = p1 then begin write( i_w := 5, s_w := 0, p1, " and ", pdSumP1, " are an amicable pair" ) end if_pdSumP1_gt_p1_and_le_MAX_NUMBER end for_p1 endend.
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 pair10744 and 10856 are an amicable pair12285 and 14595 are an amicable pair17296 and 18416 are an amicable pair
-- AMICABLE PAIRS -------------------------------------------------------------- amicablePairsUpTo :: Int -> IntonamicablePairsUpTo(max)-- amicable :: [Int] -> Int -> Int -> [Int] -> [Int]scriptamicableon|λ|(a,m,n,lstSums)if(m>n)and(m≤max)and((itemmoflstSums)=n)thena&[[n,m]]elseaendifend|λ|endscript-- divisorsSummed :: Int -> IntscriptdivisorsSummed-- sum :: Int -> Int -> Intscriptsumon|λ|(a,b)a+bend|λ|endscripton|λ|(n)foldl(sum,0,properDivisors(n))end|λ|endscriptfoldl(amicable,{},¬map(divisorsSummed,enumFromTo(1,max)))endamicablePairsUpTo-- TEST ----------------------------------------------------------------------onrunamicablePairsUpTo(20000)endrun-- PROPER DIVISORS ------------------------------------------------------------- properDivisors :: Int -> [Int]onproperDivisors(n)-- isFactor :: Int -> BoolscriptisFactoron|λ|(x)nmodx=0end|λ|endscript-- integerQuotient :: Int -> IntscriptintegerQuotienton|λ|(x)(n/x)asintegerend|λ|endscriptifn=1then{1}elsesetrealRootton^(1/2)setintRoottorealRootasintegersetblnPerfectSquaretointRoot=realRoot-- Factors up to square root of n,setlowstofilter(isFactor,enumFromTo(1,intRoot))-- and quotients of these factors beyond the square root,-- excluding n itself (last item)items1thru-2of(lows&map(integerQuotient,¬items(1+(blnPerfectSquareasinteger))thru-1ofreverseoflows))endifendproperDivisors-- GENERIC FUNCTIONS ----------------------------------------------------------- enumFromTo :: Int -> Int -> [Int]onenumFromTo(m,n)ifm>nthensetdto-1elsesetdto1endifsetlstto{}repeatwithifrommtonbydsetendoflsttoiendrepeatreturnlstendenumFromTo-- filter :: (a -> Bool) -> [a] -> [a]onfilter(f,xs)tellmReturn(f)setlstto{}setlngtolengthofxsrepeatwithifrom1tolngsetvtoitemiofxsif|λ|(v,i,xs)thensetendoflsttovendrepeatreturnlstendtellendfilter-- foldl :: (a -> b -> a) -> a -> [b] -> aonfoldl(f,startValue,xs)tellmReturn(f)setvtostartValuesetlngtolengthofxsrepeatwithifrom1tolngsetvto|λ|(v,itemiofxs,i,xs)endrepeatreturnvendtellendfoldl-- map :: (a -> b) -> [a] -> [b]onmap(f,xs)tellmReturn(f)setlngtolengthofxssetlstto{}repeatwithifrom1tolngsetendoflstto|λ|(itemiofxs,i,xs)endrepeatreturnlstendtellendmap-- Lift 2nd class handler function into 1st class script wrapper-- mReturn :: Handler -> ScriptonmReturn(f)ifclassoffisscriptthenfelsescriptproperty|λ|:fendscriptendifendmReturn
{{220,284},{1184,1210},{2620,2924},{5020,5564},{6232,6368},{10744,10856},{12285,14595},{17296,18416}}
… and about 55 times as fast as the above.
onproperDivisors(n)setoutputto{}if(n>1)thensetsqrtton^0.5setlimittosqrtdiv1if(limit=sqrt)thensetendofoutputtolimitsetlimittolimit-1endifrepeatwithifromlimitto2by-1if(nmodiis0)thensetbeginningofoutputtoisetendofoutputtondiviendifendrepeatsetbeginningofoutputto1endifreturnoutputendproperDivisorsonsumList(listOfNumbers)scriptopropertyl:listOfNumbersendscriptsetsumto0repeatwithnino'slsetsumtosum+nendrepeatreturnsumendsumListonamicablePairsBelow(limitPlus1)scriptopropertypdSums:{missing value}-- Sums of proper divisors. (Dummy item for 1's.)endscriptsetlimittolimitPlus1-1repeatwithnfrom2tolimitsetendofo'spdSumstosumList(properDivisors(n))endrepeatsetoutputto{}repeatwithn1from2to(limit-1)setn2too'spdSums'sitemn1if((n1<n2)and(n2<limitPlus1)and(o'spdSums'sitemn2=n1))then¬setendofoutputto{n1,n2}endrepeatreturnoutputendamicablePairsBelowonjoin(lst,delim)setastidtoAppleScript'stext item delimiterssetAppleScript'stext item delimiterstodelimsettxttolstastextsetAppleScript'stext item delimiterstoastidreturntxtendjoinontask()setoutputtoamicablePairsBelow(20000)repeatwiththisPairinoutputsetthisPair'scontentstojoin(thisPair," & ")endrepeatreturnjoin(output,linefeed)endtasktask()
"220 & 2841184 & 12102620 & 29245020 & 55646232 & 636810744 & 1085612285 & 1459517296 & 18416"
/* ARM assembly Raspberry PI or android with termux *//* program amicable.s */ /* REMARK 1 : this program use routines in a include file see task Include a file language arm assembly for the routine affichageMess conversion10 see at end of this program the instruction include *//* for constantes see task include a file in arm assembly *//************************************//* Constantes *//************************************/.include "../constantes.inc" .equ NMAXI, 20000.equ TABMAXI, 100 /*********************************//* Initialized data *//*********************************/.datasMessResult: .asciz " @ : @\n"szCarriageReturn: .asciz "\n"szMessErr1: .asciz "Array too small !!"/*********************************//* UnInitialized data *//*********************************/.bss sZoneConv: .skip 24tResult: .skip 4 * TABMAXI/*********************************//* code section *//*********************************/.text.global main main: @ entry of program ldr r3,iNMaxi @ load limit mov r4,#2 @ number begin1: mov r0,r4 @ number bl decFactor @ compute sum factors cmp r0,r4 @ equal ? beq 2f mov r2,r0 @ factor sum 1 bl decFactor cmp r0,r4 @ equal number ? bne 2f mov r0,r4 @ yes -> search in array mov r1,r2 @ and store sum bl searchRes cmp r0,#0 @ find ? bne 2f @ yes mov r0,r4 @ no -> display number ans sum mov r1,r2 bl displayResult2: add r4,#1 @ increment number cmp r4,r3 @ end ? ble 1b 100: @ standard end of the program mov r0, #0 @ return code mov r7, #EXIT @ request to exit program svc #0 @ perform the system calliAdrszCarriageReturn: .int szCarriageReturniNMaxi: .int NMAXI/***************************************************//* display message number *//***************************************************//* r0 contains number 1 *//* r1 contains number 2 */displayResult: push {r1-r3,lr} @ save registers mov r2,r1 ldr r1,iAdrsZoneConv bl conversion10 @ call décimal conversion ldr r0,iAdrsMessResult ldr r1,iAdrsZoneConv @ insert conversion in message bl strInsertAtCharInc mov r3,r0 mov r0,r2 ldr r1,iAdrsZoneConv bl conversion10 @ call décimal conversion mov r0,r3 ldr r1,iAdrsZoneConv @ insert conversion in message bl strInsertAtCharInc bl affichageMess @ display message pop {r1-r3,pc} @ restaur des registresiAdrsMessResult: .int sMessResultiAdrsZoneConv: .int sZoneConv/***************************************************//* compute factors sum *//***************************************************//* r0 contains the number */decFactor: push {r1-r5,lr} @ save registers mov r5,#1 @ init sum mov r4,r0 @ save number mov r1,#2 @ start factor -> divisor1: mov r0,r4 @ dividende bl division cmp r1,r2 @ divisor > quotient ? bgt 3f cmp r3,#0 @ remainder = 0 ? bne 2f add r5,r5,r1 @ add divisor to sum cmp r1,r2 @ divisor = quotient ? beq 3f @ yes -> end add r5,r5,r2 @ no -> add quotient to sum2: add r1,r1,#1 @ increment factor b 1b @ and loop3: mov r0,r5 @ return sum pop {r1-r5,pc} @ restaur registers/***************************************************//* search and store result in array *//***************************************************//* r0 contains the number *//* r1 contains factors sum *//* r0 return 1 if find 0 else -1 if error */searchRes: push {r1-r4,lr} @ save registers ldr r4,iAdrtResult @ array address mov r2,#0 @ indice begin1: ldr r3,[r4,r2,lsl #2] @ load one result array cmp r3,#0 @ if 0 store new result beq 2f cmp r3,r0 @ equal ? moveq r0,#1 @ find -> return 1 beq 100f add r2,r2,#1 @ increment indice cmp r2,#TABMAXI @ maxi array ? blt 1b ldr r0,iAdrszMessErr1 @ error bl affichageMess mov r0,#-1 b 100f2: str r1,[r4,r2,lsl #2] mov r0,#0 @ not find -> store and retun 0100: pop {r1-r4,pc} @ restaur registersiAdrtResult: .int tResultiAdrszMessErr1: .int szMessErr1/***************************************************//* ROUTINES INCLUDE *//***************************************************/.include "../affichage.inc"220 : 284 1184 : 1210 2620 : 2924 5020 : 5564 6232 : 6368 10744 : 10856 12285 : 14595 17296 : 18416
properDivs:function[x]->(factors x)--xamicable:function[x][y:sumproperDivsxifand?x=sumproperDivsyx<>y->return@[x,y]returnø]amicables:[]loop1..20000'n[am:amicablenifam<>ø->'amicables++@[sortam]]printuniqueamicables
[220 284] [1184 1210] [2620 2924] [5020 5564] [6232 6368] [10744 10856] [12285 14595] [17296 18416]
(* ****** ****** *)//#include"share/atspre_staload.hats"#include"share/HATS/atspre_staload_libats_ML.hats"//(* ****** ****** *)//funsum_list_vt (xs: List_vt(int)): int =( case+ xs of | ~list_vt_nil() => 0 | ~list_vt_cons(x, xs) => x + sum_list_vt(xs))//(* ****** ****** *)funpropDivs( x0: int) : List0_vt(int) = loop(x0, 2, list_vt_sing(1)) where{//funloop(x0: int, i: int, res: List0_vt(int)) : List0_vt(int) =(if(i * i) > x0then list_vt_reverse(res)else( if x0 % i != 0 then loop(x0, i+1, res) // end of [then] else let val res = cons_vt(i, res) // end of [val] val res = ( if i * i = x0 then res else cons_vt(x0 / i, res) ) : List0_vt(int) // end of [val] in loop(x0, i+1, res) end // end of [else] // end of [if])) (* end of [loop] *)//} // end of [propDivs](* ****** ****** *)funsum_propDivs(x: int): int = sum_list_vt(propDivs(x))(* ****** ****** *)valtheNat2 = auxmain(2) where{funauxmain( n: int) : stream_vt(int) = $ldelay(stream_vt_cons(n, auxmain(n+1)))}(* ****** ****** *)//valtheAmicable =(stream_vt_takeLte(theNat2, 20000)).filter()(lam x =>let val x2 = sum_propDivs(x)in x < x2 && x = sum_propDivs(x2) end)//(* ****** ****** *)val () =theAmicable.foreach()( lam x => println! ("(", x, ", ", sum_propDivs(x), ")"))(* ****** ****** *)implement main0 () = ()(* ****** ****** *)(220, 284)(1184, 1210)(2620, 2924)(5020, 5564)(6232, 6368)(10744, 10856)(12285, 14595)(17296, 18416)
SetBatchLines-1Loop,20000{m:=A_index;Gettingfactorsloop%floor(sqrt(m)){if(mod(m,A_index)=0){if(A_index**2==m){sum+=A_indexcontinue}elseif(A_index!=1){sum+=A_index+m//A_index}elseif(A_index=1){sum+=A_index}}};Factorsobtained;Checkingfactorsofsumif(sum>1){loop%floor(sqrt(sum)){if(mod(sum,A_index)=0){if(A_index**2==sum){sum2+=A_indexcontinue}elseif(A_index!=1){sum2+=A_index+sum//A_index}elseif(A_index=1){sum2+=A_index}}}if(m=sum2)&&(m!=sum)&&(m<sum)final.=m.":".sum."`n"};Checkedsum:=0sum2:=0}MsgBox%finalExitApp
220:2841184:12102620:29245020:55646232:636810744:1085612285:1459517296:18416
#!/bin/awk -ffunctionsumprop(num,i,sum,root){if(num<2)return0sum=1root=sqrt(num)for(i=2;i<root;i++){if(num%i==0){sum=sum+i+num/i}}if(num%root==0){sum=sum+root}returnsum}BEGIN{limit=20000print"Amicable pairs < ",limitfor(n=1;n<limit+1;n++){m=sumprop(n)if(n==sumprop(m)&&n<m)printn,m}}}
# ./amicable Amicable pairs < 20000220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
100DECLAREEXTERNALFUNCTIONsum_proper_divisors110CLEAR120!130DIMf(20001)!sumofproperfactorsforeachn140FORi=1TO20000150LETf(i)=sum_proper_divisors(i)160NEXTi170!lookforpairs180FORi=1TO20000190FORj=i+1TO20000200IFf(i)=jANDi=f(j)THEN210PRINT"Amicable pair ";i;" ";j220ENDIF230NEXTj240NEXTi250!260PRINT270PRINT"-- found all amicable pairs"280END290!300!Computethesumofproperdivisorsofgivennumber310!320EXTERNALFUNCTIONsum_proper_divisors(n)330!340IFn>1THEN!nmustbe2orlarger350LETsum=1!startwith1360LETroot=SQR(n)!notethatrootisaninteger370!checkpossiblefactors,uptosqrt380FORi=2TOroot390IFMOD(n,i)=0THEN400LETsum=sum+i!iisafactor410IFi*i<>nTHEN!checkiisnotactualsquarerootofn420LETsum=sum+n/i!son/iwillalsobeafactor430ENDIF440ENDIF450NEXTi460ENDIF470LETsum_proper_divisors=sum480ENDFUNCTION
Amicable pair 220 284 Amicable pair 1184 1210 Amicable pair 2620 2924 Amicable pair 5020 5564 Amicable pair 6232 6368 Amicable pair 10744 10856 Amicable pair 12285 14595 Amicable pair 17296 18416 -- found all amicable pairs
function SumProperDivisors(number)if number < 2 then return 0sum = 0for i = 1 to number \ 2if number mod i = 0 then sum += inext ireturn sumend functiondim sum(20000)for n = 1 to 19999sum[n] = SumProperDivisors(n)next nprint "The pairs of amicable numbers below 20,000 are :"printfor n = 1 to 19998f = sum[n]if f <= n or f < 1 or f > 19999 then continue forif f = sum[n] and n = sum[f] thenprint rjust(string(n), 5); " and "; sum[n]end ifnext nend
The pairs of amicable numbers below 20,000 are : 220 and 284 1184 and 1210 2620 and 2924 5020 and 5564 6232 and 636810744 and 1085612285 and 1459517296 and 18416
search.limit%=20000dimsumf%(search.limit%)print"Searching up to";search.limit%;"for amicable pairs:"rem-createatableofproperdivisorsumsfora%=1tosearch.limit%sumf%(a%)=1nexta%fora%=2tosearch.limit%b%=a%+a%while(b%>0)and(b%<=search.limit%)sumf%(b%)=sumf%(b%)+a%b%=b%+a%wendnexta%rem-searchforpairscount%=0a%=2whilea%<search.limit%b%=sumf%(a%)rem-protectagainstinvalidarrayindexifb%<=0orb%>=search.limit%thenb%=2rem-otherwise,we're good to goif(b%>a%)and(a%=sumf%(b%))then\printusing"##### #####";a%;b%:\count%=count%+1a%=a%+1wendprintcount%;"pairs were found"end
Searching up to 20000 for amicable pairs: 220 284 1184 1210 2620 2924 5020 5564 6232 636810744 1085612285 1459517296 18416 8 pairs were found
100cls:rem10HOMEforApplesoftBASIC110print"The pairs of amicable numbers below 20,000 are :"120print130size=18500140forn=1tosize150m=amicable(n)160ifm>nandamicable(m)=nthen170printusing"#####";n;180print" and ";190printusing"#####";m200endif210next220end230functionamicable(nr)240suma=1250ford=2tosqr(nr)260ifnrmodd=0thensuma=suma+d+nr/d270next280amicable=suma290endfunction
Same as FreeBASIC entry.
Publicsum[19999]AsIntegerPublicSubMain()DimnAsInteger,fAsIntegerForn=0To19998sum[n]=SumProperDivisors(n)NextPrint"The pairs of amicable numbers below 20,000 are :\n"Forn=0To19998' f = SumProperDivisors(n)f=sum[n]Iff<=nOrf<1Orf>19999ThenContinueIff=sum[n]Andn=sum[f]ThenPrintFormat$(Str$(n),"#####");" And ";Format$(Str$(sum[n]),"#####")EndIfNextEndFunctionSumProperDivisors(numberAsInteger)AsIntegerIfnumber<2ThenReturn0DimsumAsInteger=0ForiAsInteger=1Tonumber\2IfnumberModi=0Thensum+=iNextReturnsumEndFunction
Same as FreeBASIC entry.
FUNCTIONamicable(nr)suma=1FORd=2TOSQR(nr)IFnrMODd=0THENsuma=suma+d+nr/dNEXTamicable=sumaENDFUNCTIONPRINT"The pairs of amicable numbers below 20,000 are :"PRINTsize=18500FORn=1TOsizem=amicable(n)IFm>nANDamicable(m)=nTHENPRINTUSING"##### and #####";n;mENDIFNEXTEND
Same as FreeBASIC entry.
FUNCTIONamicable(nr)LETsuma=1FORd=2TOSQR(nr)IFREMAINDER(nr,d)=0THENLETsuma=suma+d+nr/dENDIFNEXTdLETamicable=sumaENDFUNCTIONPRINT"The pairs of amicable numbers below 20,000 are :"PRINTLETsize=18500FORn=1TOsizeLETm=amicable(n)IFm>nANDamicable(m)=nTHENPRINTUSING"##### and #####":n,mNEXTnEND
Same as FreeBASIC entry.
get "libhdr"manifest $( MAXIMUM = 20000$)// Calculate proper divisors for 1..Nlet propDivSums(n) = valof$( let v = getvec(n) for i = 1 to n do v!i := 1 for i = 2 to n/2 do $( let j = i*2 while j < n do $( v!j := v!j + i j := j + i $) $) resultis v$)// Are A and B an amicable pair, given the list of sums of proper divisors?let amicable(pdiv, a, b) = a = pdiv!b & b = pdiv!alet start() be$( let pds = propDivSums(MAXIMUM) for x = 1 to MAXIMUM do for y = x+1 to MAXIMUM do if amicable(pds, x, y) do writef("%N, %N*N", x, y)$)220, 2841184, 12102620, 29245020, 55646232, 636810744, 1085612285, 1459517296, 18416
v_@#-*8*:"2":$_:#!2#*8#g*#6:#0*#!:#-*#<v>*/.55+,1>$$:28*:*:*%\28*:*:*/`06p28*:*:*/\2v%%^:*:<>*v+|!:-1g60/*:*:*82::+**:*:<<>:#**#8:#<*^>.28*^8::v>>*:*%/\28*:*:*%+\v>8+#$^#_+#`\:#0<:\`1/*:*2#<2v^:*82\/*:*:*82:::_v#!%%*:*:*82\/*:*:*82::<_^#<>>06p:28*:*:**1+01-\>1+::28*:*:*/\28*:*:*%:*\`!^
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
Remark:Look at Pascal Alternative [[2]].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)
#include<stdio.h>#include<stdlib.h>typedefunsignedintuint;intmain(intargc,char**argv){uinttop=atoi(argv[1]);uint*divsum=malloc((top+1)*sizeof(*divsum));uintpows[32]={1,0};for(uinti=0;i<=top;i++)divsum[i]=1;// sieve// only sieve within lower half , the modification starts at 2*pfor(uintp=2;p+p<=top;p++){if(divsum[p]>1){divsum[p]-=p;// subtract number itself from divisor sum ('proper')continue;}// p not primeuintx;// highest power of p we need//checking x <= top/y instead of x*y <= top to avoid overflowfor(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-1uintk=p-1;for(uintn=p+p;n<=top;n+=p){uints=1+pows[1];k--;// search the right power only if neededif(k==0){for(uinti=2;i<x&&!(n%pows[i]);s+=pows[i++]);k=p;}divsum[n]*=s;}}//now correct the upper halffor(uintp=(top>>1)+1;p<=top;p++){if(divsum[p]>1){divsum[p]-=p;}}uintcnt=0;for(uinta=1;a<=top;a++){uintb=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);return0;}
% ./a.out 20000220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416Top 20000 count : 8% ./a.out 524000000..475838415 514823985491373104 511419856509379344 523679536Top 524000000 count : 442real 0m16.285suser 0m16.156s
usingSystem;usingSystem.Collections.Generic;usingSystem.Linq;namespaceRosettaCode.AmicablePairs{internalstaticclassProgram{privateconstintLimit=20000;privatestaticvoidMain(){foreach(varpairinGetPairs(Limit)){Console.WriteLine("{0} {1}",pair.Item1,pair.Item2);}}privatestaticIEnumerable<Tuple<int,int>>GetPairs(intmax){List<int>divsums=Enumerable.Range(0,max+1).Select(i=>ProperDivisors(i).Sum()).ToList();for(inti=1;i<divsums.Count;i++){intsum=divsums[i];if(i<sum&&sum<=divsums.Count&&divsums[sum]==i){yieldreturnnewTuple<int,int>(i,sum);}}}privatestaticIEnumerable<int>ProperDivisors(intnumber){returnEnumerable.Range(1,number/2).Where(divisor=>number%divisor==0);}}}
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
#include<vector>#include<unordered_map>#include<iostream>intmain(){std::vector<int>alreadyDiscovered;std::unordered_map<int,int>divsumMap;intcount=0;for(intN=1;N<=20000;++N){intdivSumN=0;for(inti=1;i<=N/2;++i){if(fmod(N,i)==0){divSumN+=i;}}// populate map of integers to the sum of their proper divisorsif(divSumN!=1)// do not include primesdivsumMap[N]=divSumN;for(std::unordered_map<int,int>::iteratorit=divsumMap.begin();it!=divsumMap.end();++it){intM=it->first;intdivSumM=it->second;intdivSumN=divsumMap[N];if(N!=M&&divSumM==N&&divSumN==M){// do not print duplicate pairsif(std::find(alreadyDiscovered.begin(),alreadyDiscovered.end(),N)!=alreadyDiscovered.end())break;std::cout<<"["<<M<<", "<<N<<"]"<<std::endl;alreadyDiscovered.push_back(M);alreadyDiscovered.push_back(N);count++;}}}std::cout<<count<<" amicable pairs discovered"<<std::endl;}
[220, 284][1184, 1210][2620, 2924][5020, 5564][6232, 6368][10744, 10856][12285, 14595][17296, 18416]8 amicable pairs discovered
(nsexample(:gen-class))(defnfactors[n]" Find the proper factors of a number "(into(sorted-set)(mapcat(fn[x](if(=x1)[x][x(/nx)]))(filter#(zero?(remn%))(range1(inc(Math/sqrtn)))))))(deffind-pairs(into#{}(for[n(range220000):let[f(factorsn); Factors of nM(apply+f); Sum of factorsg(factorsM); Factors of sumN(apply+g)]; Sum of Factors of sum:when(=nN); (sum(proDivs(N)) = M and sum(propDivs(M)) = N:when(not=MN)]; N not-equal M(sorted-setnM)))); Found pair;; Output Results(doseq[qfind-pairs](printlnq))
#{220 284}#{6232 6368}#{1184 1210}#{5020 5564}#{2620 2924}#{12285 14595}#{17296 18416}#{10744 10856}% Generate proper divisors from 1 to maxproper_divisors = proc (max: int) returns (array[int]) divs: array[int] := array[int]$fill(1, max, 0) for i: int in int$from_to(1, max/2) do for j: int in int$from_to_by(i*2, max, i) do divs[j] := divs[j] + i end end return(divs)end proper_divisors% Are A and B and amicable pair, given the proper divisors?amicable = proc (divs: array[int], a, b: int) returns (bool) return(divs[a] = b & divs[b] = a)end amicable% Find all amicable pairs up to 20 000start_up = proc () max = 20000 po: stream := stream$primary_output() divs: array[int] := proper_divisors(max) for a: int in int$from_to(1, max) do for b: int in int$from_to(a+1, max) do if amicable(divs, a, b) then stream$putl(po, int$unparse(a) || ", " || int$unparse(b)) end end endend start_up
220, 2841184, 12102620, 29245020, 55646232, 636810744, 1085612285, 1459517296, 18416
(let((cache(make-hash-table)))(defunsum-proper-divisors(n)(or(gethashncache)(setf(gethashncache)(loopforxfrom1to(/n2)when(zerop(remnx))sumx)))))(defunamicable-pairs-up-to(n)(loopforxfrom1tonforsum-divs=(sum-proper-divisorsx)when(and(<xsum-divs)(=x(sum-proper-divisorssum-divs)))collect(listxsum-divs)))(amicable-pairs-up-to20000)
((220 284) (1184 1210) (2620 2924) (5020 5564) (6232 6368) (10744 10856) (12285 14595) (17296 18416))
include "cowgol.coh";const LIMIT := 20000;# Calculate sums of proper divisorsvar divSum: uint16[LIMIT + 1];var i: @indexof divSum;var j: @indexof divSum;i := 2;while i <= LIMIT loop divSum[i] := 1; i := i + 1;end loop;i := 2;while i <= LIMIT/2 loop j := i * 2; while j <= LIMIT loop divSum[j] := divSum[j] + i; j := j + i; end loop; i := i + 1;end loop;# Test each pairi := 2;while i <= LIMIT loop j := i + 1; while j <= LIMIT loop if divSum[i] == j and divSum[j] == i then print_i32(i as uint32); print(", "); print_i32(j as uint32); print_nl(); end if; j := j + 1; end loop; i := i + 1;end loop;220, 2841184, 12102620, 29245020, 55646232, 636810744, 1085612285, 1459517296, 18416
MX=524_000_000N=Math.sqrt(MX).to_u32x=Array(Int32).new(MX+1,1)(2..N).each{|i|p=i*ix[p]+=ik=i+i+1(p+i..MX).step(i){|j|x[j]+=kk+=1}}(4..MX).each{|m|n=x[m]ifn<m&&n!=0&&m==x[n]puts"#{n}#{m}"end}
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416...... ............. .......426191535 514780497475838415 514823985509379344 523679536
voidmain()@safe/*@nogc*/{importstd.stdio,std.algorithm,std.range,std.typecons,std.array;immutableproperDivs=(inuintn)purenothrow@safe/*@nogc*/=>iota(1,(n+1)/2+1).filter!(x=>n%x==0);enumrangeMax=20_000;auton2d=iota(1,rangeMax+1).map!(n=>properDivs(n).sum);foreach(immutablen,immutabledivSum;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));}
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]
SeePascal.
/* Fill a given array such that for each N, * P[n] is the sum of proper divisors of N */proc nonrec propdivs([*] word p) void: word i, j, max; max := dim(p,1)-1; for i from 0 upto max do p[i] := 0 od; for i from 1 upto max/2 do for j from i*2 by i upto max do p[j] := p[j] + i od odcorp/* Find all amicable pairs between 0 and 20,000 */proc nonrec main() void: word MAX = 20000; word i, j; [MAX] word p; propdivs(p); for i from 1 upto MAX-1 do for j from i+1 upto MAX-1 do if p[i]=j and p[j]=i then writeln(i:5, ", ", j:5) fi od odcorp
220, 284 1184, 1210 2620, 2924 5020, 5564 6232, 636810744, 1085612285, 1459517296, 18416
func sumdivs n . sum = 1 for d = 2 to sqrt n if n mod d = 0 sum += d + n div d . . return sum.for n = 1 to 20000 m = sumdivs n if m > n if sumdivs m = n print n & " " & m . ..
;; using (sum-divisors) from math.lib(lib'math)(define(amicableN)(definen0)(for/list((m(in-range2N)))(set!n(sum-divisorsm))#:continue(>=n(*1.5m));; assume n/m < 1.5#:continue(<=nm);; prevent perfect numbers#:continue(!=(sum-divisorsn)m)(consmn)))(amicable20000)→((220.284)(1184.1210)(2620.2924)(5020.5564)(6232.6368)(10744.10856)(12285.14595)(17296.18416))(amicable1_000_000);; 42 pairs→(...(802725.863835)(879712.901424)(898216.980984)(947835.1125765)(998104.1043096))
open monad io number listdivisors n = filter ((0 ==) << (n `mod`)) [1..(n `div` 2)]range = [1 .. 20000]divs = zip range $ map (sum << divisors) rangepairs = [(n, m) \\ (n, nd) <- divs, (m, md) <- divs | n < m && nd == m && md == n]do putLn pairs ::: IO
[(220,284),(1184,1210),(2620,2924),(5020,5564),(6232,6368),(10744,10856),(12285,14595),(17296,18416)]
ELENA 6.x :
import extensions;import system'routines; const int N = 20000; extension op{ ProperDivisors = Range.new(1,self / 2).filterBy::(n => self.mod(n) == 0); get AmicablePairs() { var divsums := Range .new(0, self + 1) .selectBy::(i => i.ProperDivisors.summarize(Integer.new())) .toArray(); ^ 1.repeatTill(divsums.Length) .filterBy::(i) { var ii := i; var sum := divsums[i]; ^ (i < sum) && (sum < divsums.Length) && (divsums[sum] == i) } .selectBy::(i => new { Item1 = i; Item2 = divsums[i]; }) }} public program(){ N.AmicablePairs.forEach::(pair) { console.printLine(pair.Item1, " ", pair.Item2) }}220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
import extensions;import system'routines'stex;import system'collections; const int N = 20000; extension op : IntNumber{ Enumerator<int> ProperDivisors = new Range(1,self / 2).filterBy::(int n => self.mod(n) == 0); get AmicablePairs() { auto divsums := new List<int>( cast Enumerator<int>( new Range(0, self).selectBy::(int i => i.ProperDivisors.summarize(0)))); ^ new Range(0, divsums.Length) .filterBy::(int i) { auto sum := divsums[i]; ^ (i < sum) && (sum < divsums.Length) && (divsums[sum] == i) } .selectBy::(int i => new Tuple<int,int>(i,divsums[i])); }} public program(){ N.AmicablePairs.forEach::(var Tuple<int,int> pair) { console.printLine(pair.Item1, " ", pair.Item2) }}import extensions;import system'routines'stex;import system'collections; const int Limit = 20000;singleton ProperDivisors{ Enumerator<int> function(int number) = Range.new(1, number / 2).filterBy::(int n => number.mod(n) == 0);} public sealed AmicablePairs{ int max; constructor(int max) { this max := max } yieldable Tuple<int, int> next() { List<int> divsums := Range.new(0, max + 1).selectBy::(int i => ProperDivisors(i).summarize(0)); for (int i := 1; i < divsums.Length; i += 1) { int sum := divsums[i]; if(i < sum && sum <= divsums.Length && divsums[sum] == i) { $yield new Tuple<int, int>(i, sum); } }; ^ nil } } public program(){ auto e := new AmicablePairs(Limit); for(auto pair := e.next(); pair != nil) { console.printLine(pair.Item1, " ", pair.Item2) }}220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
Withproper_divisors#Elixir in place:
defmoduleProperdodefdivisors(1),do:[]defdivisors(n),do:[1|divisors(2,n,:math.sqrt(n))]|>Enum.sortdefpdivisors(k,_n,q)whenk>q,do:[]defpdivisors(k,n,q)whenrem(n,k)>0,do:divisors(k+1,n,q)defpdivisors(k,n,q)whenk*k==n,do:[k|divisors(k+1,n,q)]defpdivisors(k,n,q),do:[k,div(n,k)|divisors(k+1,n,q)]endmap=Map.new(1..20000,fnn->{n,Proper.divisors(n)|>Enum.sum}end)Enum.filter(map,fn{n,sum}->map[sum]==nandn<sumend)|>Enum.sort|>Enum.each(fn{i,j}->IO.puts"#{i} and#{j}"end)
220 and 2841184 and 12102620 and 29245020 and 55646232 and 636810744 and 1085612285 and 1459517296 and 18416
Very slow solution. Same functions by and large as in proper divisors and co.
-module(properdivs).-export([amicable/1,divs/1,sumdivs/1]).amicable(Limit)->amicable(Limit,[],3,2).amicable(Limit,List,_Current,Acc)whenAcc>=Limit->List;amicable(Limit,List,Current,Acc)whenCurrent=<Acc/2->amicable(Limit,List,Acc,Acc+1);amicable(Limit,List,Current,Acc)->CS=sumdivs(Current),AS=sumdivs(Acc),ifCS==AccandalsoAS==CurrentandalsoAcc=/=Current->io:format("A:~w, B:~w,~nL:~w~w~n",[Current,Acc,divs(Current),divs(Acc)]),NL=List++[{Current,Acc}],amicable(Limit,NL,Acc+1,Acc+1);true->amicable(Limit,List,Current-1,Acc)end.divs(0)->[];divs(1)->[];divs(N)->lists:sort(divisors(1,N)).divisors(1,N)->[1]++divisors(2,N,math:sqrt(N)).divisors(K,_N,Q)whenK>Q->[];divisors(K,N,_Q)whenNremK=/=0->[]++divisors(K+1,N,math:sqrt(N));divisors(K,N,_Q)whenK*K==N->[K]++divisors(K+1,N,math:sqrt(N));divisors(K,N,_Q)->[K,NdivK]++divisors(K+1,N,math:sqrt(N)).sumdivs(N)->lists:sum(divs(N)).
3> properdivs:amicable(20000).A: 220, B: 284, L: [1,2,4,5,10,11,20,22,44,55,110][1,2,4,71,142]A: 1184, B: 1210, L: [1,2,4,8,16,32,37,74,148,296,592][1,2,5,10,11,22,55,110,121,242,605]A: 2620, B: 2924, L: [1,2,4,5,10,20,131,262,524,655,1310][1,2,4,17,34,43,68,86,172,731,1462]A: 5020, B: 5564, L: [1,2,4,5,10,20,251,502,1004,1255,2510][1,2,4,13,26,52,107,214,428,1391,2782]A: 6232, B: 6368, L: [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]A: 10744, B: 10856, L: [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]A: 12285, B: 14595, L: [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]A: 17296, B: 18416, L: [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][{220,284}, {1184,1210}, {2620,2924}, {5020,5564}, {6232,6368}, {10744,10856}, {12285,14595}, {17296,18416}]This is lazy AND depends on the fun fact that we're not really identifying pairs. They just happen to order.Probably, this answer is false in some sense. But a good deal faster :) As above with the additional function.
[See the talk section amicable pairs, out of order for this Rosetta Code task.]
friendly(Limit)->List=[{X,properdivs:sumdivs(X)}||X<-lists:seq(3,Limit)],Final=[X||X<-lists:seq(3,Limit),X==properdivs:sumdivs(proplists:get_value(X,List))andalsoX=/=proplists:get_value(X,List)],io:format("L:~w~n",[Final]).
45> properdivs:friendly(20000).L: [220,284,1184,1210,2620,2924,5020,5564,6232,6368,10744,10856,12285,14595,17296,18416]ok
We might answer a challenge by saying:
friendly(Limit)->List=[{X,properdivs:sumdivs(X)}||X<-lists:seq(3,Limit)],Final=[X||X<-lists:seq(3,Limit),X==properdivs:sumdivs(proplists:get_value(X,List))andalsoX=/=proplists:get_value(X,List)],findfriendlies(Final,[]).findfriendlies(List,Acc)whenlength(List)=<0->Acc;findfriendlies(List,Acc)->A=lists:nth(1,List),AS=sumdivs(A),B=lists:nth(2,List),BS=sumdivs(B),ifAS==BandalsoBS==A->{_,BL}=lists:split(2,List),findfriendlies(BL,Acc++[{A,B}]);true->falseend.
94> properdivs:friendly(20000).[{220,284}, {1184,1210}, {2620,2924}, {5020,5564}, {6232,6368}, {10744,10856}, {12285,14595}, {17296,18416}]In either case, it's a lot faster than the recursion in my first example.
PROGRAM AMICABLECONST LIMIT=20000PROCEDURE SUMPROP(NUM->M) IF NUM<2 THEN M=0 EXIT PROCEDURE SUM=1 ROOT=SQR(NUM) FOR I=2 TO ROOT-1 DO IF (NUM=I*INT(NUM/I)) THEN SUM=SUM+I+NUM/I END IF IF (NUM=ROOT*INT(NUM/ROOT)) THEN SUM=SUM+ROOT END IF END FOR M=SUMEND PROCEDUREBEGIN PRINT(CHR$(12);) ! CLS PRINT("Amicable pairs < ";LIMIT) FOR N=1 TO LIMIT DO SUMPROP(N->M1) SUMPROP(M1->M2) IF (N=M2 AND N<M1) THEN PRINT(N,M1) END FOREND PROGRAMAmicable pairs < 20000 220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416
[2..20000-1]|>List.map(funn->n,([1..n/2]|>List.filter(funx->n%x=0)|>List.sum))|>List.map(fun(a,b)->ifa<bthen(a,b)else(b,a))|>List.groupByid|>List.mapsnd|>List.filter(List.length>>((=)2))|>List.mapList.head|>List.iter(printfn"%A")
(220, 284)(1184, 1210)(2620, 2924)(5020, 5564)(6232, 6368)(10744, 10856)(12285, 14595)(17296, 18416)
This solution focuses on the language's namesake: factoring code into small words which are subsequently composed to form more powerful — yet just as simple — words. Using this approach, the final word naturally arrives at the solution. This is often referred to as the bottom-up approach, which is a way in which Factor (and other concatenative languages) commonly differs from other languages.
USING:groupingmath.primes.factorsmath.ranges;:pdivs(n--seq)divisorsbut-last;:dsum(n--sum)pdivssum;:dsum=(nm--?)dsum=;:both-dsum=(nm--?)[dsum=][swapdsum=]2biand;:amicable?(nm--?)[both-dsum=][=not]2biand;:drange(--seq)2 20000[a,b);:dsums(--seq)drange[dsum]map;:is-am?-seq(--seq)dsumsdrange[amicable?]2map;:am-nums(--seq)tis-am?-seqindices;:am-nums-c(--seq)am-nums[2+]map;:am-pairs(--seq)am-nums-c2group;:print-am(--)am-pairs[>array.]each;print-am
{ 220 284 }{ 1184 1210 }{ 2620 2924 }{ 5020 5564 }{ 6232 6368 }{ 10744 10856 }{ 12285 14595 }{ 17296 18416 }Calculate many times the divisors.
:proper-divisors( n -- 1..n )dup2/1+1?dodupimod0=ifiswapthenloopdrop;:divisors-sum( 1..n -- n )dup1=ifexitthenbeginover+swap1=until;:pair( n -- n )dup1=ifexitthenproper-divisorsdivisors-sum;:?paired( n -- t | f )duppair2duppair=>r<r>and;:amicable-list1+1doi?pairedifcri.ipair.thenloop;20000amicable-list
220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416 ok
Use memory to store sum of divisors, a little quicker.
variableamicable-table:proper-divisors( n -- 1..n )dup1=ifexitthen( not really but useful )dup2/1+1?dodupimod0=ifiswapthenloopdrop;:divisors-sum( 1..n -- n )dup1=ifexitthenbeginover+swap1=until;:build-amicable-tablehereamicable-table!1+dup,1doiproper-divisorsdivisors-sum,loop;:pairedcellsamicable-table@+@;:.amicablesamicable-table@@1doipairedpairedi=ipairedi>andifcri.ipaired.thenloop;:amicable-listbuild-amicable-table.amicables;20000amicable-list
220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416 ok
This version uses some latter-day facilities such as array assignment that could be replaced by an ordinary DO-loop, as could the FOR ALL statement that for two adds two to every second element, for three adds three to every third, etc. Each FORALL statement applies its DO-given increment to all the selected array elements potentially in any order or even simultaneously. Likewise, the "MODULE" protocol could be abandoned, which would mean that the KNOWNSUM array would have to be declared COMMON for access across routines - or the whole re-written as a single mainline. And if the PARAMETER statements were replaced appropriately, this source could be compiled using Fortran 77.
Output:
Perfect!! 6Perfect!! 28Amicable! 220 284Perfect!! 496Amicable! 1184 1210Amicable! 2620 2924Amicable! 5020 5564Amicable! 6232 6368Perfect!! 8128Amicable! 10744 10856Amicable! 12285 14595Amicable! 17296 18416
MODULEFACTORSTUFF!This protocol evades the need for multiple parameters, or COMMON, or one shapeless main line...ConcoctedbyR.N.McLean,MMXV.INTEGERLOTS,ILIMIT!Some bounds.PARAMETER(ILIMIT=2147483647)!Computer arithmetic is not with real numbers.PARAMETER(LOTS=22000)!Nor is computer storage infinite.INTEGERKNOWNSUM(LOTS)!Calculate these once as multiple references are expected.CONTAINS!Assistants.INTEGERFUNCTIONSUMF(N)!Sum of the proper divisors of N.INTEGERN!The number in question.INTEGERS,F,F2,INC,BOOST!Assistants.IF(N.LE.LOTS)THEN!If we're within reach,SUMF=KNOWNSUM(N)!The result is to hand.ELSE!Otherwise, some on-the-spot effort ensues.CoulduseSUMFinplaceofS,butsomecompilershavebeenconfusedbysuchusage.S=1!1 is always a factor of N, but N is deemed not.F=1!Prepare a crude search for factors.INC=1!One by plodding one.IF(MOD(N,2).EQ.1)INC=2!Ah, but an odd number cannot have an even number as a divisor.1F=F+INC!So half the time we can doubleplod.F2=F*F!Up to F2 < N rather than F < SQRT(N) and worries over inexact arithmetic.IF(F2.LT.N)THEN!F2 = N handled below.IF(MOD(N,F).EQ.0)THEN!Does F divide N?BOOST=F+N/F!Yes. The divisor and its counterpart.IF(S.GT.ILIMIT-BOOST)GOTO666!Would their augmentation cause an overflow?S=S+BOOST!No, so count in the two divisors just discovered.END IF!So much for a divisor discovered.GOTO1!Try for another.END IF!So much for the horde.IF(F2.EQ.N)THEN!Special case: N may be a perfect square, not necessarily of a prime number.IF(S.GT.ILIMIT-F)GOTO666!It is. And it too might cause overflow.S=S+F!But if not, count F once only.END IF!All done.SUMF=S!This is the result.END IF!Whichever way obtained,RETURN!Done.Cannotcalculatethesum,becauseitexceedstheintegerlimit.666SUMF=-666!An expression of dismay that the caller will notice.END FUNCTIONSUMF!Alternatively, find the prime factors, and combine them...SUBROUTINEPREPARESUMF!Initialise the KNOWNSUM array.ConverttheSieveofEratoshenestohaveeachslotcontainthesumoftheproperdivisorsofitsslotnumber.Changestoinsteadcountthenumberoffactors,orprimefactors,etc.wouldbesimpleenough.INTEGERF!A factor for numbers such as 2F, 3F, 4F, 5F, ...KNOWNSUM(1)=0!Proper divisors of N do not include N.KNOWNSUM(2:LOTS)=1!So, although 1 is a proper divisor of all N, 1 is excluded for itself.DOF=2,LOTS/2!Step through all the possible divisors of numbers not exceeding LOTS.FORALL(I=F+F:LOTS:F)KNOWNSUM(I)=KNOWNSUM(I)+F!And augment each corresponding slot.END DO!Different divisors can hit the same slot. For instance, 6 by 2 and also by 3.END SUBROUTINEPREPARESUMF!Could alternatively generate all products of prime numbers.END MODULEFACTORSTUFF!Enough assistants.PROGRAMAMICABLE!Seek N such that SumF(SumF(N)) = N, for N up to 20,000.USEFACTORSTUFF!This should help.INTEGERI,N!Steppers.INTEGERS1,S2!Sums of factors.CALLPREPARESUMF!Values for every N up to the search limit will be called for at least once.cWRITE(6,66)(I,KNOWNSUM(I),I=1,48)c66FORMAT(10(I3,":",I5,"|"))DON=2,20000!Step through the specified search space.S1=SUMF(N)!Only even numbers appear in the results, but check every one anyway.IF(S1.EQ.N)THEN!Catch a tight loop.WRITE(6,*)"Perfect!!",N!Self amicable! Would otherwise appear as Amicable! n,n.ELSE IF(S1.GT.N)THEN!Look for a pair going upwards only.S2=SUMF(S1)!Since otherwise each would appear twice.IF(S2.EQ.N)WRITE(6,*)"Amicable!",N,S1!Aha!END IF!So much for that candidate.END DO!On to the next.END!Done.
' FreeBASIC v1.05.0 win64FunctionSumProperDivisors(numberAsInteger)AsIntegerIfnumber<2ThenReturn0DimsumAsInteger=0ForiAsInteger=1Tonumber\2IfnumberModi=0Thensum+=iNextReturnsumEndFunctionDimAsIntegern,fDimAsIntegersum(19999)Forn=1To19999sum(n)=SumProperDivisors(n)NextPrint"The pairs of amicable numbers below 20,000 are :"PrintForn=1To19998' f = SumProperDivisors(n)f=sum(n)Iff<=nOrElsef<1OrElsef>19999ThenContinueForIff=sum(n)AndAlson=sum(f)ThenPrintUsing"#####";n;Print" and ";Using"#####";sum(n)EndIfNextPrintPrint"Press any key to exit the program"SleepEnd
The pairs of amicable numbers below 20,000 are : 220 and 284 1184 and 1210 2620 and 2924 5020 and 5564 6232 and 636810744 and 1085612285 and 1459517296 and 18416
' version 04-10-2016' compile with: fbc -s console' replaced the function with 2 FOR NEXT loops#Definemax20000' test for pairs below max#Definemax_1max-1DimAsStringu_str=String(Len(Str(max))+1,"#")DimAsUIntegern,fDimSharedAsUIntegersum(max_1)Forn=2Tomax_1sum(n)=1NextForn=2Tomax_1\2Forf=n*2Tomax_1Stepnsum(f)+=nNextNextPrintPrintUsing" The pairs of amicable numbers below"&u_str&", are :";maxPrintForn=1Tomax_1-1f=Sum(n)Iff<=nOrElsef>maxThenContinueForIff=sum(n)AndAlson=sum(f)ThenPrintUsingu_str&" and"&u_str;n;fEndIfNext' empty keyboard bufferWhileInkey<>"":WendPrint:Print:Print" Hit any key to end program"SleepEnd
The pairs of amicable numbers below 20,000 are : 220 and 284 1184 and 1210 2620 and 2924 5020 and 5564 6232 and 6368 10744 and 10856 12285 and 14595 17296 and 18416
This example uses Frink's built-in efficient factorization algorithms. It can work for arbitrarily large numbers.
n = 1seen = new setdo{ n = n + 1 if seen.contains[n] next sum = sum[allFactors[n, true, false, false]] if sum != n and sum[allFactors[sum, true, false, false]] == n { println["$n, $sum"] seen.put[sum] }} while n <= 20000220, 2841184, 12102620, 29245020, 55646232, 636810744, 1085612285, 1459517296, 18416
| This example does not show the output mentioned in the task descriptionon this page (or a page linked to from here). Please ensure that it meets all task requirements and remove this message. Note that phrases in task descriptions such as"print and display" and"print and show" for example, indicate that (reasonable length) output be a part of a language's solution. |
This program is much too parallel and manifests all the pairs, which requires a giant amount of memory.
fun divisors(n: int): []int = filter (fn x => n%x == 0) (map (1+) (iota (n/2)))fun amicable((n: int, nd: int), (m: int, md: int)): bool = n < m && nd == m && md == nfun getPair (divs: [upper](int, int)) (flat_i: int): ((int,int), (int,int)) = let i = flat_i / upper let j = flat_i % upper in unsafe (divs[i], divs[j])fun main(upper: int): [][2]int = let range = map (1+) (iota upper) let divs = zip range (map (fn n => reduce (+) 0 (divisors n)) range) let amicable = filter amicable (map (getPair divs) (iota (upper*upper))) in map (fn (np,mp) => [#1 np, #1 mp]) amicable
local fn Sigma( n as long ) as long long i, root, sum = 1 if n == 1 then exit fn = 0 root = sqr(n) for i = 2 to root if ( n mod i == 0 ) then sum += i + n/i next if root * root == n then sum -= rootend fn = sumvoid local fn CalculateAmicablePairs( limit as long ) long i, m printf @"\nAmicable pairs through %ld are:\n", limit for i = 2 to limit m = fn Sigma(i) if ( m > i ) if ( fn Sigma(m) == i ) then printf @"%6ld and %ld", i, m end if nextend fnCFTimeInterval tt = fn CACurrentMediaTimefn CalculateAmicablePairs( 20000 )printf @"\nCompute time: %.3f ms",(fn CACurrentMediaTime-t)*1000HandleEvents
Amicable pairs through 20000 are: 220 and 284 1184 and 1210 2620 and 2924 5020 and 5564 6232 and 6368 10744 and 10856 12285 and 14595 17296 and 18416Compute time: 28.701 ms
OPENW 1CLEARW 1'DIM f%(20001) ! sum of proper factors for each nFOR i%=1 TO 20000 f%(i%)=@sum_proper_divisors(i%)NEXT i%' look for pairsFOR i%=1 TO 20000 FOR j%=i%+1 TO 20000 IF f%(i%)=j% AND i%=f%(j%) PRINT "Amicable pair ";i%;" ";j% ENDIF NEXT j%NEXT i%'PRINTPRINT "-- found all amicable pairs"~INP(2)CLOSEW 1'' Compute the sum of proper divisors of given number'FUNCTION sum_proper_divisors(n%) LOCAL i%,sum%,root% ' IF n%>1 ! n% must be 2 or larger sum%=1 ! start with 1 root%=SQR(n%) ! note that root% is an integer ' check possible factors, up to sqrt FOR i%=2 TO root% IF n% MOD i%=0 sum%=sum%+i% ! i% is a factor IF i%*i%<>n% ! check i% is not actual square root of n% sum%=sum%+n%/i% ! so n%/i% will also be a factor ENDIF ENDIF NEXT i% ENDIF RETURN sum%ENDFUNC
Output is:
Amicable pair: 220 284Amicable pair: 1184 1210Amicable pair: 2620 2924Amicable pair: 5020 5564Amicable pair: 6232 6368Amicable pair: 10744 10856Amicable pair: 12285 14595Amicable pair: 17296 18416-- found all amicable pairs
packagemainimport"fmt"funcpfacSum(iint)int{sum:=0forp:=1;p<=i/2;p++{ifi%p==0{sum+=p}}returnsum}funcmain(){vara[20000]intfori:=1;i<20000;i++{a[i]=pfacSum(i)}fmt.Println("The amicable pairs below 20,000 are:")forn:=2;n<19999;n++{m:=a[n]ifm>n&&m<20000&&n==a[m]{fmt.Printf(" %5d and %5d\n",n,m)}}}
The amicable pairs below 20,000 are: 220 and 284 1184 and 1210 2620 and 2924 5020 and 5564 6232 and 6368 10744 and 10856 12285 and 14595 17296 and 18416
divisors::(Integrala)=>a->[a]divisorsn=filter((0==).(n`mod`))[1..(n`div`2)]main::IO()main=doletrange=[1..20000::Int]divs=ziprange$map(sum.divisors)rangepairs=[(n,m)|(n,nd)<-divs,(m,md)<-divs,n<m,nd==m,md==n]printpairs
[(220,284),(1184,1210),(2620,2924),(5020,5564),(6232,6368),(10744,10856),(12285,14595),(17296,18416)]
Or, deriving proper divisors above the square root as cofactors (for better performance)
importData.Bool(bool)amicablePairsUpTo::Int->[(Int,Int)]amicablePairsUpTon=letsigma=sum.properDivisorsin[1..n]>>=(\x->lety=sigmaxinbool[][(x,y)](x<y&&x==sigmay))properDivisors::Integrala=>a->[a]properDivisorsn=letroot=(floor.sqrt)(fromIntegraln::Double)lows=filter((0==).remn)[1..root]ininit$lows++drop(bool01(root*root==n))(reverse(quotn<$>lows))main::IO()main=mapM_print$amicablePairsUpTo20000
(220,284)(1184,1210)(2620,2924)(5020,5564)(6232,6368)(10744,10856)(12285,14595)(17296,18416)
Proper Divisor implementation:
factors=:[:/:~@,*/&>@{@((^i.@>:)&.>/)@q:~&__properDivisors=:factors-.-.&1
(properDivisors is all factors except "the number itself when that number is greater than 1".)
Amicable pairs:
1+($#:I.@,)(</~@i.@#*(*|:))(=/+/@properDivisors@>)1+i.2000022028411841210262029245020556462326368107441085612285145951729618416
Explanation: We generate sequence of positive integers, starting from 1, and compare each of them to the sum of proper divisors of each of them. Then we fold this comparison diagonally, keeping only the values where the comparison was equal both ways and the smaller value appears before the larger value. Finally, indices into true values give us our amicable pairs.
importjava.util.Map;importjava.util.function.Function;importjava.util.stream.Collectors;importjava.util.stream.LongStream;publicclassAmicablePairs{publicstaticvoidmain(String[]args){intlimit=20_000;Map<Long,Long>map=LongStream.rangeClosed(1,limit).parallel().boxed().collect(Collectors.toMap(Function.identity(),AmicablePairs::properDivsSum));LongStream.rangeClosed(1,limit).forEach(n->{longm=map.get(n);if(m>n&&m<=limit&&map.get(m)==n)System.out.printf("%s %s %n",n,m);});}publicstaticLongproperDivsSum(longn){returnLongStream.rangeClosed(1,(n+1)/2).filter(i->n%i==0).sum();}}
220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416
(function(max){// Proper divisorsfunctionproperDivisors(n){if(n<2)return[];else{varrRoot=Math.sqrt(n),intRoot=Math.floor(rRoot),lows=range(1,intRoot).filter(function(x){return(n%x)===0;});returnlows.concat(lows.slice(1).map(function(x){returnn/x;}).reverse().slice((rRoot===intRoot)|0));}}// [m..n]functionrange(m,n){vara=Array(n-m+1),i=n+1;while(i--)a[i-1]=i;returna;}// Filter an array of proper divisor sums,// reading the array index as a function of N (N-1)// and the sum of proper divisors as a potential Mvarpairs=range(1,max).map(function(x){returnproperDivisors(x).reduce(function(a,d){returna+d;},0)}).reduce(function(a,m,i,lst){varn=i+1;return(m>n)&&lst[m-1]===n?a.concat([[n,m]]):a;},[]);// [[a]] -> bool -> s -> sfunctionwikiTable(lstRows,blnHeaderRow,strStyle){return'{| class="wikitable" '+(strStyle?'style="'+strStyle+'"':'')+lstRows.map(function(lstRow,iRow){varstrDelim=((blnHeaderRow&&!iRow)?'!':'|');return'\n|-\n'+strDelim+' '+lstRow.map(function(v){returntypeofv==='undefined'?' ':v;}).join(' '+strDelim+strDelim+' ');}).join('')+'\n|}';}returnwikiTable([['N','M']].concat(pairs),true,'text-align:center')+'\n\n'+JSON.stringify(pairs);})(20000);
| N | M |
|---|---|
| 220 | 284 |
| 1184 | 1210 |
| 2620 | 2924 |
| 5020 | 5564 |
| 6232 | 6368 |
| 10744 | 10856 |
| 12285 | 14595 |
| 17296 | 18416 |
[[220,284],[1184,1210],[2620,2924],[5020,5564],[6232,6368],[10744,10856],[12285,14595],[17296,18416]]
(()=>{'use strict';// amicablePairsUpTo :: Int -> [(Int, Int)]constamicablePairsUpTo=n=>{constsigma=compose(sum,properDivisors);returnenumFromTo(1)(n).flatMap(x=>{consty=sigma(x);returnx<y&&x===sigma(y)?([[x,y]]):[];});};// properDivisors :: Int -> [Int]constproperDivisors=n=>{constrRoot=Math.sqrt(n),intRoot=Math.floor(rRoot),lows=enumFromTo(1)(intRoot).filter(x=>0===(n%x));returnlows.concat(lows.map(x=>n/x).reverse().slice((rRoot===intRoot)|0,-1));};// TEST -----------------------------------------------// main :: IO ()constmain=()=>console.log(unlines(amicablePairsUpTo(20000).map(JSON.stringify)));// GENERIC FUNCTIONS ----------------------------------// compose (<<<) :: (b -> c) -> (a -> b) -> a -> cconstcompose=(...fs)=>x=>fs.reduceRight((a,f)=>f(a),x);// enumFromTo :: Int -> Int -> [Int]constenumFromTo=m=>n=>Array.from({length:1+n-m},(_,i)=>m+i);// sum :: [Num] -> Numconstsum=xs=>xs.reduce((a,x)=>a+x,0);// unlines :: [String] -> Stringconstunlines=xs=>xs.join('\n');// MAIN ---returnmain();})();
[220,284][1184,1210][2620,2924][5020,5564][6232,6368][10744,10856][12285,14595][17296,18416]
# unordereddef 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 emptyend) 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)
$jq-c-n-famicable_pairs.jq220and284areamicable1184and1210areamicable2620and2924areamicable5020and5564areamicable6232and6368areamicable10744and10856areamicable12285and14595areamicable17296and18416areamicable
Givenfactor, it is not necessary to calculate the individual divisors to compute their sum. SeeAbundant, deficient and perfect number classifications for the details.It is safe to exclude primes from consideration; their proper divisor sum is always 1. This code also 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.
usingPrimes,Printffunctionpcontrib(p::Int64,a::Int64)n=one(p)pcon=one(p)foriin1:an*=ppcon+=nendreturnpconendfunctiondivisorsum(n::Int64)dsum=one(n)for(p,a)infactor(n)dsum*=pcontrib(p,a)enddsum-=nendfunctionamicables(L=2*10^7)acnt=0println("Amicable pairs not greater than ",L)foriin2:L!isprime(i)||continuej=divisorsum(i)j<i&&divisorsum(j)==i||continueacnt+=1println(@sprintf("%4d",acnt)," => ",j,", ",i)endendamicables()
Note, the output isnot ordered by the first figure, see e.g. counters 11, 15, ..., 139, 141, etc.
Amicable pairs not greater than 20000000 1 => 220, 284 2 => 1184, 1210 3 => 2620, 2924 4 => 5020, 5564 5 => 6232, 6368 6 => 10744, 10856 7 => 12285, 14595 8 => 17296, 18416 9 => 66928, 66992 10 => 67095, 71145 11 => 63020, 76084 12 => 69615, 87633 13 => 79750, 88730 14 => 122368, 123152 15 => 100485, 124155 16 => 122265, 139815[...] 138 => 18655744, 19154336 139 => 16871582, 19325698 140 => 17844255, 19895265 141 => 17754165, 19985355
Using thefactor() function from thePrimes package allows for a quicker calculation, especially when it comes to big numbers. Here we use a busy one-liner with an iterator. The following code prints the amicable pairs inascending order and also prints the sum of the amicable pair and the cumulative sum of all pairs found so far; this allows to check results, when solvingProject Euler problem #21.
usingPrimesfunctionamicable_numbers(max::Integer=200_000_000)functionsum_proper_divisors(n::Integer)sum(vec(map(prod,Iterators.product((p.^(0:m)for(p,m)infactor(n))...))))-nendcount=0cumsum=0println("count, a, b, a+b, Sum(a+b)")forain2:maxisprime(a)&&continueb=sum_proper_divisors(a)ifa<b&&sum_proper_divisors(b)==acount+=1sumab=a+bcumsum+=sumabprintln("$count,$a,$b,$sumab,$cumsum")endendendamicable_numbers()
count, a, b, a+b, Sum(a+b)1, 220, 284, 504, 5042, 1184, 1210, 2394, 28983, 2620, 2924, 5544, 84424, 5020, 5564, 10584, 190265, 6232, 6368, 12600, 316266, 10744, 10856, 21600, 532267, 12285, 14595, 26880, 801068, 17296, 18416, 35712, 1158189, 63020, 76084, 139104, 25492210, 66928, 66992, 133920, 38884211, 67095, 71145, 138240, 52708212, 69615, 87633, 157248, 68433013, 79750, 88730, 168480, 85281014, 100485, 124155, 224640, 107745015, 122265, 139815, 262080, 133953016, 122368, 123152, 245520, 1585050[...]300, 189406984, 203592056, 392999040, 31530421032301, 190888155, 194594085, 385482240, 31915903272302, 195857415, 196214265, 392071680, 32307974952303, 196421715, 224703405, 421125120, 32729100072304, 199432948, 213484172, 412917120, 33142017192
propdivs:{1+&0=x!'1+!x%2}(8,2)#v@&{(x=+/propdivs[a])&~x=a:+/propdivs[x]}'v:1+!20000(22028411841210262029245020556462326368107441085612285145951729618416)
// version 1.1funsumProperDivisors(n:Int):Int{if(n<2)return0return(1..n/2).filter{(n%it)==0}.sum()}funmain(args:Array<String>){valsum=IntArray(20000,{sumProperDivisors(it)})println("The pairs of amicable numbers below 20,000 are:\n")for(nin2..19998){valm=sum[n]if(m>n&&m<20000&&n==sum[m]){println(n.toString().padStart(5)+" and "+m.toString().padStart(5))}}}
The pairs of amicable numbers below 20,000 are: 220 and 284 1184 and 1210 2620 and 2924 5020 and 5564 6232 and 636810744 and 1085612285 and 1459517296 and 18416
Avoids unnecessary divisor sum calculations.
Runs in around 0.11 seconds on TIO.RUN.
functionsumDivs(n)localsum=1ford=2,math.sqrt(n)doifn%d==0thensum=sum+dsum=sum+n/dendendreturnsumendforn=2,20000dom=sumDivs(n)ifm>nthenifsumDivs(m)==nthenprint(n,m)endendend
22028411841210262029245020556462326368107441085612285145951729618416
Alternative version using a table of proper divisors, constructed without division/modulo.
Runs is around 0.02 seconds on TIO.RUN.
MAX_NUMBER=20000sumDivs={}-- table of proper divisorsfori=1,MAX_NUMBERdosumDivs[i]=1endfori=2,MAX_NUMBERdoforj=i+i,MAX_NUMBER,idosumDivs[j]=sumDivs[j]+iendendforn=2,MAX_NUMBERdom=sumDivs[n]ifm>nthenifsumDivs[m]==nthenprint(n,m)endendend
22028411841210262029245020556462326368107441085612285145951729618416
// find amicable pairs p1, p2 where each is equal to the other's proper divisor sumMAX_NUMBER=20000pdSum=[1]*(MAX_NUMBER+1)// table of proper divisorsforiinrange(2,MAX_NUMBER)forjinrange(i+i,MAX_NUMBER,i)pdSum[j]+=iendforendfor// find the amicable pairs up to 20 000ap=[]forp1inrange(1,MAX_NUMBER-1)pdSumP1=pdSum[p1]ifpdSumP1>p1andpdSumP1<=MAX_NUMBERandpdSum[pdSumP1]==p1thenprintstr(p1)+" and "+str(pdSumP1)+" are an amicable pair"endifendfor
220 and 284 are an amicable pair1184 and 1210 are an amicable pair2620 and 2924 are an amicable pair5020 and 5564 are an amicable pair6232 and 6368 are an amicable pair10744 and 10856 are an amicable pair12285 and 14595 are an amicable pair17296 and 18416 are an amicable pair
NORMAL MODE IS INTEGER DIMENSION DIVS(20000) PRINT COMMENT $ AMICABLE PAIRS$ R CALCULATE SUM OF DIVISORS OF N INTERNAL FUNCTION(N) ENTRY TO DIVSUM. DS = 0 THROUGH SUMMAT, FOR DIVC=1, 1, DIVC.GE.NSUMMAT WHENEVER N/DIVC*DIVC.E.N, DS = DS+DIVC FUNCTION RETURN DS END OF FUNCTION R CALCULATE SUM OF DIVISORS FOR ALL NUMBERS 1..20000 THROUGH MEMO, FOR I=1, 1, I.GE.20000MEMO DIVS(I) = DIVSUM.(I) R FIND ALL MATCHING PAIRS THROUGH CHECK, FOR I=1, 1, I.GE.20000 THROUGH CHECK, FOR J=1, 1, J.GE.ICHECK WHENEVER DIVS(I).E.J .AND. DIVS(J).E.I, 0 PRINT FORMAT AMI,I,J VECTOR VALUES AMI = $I6,I6*$ END OF PROGRAM
AMICABLE PAIRS 284 220 1210 1184 2924 2620 5564 5020 6368 6232 10856 10744 14595 12285 18416 17296
| This example does not show the output mentioned in the task descriptionon this page (or a page linked to from here). Please ensure that it meets all task requirements and remove this message. Note that phrases in task descriptions such as"print and display" and"print and show" for example, indicate that (reasonable length) output be a part of a language's solution. |
with(NumberTheory):pairs:=[];for i from 1 to 20000 dofor j from i+1 to 20000 dosum1:=SumOfDivisors(j)-j;sum2:=SumOfDivisors(i)-i;if sum1=i and sum2=j and i<>j thenpairs:=[op(pairs),[i,j]];printf("%a", pairs);end if;end do;end do;pairs;amicableQ[n_]:=Module[{sum=Total[Most@Divisors@n]},sum!=n&&n==Total[Most@Divisors@sum]]Grid@Partition[Cases[Range[4,20000],_?(amicableQ@#&)],2]
22028411841210262029245020556462326368107441085612285145951729618416
functionamicableticN=2:1:20000;aN=[];N(isprime(N))=[];%erase prime numbersI=1;a=N(1);b=sum(pd(a));whilelength(N)>1ifa==b%erase perfect numbers;N(N==a)=[];a=N(1);b=sum(pd(a));elseifb<a%the first member of an amicable pair is abundant not defectiveN(N==a)=[];a=N(1);b=sum(pd(a));elseif~ismember(b,N)%the other member was previously erasedN(N==a)=[];a=N(1);b=sum(pd(a));elsec=sum(pd(b));ifa==caN(I,:)=[Iab];I=I+1;N(N==b)=[];elseif~ismember(c,N)%the other member was previously erasedN(N==b)=[];endendN(N==a)=[];a=N(1);b=sum(pd(a));clearcendenddisp(array2table(aN,'Variablenames',{'N','Amicable1','Amicable2'}))tocendfunctionD=pd(x)K=1:ceil(x/2);D=K(~(rem(x,K)));end
N Amicable1 Amicable2 _ _________ _________ 1 220 284 2 1184 1210 3 2620 2924 4 5020 5564 5 6232 6368 6 10744 10856 7 12285 14595 8 17296 18416 Elapsed time is 8.958720 seconds.
Being a novice, I submitted my code to the Nim community for review and received much feedback and advice.They were instrumental in fine-tuning this code for style and readability, I can't thank them enough.
frommathimportsqrtconstN=524_000_000.int32procsumProperDivisors(someNum:int32,chk4less:bool):int32=result=1letmaxPD=sqrt(someNum.float).int32letoffset=someNummod2fordivNumincountup(2+offset,maxPD,1+offset):ifsomeNummoddivNum==0:result+=divNum+someNumdivdivNumifchk4lessandresult>=someNum:return0fornincountdown(N,2):letm=sumProperDivisors(n,true)ifm!=0andn==sumProperDivisors(m,false):echo$n," ",$m
523679536 509379344511419856 491373104514823985 475838415...... ........... .....18416 1729614595 1228510856 107446368 62325564 50202924 26201210 1184284 220
Total number of pairs is 442, on my machine the code takes ~389 minutes to run.
Here's a second version that uses a large amount of memory but runs in 2m32seconds.Again, thanks to the Nim community
frommathimportsqrtconstN=524_000_000.int32varx=newSeq[int32](N+1)foriin2..sqrt(N.float).int32:varp=i*ix[p]+=ivarj=i+iwhile(p+=i;p<=N):j.incx[p]+=jformin4..N:letn=x[m]+1ifn<mandn!=0andm==x[n]+1:echon," ",m
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416..... ........... ......426191535 514780497475838415 514823985509379344 523679536
MODULEAmicablePairs;IMPORTOut;CONSTmax=20000;VARi,j:INTEGER;pd:ARRAYmax+1OFLONGINT;PROCEDUREProperDivisorsSum(n:LONGINT):LONGINT;VARi,sum:LONGINT;BEGINsum:=0;IFn>1THENINC(sum,1);i:=2;WHILE(i<n)DOIF(nMODi)=0THENINC(sum,i)END;INC(i)ENDEND;RETURNsumENDProperDivisorsSum;BEGINFORi:=0TOmaxDOpd[i]:=ProperDivisorsSum(i)END;FORi:=2TOmaxDOFORj:=i+1TOmaxDOIF(pd[i]=j)&(pd[j]=i)THENOut.Char('[');Out.Int(i,0);Out.Char(',');Out.Int(j,0);Out.Char("]");Out.LnENDENDENDENDAmicablePairs.
[220,284][1184,1210][2620,2924][5020,5564][6232,6368][10744,10856][12285,14595][17296,18416]
Using properDivs implementation tasks without optimization (calculating proper divisors sum without returning a list for instance) :
import: mappingInteger method: properDivs -- [] #[ self swap mod 0 == ] self 2 / seq filter ; : amicables| i j | Array new 20000 loop: i [ i properDivs sum dup ->j i <= if continue then j properDivs sum i <> if continue then [ i, j ] over add ];
amicables .[[220, 284], [1184, 1210], [2620, 2924], [5020, 5564], [6232, 6368], [10744, 10856], [12285, 14595], [17296, 18416]]
letrecisqrtn=ifn=1then1elselet_n=isqrt(n-1)in(_n+(n/_n))/2letsum_divsn=letsum=ref1inford=2toisqrtndoif(nmodd)=0thensum:=!sum+(n/d+d);done;!sumlet()=forn=2to20000doletm=sum_divsninif(m>n)thenif(sum_divsm)=nthenPrintf.printf"%d %d\n"nm;done
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
for(x=1,20000,my(y=sigma(x)-x); if(y>x && x == sigma(y)-y,print(x" "y)))
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
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:
ProgramSumOfFactors;usescrt;{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}varoutf:text;constLimit=2147483647;constlots=20000;{This should be much bigger, but problems apply.}varKnownSum:array[1..lots]oflongint;FunctionSumF(N:Longint):Longint;varf,f2,s,ulp:longint;Beginifn<=lotsthenSumF:=KnownSum[N]{Hurrah!}elsebegin{This is really crude...}s:=1;{1 is always a factor, but N is not.}f:=2;f2:=f*f;whilef2<NdobeginifNmodf=0thenbegin{We have a divisor, and its friend.}ulp:=f+(Ndivf);ifs>Limit-ulpthenbeginSumF:=-666;exit;end;s:=s+ulp;end;f:=f+1;f2:=f*f;end;iff2=Nthen{A perfect square gets its factor in once only.}ifs<=Limit-fthens:=s+felsebeginSumF:=-667;exit;end;SumF:=s;end;End;vari,j,l,sf,fs:LongInt;constenuff=666;{Only so much sociability.}vartrail:array[0..enuff]oflongint;BEGINClrScr;WriteLn('Chasing Chains of Sums of Factors of Numbers.');fori:=1tolotsdoKnownSum[i]:=1;{Sigh. KnownSum:=1;}{start summing every divisor }fori:=2tolotsdobeginj:=i+i;Whilej<=lotsdo{Sigh. For j:=i + i:Lots:i do KnownSum[j]:=KnownSum[j] + i;}beginKnownSum[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.');fori:=2tolotsdo{Search.}beginl:=0;sf:=SumF(i);while(sf>i)and(l<enuff)dobeginl:=l+1;trail[l]:=sf;sf:=SumF(sf);end;ifl>=enuffthenwriteln('Rope ran out! ',i);{$IFDEF ShowOverflow}ifsf<0thenwriteln('Overflow with ',i);{$ENDIF}ifi=sfthen{A loop?}begin{Yes. Reveal its members.}trail[0]:=i;{The first.}ifl=0thenwrite('Perfect!! ')elseifl=1thenwrite('Amicable! ')elsewrite('Sociable: ');forj:=0toldoWrite(Trail[j],',');WriteLn;ifl=0thenwrite(outf,'Perfect!! ')elseifl=1thenwrite(outf,'Amicable! ')elsewrite(outf,'Sociable: ');forj:=0toldowrite(outf,Trail[j],',');WriteLn(outf);end;end;Close(outf);END.
a "normal" Version. Nearly fast as perl using nTheory.
programAmicablePairs;{$IFDEF FPC}{$MODE DELPHI}{$H+}{$ELSE}{$APPTYPE CONSOLE}{$ENDIF}usessysutils;constMAX=20000;//MAX = 20*1000*1000;typetValue=LongWord;tpValue=^tValue;tPower=array[0..31]oftValue;tIndex=recordidxI,idxS:Uint64;end;varIndices:array[0..511]oftIndex;//primes up to 65536 enough until 2^32primes:array[0..6542]oftValue;procedureInitPrimes;// sieve of erathosthenes without multiples of 2typetSieve=array[0..(65536-1)div2]ofansichar;varESieve:^tSieve;idx,i,j,p:LongINt;Beginnew(ESieve);fillchar(ESieve^[0],SizeOF(tSieve),#1);primes[0]:=2;idx:=1;//sievingj:=1;p:=2*j+1;repeatifEsieve^[j]=#1thenbegini:=(2*j+2)*j;// i := (sqr(p) -1) div 2;ifi>High(tSieve)thenBREAK;repeatESIeve^[i]:=#0;inc(i,p);untili>High(tSieve);end;inc(j);inc(p,2);untilj>High(tSieve);//collectingFori:=1toHigh(tSieve)doIFEsieve^[i]=#1thenBeginprimes[idx]:=2*i+1;inc(idx);IFidx>High(primes)thenBREAK;end;dispose(Esieve);end;procedureSu_append(n,factor:tValue;varsu:string);varq,p:tValue;beginp:=0;repeatq:=ndivfactor;IFq*factor<>nthenBreak;inc(p);n:=q;untilfalse;IFp>0thenIFp=1thensu:=su+IntToStr(factor)+'*'elsesu:=su+IntToStr(factor)+'^'+IntToStr(p)+'*';end;procedureProperDivs(n:Uint64);//output of prime factorizationvarsu:string;primNo:tValue;p:tValue;beginstr(n:8,su);su:=su+' [';primNo:=0;p:=primes[0];repeatSu_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;procedureAmPairOutput(cnt:tValue);vari:tValue;r_max,r_min,r:double;beginr_max:=1.0;r_min:=16.0;Fori:=0tocnt-1dowithIndices[i]dobeginr:=IdxS/IDxI;writeln(i+1:4,IdxI:16,IDxS:16,' ratio ',r:10:7);IFr<1thenbeginwriteln(i);readln;halt;end;ifr_max<rthenr_max:=relseifr_min>rthenr_min:=r;IFcnt<20thenbeginProperDivs(IdxI);ProperDivs(IdxS);end;end;writeln(' min ratio ',r_min:12:10);writeln(' max ratio ',r_max:12:10);end;procedureSumOFProperDiv(n:tValue;varSumOfProperDivs:tValue);// calculated by prime factorizationvari,q,primNo,Prime,pot:tValue;SumOfDivs:tValue;begini:=N;SumOfDivs:=1;primNo:=0;Prime:=Primes[0];q:=iDIVPrime;repeatifq*Prime=ithenBeginpot:=1;repeati:=q;q:=idivPrime;Pot:=Pot*Prime+1;untilq*Prime<>i;SumOfDivs:=SumOfDivs*pot;end;Inc(primNo);Prime:=Primes[primNo];q:=iDIVPrime;{check if i already prime}ifPrime>qthenbeginprime:=i;q:=1;end;untili=1;SumOfProperDivs:=SumOfDivs-N;end;functionCheck:tValue;const//going backwardsDIV23:array[0..5]ofbyte=//== 5,4,3,2,1,0(1,0,0,0,1,0);vari,s,k,n:tValue;idx:nativeInt;beginn:=0;idx:=3;Fori:=2toMAXdobegin//must be divisble by 2 or 3 ( n < High(tValue) < 1e14 )IFDIV23[idx]=0thenbeginSumOFProperDiv(i,s);//only 24.7...%IFs>ithenBeginSumOFProperDiv(s,k);IFk=ithenbeginWithindices[n]dobeginidxI:=i;idxS:=s;end;inc(n);end;end;end;dec(idx);IFidx<0thenidx:=high(DIV23);end;result:=n;end;varT2,T1:TDatetime;APcnt:tValue;beginInitPrimes;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.
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]
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
programAmicPair;{find amicable pairs in a limited region 2..MAXbeware that >both< numbers must be smaller than MAXthere are 455 amicable pairs up to 524*1000*1000correct 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}usessysutils;typetValue=LongWord;tpValue=^tValue;tDivSum=array[0..0]oftValue;// evil, but dynamic arrays are slowertpDivSum=^tDivSum;tPower=array[0..31]oftValue;tIndex=recordidxI,idxS:tValue;end;varpower,PowerFac:tPower;ds:arrayoftValue;Indices:array[0..511]oftIndex;DivSumField:tpDivSum;MAX:tValue;procedureInit;vari:LongInt;beginDivSumField[0]:=0;Fori:=1toMAXdoDivSumField[i]:=1;end;procedureProperDivs(n:tValue);//Only for output, normally a factorication would dovarsu,so:string;i,q:tValue;beginsu:='1';so:='';i:=2;whilei*i<=ndobeginq:=ndivi;IFq*i-n=0thenbeginsu:=su+','+IntToStr(i);IFq<>ithenso:=','+IntToStr(q)+so;end;inc(i);end;writeln(' [',su+so,']');end;procedureAmPairOutput(cnt:tValue);vari:tValue;r:double;beginr:=1.0;Fori:=0tocnt-1dowithIndices[i]dobeginwriteln(i+1:4,IdxI:12,IDxS:12,' ratio ',IdxS/IDxI:10:7);ifr<IdxS/IDxIthenr:=IdxS/IDxI;IFcnt<20thenbeginProperDivs(IdxI);ProperDivs(IdxS);end;end;writeln(' max ratio ',r:10:4);end;functionCheck:tValue;vari,s,n:tValue;beginn:=0;Fori:=1toMAXdobegin//s = sum of proper divs (I) == sum of divs (I) - Is:=DivSumField^[i];IF(s<=MAX)AND(s>i)AND(DivSumField^[s]=i)thenbeginWithindices[n]dobeginidxI:=i;idxS:=s;end;inc(n);end;end;result:=n;end;ProcedureCalcPotfactor(prim:tValue);//PowerFac[k] = (prim^(k+1)-1)/(prim-1) == Sum (i=0..k) prim^ivark:tValue;Pot,//== prim^kPFac:Int64;beginPot:=prim;PFac:=1;Fork:=0toHigh(PowerFac)dobeginPFac:=PFac+Pot;IF(POT>MAX)thenBREAK;PowerFac[k]:=PFac;Pot:=Pot*prim;end;end;procedureInitPW(prim:tValue);beginfillchar(power,SizeOf(power),#0);CalcPotfactor(prim);end;functionNextPotCnt(p:tValue):tValue;//return the first power <> 0//power == n to base primvari:tValue;beginresult:=0;repeati:=power[result];Inc(i);IFi<pthenBREAKelsebegini:=0;power[result]:=0;inc(result);end;untilfalse;power[result]:=i;end;procedureSieve(prim:tValue);varactNumber,idx:tValue;begin//sieve with "small" primeswhileprim*prim<=MAXdobeginInitPW(prim);Begin//actNumber = actual number = n*primactNumber:=prim;idx:=prim;whileactNumber<=MAXdobegindec(idx);IFidx>0thenDivSumField^[actNumber]*=PowerFac[0]elseBeginDivSumField^[actNumber]*=PowerFac[NextPotCnt(prim)+1];idx:=Prim;end;inc(actNumber,prim);end;end;//next primerepeatinc(prim);untilDivSumField^[prim]=1;//(DivSumField[prim] = 1);end;//sieve with "big" primes, only one factor is possiblewhile2*prim<=MAXdobeginInitPW(prim);BeginactNumber:=prim;idx:=PowerFac[0];whileactNumber<=MAXdobeginDivSumField^[actNumber]*=idx;inc(actNumber,prim);end;end;repeatinc(prim);untilDivSumField^[prim]=1;end;Foridx:=2toMAXdodec(DivSumField^[idx],idx);end;varT2,T1,T0:TDatetime;APcnt:tValue;i:NativeInt;beginMAX:=20000;IFParamCount>0thenMAX:=StrToInt(ParamStr(1));setlength(ds,MAX);DivSumField:=@ds[0];T0:=time;Fori:=1to1doBeginInit;Sieve(2);end;T1:=time;APCnt:=Check;T2:=time;AmPairOutput(APCnt);writeln(APCnt,' amicable pairs til ',MAX);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));setlength(ds,0);{$IFNDEF UNIX}readln;{$ENDIF}end.
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 2000000:00:00.000.... Test with 524*1000*1000 Linux32, FPC 3.0.1, i4330 3.5 Ghz //Win32 swaps first to allocate 2 GB ) 440 475838415 514823985 ratio 1.0819303 441 491373104 511419856 ratio 1.0407974 442 509379344 523679536 ratio 1.0280738 max ratio 1.3537442 amicable pairs til 524000000Time to calc sum of divs 00:00:12.601Time to find amicable pairs 00:00:02.557
PascalABC.NET 3.10.3.3614
functionk(n:Integer):=1.to(ndiv2).where(d->nmodd=0).sum;begin1.to(20000).Where(n->n=k(k(n))).Select(n->beginvarm:=k(n);Result:=(min(n,m),max(n,m));end).where(v->v[0]<v[1]).Distinct.Select(v->$'{v[0]}-{v[1]}').Printlines;end.
220-2841184-12102620-29245020-55646232-636810744-1085612285-1459517296-18416
Not particularly clever, but instant for this example, and does up to 20 million in 11 seconds.
usentheoryqw/divisor_sum/;formy$x(1..20000){my$y=divisor_sum($x)-$x;say"$x $y"if$y>$x&&$x==divisor_sum($y)-$y;}
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
withjavascript_semanticsform=1to20000dointegern=sum(factors(m,-1))ifm<nandm=sum(factors(n,-1))then?{m,n}endifendfor
{220,284}{1184,1210}{2620,2924}{5020,5564}{6232,6368}{10744,10856}{12285,14595}{17296,18416}def sumDivs var n 1 var sum n sqrt 2 swap 2 tolist for var d n d mod not if sum d + n d / + var sum endif endfor sumenddef2 20000 2 tolist for var i i sumDivs var m m i > if m sumDivs i == if i print "\t" print m print nl endif endifendfornl msec print " s" print
<?phpfunctionsumDivs($n){$sum=1;for($d=2;$d<=sqrt($n);$d++){if($n%$d==0)$sum+=$n/$d+$d;}return$sum;}for($n=2;$n<20000;$n++){$m=sumDivs($n);if($m>$n){if(sumDivs($m)==$n)echo$n." ".$m."<br />";}}?>
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
Different approaches to solve this task:
* foreach loop (two variants) * list comprehension * while loop.
Also, the calculation of the sum of divisors is tabled (the table is cleared between each run).
go => N = 20000, println(amicable1), time(amicable1(N)), % initialize_table is needed to clear the table cache % of sum_divisors/1 between each run. initialize_table, println(amicable2), time(amicable2(N)), initialize_table, println(amicable3), time(amicable3(N)), initialize_table, println(amicable4), time(amicable4(N)), nl.% Foreach loop and a map (hash table)amicable1(N) => Pairs = new_map(), foreach(A in 1..N) B = sum_divisors(A), C = sum_divisors(B), if A != B, A == C then Pairs.put([A,B].sort(),1) end end, println(Pairs.keys().sort()).% List comprehensionamicable2(N) => println([[A,B].sort() : A in 1..N, B = sum_divisors(A), C = sum_divisors(B), A != B, A == C].remove_dups()).% While loopamicable3(N) => A = 1, while(A <= N) B = sum_divisors(A), if A < B, A == sum_divisors(B) then print([A,B]), print(" ") end, A := A + 1 end, nl.% Foreach loop, everything in the conditionamicable4(N) => foreach(A in 1..N, B = sum_divisors(A), A < B, A == sum_divisors(B)) print([A,B]), print(" ") end, nl.%% Sum of divisors of N%tablesum_divisors(N) = Sum => sum_divisors(2,N,1,Sum).% Base case: exceeding the limitsum_divisors(I,N,Sum0,Sum), I > floor(sqrt(N)) => Sum = Sum0.% I is a divisor of Nsum_divisors(I,N,Sum0,Sum), N mod I == 0 => Sum1 = Sum0 + I, (I != N div I -> Sum2 = Sum1 + N div I ; Sum2 = Sum1 ), sum_divisors(I+1,N,Sum2,Sum).% I is not a divisor of N.sum_divisors(I,N,Sum0,Sum) => sum_divisors(I+1,N,Sum0,Sum).amicable1[[220,284],[1184,1210],[2620,2924],[5020,5564],[6232,6368],[10744,10856],[12285,14595],[17296,18416]]CPU time 0.114 seconds.amicable2[[220,284],[1184,1210],[2620,2924],[5020,5564],[6232,6368],[10744,10856],[12285,14595],[17296,18416]]CPU time 0.106 seconds.amicable3[220,284] [1184,1210] [2620,2924] [5020,5564] [6232,6368] [10744,10856] [12285,14595] [17296,18416] CPU time 0.111 seconds.amicable4[220,284] [1184,1210] [2620,2924] [5020,5564] [6232,6368] [10744,10856] [12285,14595] [17296,18416] CPU time 0.107 seconds.
(de accud (Var Key) (if (assoc Key (val Var)) (con @ (inc (cdr @))) (push Var (cons Key 1)) ) Key )(de **sum (L) (let S 1 (for I (cdr L) (inc 'S (** (car L) I)) ) S ) )(de factor-sum (N) (if (=1 N) 0 (let (R NIL D 2 L (1 2 2 . (4 2 4 2 4 6 2 6 .)) M (sqrt N) N1 N S 1 ) (while (>= M D) (if (=0 (% N1 D)) (setq M (sqrt (setq N1 (/ N1 (accud 'R D)))) ) (inc 'D (pop 'L)) ) ) (accud 'R N1) (for I R (setq S (* S (**sum I))) ) (- S N) ) ) )(bench (for I 20000 (let X (factor-sum I) (and (< I X) (= I (factor-sum X)) (println I X) ) ) ) )
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 184160.101 sec
*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;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
Rather than populating an array with the sum of the proper divisors and then searching, the approach here calculates them on the fly as needed, saving memory, and avoiding a noticeable lag on program startup on the ancient systems that hosted PL/I-80.
amicable: procedure options (main); %replace search_limit by 20000; dcl (a, b, found) fixed bin; put skip list ('Searching for amicable pairs up to '); put edit (search_limit) (f(5)); found = 0; do a = 2 to search_limit; b = sumf(a); if (b > a) then do; if (sumf(b) = a) then do; found = found + 1; put skip edit (a,b) (f(7)); end; end; end; put skip list (found, ' pairs were found'); stop;/* return sum of the proper divisors of n */sumf: procedure(n) returns (fixed bin); dcl (n, sum, f1, f2) fixed bin; sum = 1; /* 1 is a proper divisor of every number */ f1 = 2; do while ((f1 * f1) < n); if mod(n, f1) = 0 then do; sum = sum + f1; f2 = n / f1; /* don't double count identical co-factors! */ if f2 > f1 then sum = sum + f2; end; f1 = f1 + 1; end; return (sum);end sumf;end amicable;Searching for amicable pairs up to 20000 220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416 8 pairs were found
amicable: procedure options (main); %replace search_limit by 20000; dcl sumf( 1 : search_limit ) fixed bin; dcl (a, b, found) fixed bin; put skip list ('Searching for amicable pairs up to '); put edit (search_limit) (f(5)); do a = 1 to search_limit; sumf( a ) = 1; end; do a = 2 to search_limit; do b = a + a to search_limit by a; sumf( b ) = sumf( b ) + a; end; end; found = 0; do a = 2 to search_limit; b = sumf(a); if (b > a) then do; if (sumf(b) = a) then do; found = found + 1; put skip edit (a,b) (f(7)); end; end; end; put skip list (found, ' pairs were found'); stop;end amicable;Same as the other PLI-80 sample.
100H:/* CP/M CALLS */BDOS: PROCEDURE (FN, ARG); DECLARE FN BYTE, ARG ADDRESS; GO TO 5; END BDOS;EXIT: PROCEDURE; CALL BDOS(0,0); END EXIT;PRINT: PROCEDURE (S); DECLARE S ADDRESS; CALL BDOS(9,S); END PRINT;/* PRINT A NUMBER */PRINT$NUMBER: PROCEDURE (N); DECLARE S (6) BYTE INITIAL ('.....$'); DECLARE (N, P) ADDRESS, C BASED P BYTE; P = .S(5);DIGIT: P = P - 1; C = N MOD 10 + '0'; N = N / 10; IF N > 0 THEN GO TO DIGIT; CALL PRINT(P);END PRINT$NUMBER;/* CALCULATE SUMS OF PROPER DIVISORS */DECLARE DIV$SUM (20$001) ADDRESS;DECLARE (I, J) ADDRESS;DO I=2 TO 20$000; DIV$SUM(I) = 1; END;DO I=2 TO 10$000; DO J = I*2 TO 20$000 BY I; DIV$SUM(J) = DIV$SUM(J) + I; END;END;/* TEST EACH PAIR */DO I=2 TO 20$000; J = DIVSUM(I); IF J > I AND J <= 20$000 THEN DO; IF DIV$SUM(J) = I THEN DO; CALL PRINT$NUMBER(I); CALL PRINT(.', $'); CALL PRINT$NUMBER(J); CALL PRINT(.(13,10,'$')); END; END;END;CALL EXIT;EOF220, 2841184, 12102620, 29245020, 55646232, 636810744, 1085612285, 1459517296, 18416
functionGet-ProperDivisorSum([int]$N){$Sum=1If($N-gt3){$SqrtN=[math]::Sqrt($N)ForEach($Divisor1in2..$SqrtN){$Divisor2=$N/$Divisor1If($Divisor2-is[int]){$Sum+=$Divisor1+$Divisor2}}If($SqrtN-is[int]){$Sum-=$SqrtN}}return$Sum}functionGet-AmicablePairs($N=300){ForEach($Xin1..$N){$Sum=Get-ProperDivisorSum$XIf($Sum-gt$X-and$X-eq(Get-ProperDivisorSum$Sum)){"$X, $Sum"}}}Get-AmicablePairs20000
220, 2841184, 12102620, 29245020, 55646232, 636810744, 1085612285, 1459517296, 18416
With some guidance from other solutions here:
divisor(N,Divisor):-UpperBoundisround(sqrt(N)),between(1,UpperBound,D),0isNmodD,(Divisor=D;LargerDivisorisN/D,LargerDivisor=\=D,Divisor=LargerDivisor).proper_divisor(N,D):-divisor(N,D),D=\=N.assoc_num_divsSum_in_range(Low,High,Assoc):-findall(Num-DivSum,(between(Low,High,Num),aggregate_all(sum(D),proper_divisor(Num,D),DivSum)),Pairs),list_to_assoc(Pairs,Assoc).get_amicable_pair(Assoc,M-N):-gen_assoc(M,Assoc,N),M<N,get_assoc(N,Assoc,M).amicable_pairs_under_20000(Pairs):-assoc_num_divsSum_in_range(1,20000,Assoc),findall(P,get_amicable_pair(Assoc,P),Pairs).
Output:
?-amicable_pairs_under_20000(R).R=[220-284,1184-1210,2620-2924,5020-5564,6232-6368,10744-10856,12285-14595,17296-18416].
EnableExplicitProcedure.iSumProperDivisors(Number)IfNumber<2:ProcedureReturn0:EndIfProtectedi,sum=0Fori=1ToNumber/2IfNumber%i=0sum+iEndIfNextProcedureReturnsumEndProcedureDefinen,fDefineDimsum(19999)IfOpenConsole()Forn=1To19999sum(n)=SumProperDivisors(n)NextPrintN("The pairs of amicable numbers below 20,000 are : ")PrintN("")Forn=1To19998f=sum(n)Iff<=nOrf<1Orf>19999:Continue:EndIfIff=sum(n)Andn=sum(f)PrintN(RSet(Str(n),5)+" and "+RSet(Str(sum(n)),5))EndIfNextPrintN("")PrintN("Press any key to close the console")Repeat:Delay(10):UntilInkey()<>""CloseConsole()EndIf
The pairs of amicable numbers below 20,000 are : 220 and 284 1184 and 1210 2620 and 2924 5020 and 5564 6232 and 636810744 and 1085612285 and 1459517296 and 18416
ImportingProper divisors from prime factors:
fromproper_divisorsimportproper_divsdefamicable(rangemax=20000):n2divsum={n:sum(proper_divs(n))forninrange(1,rangemax+1)}fornum,divsuminn2divsum.items():ifnum<divsumanddivsum<=rangemaxandn2divsum[divsum]==num:yieldnum,divsumif__name__=='__main__':fornum,divsuminamicable():print('Amicable pair:%i and%i With proper divisors:\n%r\n%r'%(num,divsum,sorted(proper_divs(num)),sorted(proper_divs(divsum))))
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]
Or, supplying our ownproperDivisors function, and defining the harvest in terms of a genericconcatMap:
'''Amicable pairs'''fromitertoolsimportchainfrommathimportsqrt# amicablePairsUpTo :: Int -> [(Int, Int)]defamicablePairsUpTo(n):'''List of all amicable pairs of integers below n. '''sigma=compose(sum)(properDivisors)defamicable(x):y=sigma(x)return[(x,y)]if(x<yandx==sigma(y))else[]returnconcatMap(amicable)(enumFromTo(1)(n))# TEST ----------------------------------------------------# main :: IO ()defmain():'''Amicable pairs of integers up to 20000'''forxinamicablePairsUpTo(20000):print(x)# GENERIC -------------------------------------------------# compose (<<<) :: (b -> c) -> (a -> b) -> a -> cdefcompose(g):'''Right to left function composition.'''returnlambdaf:lambdax:g(f(x))# concatMap :: (a -> [b]) -> [a] -> [b]defconcatMap(f):'''A concatenated list or string over which a function f has been mapped. The list monad can be derived by using an (a -> [b]) function which wraps its output in a list (using an empty list to represent computational failure). '''returnlambdaxs:(''.joinifisinstance(xs,str)elselist)(chain.from_iterable(map(f,xs)))# enumFromTo :: Int -> Int -> [Int]defenumFromTo(m):'''Enumeration of integer values [m..n]'''defgo(n):returnlist(range(m,1+n))returnlambdan:go(n)# properDivisors :: Int -> [Int]defproperDivisors(n):'''Positive divisors of n, excluding n itself'''root_=sqrt(n)intRoot=int(root_)blnSqr=root_==intRootlows=[xforxinrange(1,1+intRoot)if0==n%x]returnlows+[n//xforxinreversed(lows[1:-1]ifblnSqrelselows[1:])]# MAIN ---if__name__=='__main__':main()
(220, 284)(1184, 1210)(2620, 2924)(5020, 5564)(6232, 6368)(10744, 10856)(12285, 14595)(17296, 18416)
properdivisors is defined atProper divisors#Quackery.
[ properdivisors dup size 0 = iff [ drop 0 ] done behead swap witheach + ] is spd ( n --> n ) [ dup dup spd dup spd rot = unrot > and ] is largeamicable ( n --> b ) [ [] swap times [ i^ largeamicable if [ i^ dup spd swap join nested join ] ] ] is amicables ( n --> [ ) 20000 amicables witheach [ witheach [ echo sp ] cr ]
220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416
divisors<-function(n){Filter(function(m)0==n%%m,1:(n/2))}table=sapply(1:19999,function(n)sum(divisors(n)))for(nin1:19999){m=table[n]if((m>n)&&(m<20000)&&(n==table[m]))cat(n," ",m,"\n")}
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
WithProper_divisors#Racket in place:
#langracket(require"proper-divisors.rkt")(defineSCOPE20000)(defineP(let((P-v(vector)))(λ(n)(set!P-v(fold-divisorsP-vn0+))(vector-refP-vn))));; returns #f if not an amicable number, amicable pairing otherwise(define(amicable?n)(definem(Pn))(definem-sod(Pm))(and(=m-sodn)(<mn); each pair exactly once, also eliminates perfect numbersm))(void(amicable?SCOPE)); prime the memoisation(for*((n(in-range1(add1SCOPE)))(m(in-value(amicable?n)))#:whenm)(printf#<<EOSamicable pair: ~a, ~a ~a: divisors: ~a ~a: divisors: ~aEOSnmn(proper-divisorsn)m(proper-divisorsm)))
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)
(formerly Perl 6)
subpropdivsum (\x) {my@l =1ifx >1; (2 ..x.sqrt.floor).map: -> \d {unlessx %d {@l.push:d;my \y =xdivd;@l.push:yify !=d } }sum@l}(1..20000).race.map: ->$i {my$j =propdivsum($i);say"$i $j"if$j >$iand$i ==propdivsum($j);}
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
;- based on Lua code ;-)sum-of-divisors:func[n/localsum][sum:1; using `to-integer` for compatibility with Rebol2ford2(to-integersquare-root n)1[if0=remaindernd[sum:n/d+sum+d]]sum]forn2200001[ifn<m:sum-of-divisorsn[ifn=sum-of-divisorsm[print[ntabm]]]]
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
let isqrt = (v) => { Belt.Float.toInt( sqrt(Belt.Int.toFloat(v)))}let sum_divs = (n) => { let sum = ref(1) for d in 2 to isqrt(n) { if mod(n, d) == 0 { sum.contents = sum.contents + (n / d + d) } } sum.contents}{ for n in 2 to 20000 { let m = sum_divs(n) if (m > n) { if sum_divs(m) == n { Printf.printf("%d %d\n", n, m) } } }}$ bsc ampairs.res > ampairs.bs.js$ node ampairs.bs.js220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
As given by the OEIS A259180
1 220 284 2 1184 1210 3 2620 2924 4 5020 5564 5 6232 6368 6 10744 10856 7 12285 14595 8 17296 18416 9 63020 7608410 66928 6699211 67095 7114512 69615 8763313 79750 8873014 100485 12415515 122265 13981516 122368 12315217 141664 15317618 142310 16873019 171856 17633620 176272 18084821 185368 20343222 196724 20244423 280540 365084
/*REXX*/Calltime'R'Dox=1To20000pd=proper_divisors(x)sumpd.x=sum(pd)EndSay'sum(pd) computed in'time('E')'seconds'Calltime'R'Dox=1To20000/* If x//1000=0 Then Say x time() */Doy=x+1To20000Ify=sumpd.x&,x=sumpd.yThenSayxy'found after'time('E')'seconds'EndEndSaytime('E')'seconds total search time'Exitproper_divisors:ProcedureParseArgnPd=''Ifn=1ThenReturn''Ifn//2=1Then/* odd number */delta=2Else/* even number */delta=1Dod=1Ton%2BydeltaIfn//d=0Thenpd=pddEndReturnspace(pd)sum:ProcedureParseArglistsum=0Doi=1Towords(list)sum=sum+word(list,i)EndReturnsum
sum(pd) computed in 22.449000 seconds220 284 found after 0.824000 seconds1184 1210 found after 4.447000 seconds2620 2924 found after 9.595000 seconds5020 5564 found after 17.367000 seconds6232 6368 found after 20.854000 seconds10744 10856 found after 31.320000 seconds12285 14595 found after 34.160000 seconds17296 18416 found after 39.598000 seconds40.354000 seconds total search time
Modules:How to use
Modules:Source code
Function Amicables(n) in module Sequences generates all amicable pairs up to n and returns the number of pairs.
--22Mar2025includeSettingssay'AMICABLE PAIRS'sayversionsayargnnumericdigits16ifn=''thenn=20000show=(n>0);n=Abs(n)a=Amicables(n)saytime('e')a'amicable pairs collected'ifshowthendodoi=1toasaytime('e')amic.1.i'and'amic.2.i'are an amicable pair'endendsaytime('e')'seconds'exitincludeSequencesincludeNumbersincludeFunctionsincludeAbend
AMICABLE PAIRSREXX-Regina_3.9.6(MT) 5.00 29 Apr 20240.003000 220 284 is an amicable pair0.023000 1184 1210 is an amicable pair0.072000 2620 2924 is an amicable pair0.161000 5020 5564 is an amicable pair0.194000 6232 6368 is an amicable pair0.406000 10744 10856 is an amicable pair0.625000 12285 14595 is an amicable pair0.883000 17296 18416 is an amicable pair1.000000 seconds
Up to 20,000Version 1 62.803sVersion 2 1.000sUp to 200,000Version 1 InfVersion 2 36.109s
size = 18500for n = 1 to size m = amicable(n) if m>n and amicable(m)=n see "" + n + " and " + m + nl oknextsee "OK" + nlfunc amicable nr sum = 1 for d = 2 to sqrt(nr) if nr % d = 0 sum = sum + d sum = sum + nr / d ok next return sum
≪ {} 2 ROTFOR jIF DUP j POS NOTTHEN@ avoiding duplicates j DIVIS ∑LIST j - DUPIF DUP 1 ≠ OVER j ≠ ANDTHENIF DUP DIVIS ∑LIST SWAP - j ==THEN + j +ELSE DROPENDELSE DROP2ENDENDNEXT {} 1 3 PICK SIZEFOR j@ formatting the list for output OVER j DUP 1 + SUB REVLIST 1 →LIST + 2STEP NIP≫ 'TASK' STO200000TASK1: {{220 284} {1184 1210} {2620 2924} {5020 5564} {6232 6368} {10744 10856} {12285 14595} {17296 18416}}Withproper_divisors#Ruby in place:
h={}(1..20_000).each{|n|h[n]=n.proper_divisors.sum}h.select{|k,v|h[v]==k&&k<v}.eachdo|key,val|# k<v filters out doubles and perfectsputs"#{key} and#{val}"end
220 and 2841184 and 12102620 and 29245020 and 55646232 and 636810744 and 1085612285 and 1459517296 and 18416
size = 18500for n = 1 to size m = amicable(n) if m > n and amicable(m) = n then print n ; " and " ; mnext function amicable(nr) amicable = 1 for d = 2 to sqr(nr) if nr mod d = 0 then amicable = amicable + d + nr / d next end function
220 and 2841184 and 12102620 and 29245020 and 55646232 and 636810744 and 1085612285 and 1459517296 and 18416
fnsum_of_divisors(val:u32)->u32{(1..val/2+1).filter(|n|val%n==0).fold(0,|sum,n|sum+n)}fnmain(){letiter=(1..20_000).map(|i|(i,sum_of_divisors(i))).filter(|&(i,div_sum)|i>div_sum);for(i,sum1)initer{ifsum_of_divisors(sum1)==i{println!("{} {}",i,sum1);}}}
284 2201210 11842924 26205564 50206368 623210856 1074414595 1228518416 17296
fnmain(){constRANGE_MAX:u32=20_000;letproper_divs=|n:u32|->Vec<u32>{(1..=(n+1)/2).filter(|&x|n%x==0).collect()};letn2d:Vec<u32>=(1..=RANGE_MAX).map(|n|proper_divs(n).iter().sum()).collect();for(n,&div_sum)inn2d.iter().enumerate(){letn=nasu32+1;ifn<div_sum&&div_sum<=RANGE_MAX&&n2d[div_sumasusize-1]==n{println!("Amicable pair: {} and {} with proper divisors:",n,div_sum);println!(" {:?}",proper_divs(n));println!(" {:?}",proper_divs(div_sum));}}}
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]
# Define the sum of proper divisors functiondefsum_of_proper_divisors(n):returnsum(divisors(n))-n# Iterate over the desired rangeforxinrange(1,20001):y=sum_of_proper_divisors(x)ify>x:ifx==sum_of_proper_divisors(y):print(f"{x}{y}")
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
$lines$constantsearch_limit=20000vara,b,count=integerdimintegersumf(search_limit)print"Searching up to";search_limit;" for amicable pairs:"rem-setupthetableofproperdivisorsumsfora=1tosearch_limitsumf(a)=1nextafora=2tosearch_limitb=a+awhile(b>0)and(b<=search_limit)dobeginsumf(b)=sumf(b)+ab=b+aendnextarem-searchforpairsusingthetablecount=0fora=2tosearch_limitb=sumf(a)if(b>a)and(b<search_limit)thenifa=sumf(b)thenbeginprintusing"##### #####";a;bcount=count+1endnextaprintcount;" pairs were found"end
Searching up to 20000 for amicable pairs: 220 284 1184 1210 2620 2924 5020 5564 6232 636810744 1085612285 1459517296 18416 8 pairs were found
defproperDivisors(n:Int)=(1ton/2).filter(i=>n%i==0)valdivisorsSum=(1to20000).map(i=>i->properDivisors(i).sum).toMapvalresult=divisorsSum.filter(v=>v._1<v._2&&divisorsSum.get(v._2)==Some(v._1))println(resultmkString", ")
5020 -> 5564, 220 -> 284, 6232 -> 6368, 17296 -> 18416, 2620 -> 2924, 10744 -> 10856, 12285 -> 14595, 1184 -> 1210
(import(schemebase)(schemeinexact)(schemewrite)(only(srfi1)fold));; return a list of the proper-divisors of n(define(proper-divisorsn)(let((root(sqrtn)))(letloop((divisors(list1))(i2))(if(>iroot)divisors(loop(if(zero?(moduloni))(if(=(squarei)n)(considivisors)(append(listi(quotientni))divisors))divisors)(+1i))))))(define(sum-proper-divisorsn)(if(<n2)0(fold+0(proper-divisorsn))))(define*max-n*20000);; hold sums of proper divisors in a cache, to avoid recalculating(define*cache*(make-vector(+1*max-n*)))(for-each(lambda(i)(vector-set!*cache*i(sum-proper-divisorsi)))(iota*max-n*1))(define(amicable-pair?ij)(and(not(=ij))(=i(vector-ref*cache*j))(=j(vector-ref*cache*i))));; double loop to *max-n*, displaying all amicable pairs(letloop-i((i1))(when(<=i*max-n*)(letloop-j((ji))(when(<=j*max-n*)(when(amicable-pair?ij)(display(string-append"Amicable pair: "(number->stringi)" "(number->stringj)))(newline))(loop-j(+1j))))(loop-i(+1i))))
Amicable pair: 220 284Amicable pair: 1184 1210Amicable pair: 2620 2924Amicable pair: 5020 5564Amicable pair: 6232 6368Amicable pair: 10744 10856Amicable pair: 12285 14595Amicable pair: 17296 18416
program amicable_pairs; p := propDivSums(20000); loop for [n,m] in p | n = p(p(n)) and n<m do print([n,m]); end loop; proc propDivSums(n); divs := {}; loop for i in [1..n] do loop for j in [i*2, i*3..n] do divs(j) +:= i; end loop; end loop; return divs; end proc;end program;[220 284][1184 1210][2620 2924][5020 5564][6232 6368][10744 10856][12285 14595][17296 18416]
funcpropdivsum(n){n.sigma-n}foriin(1..20000){varj=propdivsum(i)say"#{i}#{j}"if(j>i&&i==propdivsum(j))}
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
importfuncDarwin.sqrtfuncsqrt(x:Int)->Int{returnInt(sqrt(Double(x)))}funcproperDivs(n:Int)->[Int]{ifn==1{return[]}varresult=[Int]()fordivinfilter(1...sqrt(n),{n%$0==0}){result.append(div)ifn/div!=div&&n/div!=n{result.append(n/div)}}returnsorted(result)}funcsumDivs(n:Int)->Int{structCache{staticvarsum=[Int:Int]()}ifletsum=Cache.sum[n]{returnsum}letsum=properDivs(n).reduce(0){$0+$1}Cache.sum[n]=sumreturnsum}funcamicable(n:Int,m:Int)->Bool{ifn==m{returnfalse}ifsumDivs(n)!=m||sumDivs(m)!=n{returnfalse}returntrue}varpairs=[(Int,Int)]()fornin1..<20_000{forminn+1...20_000{ifamicable(n,m){pairs.append(n,m)println("\(n,m)")}}}
about 800 times faster.
importfuncDarwin.sqrtfuncsqrt(x:Int)->Int{returnInt(sqrt(Double(x)))}funcsigma(n:Int)->Int{ifn==1{return0}// definition of aliquot sumvarresult=1letroot=sqrt(n)forvardiv=2;div<=root;++div{ifn%div==0{result+=div+n/div}}ifroot*root==n{result-=root}return(result)}funcamicables(upTo:Int)->(){varaliquot=Array(count:upTo+1,repeatedValue:0)foriin1...upTo{// fill lookup arrayaliquot[i]=sigma(i)}foriin1...upTo{leta=aliquot[i]ifa>upTo{continue}//second part of pair out-of-boundsifa==i{continue}//skip perfect numbersifi==aliquot[a]{print("\(i,a)")aliquot[a]=upTo+1//prevent second display of pair}}}amicables(20_000)
(220, 284)(1184, 1210)(2620, 2924)(5020, 5564)(6232, 6368)(10744, 10856)(12285, 14595)(17296, 18416)
dimsums(20000)subsum_proper_divisors(n)dimsum=0dimiifn>1thenfori=1to(n\2)ifn%%i=0thensum=sum+iendifnextendifreturnsumendsubdimi,jfori=1to20000sums(i)=sum_proper_divisors(i)forj=i-1to2by-1ifsums(i)=jandsums(j)=ithenprint"Amicable pair:";sums(i);"-";sums(j)exitforendifnextnext
>tbas amicable_pairs.basAmicable pair: 220 - 284Amicable pair: 1184 - 1210Amicable pair: 2620 - 2924Amicable pair: 5020 - 5564Amicable pair: 6232 - 6368Amicable pair: 10744 - 10856Amicable pair: 12285 - 14595Amicable pair: 17296 - 18416
procproperDivisors{n}{if{$n==1}returnsetdivs1setsum1for{seti2}{$i*$i<=$n}{incri}{if{!($n%$i)}{lappenddivs$iincrsum$iif{$i*$i<$n}{lappenddivs[setd[expr{$n/$i}]]incrsum$d}}}return[list$sum$divs]}procamicablePairs{limit}{setresult{}setsums[setdivs{{}}]for{setn1}{$n<$limit}{incrn}{lassign[properDivisors$n]sumdlappendsums$sumlappenddivs[lsort-integer$d]}for{setn1}{$n<$limit}{incrn}{setnsum[lindex$sums$n]for{setm1}{$m<$n}{incrm}{if{$n==[lindex$sums$m]&&$m==$nsum}{lappendresult$m$n[lindex$divs$m][lindex$divs$n]}}}return$result}foreach{mnmdnd}[amicablePairs20000]{puts"$m and $n are an amicable pair with these proper divisors"puts"\t$m : $md"puts"\t$n : $nd"}
220 and 284 are an amicable pair with these proper divisors220 : 1 2 4 5 10 11 20 22 44 55 110284 : 1 2 4 71 1421184 and 1210 are an amicable pair with these proper divisors1184 : 1 2 4 8 16 32 37 74 148 296 5921210 : 1 2 5 10 11 22 55 110 121 242 6052620 and 2924 are an amicable pair with these proper divisors2620 : 1 2 4 5 10 20 131 262 524 655 13102924 : 1 2 4 17 34 43 68 86 172 731 14625020 and 5564 are an amicable pair with these proper divisors5020 : 1 2 4 5 10 20 251 502 1004 1255 25105564 : 1 2 4 13 26 52 107 214 428 1391 27826232 and 6368 are an amicable pair with these proper divisors6232 : 1 2 4 8 19 38 41 76 82 152 164 328 779 1558 31166368 : 1 2 4 8 16 32 199 398 796 1592 318410744 and 10856 are an amicable pair with these proper divisors10744 : 1 2 4 8 17 34 68 79 136 158 316 632 1343 2686 537210856 : 1 2 4 8 23 46 59 92 118 184 236 472 1357 2714 542812285 and 14595 are an amicable pair with these proper divisors12285 : 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 409514595 : 1 3 5 7 15 21 35 105 139 417 695 973 2085 2919 486517296 and 18416 are an amicable pair with these proper divisors17296 : 1 2 4 8 16 23 46 47 92 94 184 188 368 376 752 1081 2162 4324 864818416 : 1 2 4 8 16 1151 2302 4604 9208
#langtransdMainModule:{_start:(lambda(withsum0d0fFilter(from:1to:20000apply:(lambda(=sum1)(foriinRange(2(to-Int(sqrt@it)))do(if(not(mod@iti))(=d(/@iti))(+=sumi)(if(neqdi)(+=sumd))))(retsum)))(withv(to-vectorf)(foriinvdo(if(and(<i(sizev))(eq(+@idx1)(getv(-i1)))(<i(getv(-i1))))(textout(+@idx1)", "i"\n"))))))}
284, 2201210, 11842924, 26205564, 50206368, 623210856, 1074414595, 1228518416, 17296
Input "Limit: ";lPrint "Amicable pairs < ";lFor n = 1 To l m = FUNC(_SumDivisors (n))-n If m = 0 Then Continue ' No division by zero, please p = FUNC(_SumDivisors (m))-m If (n=p) * (n<m) Then Print n;" and ";mNextEnd_LeastPower Param(2) Local(1) c@ = a@ Do While (b@ % c@) = 0 c@ = c@ * a@ LoopReturn (c@)' Return the sum of the proper divisors of a@_SumDivisors Param(1) Local(4) b@ = a@ c@ = 1 ' Handle two specially d@ = FUNC(_LeastPower (2,b@)) c@ = c@ * (d@ - 1) b@ = b@ / (d@ / 2) ' Handle odd factors For e@ = 3 Step 2 While (e@*e@) < (b@+1) d@ = FUNC(_LeastPower (e@,b@)) c@ = c@ * ((d@ - 1) / (e@ - 1)) b@ = b@ / (d@ / e@) Loop ' At this point, t must be one or prime If (b@ > 1) c@ = c@ * (b@+1)Return (c@)
Limit: 20000Amicable pairs < 20000220 and 2841184 and 12102620 and 29245020 and 55646232 and 636810744 and 1085612285 and 1459517296 and 184160 OK, 0:238
···http://rosettacode.org/wiki/Amicable_pairs···■ AmicablePairs § static ▶ main • args⦂ String[] ∀ n ∈ 1…20000 m⦂ int: sumPropDivs n if m < n = sumPropDivs m System.out.println "⸨m⸩ ; ⸨n⸩" ▶ sumPropDivs⦂ int • n⦂ int m⦂ int: 1 ∀ i ∈ √n ⋯> 1 m +: n \ i = 0 ? i + (i = n / i ? 0 ! n / i) ! 0 ⏎ m
OptionExplicitPublicSubAmicablePairs()Dima(2To20000)AsLong,cAsNewCollection,iAsLong,jAsLong,t#t=TimerFori=LBound(a)ToUBound(a)'collect the sum of the proper divisors'of each numbers between 2 and 20000a(i)=S(i)Next'Double Loops to test the amicableFori=LBound(a)ToUBound(a)Forj=i+1ToUBound(a)Ifi=a(j)ThenIfa(i)=jThenOnErrorResumeNextc.Addi&" : "&j,CStr(i*j)OnErrorGoTo0ExitForEndIfEndIfNextNext'End. Return :Debug.Print"Execution Time : "&Timer-t&" seconds."Debug.Print"Amicable pairs below 20 000 are : "Fori=1Toc.CountDebug.Printc.Item(i)NextiEndSubPrivateFunctionS(nAsLong)AsLong'returns the sum of the proper divisors of nDimjAsLongForj=1Ton\2IfnModj=0ThenS=j+SNextEndFunction
Execution Time : 7,95703125 seconds.Amicable pairs below 20 000 are : 220 : 2841184 : 12102620 : 29245020 : 55646232 : 636810744 : 1085612285 : 1459517296 : 18416
Not at all optimal. :-(
start=NowSetnlookup=CreateObject("Scripting.Dictionary")Setuniquepair=CreateObject("Scripting.Dictionary")Fori=1To20000sum=0Forn=1To20000Ifn<iThenIfiModn=0Thensum=sum+nEndIfEndIfNextnlookup.Addi,sumNextForj=1To20000sum=0Form=1To20000Ifm<jThenIfjModm=0Thensum=sum+mEndIfEndIfNextIfnlookup.Exists(sum)Andnlookup.Item(sum)=jAndj<>sum_Anduniquepair.Exists(sum)=FalseThenuniquepair.Addj,sumEndIfNextForEachkeyInuniquepair.KeysWScript.Echokey&":"&uniquepair.Item(key)NextWScript.Echo"Execution Time: "&DateDiff("s",Start,Now)&" seconds"
220:2841184:12102620:29245020:55646232:636810744:1085612285:1459517296:18416Execution Time: 162 seconds
fnpfac_sum(iint)int{mutsum:=0forp:=1;p<=i/2;p++{ifi%p==0{sum+=p}}returnsum}fnmain(){a:=[]int{len:20000,init:pfac_sum(it)}println('Theamicablepairsbelow20,000are:')fornin2..a.len{m:=a[n]ifm>n&&m<20000&&n==a[m]{println('${n:5}and${m:5}')}}}
The amicable pairs below 20,000 are: 220 and 284 1184 and 1210 2620 and 2924 5020 and 5564 6232 and 636810744 and 1085612285 and 1459517296 and 18416
10 M=2000020 I=130 :I)=140 I=I+150 #=M>I*3060 I=270 J=I+I80 :J)=:J)+I90 J=J+I100 #=M>J*80110 I=I+1120 #=(M/2)>I*70130 I=1140 J=:I)150 #=(I<J)*I=:J)*190160 I=I+1170 #=M>I*140180 #=999190 ?=I200 $=9210 ?=J220 ?=""230 #=!
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
import"./fmt"forFmtimport"./math"forInt,Numsvara=List.filled(20000,0)for(iin1...20000)a[i]=Nums.sum(Int.properDivisors(i))System.print("The amicable pairs below 20,000 are:")for(nin2...19999){varm=a[n]if(m>n&&m<20000&&n==a[m]){Fmt.print(" $5d and $5d",n,m)}}
The amicable pairs below 20,000 are: 220 and 284 1184 and 1210 2620 and 2924 5020 and 5564 6232 and 6368 10744 and 10856 12285 and 14595 17296 and 18416
func SumDiv(Num); \Return sum of proper divisors of Numint Num, Div, Sum, Quot;[Div:= 2;Sum:= 0;loop [Quot:= Num/Div; if Div > Quot then quit; if rem(0) = 0 then [Sum:= Sum + Div; if Div # Quot then Sum:= Sum + Quot; ]; Div:= Div+1; ];return Sum+1;];def Limit = 20000;int Tbl(Limit), N, M;[for N:= 0 to Limit-1 do Tbl(N):= SumDiv(N);for N:= 1 to Limit-1 do [M:= Tbl(N); if M<Limit & N=Tbl(M) & M>N then [IntOut(0, N); ChOut(0, 9\tab\); IntOut(0, M); CrLf(0); ]; ];]
220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416
sub sumDivs(n) local sum, d sum = 1 for d = 2 to sqrt(n) if not mod(n, d) then sum = sum + d sum = sum + n / d end if next return sumend subfor n = 2 to 20000 m = sumDivs(n) if m > n then if sumDivs(m) = n print n, "\t", m end ifnextprint : print peek("millisrunning"), " ms"constMAXIMUM:u32=20_000;// Fill up a given array with arr[n] = sum(propDivs(n))pubfncalcPropDivs(divs:[]u32)void{for(divs)|*d|d.*=1;vari:u32=2;while(i<=divs.len/2):(i+=1){varj=i*2;while(j<divs.len):(j+=i)divs[j]+=i;}}// Are (A, B) an amicable pair?pubfnamicable(divs:[]constu32,a:u32,b:u32)bool{returndivs[a]==banda==divs[b];}pubfnmain()!void{conststdout=@import("std").io.getStdOut().writer();vardivs:[MAXIMUM+1]u32=undefined;calcPropDivs(divs[0..]);vara:u32=1;while(a<divs.len):(a+=1){varb=a+1;while(b<divs.len):(b+=1){if(amicable(divs[0..],a,b))trystdout.print("{d}, {d}\n",.{a,b});}}}
220, 2841184, 12102620, 29245020, 55646232, 636810744, 1085612285, 1459517296, 18416
Slooooow
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();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))
10LETlimit=2000020PRINT"Amicable pairs < ";limit30FORn=1TOlimit40LETnum=n:GOSUB100050LETm=num60GOSUB100070IFn=numANDn<mTHENPRINTn;" ";m80NEXTn90STOP1000REM sumprop1010IFnum<2THENLETnum=0:RETURN1020LETsum=11030LETroot=SQRnum1040FORi=2TOroot-.011050IFnum/i=INT(num/i)THENLETsum=sum+i+num/i1060NEXTi1070IFnum/root=INT(num/root)THENLETsum=sum+root1080LETnum=sum1090RETURN
Amicable pairs < 20000220 2841184 12102620 29245020 55646232 636810744 1085612285 1459517296 18416