# python - What is the best algorithm to solve this puzzle?

So I came across this question:

How many numbers are there from 1 to 1000 which are not divisible by the digits 2, 3 and 5?

It seems pretty easy at first, so I wrote a quick python program to solve it:

``count = 0for number in range(1,1000):if number % 2 != 0 and number % 3 != 0 and number % 5 != 0:count += 1print(count)``

I got the correct answer (266), but I thought that doing it that way was a lot of typing if I ever wanted to check more than just 3 values. I also wanted to do a mathematical solution so I came across this:

`1000 - ((1000/2 +1000/3 +1000/5) -(1000/2x3 +1000/2x5 + 1000/3x5)+ (1000/2x3x5)) = 1000-((500+333+200) - (166 +100 + 66) + 33) = 1000- 734 = 266`

I thought it was a good approach so I implemented it in code:

``def foo(ln = 1000), numbers = [2,3,5]:div = 0muldiv = 0totdiv = 1for n in numbers:div += ln/nfor i in numbers:for n in range(numbers.index(i)+1, len(numbers)):muldiv += ln/(i * numbers[n])for n in numbers:totdiv *= nanswer = ln - (div - muldiv + ln/totdiv)print("answer is ", math.floor(answer))``

Now I am pretty sure I messed up somewhere in my second function because it doesn't seem to work for more numbers. For example, if I were to try to find

How many numbers are there from 1 to 1000 which are not divisible by the digits 2, 3, 5 and 7?

the first method returns `228` and `foo(numbers = [2,3,5,7])` returns `300`... I'm pretty sure `228` is the correct answer since one more number would mean that there are LESS factors instead of more, but where did I go wrong? and is there a better way to solve this problem?

You do not need an algorithm for that, simple mathematics is enough:

Say you want to count the amount of numbers from 1 to N (inclusive) dividable by k, that is simply equivalent to:

floor(N/k).

So the amount of numbers dividable by 3 in this case is 333.

Now you can't however simply use calculate the amount of numbers dividable by 2, 3 and 5; and sum them up, because there are common ones. Indeed: for instance 15 is dividable by both 3 and 5.

You can solve this however using the inclusion-exclusion principle:

the amount of numbers dividable by 2, 3 and 5 is the same as

• the amount numbers dividable by 2
• plus the amount of numbers dividable by 3
• plus the amount of numbers dividable by 5
• minus the amount of numbers dividable by 2 and 3
• minus the amount of numbers dividable by 2 and 5
• minus the amount of numbers dividable by 3 and 5
• plus the amount of numbers dividable by 2, 3 and 5.

So in order to solve your first problem, you can simply state:

``````def n_div(N,k):
return N//k

def n_div235(N):
return n_div(N,2)+n_div(N,3)+n_div(N,5)-n_div(N,2*3)-n_div(N,2*5)-n_div(N,3*5)+n_div(N,2*3*5)

def not_div235(N):
return N-n_div235(N)
``````

As you can see it generates the correct result:

``````>>> not_div235(1000)
266
``````

As long as N is very large compared to the number of divisors, you better use the inclusion-exclusion approach:

you can do this like:

``````import itertools
from functools import reduce
import operator

def gcd(a, b):
while b:
a, b = b, a % b
return a

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

def lcm_list(ks):
res = 1
for k in ks:
res = lcm(res,k)
return res

def n_div_gen(N,ks):
nk = len(ks)
sum = 0
factor = 1
for i in range(1,nk+1):
subsum = 0
for comb in itertools.combinations(ks,i):
subsum += n_div(N,lcm_list(comb))
sum += factor * subsum
factor = -factor
return sum

def not_div_gen(N,ks):
return N-n_div_gen(N,ks)
``````

For small N, this will not pay off, but say to want to calculate the amount of numbers dividable by 3, 5 and 7 from 1 to 1 000 000 000 is:

``````>>> not_div_gen(1000000000,[3,5,7])
457142857
``````

You can do this with:

``````>>> sum(i%3!=0 and i%5!=0 and i%7!=0 for i in range(1,1000000001))
457142857
``````

But it takes minutes to calculate that whereas our own approach uses milliseconds. Note that this only works for a huge N.

Use the built-in functions `sum` and `all` with a nested generator:

``````def f(r=1000, nums=(2,3,5)):
return sum(all(x%n for n in nums) for x in range(1, r+1))
``````

