Sum of numbers which are not sum of two abundant numbers

Question

The following program is a solution to Project Euler Problem 23. That problem defines an abundant number to be one whose sum of proper divisors is greater than it, tells us that all integers at least 28124 can be written as a sum of two abundant numbers and asks for the sum of all numbers that cannot be written as a sum of two abundant numbers.

I've tried to to follow PEP8, break my code up into functions so that I can unit test, and generally write "production-quality" code. I'd be very interested in any advice on that front. I wasn't really concentrating on making it the fastest solution but comments on speed/algorithm improvements would be welcome too.

ProjectEuler23.py:

""" Project Euler Problem 23: find the sum of all numbers which cannot be
written as the sum of two abundant numbers. An abundant number is one whose
proper divisors sum to be greater than the number. We are given that all
integers at least 28124 can be written as the sum of two abundant numbers.

Creator: Martin Leslie
Date created: 10/18/14
Date modified: 10/21/14
"""

def sum_of_divisors(input_integer):
    """ return sum of divisors of positive integer.
    For example, sum_of_divisors(12) == 28 because divisors of 12 are
    1, 2, 3, 4, 6, 12
    """
    sum_so_far = 0
    low_divisor = 1
    while low_divisor**2 <= input_integer:
        if input_integer % low_divisor == 0:
            high_divisor = input_integer / low_divisor
            sum_so_far += low_divisor
            if high_divisor != low_divisor:
                sum_so_far += high_divisor
        low_divisor += 1
    return sum_so_far

def list_abundant_numbers(upper_limit):
    """ return list of all abundant numbers < upper_limit"""
    abundants = []
    for n in xrange(1, upper_limit):
        if sum_of_divisors(n) > 2*n:
        # abundant if sum of proper divisors > n,
        # this is equivalent to sum of divisors > 2*n
            abundants.append(n)
    return abundants

def set_of_sums_of_pairs(input_list, max_sum=float('inf')):
    """ return set of all distinct n+m for m,n in input_list with m+n < max_sum
    """
    set_of_sums = set()
    for n in input_list:
        for m in input_list:
            if m+n < max_sum:
                set_of_sums.add(m+n)
    return set_of_sums

def sum_of_one_up_to_n(n):
    """ return 1+2+...+n = n(n+1)/2"""
    return n*(n+1)/2

UPPER_LIMIT = 28124
# integers at least 28124 can be written as the sum of two abundant numbers

def main():
    """ compute total of all numbers < UPPER_LIMIT and subtract the
    numbers < UPPER_LIMIT which are sums of two abundant numbers
    """
    total_of_all_numbers = sum_of_one_up_to_n(UPPER_LIMIT-1)
    abundants = list_abundant_numbers(UPPER_LIMIT)
    set_of_sums_of_abundants = set_of_sums_of_pairs(abundants, UPPER_LIMIT)
    total_of_not_abundunt_sums = (total_of_all_numbers -
                                  sum(list(set_of_sums_of_abundants)))
    print ("Sum of numbers not a sum of two abundant numbers is " +
           str(total_of_not_abundunt_sums))

if __name__ == "__main__":
    main()

ProjectEuler23Test.py:

""" Unit tests for Project Euler problem 23.

Creator: Martin Leslie
Date created: 10/21/14
Date modified: 10/21/14
"""

from ProjectEuler23 import (sum_of_divisors, list_abundant_numbers,
                            set_of_sums_of_pairs, sum_of_one_up_to_n)
import unittest

