Amicable Numbers¶
The sum of the proper divisors of a number is called the aliquot sum.
SageMath provides the sigma function, which can compute the sum of the divisors of $n$. The aliquot sum is then just sigma(n) - n.
def aliquot_sum(n): return sigma(n) - n
limit = 10000
If a number is equal to its own aliquot sum, it's called a perfect number, and we exclude those numbers from our total.
amicables = set()
for a in range(2, limit):
if a in amicables:
continue
b = aliquot_sum(a)
if a == b:
continue
if a == aliquot_sum(b):
amicables.update({a, b})
sum(amicables)
31626
Funny enough, there's only five pairs of amicable numbers below 10,000. If you looked up amicable numbers, you may have stumbled on all the numbers you need to answer the problem!
amicables
{220, 284, 1184, 1210, 2620, 2924, 5020, 5564, 6232, 6368}
Sum of divisors¶
Of course, you could implement your own divisor sum function. In problem 12, we implemented a divisor counting function, which is related.
One important property of the divisor sum function $\sigma(n)$ is that it is multiplicative. This means that if $a$ and $b$ are coprime, $\sigma(ab) = \sigma(a)\sigma(b)$. It follows that if $n = 2^{r_1}3^{r_2}5^{r_3}7^{r_4}\cdots$, then $$\sigma(n) = \sigma(2^{r_1})\sigma(3^{r_2})\sigma(5^{r_3})\sigma(7^{r_4})\cdots$$
Furthermore, for a prime $p$, $\sigma(p^k) = 1 + p + p^2 + \cdots + p^k$. This expression is actually a partial sum of a geometric series, which has a closed formula: $$1 + p + p^2 + \cdots + p^k = \frac{p^{k+1}-1}{p-1}$$
Putting it all together, we can compute $\sigma(n)$ as $$\sigma(n) = \frac{2^{r_1+1}-1}{2-1} \cdot \frac{3^{r_2+1}-1}{3-1} \cdot \frac{5^{r_3+1}-1}{5-1} \cdot \frac{7^{r_4+1}-1}{7-1} \cdot \cdots$$
Therefore, if you have the number's factorization (see problem 3), you can use it to compute the sum of its divisors.
Avoiding factoring¶
Since we need all the sums of divisors up to 10000, instead of factoring each number individually, we can compute each value of $\sigma$ in a loop, once again taking advantage of $\sigma$ being multiplicative. (In some ways, this method is similar to the sieve of Eratosthenes - see problem 10.)
def update_multiples(dsum, p, limit):
q = p
while True:
# sigma(a*b) = sigma(a) * sigma(b) if gcd(a, b) = 1
for k in range(2 * q, limit, q):
if k % (p*q) != 0:
dsum[k] *= dsum[q]
if p * q >= limit:
break
# sigma(p^k) = p^k + sigma(p^(k-1))
dsum[p*q] = p * q + dsum[q]
q *= p
def sum_of_divisors_range(limit):
dsum = [1 for n in range(0, limit)]
dsum[0] = 0
for n in range(0, limit):
if n == 0 or n == 1 or dsum[n] != 1:
# n is 0, 1, or composite
yield dsum[n]
continue
# n is prime
dsum[n] = n + 1
update_multiples(dsum, n, limit)
yield dsum[n]
With this method, we can redefine our aliquot sum function:
divisor_sums = list(sum_of_divisors_range(limit))
def aliquot_sum(n): return divisor_sums[n] - n
Then we compute the amicable numbers the same as before.
amicables = set()
for a in range(2, limit):
if a in amicables:
continue
b = aliquot_sum(a)
if a == b:
continue
# if b is greater than limit, it will cause an IndexError
if b < limit and a == aliquot_sum(b):
amicables.update({a, b})
sum(amicables)
31626
Relevant sequences¶
Copyright (C) 2025 filifa¶
This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International license and the BSD Zero Clause license.