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)
