# Generalized Project Euler 1: A sledgehammer to crack a nut

## The problem

Project Euler 1 is one of the most asked questions on site. However I wanted to solve the more general problem of division.

Multiples of a list

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

Find the sum of all the multiples in a list from start to stop, where start and stop are natural numbers.

In the classical problem, this list would be [3, 5]. However I wanted to create a code that could solve [3, 2], [3, 2, 7] and really any combination of numbers. The naive solution is shown below

def sum_divisible_naive(divisors=[3, 5], start=0, stop=100):
count = 0
for num in xrange(start, stop):
for d in divisors:
if num % d == 0:
count += num
break
return count


I do not want comment on this code, but include it as a way to test the correctness of my implementation. The above code is correct, but has a linear $\mathcal{O}(n)$ running time. I wanted to create something faster. I did however run into some problems which made my code longer than wanted.

## Improved algorithm

To find the number of numbers divisible by some number k between start and stop, I use

def divisible_by_k(k, start=0, limit=100):
stop = int((limit-1)/float(k))
return k*(stop+start)*(stop-start+1)/2


This is just a slightly modified version of the sum of the first n natural numbers. Take [2, 3] as an example. If I were to do

total = divisible_by_k(3) + divisible_by_k(5)


I would count too many numbers since there are some numbers which occur in both of the lists: [3, 6, 9, 12, 15] and [5, 10, 15]. To correct this we have to subtract all the multiples of $15$.

Another problem is with [6, 9], if we did it like above we would remove all multiples of $6*9 = 54$, however $18$ is the first number that occurs in both lists. A solution to this is to divide the product by the gcd. Eg $6\cdot9/\gcd(6,9)=56/3 = 18$ as wanted.

total = divisible_by_k(a) + divisible_by_k(b) - sum_divisible_by_k(a*b/gcd(a, b))


I wanted to generalize this further to allow for a list of such numbers. Say we have [2, 3, 8], there is no need to include $8$ because all multiples of $8$ are already counted in $2$. I remove the unnecessary multiples with the following code

def remove_multiples(divisors):
i = 0
divisors = sorted(set(divisors))
divisors_len = len(divisors)
while i < divisors_len:
j = i + 1
while j < divisors_len:
if divisors[j] % divisors[i] == 0:
divisors.remove(divisors[j])
divisors_len -= 1
else:
j += 1
i += 1
return divisors


