3
$\begingroup$

The following text is from Problem 21 - Project Euler :

Let $d(n)$ be defined as the sum of proper divisors of $n$ (numbers less than n which divide evenly into $n$). If $d(a) = b$ and $d(b) = a$, where $a ≠ b$, then $a$ and $b$ are an amicable pair and each of $a$ and $b$ are called amicable numbers.

For example, the proper divisors of $220$ are $1$, $2$, $4$, $5$, $10$, $11$, $20$, $22$, $44$, $55$ and $110$; therefore $d(220) = 284$. The proper divisors of $284$ are $1$, $2$, $4$, $71$ and $142$; so $d(284) = 220$.

Evaluate the sum of all the amicable numbers under $10000$.

The fastest Mathematica snippet I came up with uses the characterization of amicable numbers as proper involutive points for the divisor sum function:

DSum[n_]:=DivisorSigma[1,n]-n
AmicableCandidates=Select[Range[10000], DSum[DSum[#]]==# &];
AmicableNumbers=Select[AmicableCandidates, DSum[#]!=# &];
Sum[AmicableNumbers[[i]], {i,1,Length[AmicableNumbers]}]

Timing this yields 0.09375s on my machine.

  1. Can this code be improved?
  2. What are some other efficient ways to compute the sum of all amicable numbers under $10000$?
$\endgroup$
2
$\begingroup$

Since DivisorSigma is Listable, one can use a vectorized approach:

amicable[list_] := Module[{f1, f2},
    f1 = DivisorSigma[1, list] - list;
    f2 = DivisorSigma[1, f1] - f1;
    Pick[list, Unitize[f1-list]BitXor[Unitize[f2-list], 1], 1]
]
amicable[Range[2,10^4]] //Total //RepeatedTiming

{0.045, 31626}

I don't include 1 because DivisorSigma[1, 1] - 1 is 0, and DivisorSigma[1, 0] is not defined. Hence the output of f2 becomes unpacked when 1 is included in the list.

$\endgroup$
2
$\begingroup$

Instant speedup: Change the last line to:

Total[AmicableNumbers]

Also try this:

DSum[n_]:=DivisorSigma[1,n]-n;
Total@Select[Range[10000], (DSum[DSum[#]]==# && DSum[#]!=#) &]
$\endgroup$
1
$\begingroup$

Just for another way:

f[n_] := DivisorSigma[1, n] - n;
g[n_] := NestList[f, n, 2]
Total[Cases[g /@ Range[2, 10000], {i_, j_, i_} /; i!= j]][[1]] // 
RepeatedTiming

yields:

{0.11, 31626}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.