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 \ne 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\).
If \(n = p_1^{a_1} p_2^{a_2} \cdots p_m^{a_m}\) for distinct primes \(p_1,p_2,\ldots,p_m\), then the sum of all the divisors of \(n\) is
\[ \begin{aligned} \sigma(n) &= (1 + p_1 + p_1^2 + \cdots + p_1^{a_1}) \times \cdots \times (1 + p_m + p_m^2 + \cdots + p_m^{a_m})\\ &= \frac{p_1^{a_1+1} - 1}{p_1 - 1} \times \cdots \times \frac{p_m^{a_m+1} - 1}{p_m - 1} \end{aligned} \]and \(d(n) = \sigma(n) - n\).
My code calculates \(\sigma(n)\) for each \(n\) as follows. First, it sets \(\sigma[n] = 1\) for every \(n\). Then for every prime \(p\), it multiplies \(\sigma[kp]\) by \(\frac{p^2 - 1}{p-1}\) for all \(k\), then multiplies \(\sigma[kp^2]\) by \(\frac{p^3 - 1}{p^2 - 1}\) for all \(k\), then multiplies \(\sigma[kp^3]\) by \(\frac{p^4 - 1}{p^3 - 1}\) for all \(k\), until each multiple of \(p\) has been multiplied by the correct factor. After this has been repeated for all primes, we will have \(\sigma[n] = \sigma(n)\) for all \(n\).
From there, it's just a matter of setting \(d(n) = \sigma(n) - n\), iterating \(a\) from \(1\) up to \(10000\), and checking whether \(d(a) \neq a\) and \(d(d(a)) = a\).
from sympy import primerange
def p21(cap = 10000):
sigma = [1 for i in range(cap)]
for p in primerange(cap):
px = p
while px < cap:
for m in range(px, cap, px):
sigma[m] = sigma[m] * (px * p - 1) // (px - 1)
px *= p
res = 0
for a in range(2, cap):
b = sigma[a] - a
if b != a and b < cap and sigma[b] - b == a:
res += a
return res