This goes through the range of numbers, check whether each of those numbers has a nonzero modulus with each of the specified numbers, and sums those boolean values (`False` is 0 and `True` is 1). A `nums` of `(2,3,5,7)` produces a result of 228, which is in agreement with your shorter, simpler code (which, reassuringly, doesn't use any floating-point arithmetic, as your second code block does).

The number of integers up to N not divisible by n1,n2,...,nt (assumed to be pairwise-coprime) is

the number of integers up to N minus
( SUMi in 1..t ( the number of integers up to N divisible by ni)) plus
( SUMi,j in 1..t, i<j ( the number of integers up to N divisible by ninj)) minus
( SUMi,j,k in 1..t, i<j<k ( the number of integers up to N divisible by ninjnk)) plus
( SUMi,j,k,l in 1..t, i<j<k<l ( the number of integers up to N divisible by ninjnknl)) minus
... ... ... ...
( SUMi,j,k,l,...q in 1..t, i<j<k<l<...<q ( the number of integers up to N divisible by ninjnknl...nq))

The series continues until the subscript contains all t integers from the original list.

For numbers that are not known to be pairwise-coprime, replace their product by the least common multiple.

This is why your method works only for 3 numbers. You only compute the first four members of the series.

Here's another implementation that uses inclusion-exclusion. It's simpler than the code in Willem Van Onsem's excellent answer (which I didn't see before I wrote this code), but this one only works if the numbers in the list of divisors are all coprime to each other. For the the more general case, you need to use Willem's approach.

``````from itertools import combinations
from functools import reduce

def prod(seq, mul=int.__mul__):
return reduce(mul, seq, 1)

def count_coprimes(n, divisors):
total = n
sign = -1
for i in range(1, len(divisors) + 1):
for k in combinations(divisors, i):
total += n // prod(k) * sign
sign = -sign
return total

print(count_coprimes(1000, [2, 3, 5]))
``````

output

``````266
``````

FWIW, here's the same algorithm as a "one-liner" (split over several lines to improve readability). It's a little less efficient due to the `(-1)**i` in the inner loop.

``````def count_coprimes(n, divisors):
return n + sum(n // prod(k) * (-1)**i
for i in range(1, len(divisors) + 1)
for k in combinations(divisors, i))

print(count_coprimes(1000000000, [3, 5, 7]))
``````

output

``````457142857
``````

We can get rid of that `(-1)**i` by negating the divisors and using a modified integer division function:

``````def div(p, q):
return p // q if q > 0 else -(p // -q)

def count_coprimes(n, divisors):
return sum(div(n, prod(k))
for i in range(len(divisors) + 1)
for k in combinations([-u for u in divisors], i))
``````

A very small change you can make to roughly halve the amount of time that it takes is rather than generating all numbers from 1 to 1000, generate all odd numbers from 1 to 1000:

``````count = 0
for number in range(1,1001,2):
if number % 3 != 0 and number % 5 != 0:
count += 1
print(count)
``````

While this is not a huge change and is not a mathematical solution it makes the code no less readable and slightly more efficient.

And to take into account your other code, you can use a list comprehension in the if statement to check other numbers(note that I also use the first number to generate the initial list of numbers, rather than performing the modulo operation for all 1000):

``````def foo(lst):
count = 0
for number in range(1,1001,lst[0]):
if not any([number % i == 0 for i in lst[1:]]):
count += 1
return count

>>> foo([2,3,5])
266
>>> foo([2,3,5,7])
228
``````

There are many ways to solve this little problem iteratively, all of them with pretty much similar performance, here's few examples:

``````import timeit

def f1(l, h):
count = 0
for number in range(l, h):
if number % 2 != 0 and number % 3 != 0 and number % 5 != 0:
count += 1

return count

def f2(l, h):
return len(filter(lambda x: x % 2 != 0 and x % 3 != 0 and x % 5 != 0, range(l, h)))

def f3(l, h):
count = 0
for number in range(l, h):
if number % 2 == 0:
continue
if number % 3 == 0:
continue
if number % 5 == 0:
continue

count += 1

return count

def f4(l, h):
return len([x for x in range(l, h) if x % 2 != 0 and x % 3 != 0 and x % 5 != 0])

a, b, N = 1, 1000, 10000

print timeit.timeit('f1(a,b)', setup='from __main__ import f1, a, b', number=N)
print timeit.timeit('f2(a,b)', setup='from __main__ import f2, a, b', number=N)
print timeit.timeit('f3(a,b)', setup='from __main__ import f3, a, b', number=N)
print timeit.timeit('f4(a,b)', setup='from __main__ import f4, a, b', number=N)
``````