So remove_multiples([2, 3, 6, 8) return [2, 3].

Take [2, 3, 5] as the next example. We now have to subtract $2 \cdot 3$, $2 \cdot5$ and $3 \cdot 5$ to avoid the double counting. However we now remove too many numbers and have to include $2 \cdot 3 \cdot 5$. Hopefully this becomes clearer when studying the following image

The general case of n numbers can be found by studying the previous examples. We can get the combinations of k<=n numbers using the combinations module from collections. Then I just have to switch between adding and subtracting.

This means I also had to find the gcd of a list and not just a pair of numbers. The gcd of a list is computed many times because we go through every combination. I thought making a dict would be a good idea to lessen the work. Since perm is sorted a will always be less than b

gcd_dict = defaultdict(int)

def gcd_list(perm):
perm_len = len(perm)
a = perm[0]
for i in range(1, perm_len):
b = perm[i]
if gcd_dict[(a, b)] == 0:
gcd_dict[(a, b)] = gcd(a, b)
a = gcd_dict[(a, b)]
return a


I have not put much though into the naming nor adding docstrings or comments, because the code is not quite complete.

## Questions:

• Is there some way to write this code clearer, and shorter?

## The code

from collections import defaultdict
from itertools import combinations
from operator import mul
from fractions import gcd

gcd_dict = defaultdict(int)

def gcd_list(perm):
perm_len = len(perm)
a = perm[0]
for i in range(1, perm_len):
b = perm[i]
if gcd_dict[(a, b)] == 0:
gcd_dict[(a, b)] = gcd(a, b)
a = gcd_dict[(a, b)]
return a

def remove_multiples(divisors):
i = 0
divisors = sorted(set(divisors))
divisors_len = len(divisors)
while i < divisors_len:
j = i + 1
while j < divisors_len:
if divisors[j] % divisors[i] == 0:
divisors.remove(divisors[j])
divisors_len -= 1
else:
j += 1
i += 1
return divisors

def sum_divisible_by_k(k, start, limit):
stop = int((limit-1)/float(k))
return k*(stop+start)*(stop-start+1)/2

def sum_divisible_fast(div, start, stop):
divisors = remove_multiples(div)

total = 0
for divisor in divisors:
total += sum_divisible_by_k(divisor, start, stop)

k = -1
for i in range(2, len(divisors)+1):
for perm in combinations(divisors, i):
product = reduce(mul, perm)/gcd_list(perm)
total += k*sum_divisible_by_k(product, start, stop)
k *= -1
return int(total)

def sum_divisible_naive(divisors=[3, 5], start=0, stop=100):
count = 0
for num in xrange(start, stop):
for d in divisors:
if num % d == 0:
count += num
break
return count

if __name__ == '__main__':

list_to_test = [2, 6, 9, 8]
start = 1
stop = 10**5
print sum_divisible_naive(list_to_test, start, stop)
print sum_divisible_fast(list_to_test, start, stop)


# Background

Let $N$ be the size of the range [start, stop). $N = \textrm{stop} - \textrm{start}$.

Let $M$ be the size of the divisor list. $M$ = $\lvert \textrm{divisors} \rvert$

The naive method is $\Omega(N)$ and your method is $\mathcal{O}(2^M)$ since you must iterate through the powerset of the divisor list.

# Variable naming

I am aware that k is just a generic iterator name but the following line is confusing to me:

total += k*sum_divisible_by_k(product, start, stop)


Use a different generic name for your sign-switching iterator such as j, s, or basically anything other than k.

# Errors

In sum_divisible_by_k using start directly is incorrect. Instead you should be computing $\lceil \frac{start}{k} \rceil$. Also, for stop you can just use integer division $\lfloor \frac{limit-1}{k} \rfloor$.

# Anti-patterns

Iterating through a list while removing values is an anti-pattern in most languages I know. As mentioned in another answer here I would recommend building a new list.

# Memoization

Your dict of previously computed gcd's is suboptimal on two fronts.

• You are only storing the previous value of the gcd so you still have to perform the multiplications for each combination of divisors. For example in the divisor set $\{a,b,c,d,e\}$ you will compute $a \cdot b\textrm{, } a \cdot b \cdot c\textrm{, } a \cdot b \cdot c \cdot d\textrm{, } a \cdot b \cdot c \cdot d \cdot e$. You should be storing the actual lcm not just the intermediary step of calculating the gcd.

• $\textrm{lcm}(a,b,c) = \textrm{lcm}(\textrm{lcm}(a,b),c)$ as I am sure you are aware. However, when you are computing the $\textrm{lcm}(a,b,c,d,e)$ you are performing a lookup of $\textrm{lcm}(a,b)$ then $\textrm{lcm}(\textrm{lcm}(a,b),c)$ then $\textrm{lcm}(\textrm{lcm}(\textrm{lcm}(a,b),c), d)$, etc. Instead of performing $L-1$ lookups for the lcm of a list of length $L$ you can perform $1$ lookup.

From itertools:

Combinations are emitted in lexicographic sort order. So, if the input iterable is sorted, the combination tuples will be produced in sorted order.

Assuming you start with the smallest combinations first (which you currently are) you can be ensured that the current combination with the last element removed is already in the dict. That is, when determining $\textrm{lcm}(a,b,c,d)$, we can always grab $\textrm{lcm}(a,b,c)$ from the dict.

# Combining logic

Often in mathematics values are given to certain expressions for convenience. For instance, $0!=1$. A convenient value for you would be to define $\textrm{lcm}(\textrm{ () })=1$ so that $\textrm{lcm}(\textrm{ () }, a) = a$. You could then combine your two outer loops in sum_divisors_fast into one loop.

# Updated Code

from itertools import combinations
from fractions import gcd

lcm_dict = {():1}

def lcm(a, b):
return a * b // gcd(a, b)

def lcm_list(perm):
l = lcm(lcm_dict[perm[:-1]], perm[-1])
lcm_dict[perm] = l
return l

def remove_multiples(divisors):
# Directly from Joe Wallis's answer
new_divisors = []
divisors = sorted(set(divisors))
for divisor in divisors:
if not any(divisor % d == 0 for d in new_divisors):
new_divisors.append(divisor)
return new_divisors

def ceil_div(a, b):
return -(-a // b)

def summation_from_zero(n):
return n * (n+1) // 2

def sum_divisible_by_k(k, start, stop):
sum_start = ceil_div(start, k)
sum_stop = (stop - 1) // k
return k * (summation_from_zero(sum_stop) - summation_from_zero(sum_start-1))

def sum_divisible_fast(div=[3,5], start=0, stop=100):
divisors = remove_multiples(div)
total = 0

j = 1
for i in xrange(1, len(divisors)+1):
for perm in combinations(divisors, i):
product = lcm_list(perm)
total += j * sum_divisible_by_k(product, start, stop)
j = -j

def sum_divisible_naive(divisors=[3, 5], start=0, stop=100):
count = 0
for num in xrange(start, stop):
for d in divisors:
if num % d == 0:
count += num
break
return count

if __name__ == '__main__':

list_to_test = [2, 6, 9, 8]
start = 2
stop = 10**5
print sum_divisible_naive(list_to_test, start, stop)
print sum_divisible_fast(list_to_test, start, stop)

print sum_divisible_naive()
print sum_divisible_fast()


# Timings

• Original

timeit.timeit('orig.sum_divisible_naive([3,5,7,11,13,17,19,23,29,31,37,39,41], 0, 10**7)', 'from __main__ import orig', number=1)
5.059425115585327
timeit.timeit('orig.sum_divisible_fast([3,5,7,11,13,17,19,23,29,31,37,39,41], 0, 10**7)', 'from __main__ import orig', number=10)
0.2042551040649414
timeit.timeit('orig.sum_divisible_fast([3,5,7,11,13,17,19,23,29,31,37,39,41], 0, 10**7)', 'from __main__ import orig', number=100)
1.8143179416656494

• Updated

timeit.timeit('updated_orig.sum_divisible_fast([3,5,7,11,13,17,19,23,29,31,37,39,41], 0, 10**7)', 'from __main__ import updated_orig', number=10)
0.14538812637329102

timeit.timeit('updated_orig.sum_divisible_fast([3,5,7,11,13,17,19,23,29,31,37,39,41], 0, 10**7)', 'from __main__ import updated_orig', number=100)
1.0789260864257812

• I really do not think the naive method is as slow as $\mathcal{O}(MN)$. That is because we use a if loop and stop when we reach a divisble number. So half of every number is divisible by 2, and hence no more tests is done on those numbers. I think one has to look at the average number of loops the function has to go through to find the excact running time. Other than that fantastic answer! – N3buchadnezzar Jul 5 '16 at 9:11
• @N3buchadnezzar You are absolutely right. And to make matters worse I said "of course" its $\mathcal{O}(MN)$. Unfortunately I am not sure right now what we can say about the runtime of this algorithm without knowing the exact distribution of the divisors in $M$. – twohundredping Jul 5 '16 at 12:46
• I think $\mathcal{O}(2^M)$ versus $\mathcal{O}(N)$ is a fair comparison. One is dominated by the number of divisors while the other one is dominated by the length of our interval. I think that is the most important part, and the rest is nitpicking ^^ – N3buchadnezzar Jul 5 '16 at 13:39
• I would prefer to say $\omega(\textrm{max}(N, \frac{NM}{logN}))$ but yes I agree that specifying a better bound is nitpicking, particularly for code review. I will save the pedantry for more suitable networks. – twohundredping Jul 5 '16 at 14:40

Your maths is much better than mine, and so I can't comment on it.

However your function remove_multiples can be simplified if you come at the problem in a different way. It's normally easier to add to a new list, than to remove from an old one.

The way this works is rather than checking if all the numbers after the number you are adding are ok, you instead check if the number is ok, by looking at the numbers before it. Then you can achieve a simpler function:

def remove_multiples(divisors):
new_divisors = []
divisors = sorted(set(divisors))
for divisor in divisors:
if not any(divisor % d == 0 for d in new_divisors):
new_divisors.append(divisor)
return new_divisors


The sorted(set(divisors)) from your old code allows us to say that any new_divisors is smaller than any divisors. Unless you pass negatives to it. From this we can go through the divisors, and if it's divisible by any of the new_divisors, then it's to be removed.

gcd_dict / gcd_list stuff could be more expressive.

For example, since you are touching memoization, you could say so explicitly. Simple decorator that does the job would take less space than your gcd_list and it will be reusable.

Take a look at this one - sure, it's an overkill for one task, but you'll find some uses for memoizing decorator later on if you continue with Project Euler tasks so investing time in good one would pay off.

remove_multiples was already covered by another poster. I'll add that it could be shortened with some crafty reduce instead of iteration but I can't say that it definitely would be an improvement.

You have 2 pieces of code that increase/decrease total with sum_divisible_by_k. Sure, it would be simpler for reader to understand what first part does but reader would have to make sense of second part anyway. Thus, you can replace range(2, len(divisors)+1) with range(1, len(divisors)+1) and set initial k value to 1 - then you can get rid of for divisor in divisors: cycle altogether.

Also, meaning of k is not obvious. If it means '-1 for even-sized permutations, 1 for odd-sized', it should say so. For example 1 if (len(perm) % 2) else -1. Don't worry, len() is fast.

product could be moved to a separate function if you don't find reduce's an eye candy, but that's just a matter of personal preference.

I think there is an additional speed up. Let N be stop - start, and k be the number of divisors. Consider the case where k << N, say divisors = {3,5}, start = 1, stop = 1Million.

We don't need to perform $O(N)$ lookups, we can in $O(1)$ calculate the sum of multiples of 3. It's 3*(1+2+3+4...floor(1 Million/3)) which reduces to a well known 3*n(n-1)/2 formula.

Now your problem takes only $O(k^2)$ instead of $O(N)$ which is likely much much faster.

• I already does this. I also improved it such that I can use this technique on a list longer than two digits. The question is not about this trick, but rather on my solution of it. – N3buchadnezzar Jul 1 '16 at 18:03