class TestProjectEuler23Functions(unittest.TestCase):

    def test_sum_of_divisors(self):
        """ Test sum_of_divisors agrees with https://oeis.org/A000203"""
        self.assertEqual(sum_of_divisors(1), 1)
        self.assertEqual(sum_of_divisors(2), 3)
        self.assertEqual(sum_of_divisors(4), 7)
        self.assertEqual(sum_of_divisors(12), 28)
        self.assertEqual(sum_of_divisors(24), 60)

    def test_list_abundant_numbers(self):
        """ Test list_abundant_numbers follows https://oeis.org/A005101"""
        self.assertEqual(list_abundant_numbers(60),
                         [12, 18, 20, 24, 30, 36, 40, 42, 48, 54, 56])
        self.assertEqual(list_abundant_numbers(61),
                         [12, 18, 20, 24, 30, 36, 40, 42, 48, 54, 56, 60])

    def test_set_of_sums_of_pairs(self):
        self.assertEqual(set_of_sums_of_pairs([]), set())
        self.assertEqual(set_of_sums_of_pairs([1]), {2})
        self.assertEqual(set_of_sums_of_pairs([1, 2, 4]), {2, 3, 4, 5, 6, 8})
        self.assertEqual(set_of_sums_of_pairs([1, 2, 4], 6), {2, 3, 4, 5})

    def test_sum_of_one_up_to_n(self):
        self.assertEqual(sum_of_one_up_to_n(0), 0) # empty sum
        self.assertEqual(sum_of_one_up_to_n(1), 1)
        self.assertEqual(sum_of_one_up_to_n(2), 3)
        self.assertEqual(sum_of_one_up_to_n(50), 1275)

if __name__ == '__main__':
    unittest.main()

Answer 1

Your code looks good and is properly tested but you can easily make it better.

Also, your trick about sums of numbers is a really nice touch.

Separation of concerns and reusable functions

You should try to make your code as easy as possible to re-use. For instance, having a sum_of_divisors(n) is cool but if you had defined a divisors(n) function yielding divisors then you'd just have to call sum(divisors(n)).

def divisors(input_integer):
    """Please comment me."""
    low_divisor = 1
    while low_divisor**2 <= input_integer:
        if input_integer % low_divisor == 0:
            high_divisor = input_integer / low_divisor
            yield low_divisor
            if high_divisor != low_divisor:
                yield high_divisor
        low_divisor += 1


def sum_of_divisors(input_integer):
    """ return sum of divisors of positive integer.
    For example, sum_of_divisors(12) == 28 because divisors of 12 are
    1, 2, 3, 4, 6, 12
    """
    return sum(divisors(input_integer))

Then you can compute the limit for your loop easily :

def divisors(input_integer):
    """Please comment me."""
    for low_divisor in range(1, int(math.sqrt(input_integer)) + 1):
        if input_integer % low_divisor == 0:
            high_divisor = input_integer / low_divisor
            yield low_divisor
            if high_divisor != low_divisor:
                yield high_divisor

(there might be off-by-one errors in this code).

By the way, sum works on any iterable so you don't have to convert set_of_sums_of_pairs to a list first. It's a pretty bad idea from a performance point of view.

A better algorithm

Because you already know an upper limit for the number you'll consider, instead of going through numbers one at a time and try to compute its decomposition, you might as well perform some variation around the sieve of Eratosthene :

def divisors_sieve(lim):
    """Computes the list of divisors for values up to lim included."""
    div = [[1] for i in range(lim + 1)]
    for i in range(2, 1 + lim // 2):
        for j in range(2 * i, lim + 1, i):
            div[j].append(i)
    return div

I've not included n as a divisor of n as it didn't seem too relevant.

Then you can easily get the abundant numbers by doing something like :

abun = [i for i, l in enumerate(divisors_sieve(lim)) if i > 0 and sum(l) > i]

Using Python good stuff

I have to go but you could easily write set_of_sums_of_pairs with a set comprehension.

However, something else is to be improved : at the moment, you are iterating over all pairs (n,m) but iterating over pairs such that n<=m will be enough. This can be provided by itertools.combinations_with_replacement :

set_sum = { n+m for n,m in itertools.combinations_with_replacement(abun, 2) if n+m <lim}

Alternatively, if you were to keep this in a function with a conventional loop, you could add break when you know that the max value is reached and the current loop will not bring any potential solution.