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.

<!--
{-# 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)))

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.

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

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))

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.

References

1. #21 Amicable numbers - Project Euler. projecteuler.net . Accessed Feb 6, 2022.
2. Documentation. www.haskell.org . Accessed Feb 6, 2022.
3. Literate programming - HaskellWiki. wiki.haskell.org . Accessed Feb 6, 2022.
4. Lecture notes and assignments. Brent Yorgey. www.seas.upenn.edu . Accessed Feb 6, 2022.
5. Starting with Haskell - School of Haskell | School of Haskell. Brent Yorgey. www.schoolofhaskell.com . Accessed Feb 8, 2022.
6. Thread & Overview for 21 - Project Euler. projecteuler.net . projecteuler.net . Accessed Feb 9, 2022.
7. The Prime Glossary: amicable numbers. Chris K. Caldwell. primes.utm.edu . Accessed Feb 9, 2022.
8. formula for sum of divisors. planetmath.org . Accessed Feb 12, 2022.
9. Sum of n, n², or n³ | Brilliant Math & Science Wiki. brilliant.org . Accessed Feb 12, 2022.
10. Hoogle. hoogle.haskell.org . Accessed Feb 12, 2022.
11. ndmitchell/hlint: Haskell source code suggestions. github.com . Accessed Feb 12, 2022.