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.
<!--
{-# 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.
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
[1],
nor do I check (b > a)
[2]. 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.
[1]: I did not cache b
because I wasn’t that familiar with Haskell’s syntax
for functions with multiple expressions.
[2]: 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} $$
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.
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.