# 021. Amicable Numbers

## Problem Statement

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 \neq 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, 110$$; therefore $$d(220) = 284$$. The proper divisors of $$284$$ are $$1, 2, 4, 71, 142$$; so $$d(284) = 220$$.

Evaluate the sum of all the amicable numbers under $$10{,}000$$.

## My Solution

I don’t think this problem can be efficiently solved manually. Too much grunt work. Given a number’s prime factorization, provides a formula for summing up the divisors.

For this problem, I’ll try out Haskell. lists several resources and recommends UPenn’s for getting started. is by the same author, and seems more applicable and less academic-focused. is a Haskell API search engine that helps one avoid re-implementing functionality from the Haskell library. There are VS Code extensions that use to offer contextual suggestions; they should help me hone my Haskell.

provides lecture slides as .html and .lhs files. Investigating .lhs files introduced me to literate programming, whose main idea is to regard a program as a communication to human beings rather than as a set of instructions to a computer. Nifty! There is even multi-mode support in Emacs, which switches between haskell-mode and latex-mode depending on where the cursor is.

Pays homage to “nothing new under the sun”. I thought Jupyter notebooks were quite fancy, but Literate Haskell already existed. Granted, Jupyter has bells and whistles for interactive programming.

<!--
{-# OPTIONS_GHC -Wall #-}
-->

\begin{code}

-- Haskell functions cannot take multiple arguments, hence this syntax.
isDivisor :: Integer -> Integer -> Bool
isDivisor n f = n mod f == 0

sumOfDivisors :: Integer -> Integer
sumOfDivisors n = sum $filter (\f -> isDivisor n f) [1..(n-1)] isAmicableNumber :: Integer -> Bool isAmicableNumber n = ((sumOfDivisors n) /= n) && (sumOfDivisors (sumOfDivisors n) == n) sumOfAmicableNumbers :: Integer -> Integer sumOfAmicableNumbers n = sum(filter isAmicableNumber [1..n]) -- Every Haskell program needs to have one function called main, which is its -- entry point. -- The$ operator binds to the right, and therefore
-- (sumOfAmicableNumbers 10000) is executed first.
main = do
putStrLn "The sum of amicable numbers under 10,000 is:"
print $sumOfAmicableNumbers 9999 -- 31626 \end{code}  The 45s running time is pretty bad considering that a similar Python script takes 6s. def sum_of_divisors(n): s = 0 for i in range(1, n): if n % i == 0: s += i return s def sum_of_divisors_slow(n): def is_factor(f): return n % f == 0 return sum(filter(is_factor, range(1, n)))  I’m surprised that using sum_of_divisors_slow takes about twice the time as using sum_of_divisors. I expected sum_of_divisors to be faster (given how for-loops go brrr), but not by much. ## Learning From Others' Solutions ### Minimizing the Number of Calls to the SumOfProperDivisors Function A construction like: for a in range(1, 10000): # 1, 2, 3, ..., 9999 for b in range(a + 1, 10000): # a + 1, a + 2, ..., 9999 if sum_of_proper_divisors(a) == b and sum_of_proper_divisors(b) == a: # ...  … calls sum_of_proper_divisors at most $$(9{,}999 \times 9{,}998)$$ times. Computing $$(9{,}998 \times 9{,}997) / 2$$: The inner for-loop executes $$9{,}998 + 9{,}997 + … + 1$$ times. The closed form of $$1 + 2 + 3 + … + n$$ is $$n \cdot (n + 1) / 2$$. However, if we calculate b = sum_of_proper_divisors(a), then sum_of_proper_divisors(b) should be equal to a, and therefore we can have: s = 0 for a in range(1, 10000): # 1, 2, 3, ..., 9999 b = sum_of_proper_divisors(a) # Given a, there's at most one possible b if b > a: # Avoid double counting (a, b) and (b, a) if sum_of_proper_divisors(b) == a: s += a + b  … which calls sum_of_proper_divisors at most $$(9{,}999 \times 2)$$ times. My Haskell solution calls sum_of_proper_divisors at most $$(9{,}999 \times 3)$$ times. $$\times 3$$ because I neither cache b , nor do I check (b > a) . Still, it did not occur to me that looping through possible pairs $$(a, b)$$ was also an available, albeit inefficient, option. Applying these two optimizations halves the running time from 45s to 23s. : I did not cache b because I wasn’t that familiar with Haskell’s syntax for functions with multiple expressions. : I did not check (b > a) because it wasn’t on my radar. My approach involved deciding if each number in $$[1, 2, …, 9{,}999]$$ was amicable or not. For amicable pair $$(a, b)$$, I included $$a$$ in the sum, and disregarded $$b$$ until I got to the pair $$(b, a)$$. ### Computing SumOfProperDivisors Efficiently To find the proper divisors of $$n$$, we need not check from $$1$$ to $$n$$. At a first approximation, we only need to check till $$n/2$$. However, we can go even further and only check till $$\sqrt{n}$$, and use the fact that if $$i : i \le \lfloor \sqrt{n} \rfloor$$ is a divisor, then $$n/i$$ is also a proper divisor of $$n$$. Furthermore, odd numbers cannot have even numbers as divisors. {-# OPTIONS_GHC -Wall #-} import GHC.Float (sqrtDouble) sumOfProperDivisorPair :: Integer -> Integer -> Integer sumOfProperDivisorPair n f | n == f = 0 | f == 1 = 1 | n mod f == 0 = f + (n div f) | otherwise = 0 sumOfProperDivisors :: Integer -> Integer sumOfProperDivisors 1 = 0 sumOfProperDivisors n = s where flooredSqrt = floor$ sqrtDouble $fromIntegral n upperLimit = if flooredSqrt * flooredSqrt == n then flooredSqrt + 1 else flooredSqrt stepSize = if odd n then 2 else 1 s = sum$ map
(sumOfProperDivisorPair n)
[1, (1 + stepSize) .. (upperLimit - 1)]

contributionToSum :: Integer -> Integer
contributionToSum a = res where
b = sumOfProperDivisors a
res =
if (b > a) && (sumOfProperDivisors b == a)
then a + b
else 0

sumOfAmicableNumbers :: Integer -> Integer
sumOfAmicableNumbers n = sum(map contributionToSum [1..n])

main :: IO()
main = do
print \$ sumOfAmicableNumbers 9999 -- 31626


Applying this optimization slashes the execution time from 23s to 1s.

Suppose that $$n$$ is a positive integer whose factorization into prime factors is $$\prod_{i=1}^k p_{i}^{m_i}$$. Notice that any divisor $$d$$ must be a product of some number of each of the prime factors, i.e. $$\prod_{i=1}^k p_{i}^{\mu_i}$$, where $$0 \le \mu_i \le m_i$$. Then the sum over all divisors becomes the sum over all possible choices for the $$\mu_{i}$$’s

$$\sum_{ d \mid n} d = \sum_{0 \le \mu_i \le m_i} \prod_{i=1}^k p_{i}^{\mu_i} = \sum_{\mu_1 = 0}^{m_1} \sum_{\mu_2=0}^{m_2} … \sum_{\mu_k=0}^{m_k} \prod_{i=1}^k p_{i}^{\mu_i} = \prod_{i=1}^k \left( \sum_{\mu_i = 0}^{m_i} p_{i}^{\mu_i} \right) = \prod_{i=1}^k \frac{p_{i}^{m_i + 1} - 1}{p_i - 1}$$

uses several mathematical maneuvers: expressing sum as a multiple sum, factoring a sum of products into a product of sums, and recognizing that each sum is a geometric series.

For instance, the factorization of $$220$$ is $$2^2 \cdot 5^1 \cdot 11^1$$, and therefore the sum of its proper divisors is:

$$\left( \frac{2^3 - 1}{1} \cdot \frac{5^2 - 1}{4} \cdot \frac{11^2 - 1}{10} \right) - 220 = 284$$

#!/usr/bin/env python

def sum_of_proper_divisors(n):
if n == 0: return 0
orig_n = n
s, factor = 1, 2
while factor * factor <= n and n > 1: # We only need to check up to sqrt(n)
multiplicity = 0
while n % factor == 0:
n /= factor
multiplicity += 1

if multiplicity > 0:
s *= (factor ** (multiplicity + 1) - 1) / (factor - 1)

factor = 3 if factor == 2 else factor + 2

# Account for any remaining prime factor greater than sqrt(n). There will be
# at most one such factor. (n^2 - 1) / (n-1) simplifies to (n + 1).
if n > 1: s *= (n + 1)

return s - orig_n

def is_amicable_number(n):
s = sum_of_proper_divisors(n)
if s == n: return False
return sum_of_proper_divisors(s) == n

def sum_of_amicable_numbers(n):
return sum(filter(is_amicable_number, range(1, n)))

if __name__ == "__main__":
print(sum_of_amicable_numbers(10000))


Solve PE #021 using the sum of divisors formula in Haskell. The Python implementation runs in 0.1s (compared to the 6s Python equivalent of this naïve approach ).

### Trivia

Amicable numbers have history to them. For example, Genesis 32:14 has Jacob giving Esau 220 sheep when he was afraid that Esau was going to kill him.

1. Starting with Haskell - School of Haskell | School of Haskell. Brent Yorgey. www.schoolofhaskell.com . Accessed Feb 8, 2022.
2. The Prime Glossary: amicable numbers. Chris K. Caldwell. primes.utm.edu . Accessed Feb 9, 2022.
3. Sum of n, n², or n³ | Brilliant Math & Science Wiki. brilliant.org . Accessed Feb 12, 2022.
4. ndmitchell/hlint: Haskell source code suggestions. github.com . Accessed Feb 12, 2022.