Times on a i7-2.6ghz would be:

``````0.802361558825
1.46568073638
0.91737188946
0.846404330893
``````

Usually these times are good enough to be considered when the lower/upper bounds (1,1000) are relatively small. Now, if we're talking about really high bounds (trillions) where the computation is not feasible you could consider apply the much smarter inclusion-exclusion principle, that way you'd solve the problem analitically and you'd be granted to get constant time with your solution.

1. input

let `n` be the interval `<0,n>` to test and `d[]={2,3,5,0};` be null terminated array of divisors

2. compute LCM of `d[]`

least common multiple is the period with which the SoE will repeat itself. for `2,3,5` is `lcm=30`. Use the `max(d[])` as increment while computing it to boost speed... If LCM is too big (`LCM>=n`) then use `n` instead for speed.

3. compute SoE for `<0,LCM)`

simply create array of LCM numbers and set `a[i]=1` for non divisible `i` and `a[i]=0` for divisible `i`.

4. convert SoE to non divisible number count

simply compute `a'[i]=a[0]+a[1]+..a[i]`

5. compute the count

count is simply:

``````int(n/LCM)*a'[LCM-1] + a'[n%LCM];
``````

Here simple C++ example:

``````int non_divisibles(int n,const int *d)  // SoE
{
int i,j,cnt,lcm,m,*a;
for (m=0;d[m];m++); // compute m
for (cnt=0,i=0;i<m;i++) if (cnt<d[i]) cnt=d[i]; // cnt = max d[] (to speed up LCM)
// LCM d[]
a=new int[m]; if (a==NULL) return -1;
for (i=0;i<m;i++) a[i]=d[i];
for (lcm=cnt;lcm<=n;lcm+=cnt) // no need to test above `n` even if lcm is bigger
{
for (i=0;i<m;i++) for (;a[i]<lcm;) a[i]+=d[i];
for (i=0;i<m;i++) if (a[i]!=lcm) { i=-1; break; }
if (i>=0) break;
}
delete[] a;
// SoE <0,LCM)
a=new int[lcm]; if (a==NULL) return -1;
for (i=0;i<lcm;i++) a[i]=1;
for (j=0;d[j];j++)
for (i=0;i<lcm;i+=d[j])
a[i]=0;
// convert to cnt
for (i=1;i<lcm;i++) a[i]+=a[i-1];
// compute whole count
cnt =(n/lcm)*a[lcm-1];
cnt+=a[n%lcm];
delete[] a;
return cnt;
}
``````

Here some measurements to compare naive,SoE and this SoE(max(`n`,LCM(`d[]`))) approaches:

``````n=1000000 d[]={ 2 3 5 7 11 13 17 19 }
171021 [  27.514 ms] naive
171021 [  12.642 ms] SoE
171021 [  25.438 ms] LCM+Soe

n=1000000 d[]={ 2 3 5 7 11 13 17 }
180524 [  26.212 ms] naive
180524 [  11.781 ms] SoE
180524 [   9.807 ms] LCM+Soe

n=1000000 d[]={ 2 3 5 7 11 13 }
191808 [  24.690 ms] naive
191808 [  11.512 ms] SoE
191808 [   0.702 ms] LCM+Soe

n=1000000 d[]={ 2 3 5 }
266666 [  16.468 ms] naive
266666 [   9.744 ms] SoE
266666 [   0.006 ms] LCM+Soe

n= 1000 d[]={ 2 3 5 }
266 [   0.012 ms] naive
266 [   0.012 ms] SoE
266 [   0.001 ms] LCM+Soe

n=1000000 d[]={ 2 3 5 19 23 61 87 10001 }
237662 [  26.493 ms] naive
237662 [  10.180 ms] SoE
237662 [  19.429 ms] LCM+Soe
``````

As you can see SoE(n) is better if LCM is too big in comparison to `n`(`d[]` contains many primes or big numbers) but need `O(n)` space.

You can bail out as soon as you hit a divisible.

``````    def test(tocheck):

count = 0

for number in range(1, 1000):

for div in tocheck:
if not number % div:
break
else:
#else on a loop means it ran to completion w.o. break
count += 1

print("%d not divisible by %s" % (count, tocheck))

test([2,3,5])
test([2,3,5,7])
``````

output:

``````266 not divisible by [2, 3, 5]
228 not divisible by [2, 3, 5, 7]
``````