Problem: Find the sum of all natural numbers below 1000 which are divisible by 3 or 5.
The basic solution is to iterate over all numbers and check each:
L = 1000
total = 0
for i in range(1,L):
if i % 3 == 0 or i % 5 == 0:
total += i
print(total)
Something a little more advanced we can do is use the summations. For example, we would add the following multiples of 3:
\[3+6+9+\ldots+996+999=3\left(1+2+3+\ldots+333\right) =3\sum_{i=1}^{333}i=3\left(\frac{333\cdot334}{2}\right)\]We do the same for multiples of 5. But everything that is a multiple of \(\lcm(3,5)=15\) is counted twice so we must subtract those. This part is doable by hand calculation. Let \(S(d,L)= d\sum_{i=1}^{\lfloor L/d\rfloor}i=d\left(\frac{(\lfloor L/d\rfloor) (\lfloor L/d\rfloor+1)}{2}\right)\) be one such summation. Then the solution is:
\[S(3,999)+S(5,999)-S(15,999)\]Note: 999 is used instead of 1000 because the problem asks for numbers below 1000 so we should not count 1000.
Problem: Find the sum of all even Fibonacci numbers up to 4000000.
The Fibonacci sequence starts with \(F_0=0,F_1=1\) and then is recursively defined \(F_n=F_{n-1}+F_{n-2}\) when \(n\geq2\). The simple way to solve this is iterating this increasing sequence until going past 4000000:
L = 4000000
a,b = 0,1
total = 0
while True:
a,b = b,a+b
if b > L:
break
if b % 2 == 0:
total += b
print(total)
This is a bit tedious but doable by hand. We can actually be even faster by observing that indexes \(0,3,6,9,\ldots\) of the Fibonacci sequence are even:
\[\red{0},1,1,\red{2},3,5,\red{8},13,21,\red{34},55,89,\ldots\]Using the recursive definition, we can write \(F_n\) in terms of \(F_{n-3},F_{n-6}\):
\[\begin{align} F_n&=F_{n-1}+F_{n-2}=(F_{n-2}+F_{n-3})+(F_{n-3}+F_{n-4})\\ &=2F_{n-3}+F_{n-2}+F_{n-4}=2F_{n-3}+(F_{n-3}+F_{n-4})+(F_{n-5}+F_{n-6})\\ &=3F_{n-3}+(F_{n-4}+F_{n-5})+F_{n-6}=4F_{n-3}+F_{n-6}\\ \end{align}\]Now we can solve the problem in almost the same way, but we start with \(F_0=0,F_3=2\) and use this recurrence which generates only even terms:
L = 4000000
a,b = 0,2
total = 2
while True:
a,b = b,4*b+a
if b > L:
break
total += b
print(total)
Note: the largest term included in the sum is \(F_{33}=3524578\).
Problem: Find the largest prime factor of 600851475143.
One way to solve this problem is by using the factor command
which is included in some Linux distros. But if we challenge ourselves to
write some actual code, here is how we can do it. Let \(n\) be the number to
factor. Test divisors \(d=2,d=3,\ldots\) to find a factor. If \(n\) has a
factor, one of them will be at most \(\lfloor\sqrt{n}\rfloor\) so we can
stop as soon as \(d^2>n\). Each time we find a factor, we divide \(n\) until
it it no longer divisible by that factor. At the end, we may have \(n=1\)
or if \(n\neq1\) then it is a prime.
n = 600851475143
d = 2
dmax = 1
while d*d <= n:
while n % d == 0:
n //= d
dmax = max(dmax,d)
print(f'factor {d}')
d += 1
dmax = max(dmax,n) # if n!=1 it is prime
if n != 1:
print(f'factor {n}')
print(dmax)
Factoring methods can get a whole lot more complicated, but we'll leave it here for this problem. This method is reasonably efficient for 12 digit numbers.
Problem: Find the largest palindrome that is a product of 2 3-digit numbers.
The basic solution is to have a double loop over all possible 3 digit numbers. It is fast enough for this problem and we can check if numbers are palindrome by converting them to strings.
largest = 0
p,q = 0,0
for a in range(100,1000):
for b in range(100,1000):
n = a*b
s = str(n)
if n > largest and s == s[::-1]:
largest = n
p,q = a,b
print(f'product {p} * {q}')
print(largest)
It can be done faster by starting from the largest and terminating the loops when it is impossible to multiply to a bigger number. Additionally, we can generalize the code a bit by using a variable for the number lengths. Lastly, due to symmetry, we can require \(a\geq b\) and start the inner loop at \(a\) instead of 999.
D = 3
largest = 0
p,q = 0,0
for a in range(10**D-1,10**(D-1)-1,-1):
if a * (10**D-1) < largest:
break # stop when a*b cannot exceed largest
for b in range(a,10**(D-1)-1,-1):
n = a*b
if n < largest:
break
s = str(n)
if s == s[::-1]:
largest = n
p,q = a,b
print(f'product {p} * {q}')
print(largest)
Problem: Find the smallest number divisible by all of \(1,2,3,\ldots,20\).
It is slow, but we can have a basic loop that checks divisibility for every integer counting up. We could even go up in increments of 20, but that would still be very slow. The efficient way to solve this would be using prime factorizations, which is doable by hand. We factor all the numbers \(1,2,3,\ldots,20\) and then construct a new number that has each of the prime factors with minimum multiplicity.
Then from here, we take the highest multiplicities of each prime factor that occur in any of these factorizations and combine them to produce the smallest number that all of these divide into. This can be automated with code that finds the smallest multiplicity of each prime needed:
def is_prime(p):
return p > 1 and all(p % d != 0 for d in range(2,p))
n = 20
primes = [p for p in range(n+1) if is_prime(p)]
result = 1
for p in primes:
multiplicity = 1
while p**(multiplicity+1) <= n:
multiplicity += 1
result *= p**multiplicity
print(f'factor {p}^{multiplicity}')
print(result)
Another way is to compute \(\lcm(1,2,3,\ldots,20)\). This can be done with a number theory theorem: \(\lcm(a,b)=\frac{ab}{\gcd(a,b)}\). This theorem can be extended to more than 2 numbers. We can do this with the code here:
from math import gcd
n = 20
result = 1
for i in range(1,n+1):
result = (result*i) // gcd(result,i)
print(result)
Problem: Find the difference between the sum of squares of \(1,2,3,\ldots,100\) and the square of \(1+2+3+\ldots+100\).
This problem does not require programming, but if you really wanted to do it with code, a simple solution is:
n = 100
square_of_sum = sum(range(1,n+1))**2
sum_of_squares = sum(i**2 for i in range(1,n+1))
print(square_of_sum - sum_of_squares)
Now without code, we write this using the well known summation formulas:
\[\begin{align} \left(\sum_{i=1}^Ni\right)^2-\sum_{i=1}^Ni^2 &=\left(\frac{N(N+1)}{2}\right)^2-\frac{N(N+1)(2N+1)}{6} =\frac{N^2(N+1)^2}{4}-\frac{N(N+1)(2N+1)}{6}\\ &=N(N+1)\left(\frac{N(N+1)}{4}-\frac{2N+1}{6}\right) =N(N+1)\left(\frac{3N^2+3N}{12}-\frac{4N+2}{12}\right)\\ &=N(N+1)\frac{3N^2-N-2}{12}=\frac{N(N+1)(3N+2)(N-1)}{12} \end{align}\]With a bit of algebra, we get a nicely factored expression. Plug in \(N=100\) to get the answer.
Note: this is equivalent to the sum of cubes minus the sum of squares. The sum of cubes is the dominating part which is asymptotic to \(N^4/4\). As expected, the solution expression is also asymptotic to \(N^4/4\).
Problem: Find the 10001st prime.
This problem can be done by checking numbers for primality and counting. This is a slow method but fast enough to solve this problem in a reasonable amount of time. The prime checking function here handles 2 and even numbers separately. Then tests odd divisors up to the square root of the input. The main loop starts at 1 and counts primes until reaching the 10001st one.
from math import sqrt
def is_prime(p):
# handle 2 and even numbers first
if p < 2:
return False
if p == 2:
return True
if p % 2 == 0:
return False
# test possible odd divisors
for d in range(3,1+int(sqrt(p)),2):
if p % d == 0:
return False
return True
ind = 10001
p = 1
i = 0
while i < ind:
p += 1
if is_prime(p):
i += 1
print(p)
The alternative faster method is using a sieve. This uses more memory. The sieve of Eratosthenes works by making a list integers starting at 2:
\[2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20, 21,22,23,24,25,26,27,28,29,30,\ldots\]First, 2 is prime, and we cross off multiples of 2.
\[\red{2},3,\cancel{4},5,\cancel{6},7,\cancel{8},9,\cancel{10}, 11,\cancel{12},13,\cancel{14},15,\cancel{16},17,\cancel{18}, 19,\cancel{20},21,\cancel{22},23,\cancel{24},25,\cancel{26}, 27,\cancel{28},29,\cancel{30},\ldots\]Next 3 is prime and so are all uncrossed numbers below \(3^2\). We cross out multiples of 3 starting at 9.
\[\red{2},\red{3},\cancel{4},\red{5},\cancel{6},\red{7},\cancel{8}, \cancel{9},\cancel{10}, 11,\cancel{12},13,\cancel{14},\cancel{15},\cancel{16},17,\cancel{18}, 19,\cancel{20},\cancel{21},\cancel{22},23,\cancel{24},25,\cancel{26}, \cancel{27},\cancel{28},29,\cancel{30},\ldots\]Now we look at 5 and cross multiples starting at \(5^2\).
\[\red{2},\red{3},\cancel{4},\red{5},\cancel{6},\red{7},\cancel{8}, \cancel{9},\cancel{10},\red{11},\cancel{12},\red{13},\cancel{14},\cancel{15}, \cancel{16},\red{17},\cancel{18},\red{19},\cancel{20},\cancel{21},\cancel{22}, \red{23},\cancel{24},\cancel{25},\cancel{26},\cancel{27}, \cancel{28},29,\cancel{30},\ldots\]We can stop here because the next prime 7 is above \(\lfloor\sqrt{30}\rfloor\). 29 is uncrossed so it is also prime. This procedure allows us to create a list of primes much faster but requires more memory because we must store the sieve.
It is helpful to know how big our sieve must be before we start. Based on a table from Wikipedia, \(x/\ln(x)\) is an underestimate of the number of primes below \(x\) for any \(x\) we would care about in this problem. So if we find \(x\) such that \(x/\ln(x)=n\), then we will have enough primes below \(x\) to find the \(n\)th prime. We can use Newton's method to solve \(x/\ln(x)=n\). I'll avoid the details of Newton's method for now to focus on the sieving method. Here is the code:
from math import log,ceil
def sieve_size(ind):
# solve x/ln(x)=ind with newton's method
assert ind >= 3
f = lambda x : x/log(x) - ind
fp = lambda x : (log(x)-1)/log(x)**2
x = 3.0 # initial guess must be > e
while True:
xold = x
x -= f(xold)/fp(xold)
# terminate when new guess is close enough to previous
if abs(x-xold)/xold < 2**-44:
break
return ceil(x)
ind = 10001
siz = sieve_size(ind)
# sieve 2,3,4,...,siz
sieve = [True]*(siz+1) # initially mark all as prime
sieve[0] = sieve[1] = False
for i in range(2,siz+1):
# loop termination and ignore composite numbers
if i*i > siz:
break
if not sieve[i]:
continue
# cross off multiples of i starting at i*i
for j in range(i*i,siz+1,i):
sieve[j] = False
# count primes until reaching desired index
i = 0
result = 0
for p in range(2,siz+1):
if sieve[p]:
i += 1
if i == ind:
result = p
break
# this should always work in practical cases
# but in general, it is harder to prove
assert i == ind
print(result)
Note: the sieve size calculated to use is 116684, about 11% overestimated.
Problem: Find the greatest product formed from 13 consecutive digits in the following:
73167176531330624919225119674426574742355349194934 96983520312774506326239578318016984801869478851843 85861560789112949495459501737958331952853208805511 12540698747158523863050715693290963295227443043557 66896648950445244523161731856403098711121722383113 62229893423380308135336276614282806444486645238749 30358907296290491560440772390713810515859307960866 70172427121883998797908792274921901699720888093776 65727333001053367881220235421809751254540594752243 52584907711670556013604839586446706324415722155397 53697817977846174064955149290862569321978468622482 83972241375657056057490261407972968652414535100474 82166370484403199890008895243450658541227588666881 16427171479924442928230863465674813919123162824586 17866458359124566529476545682848912883142607690042 24219022671055626321111109370544217506941658960408 07198403850962455444362981230987879927244284909188 84580156166097919133875499200524063689912560717606 05886116467109405077541002256983155200055935729725 71636269561882670428252483600823257530420752963450
The easiest way to solve this would just be trying every product:
digits = '\
73167176531330624919225119674426574742355349194934\
96983520312774506326239578318016984801869478851843\
85861560789112949495459501737958331952853208805511\
12540698747158523863050715693290963295227443043557\
66896648950445244523161731856403098711121722383113\
62229893423380308135336276614282806444486645238749\
30358907296290491560440772390713810515859307960866\
70172427121883998797908792274921901699720888093776\
65727333001053367881220235421809751254540594752243\
52584907711670556013604839586446706324415722155397\
53697817977846174064955149290862569321978468622482\
83972241375657056057490261407972968652414535100474\
82166370484403199890008895243450658541227588666881\
16427171479924442928230863465674813919123162824586\
17866458359124566529476545682848912883142607690042\
24219022671055626321111109370544217506941658960408\
07198403850962455444362981230987879927244284909188\
84580156166097919133875499200524063689912560717606\
05886116467109405077541002256983155200055935729725\
71636269561882670428252483600823257530420752963450'
L = 13
largest = 0
for i in range(len(digits)-L+1):
product = 1
for j in range(i,i+L):
product *= int(digits[j])
if product > largest:
largest = product
print(largest)
It can be done a little faster by reusing work instead of multiplying 13 digits every time we test a product. Ignore the zeroes because a product containing them is just 0. For each block of consecutive nonzero digits, we start with the product of the first 13. Then multiply the 14th and divide the 1st. Then multiply the 15th and divide the 2nd. By repeating this, we have a sliding window of 13 digits and we shift it over 1 digit by doing 1 multiplication and 1 division. This is faster than doing 13 multiplications. Here is the code:
digits = '\
73167176531330624919225119674426574742355349194934\
96983520312774506326239578318016984801869478851843\
85861560789112949495459501737958331952853208805511\
12540698747158523863050715693290963295227443043557\
66896648950445244523161731856403098711121722383113\
62229893423380308135336276614282806444486645238749\
30358907296290491560440772390713810515859307960866\
70172427121883998797908792274921901699720888093776\
65727333001053367881220235421809751254540594752243\
52584907711670556013604839586446706324415722155397\
53697817977846174064955149290862569321978468622482\
83972241375657056057490261407972968652414535100474\
82166370484403199890008895243450658541227588666881\
16427171479924442928230863465674813919123162824586\
17866458359124566529476545682848912883142607690042\
24219022671055626321111109370544217506941658960408\
07198403850962455444362981230987879927244284909188\
84580156166097919133875499200524063689912560717606\
05886116467109405077541002256983155200055935729725\
71636269561882670428252483600823257530420752963450'
L = 13
largest = 0
i = 0
while i < len(digits):
# find the next nonzero digit
while i < len(digits) and digits[i] == '0':
i += 1
# find the end of the span of nonzero digits
j = i+1
while j < len(digits) and digits[j] != '0':
j += 1
if j-i >= L:
# with enough digits, form the initial product
product = 1
for k in range(i,i+L):
product *= int(digits[k])
largest = max(largest,product)
# slide the product over 1 digit at a time
for k in range(i+L,j):
product //= int(digits[k-L])
product *= int(digits[k])
largest = max(largest,product)
i = j # start next iteration after what this iteration found
print(largest)
Problem: find the product \(abc\) of the unique Pythagorean triple \((a,b,c)\) with \(a+b+c=1000\) and \(a<b<c\).
The first way we could do this is choosing \(a,b\) (by brute force), setting \(c=1000-a-b\), and then checking if it is a Pythagorean triple. We can speed this up a little bit with some bounds. Since \(a<b<c\), we have \(1000=a+b+c>3a \Rightarrow a\leq\lfloor1000/3\rfloor\). Similarly, \(1000=a+b+c>a+b>2b \Rightarrow b\leq\lfloor1000/2\rfloor\).
S = 1000
for a in range(1,1+S//3):
for b in range(a+1,1+S//2):
c = S - a - b
if a*a + b*b == c*c:
print(f'triple ({a},{b},{c})')
print(a*b*c)
Another method to solve this is by generating and scaling primitive Pythagorean triples, where \(a,b,c\) only have a commmon factor of 1 (if \((a,b,c)\) is a Pythagorean triple, then so is \((da,db,dc)\)).
We can generate Pythagorean triples using Euclid's formula described on Wikipedia. We choose integers \(m>n>0\). Then a Pythagorean triple is \(a=m^2-n^2,b=2mn,c=m^2+n^2\) because:
\[a^2+b^2=(m^4-2m^2n^2+n^4)+4m^2n^2=m^4+2m^2n^2+n^4=(m^2+n^2)^2=c^2\]Every primitive triple corresponds to a unique choice of \(m,n\) where they are coprime, one is even, and the other is odd. There is a proof for this on the Wikipedia page. For solving the problem, suppose \(a,b,c\) is a primitive triple and we scale it by a factor \(d\). Then:
\[d(a+b+c)=1000=d(m^2-n^2+2mn+m^2+n^2)=2dm(m+n)\]So we can choose \(m,n\) to partially factor 1000 and then choose \(d\). First, \(m\) must be a factor of \(1000/2\). Our bound for \(m\) can be determined like this:
\[1000/2=dm(m+n)\geq m(m+n)>m^2 \Rightarrow m\leq\left\lfloor\sqrt{1000/2}\right\rfloor\]Next, we can implement this idea in code:
from math import sqrt,gcd
S = 1000
assert S % 2 == 0
for m in range(1,1+int(sqrt(S//2))):
if (S//2) % m != 0:
continue # m must be a factor of S/2
n0 = 1 if m % 2 == 0 else 2 # m odd -> n even, m even -> n odd
for n in range(n0,m,2):
if (S//2//m) % (m+n) != 0:
continue # m+n must be a factor of S/(2m)
if gcd(m,n) != 1:
continue # m,n must be coprime for a primitive triple
a = m*m - n*n
b = 2*m*n
c = m*m + n*n
d = S//2//m//(m+n)
print(f'triple {d} * ({a},{b},{c})')
print(a*b*c * d**3)
Problem: Find the sum of all primes below 2000000.
For this problem, I reuse the sieve of Eratosthenes from problem 7, but make an improvement. Since we need all primes below a given limit, there is no need to estimate a sieve size. This time, to reduce the amount of memory by half, only store odd numbers in the sieve. Index \(i\) in the 0-indexed array represents \(2i+1\). It is not much different from the sieve with every number, but we have to calculate indexes a little. The inner loop goes up skipping even numbers. Finally, remember the missing even prime that is not in the sieve.
In general, the sieve of Eratosthenes can be faster and use less memory by excluding numbers divisible by more primes. For example, primesieve is a fast implementation of the sieve of Eratosthenes that skips multiples of 2, 3, 5, and 7.
L = 2000000
sieve = [True]*(L//2) # index i means number 2*i+1
sieve[0] = False
for i in range(1,L//2): # 3 is at index 1
inum = 2*i+1
if inum*inum >= L:
break
if not sieve[i]:
continue
# start at inum*inum and go in increments of
# 2*inum because this sieve skips even numbers
for jnum in range(inum*inum,L,2*inum):
sieve[jnum//2] = False
# remember to include 2 in the sum
print(2 + sum(2*i+1 for i in range(L//2) if sieve[i]))
print(f'number of primes is {1 + sum(1 for i in range(L//2) if sieve[i])}')
Note: a total of 148933 primes are summed.