#!/usr/bin/env python # coding: utf-8 # # [Project Euler](https://ProjectEuler.net) # [This Python 3 notebook](Project%20Euler%20%28Python%203%29.ipynb) contains *some* solutions for the [Project Euler](https://ProjectEuler.net) challenge. # # ### /!\ **Warning:** do not spoil yourself the pleasure of solving these problems by yourself! # # [I (Lilian Besson)](http://perso.crans.org/besson/) started in February 2015, and worked occasionally on Project Euler problems in March and April 2015. # I should try to work on it again, hence this notebook... # # ![Badge giving the number of solved problems](https://ProjectEuler.net/profile/Naereen.png?h "Badge giving the number of solved problems") # --- # ## [Problem 12 : Highly divisible triangular number](https://projecteuler.net/problem=12) # In[1]: get_ipython().run_line_magic('load_ext', 'Cython') # In[2]: from typing import List # In[3]: def divisors(n: int) -> List[int]: return [k for k in range(1, n+1) if n%k == 0] # In[4]: divisors(28) # In[6]: def number_of_divisors(n: int) -> int: return len(divisors(n)) # In[7]: number_of_divisors(28) # In[23]: def int_sqrt(n: int) -> int: if n <= 0: return 0 if n == 1: return 1 i = 0 while i*i <= n: i += 1 if i*i == n: return i return i-1 # In[24]: import math for i in range(10000): assert int_sqrt(i) == math.floor(math.sqrt(i)), f"{i}" # In[25]: def number_of_divisors2(n: int) -> int: c = 0 for k in range(1, n+1): if n%k == 0: c += 1 return c # In[26]: for n in range(1000): assert number_of_divisors(n) == number_of_divisors2(n) # In[34]: def number_of_divisors3(n: int) -> int: if n <= 0: return 0 if n == 1: return 1 c = 0 for k in range(1, int_sqrt(n)+1): # print(f"n = {n}, k = {k}, c = {c}") if n%k == 0: c += 1 if k*k == n else 2 return c # In[36]: get_ipython().run_line_magic('timeit', 'number_of_divisors(1000000)') get_ipython().run_line_magic('timeit', 'number_of_divisors2(1000000)') get_ipython().run_line_magic('timeit', 'number_of_divisors3(1000000)') # In[35]: for n in range(1000): assert number_of_divisors(n) == number_of_divisors3(n), f"{n}" # In[38]: get_ipython().run_cell_magic('cython', '-a', '\ncdef int int_sqrt(int n):\n if n <= 0: return 0\n elif n == 1: return 1\n cdef int i = 0\n while i*i <= n:\n i += 1\n if i*i == n: return i\n return i-1\n\ncdef int number_of_divisors(int n):\n if n <= 0: return 0\n elif n == 1: return 1\n cdef int c = 0\n for k in range(1, int_sqrt(n)+1):\n # print(f"n = {n}, k = {k}, c = {c}")\n if n%k == 0:\n c += 1 if k*k == n else 2\n return c\n\ndef first_highly_divisible_triangular_number(int nb_of_divisors=5):\n cdef int triangular_number = 1\n cdef int i = 2\n while number_of_divisors(triangular_number) < nb_of_divisors:\n triangular_number += i\n i += 1\n return triangular_number') # In[39]: first_highly_divisible_triangular_number(5) # In[40]: first_highly_divisible_triangular_number(150) # In[41]: first_highly_divisible_triangular_number(200) # In[42]: first_highly_divisible_triangular_number(250) # In[43]: first_highly_divisible_triangular_number(300) # In[44]: first_highly_divisible_triangular_number(350) # In[45]: first_highly_divisible_triangular_number(400) # In[46]: first_highly_divisible_triangular_number(450) # In[47]: first_highly_divisible_triangular_number(500) # That's slow if you test for divisors until n, but using this $\lfloor \sqrt{n} \rfloor$ trick speeds it up! # --- # ## [Problem 19: Counting Sundays](https://projecteuler.net/problem=19) # In[3]: from datetime import date, timedelta # In[4]: def isSunday(current_date: date) -> bool: current_iso_date = current_date.isocalendar() weekday = current_iso_date[2] return weekday == 7 # In[5]: start_date = date(year=1901, month=1, day=1) current_date = start_date tomorrow = timedelta(days=+1) max_date = date(year=2000, month=12, day=31) number_of_sundays = 0 while current_date <= max_date: if current_date.day == 1: if isSunday(current_date): number_of_sundays += 1 current_date += tomorrow print(f"Between {start_date} and {max_date}, there were {number_of_sundays} Sundays happening the first day of a month.") # ==> 171 # --- # ## [Problem 26: Reciprocal cycles](https://projecteuler.net/problem=36) TODO # In[1]: import decimal D = decimal.Decimal # Let's try with 1000 digits, and I'll just increase it if the solution is not good. # In[4]: decimal.getcontext().prec = 1000 # Example: # In[5]: D('1') / D('7') # How to detect a cycle and its length: # In[49]: def detect_cycle_length_str_offset(s: str, offset: int=0) -> int: s = s[offset:] length = 1 n = len(s) cycle = s[:length] while length < n: cut_size = n//length if s[:cut_size * length] == cycle * cut_size: return length cycle += s[length] length += 1 return 0 # In[49]: def detect_cycle_length_str(s: str, offset: int=0) -> int: s = s[offset:] length = 1 n = len(s) cycle = s[:length] while length < n: cut_size = n//length if s[:cut_size * length] == cycle * cut_size: return length cycle += s[length] length += 1 return 0 # In[ ]: # In[50]: detect_cycle_length_str('1428571428571428571428571429') # '142857' is a cycle so ==> 6 # In[48]: print(detect_cycle_length_str('1666666666666666666', offset=0)) print(detect_cycle_length_str('1666666666666666666', offset=1)) # In[20]: detect_cycle_length_str('1428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571429') # '142857' is a cycle so ==> 6 # How to detect a cycle in a decimal number of the form $\frac{1}{d}$: # In[9]: def detect_cycle_length(denominator: str, nominator: str='1') -> int: x = D(nominator) / D(denominator) s = str(x).split('.')[1] return detect_cycle_length_str(s) # In[21]: detect_cycle_length('7') # Now it's easy to find a number with maximum cycle length: # In[22]: def number_with_longest_cycle(start: int=2, end: int=1000) -> int: max_cycle_length = 0 candidate_d = None for number in range(start, end): cycle_length = detect_cycle_length(str(number)) if cycle_length > max_cycle_length: candidate_d = number max_cycle_length = cycle_length return candidate_d, max_cycle_length # In[23]: number_with_longest_cycle(start=2, end=11) # ==> 7, 6 as 1/7 has a cycle length of 6 # In[ ]: number_with_longest_cycle(start=1, end=100) # In[ ]: number_with_longest_cycle(start=1, end=1000) # I have not yet found a way to conclude. Working with strings and `decimal.Decimal` numbers is probably a wrong idea. # --- # ## [Problem 32: Pandigital products](https://projecteuler.net/problem=32) # In[223]: set(str(15234)) # In[57]: set(''.join(map(str, range(1, 5+1)))) # In[58]: from typing import Union # In[81]: def is_pandigital(n: Union[str, int], start: int=1, end: int=9): str_n = str(n) if len(str_n) != (end - start + 1): return False return set(str(n)) == set(''.join(map(str, range(start, end+1)))) # In[82]: print(is_pandigital(15234, start=1, end=5)) # True print(is_pandigital(15234, start=1, end=9)) # False # In[106]: set_1to9 = {1,2,3,4,5,6,7,8,9} def is_pandigital_1to9(n: Union[str, int]): str_n = str(n) if len(str_n) != 9: return False return ('0' not in str_n) and len(set(str_n)) == 9 # In[107]: print(is_pandigital_1to9(15234)) # False print(is_pandigital_1to9(152346987)) # False # In[108]: get_ipython().run_line_magic('timeit', 'f"{39}{186}{7254}"') get_ipython().run_line_magic('timeit', 'str(39) + str(186) + str(7254)') # In[109]: def are_pandigital(multiplicand: int, multiplier: int, product: int) -> bool: # assert multiplicand * multiplier == product n = f"{multiplicand}{multiplier}{product}" return is_pandigital_1to9(n) # In[110]: are_pandigital(39, 186, 7254) # In[111]: from typing import Tuple # In[112]: def find_pandigital_products(n: int) -> Tuple[bool, Tuple[int, int]]: x = 2 while x*x <= n: if n % x == 0: y = n // x if are_pandigital(x, y, n): return True, (x, y) x += 1 return False, (None, None) # In[113]: find_pandigital_products(7254) # In[114]: def all_pandigital_products(maxi=10000): return [n for n in range(1, maxi + 1) if find_pandigital_products(n)[0]] # In[115]: def sum_of_all_pandigital_products(maxi=10000): s = 0 for n in range(1, maxi + 1): b, (x, y) = find_pandigital_products(n) if b: print(f"n = {n} is pandigital with x = {x} and y = {y}, and current sum = {s}...") s += n return s # return sum(all_pandigital_products()) # In[116]: sum_of_all_pandigital_products() # --- # ## [Problem 33: Digit cancelling fractions](https://projecteuler.net/problem=33) # In[119]: import fractions F = fractions.Fraction # In[123]: F('49/98'), F('4/8'), F('49/98') == F('4/8') # In[133]: def all_fractions_of_n_digits(n=2): start = '1'*n end = '9'*n for i in range(int(start), int(end) + 1): for j in range(i + 1, int(end) + 1): yield F(f"{i}/{j}"), i, j # In[134]: len(list(all_fractions_of_n_digits(2))) # So there are 3916 fractions of $x/y$ with $x,y$ having exactly two numbers, such that $x/y \leq 1$. # For each of them, we will try to remove one of the digits from $x$ and from $y$, and check if the resulting "reduced" fraction is equal to the original one. If we remove $0$, we won't count it. # It's easy to remove digits when working with strings: # In[137]: si, sj = "49", "98" # In[136]: si.replace("9", "", 1) # So we can check if the fraction $i/j$ can be "wrongly" simplified/reduced by such code: # In[139]: def can_be_wrongly_simplified(i, j): si, sj = str(i), str(j) for digit in si: if digit != '0' and digit in sj: si_without_digit = int(si.replace(digit, "", 1)) sj_without_digit = int(sj.replace(digit, "", 1)) if sj_without_digit != 0 and F(f"{i}/{j}") == F(f"{si_without_digit}/{sj_without_digit}"): return True return False # Let's try it: # In[140]: for frac, i, j in all_fractions_of_n_digits(2): if can_be_wrongly_simplified(i, j): print(f"{frac} with i = {i} and j = {j}") # So there are indeed exactly four such fractions (with exactly two digits on nominator and denominator). # # Let's compute their product: # In[142]: product = F(1, 1) for frac, i, j in all_fractions_of_n_digits(2): if can_be_wrongly_simplified(i, j): print(f"{frac} with i = {i} and j = {j}") product *= frac print(f"Their product = {product}") # --- # ## [Problem 34: Digit factorials](https://projecteuler.net/problem=34) # Let's compute efficiently the factorial function, by using memoisation: # In[150]: from functools import lru_cache @lru_cache(maxsize=None) def factorial(n: int) -> int: if n <= 1: return 1 else: return n * factorial(n - 1) factorials = [factorial(n) for n in range(0, 10)] print(factorials) # In[167]: n = 100 math.ceil(math.log10(n)) # => 2 incorrect, needs + 1 n = 145 math.ceil(math.log10(n)) # => 3 correct # First, we need to bound the size of the largest integer that can be the sum of the factorials of its digits. # # Let $k$ the number of digits of $x$, for instance $k=2$. Then $x$ is between $100...0 = 10^{k-1}$ (eg $10$ for $k=2$) and $999..9 = 10^k - 1$ (eg $99 = 100-1$ for $k=2$). # The sum of the factorials of the digits of $x$ is between $1 * 1! + (k-1) * 0! = 1$ and $k * 9! = 362880 k$. # # And thus, as soon as $362880 k < 10^{k-1}$ then it is not possible for $x$ to be equal to the sum of the factorials of its digits. # # For small values of $k$, such that $k=3$ (for the example $145$), $10^{k-1} \leq 362880 k$ but we can easily find the smallest $k$ that violates this inequality. # # So the smallest $k$ that don't work is the smallest $k$ such that $362880 k < 10^{k-1}$. The right-hand side grows exponentially while the left-hand side grows linearly, so clearly there will be a solution. Intuitively, it won't be large. # Let's find out manually: # In[177]: k = 1 while 10**k <= 362880 * k: k += 1 k = max(1, k - 1) print(k) max_nb_pb34 = 10**k - 1 print(max_nb_pb34) # We can also easily test if a number is the sum of the factorials of its digits: # In[178]: def is_sum_of_factorials_of_its_digits(n: int) -> bool: digits = map(int, list(str(n))) # 145 => [1, 4, 5] sum_of_factorials_of_its_digits = sum(map(factorial, digits)) return sum_of_factorials_of_its_digits == n # In[179]: assert not is_sum_of_factorials_of_its_digits(144) assert is_sum_of_factorials_of_its_digits(145) assert not is_sum_of_factorials_of_its_digits(146) # Now we can list all the numbers up-to a certain bound, and keep the ones we want, then compute their sum: # In[180]: def list_of_such_numbers_pb34(bound=max_nb_pb34): for i in range(3, bound + 1): if is_sum_of_factorials_of_its_digits(i): yield i def sum_of_such_numbers_pb34(bound=max_nb_pb34): return sum(list_of_such_numbers_pb34(bound=bound)) # For numbers up-to $1000$, only 145 works: # In[181]: print(list(list_of_such_numbers_pb34(bound=1000))) print(sum_of_such_numbers_pb34(bound=1000)) # Now let's use the bound we computed: # In[182]: print(list(list_of_such_numbers_pb34(bound=max_nb_pb34))) print(sum_of_such_numbers_pb34(bound=max_nb_pb34)) # --- # ## [Problem 35: Circular primes](https://projecteuler.net/problem=35) # In[184]: import math # In[189]: def erathostene_sieve(n: int) -> List[int]: primes = [False, False] + [True] * (n - 1) # from 0 to n included max_divisor = math.floor(math.sqrt(n)) for divisor in range(2, max_divisor + 1): if primes[divisor]: i = 2 while divisor * i <= n: primes[divisor * i] = False i += 1 return primes # In[191]: sieve100 = erathostene_sieve(100) print(sieve100) print([i for i, b in enumerate(sieve100) if b]) # In[194]: sieve1million = erathostene_sieve(int(1e6)) print(len([i for i, b in enumerate(sieve1million) if b])) # In[196]: def is_prime(n: int) -> bool: assert 1 <= n <= len(sieve1million) - 1 return sieve1million[n] # Now we have all the prime numbers up-to one million, and a constant-time way of checking if small numbers are prime. # In[203]: s = '197' print(s) print(s[0:] + s[:0]) print(s[1:] + s[:1]) print(s[2:] + s[:2]) # In[206]: def circular_numbers(n: int) -> List[int]: str_n = str(n) return [ int(str_n[i:] + str_n[:i]) for i in range(1, len(str_n)) ] # In[207]: circular_numbers(197) # In[208]: def is_circular_prime(n: int) -> bool: if not is_prime(n): return False return all(map(is_prime, circular_numbers(n))) # In[209]: print([n for n in range(2, 100) if is_circular_prime(n)]) # In[210]: print(len([n for n in range(2, int(1e6) + 1) if is_circular_prime(n)])) # --- # ## [Problem 36: Double-base palindromes](https://projecteuler.net/problem=36) # It's very easy to check if $n$ is a palindrome in base $10$: # In[211]: def is_palindrome_base10(n: int) -> bool: str_n = str(n) return str_n == str_n[::-1] # In[212]: assert not is_palindrome_base10(584) assert is_palindrome_base10(585) assert not is_palindrome_base10(586) # It's also very easy to check if $n$ is a palindrome in base $2$: # In[215]: bin(585)[2:] # In[216]: def is_palindrome_base2(n: int) -> bool: str_n_2 = bin(n)[2:] return str_n_2 == str_n_2[::-1] # In[217]: assert not is_palindrome_base2(584) assert is_palindrome_base2(585) assert not is_palindrome_base2(586) # In[218]: def palindromes_in_base10_and_base2(start: int=1, end: int=int(1e6)) -> List[int]: return [n for n in range(start, end + 1) if is_palindrome_base10(n) and is_palindrome_base2(n)] # In[220]: palindromes_in_base10_and_base2(start=1, end=1000) # In[222]: sum(palindromes_in_base10_and_base2(start=1, end=int(1e6))) # --- # ## [Problem 37: Truncatable primes](https://projecteuler.net/problem=37) # In[224]: def truncate_from_left_to_right(n: int) -> List[int]: str_n = str(n) return [int(str_n[i:]) for i in range(0, len(str_n))] # In[225]: truncate_from_left_to_right(3797) # In[240]: def truncate_from_right_to_left(n: int) -> List[int]: str_n = str(n) return [int(str_n[:i]) for i in range(len(str_n), 0, -1)] # In[241]: truncate_from_right_to_left(3797) # In[244]: def is_truncatable_prime(n: int, sieve: List[int]) -> bool: if n <= 7 or not sieve[n]: return False return all(sieve[k] for k in truncate_from_left_to_right(n) + truncate_from_right_to_left(n)) # In[248]: assert is_truncatable_prime(37, sieve=sieve1million) assert is_truncatable_prime(3797, sieve=sieve1million) assert not is_truncatable_prime(99371, sieve=sieve1million) # 99 is not prime! # In[250]: print([i for i, b in enumerate(sieve100) if b]) # The tricky part is to bound the largest number than can be a truncatable prime. # The problem states that there are only eleven primes that are truncatable primes. # How to find an upper bound on these numbers? # # Well, one solution could just to enumerate prime numbers until we find 11 candidates, then stop. # Let's cheat and do that. # In[253]: def find_truncatable_primes(maxnumber=11, sieve=sieve1million): truncatable_primes = [] number = 0 for i, b in enumerate(sieve1million): if b: if is_truncatable_prime(i, sieve=sieve): truncatable_primes.append(i) number += 1 if number >= maxnumber: return truncatable_primes # In[264]: find_truncatable_primes(maxnumber=11) # In[265]: find_truncatable_primes(maxnumber=12) # That's a good sign, for prime numbers up-to a million we could not find the 12th truncatable prime. # In[266]: sum(find_truncatable_primes(maxnumber=11)) # --- # ## [Problem 38: Pandigital multiples](https://projecteuler.net/problem=38) # # Well I wrote the code and my notebook froze so I lost everything... Not so motivated to start again. # # The first thing to do is to mathematically bound the largest integer $x$ such that the product of $x$ by $1,\dots,n$ gives a 1-to-9 pandigital number. # # - $n$ clearly has to be smaller than $9$, otherwise the concatenated product have more than $10$ digits and cannot be 1-to-9 pandigital. # - Now, for a fixed $n$, if $x$ has $d$ digits then $nd \leq 9$, so this bounds the largest possible $x$ for a candidate $n$. # # Let's simply try them all. # In[5]: from typing import Union # In[6]: def is_pandigital_1to9(n: Union[str, int]) -> bool: str_n = str(n) if len(str_n) != 9: return False return sorted(str_n) == ['1','2','3','4','5','6','7','8','9'] # In[4]: assert is_pandigital_1to9(192384576) assert is_pandigital_1to9(918273645) # In[9]: def concatenated_product_1ton(x: int, n: int) -> str: return ''.join(str(x * i) for i in range(1, n+1)) # In[13]: assert is_pandigital_1to9(concatenated_product_1ton(192, 3)) assert is_pandigital_1to9(concatenated_product_1ton(9, 5)) assert is_pandigital_1to9(concatenated_product_1ton(1, 9)) # Now we can just try all the possible values for $x$ and $n$ and keep a list of the candidates and return the largest one. # In[26]: def all_concatenated_product_pandigital_1to9(max_n: int): candidates = [] for n in range(2, max_n + 1): print(f"\nStarting with n = {n}...") d = 1 while n*d <= 9: print(f"\tStarting with d = {d}...") for x in range(10**(d-1), 10**d): prod_x_1ton = concatenated_product_1ton(x, n) if is_pandigital_1to9(prod_x_1ton): print(f" ==> adding {prod_x_1ton} coming from x = {x}, d = {d} and n = {n}") candidates.append(int(prod_x_1ton)) d += 1 return candidates # In[27]: def largest_concatenated_product_pandigital_1to9(max_n: int): candidates = all_concatenated_product_pandigital_1to9(max_n = max_n) return max(candidates) # In[28]: largest_concatenated_product_pandigital_1to9(9) # --- # ## [Problem 39: Integer right triangles](https://projecteuler.net/problem=39) # If p is the perimeter of a right angle triangle with integral length sides, {a,b,c}, there are exactly three solutions for p = 120 : {20,48,52}, {24,45,51}, {30,40,50} # # For which value of p ≤ 1000, is the number of solutions maximised? # # Well there is no need to be smart here, we just need to try all the values $a,b,c=\sqrt{a^2+b^2}$ such that the longer slide is an integer, keep the list of values of p and the number of solutions, and find the p that maximize this number of solutions. # # We could be much smarter in narrowing down the possible values for $a,b,c$, but we don't really need to. # In[29]: import math # In[52]: def find_all_integer_right_triangle(max_p=1000): nb_of_solutions_by_p = {} solutions_by_p = {} for p in range(1, max_p + 1): nb_of_solutions_by_p[p] = 0 solutions_by_p[p] = [] max_side = p // 2 for a in range(1, max_side + 1): for b in range(a + 1, max_side + 1): c = math.sqrt(a**2 + b**2) if int(c) == c: c = int(c) if a + b + c == p: nb_of_solutions_by_p[p] += 1 solutions_by_p[p] += [(a,b,c)] max_nb_of_solutions = max(nb_solutions for nb_solutions in nb_of_solutions_by_p.values()) print(max_nb_of_solutions) print([solutions_by_p[p] for p in solutions_by_p.keys() if nb_of_solutions_by_p[p] == max_nb_of_solutions]) return [p for p in nb_of_solutions_by_p.keys() if nb_of_solutions_by_p[p] == max_nb_of_solutions] # In[53]: find_all_integer_right_triangle(max_p=120) # In[54]: find_all_integer_right_triangle(max_p=1000) # So the answer is $p=840$, which have 8 solutions. # --- # ## [Problem 40: Champernowne's constant](https://projecteuler.net/problem=40) # # Let's use a huge string and simply access its value at index $i=1,10,100,1000,10000,100000,1000000$. # It will take some time but it'll work easily. # In[55]: def construct_irrational_decimal_fraction(max_n=30): return ''.join(str(n) for n in range(1, max_n + 1)) # In[58]: assert '1' == construct_irrational_decimal_fraction(30)[11] # 12th digit is 1, ok # Now we can just build a huge string and access it. # In[59]: def product(iterator): p = 1 for value in iterator: p *= value return p # In[67]: def product_of_d1_to_dn(list_of_n=[1,10]): max_n = max(list_of_n) s = construct_irrational_decimal_fraction(max_n=1+max_n) return product(int(s[di]) for di in list_of_n) # In[68]: product_of_d1_to_dn([1, 11]) # product of 2th digit = 2 and 12th digit = 1 ==> 2 # In[70]: product_of_d1_to_dn([0, 9]) # In[71]: product_of_d1_to_dn([0, 9, 99]) # In[72]: product_of_d1_to_dn([0, 9, 99, 999]) # In[73]: product_of_d1_to_dn([0, 9, 99, 999, 9999]) # In[74]: product_of_d1_to_dn([0, 9, 99, 999, 9999, 99999]) # In[75]: product_of_d1_to_dn([0, 9, 99, 999, 9999, 99999, 999999]) # --- # ## [Problem 41: Pandigital prime](https://projecteuler.net/problem=41) # # If $p$ is a prime $n$-digit pandigital number, then $p$ is smaller than $n(n-1)..1$ (written in base $10$), eg with $n=5$, $p$ is smaller than $54321$. # For $n=9$, the largest possible $9$-digit pandigital number is $987654321$. # Let's build a Erathostene sieve of size $987654321+1$ and just check all the prime numbers for being $n$-digit pandigital, for $n=1$ to $n=9$. # It won't be fast, but it'll work. # In[1]: import math from typing import List # In[2]: get_ipython().run_line_magic('load_ext', 'Cython') # In[3]: get_ipython().run_cell_magic('cython', '', 'import math\n\ndef erathostene_sieve(int n):\n cdef list primes = [False, False] + [True] * (n - 1) #\xa0from 0 to n included\n cdef int max_divisor = math.floor(math.sqrt(n))\n cdef int i = 2\n for divisor in range(2, max_divisor + 1):\n if primes[divisor]:\n number = 2*divisor\n while number <= n:\n primes[number] = False\n number += divisor\n return primes') # For 10 million primes, this takes about 3 seconds on my machine, thanks to Cython being faster than naive Python code. # In[9]: sieve10million = erathostene_sieve(int(1e7)) primes_upto_1e7 = [p for p,b in enumerate(sieve10million) if b] # For about 1 billion primes, this takes about 30 seconds on my machine, thanks to Cython being faster than naive Python code. # In[21]: sieve1billion = erathostene_sieve(98765432+1) # For about 10 billion primes, this should take about 15 minutes on my machine, thanks to Cython being faster than naive Python code. # In[22]: sieve10billion = erathostene_sieve(987654321+1) # Now we have a LARGE bool array to check in constant time if a number is prime, we can build the list of prime numbers up-to this $max_p = 987654321$. # In[23]: primes_upto_987654321 = [p for p,b in enumerate(sieve10billion) if b] # In[24]: del sieve1billion, sieve10billion # to save so RAM # In[25]: len(primes_upto_987654321) # not so large, there is about O(log n) primes smaller than n # Let's test if a number is 1-to-n pandigital. # In[26]: from typing import Union def is_pandigital_1ton(x: Union[str, int], n: int) -> bool: str_x = str(x) if len(str_x) != n: return False # fastest version if n == 1: return sorted(str_x) == ['1'] elif n == 2: return sorted(str_x) == ['1','2'] elif n == 3: return sorted(str_x) == ['1','2','3'] elif n == 4: return sorted(str_x) == ['1','2','3','4'] elif n == 5: return sorted(str_x) == ['1','2','3','4','5'] elif n == 6: return sorted(str_x) == ['1','2','3','4','5','6'] elif n == 7: return sorted(str_x) == ['1','2','3','4','5','6','7'] elif n == 8: return sorted(str_x) == ['1','2','3','4','5','6','7','8'] elif n == 9: return sorted(str_x) == ['1','2','3','4','5','6','7','8','9'] return False # In[27]: from typing import Dict def find_pandigital_primes(list_of_primes: List[int]) -> Dict[int, List[int]]: solutions = dict() for n in range(1, 9 + 1): solutions[n] = [ prime for prime in list_of_primes if is_pandigital_1ton(prime, n) ] return solutions # In[28]: def largest_pandigital_prime(list_of_primes: List[int]) -> int: solutions = find_pandigital_primes(list_of_primes) # print("All 1-to-n pandigital primes:", {n: solutions[n] for n in range(1, 9 + 1)}) max_primes = { n: max(solutions[n]) for n in range(1, 9 + 1) if solutions[n] } print("Max 1-to-n pandigital primes:", max_primes) return max(max_primes.values()) # In[29]: largest_pandigital_prime(primes_upto_1e7) # In[30]: largest_pandigital_prime(primes_upto_987654321) # --- # ## [Problem 32: Coded triangle numbers](https://projecteuler.net/problem=42) # # - We need to be able to check (quickly) if a word is a triangle word. # - The simplest approach would be to first read all the words, compute the length $n_\max$ of the longer one, and then say that the longest word value can be $t_{up} = 26* n_\max$, so then we can compute all the triangle numbers up to $t_{up}$. # In[33]: get_ipython().system('ls p042_words.txt') get_ipython().system('wc p042_words.txt') # In[41]: with open('p042_words.txt') as f: words = f.readline() words = words.replace('"','').split(',') # In[43]: type(words), words[0], len(words) # Now we have these 1786 words. # In[45]: max_length_words = max(map(len, words)) print(f"There is {len(words)} words, of maximum length {max_length_words}") # In[46]: largest_word_value = 26 * max_length_words print(f"So the largest word is {largest_word_value}") # In[48]: def triangle_numbers_upto(max_value: int) -> List[int]: triangle_numbers = [] i = 1 t = i while t <= max_value: triangle_numbers.append(t) i += 1 t += i return triangle_numbers # In[51]: triangles = set(triangle_numbers_upto(largest_word_value)) # using a set() will make the test n in triangles faster! print(triangles) print(f"There are {len(triangles)} triangle numbers smaller than {largest_word_value}") # In[56]: def letter_to_int(letter: str) -> int: return ord(letter) - ord('A') + 1 # for letter in 'ABCDEFGHIJKLMNOPQRSTUVWXYZ': # print(f"{letter} -> {letter_to_int(letter)}") # In[57]: def word_value(word: str) -> int: return sum(letter_to_int(letter) for letter in word) # In[58]: word_value("SKY") # 19 + 11 + 25 = 55 # In[52]: def is_triangle_word(word: str) -> bool: value = word_value(word) return value in triangles # In[62]: def number_of_triangle_words(words: List[str]) -> int: return sum(1 for word in words if is_triangle_word(word)) # In[63]: number_of_triangle_words(words) # Good, it wasn't that hard! # --- # ## [Problem 43: sub-string divisibility](https://projecteuler.net/problem=43) # In[64]: from typing import Union def is_pandigital_0to9(x: Union[str, int]) -> bool: str_x = str(x) if len(str_x) != 10: return False return sorted(str_x) == ['0','1','2','3','4','5','6','7','8','9'] # In[65]: is_pandigital_0to9(1406357289) # In[67]: def has_substring_divisibility(x: Union[str, int]) -> bool: x = str(x) digits = {i+1: int(c) for i,c in enumerate(x)} make_digit = lambda i: int(f"{digits[i]}{digits[i+1]}{digits[i+2]}") divisors = [2, 3, 5, 7, 11, 13, 17] return all( make_digit(i+2) % divisor == 0 for i, divisor in enumerate(divisors) ) # In[69]: assert not has_substring_divisibility(1406357287) assert not has_substring_divisibility(1406357288) assert has_substring_divisibility(1406357289) assert not has_substring_divisibility(1406357290) assert not has_substring_divisibility(1406357291) # Now we simply have to try all the 0-to-9 pandigital numbers, from 0123456789 (smallest) to 9876543210 (largest), and count the number of such numbers that satisfy this property. # # It's easy to test all such pandigital numbers by using permutations of the string 0123456789, using [`itertools.permutations`](https://docs.python.org/3/library/itertools.html#itertools.permutations). # In[70]: import itertools # In[78]: len(list(itertools.permutations(list('0123456789')))) # In[79]: def find_all_substring_divisible_0to9_pandigital_numbers(): candidates = [] for digits in itertools.permutations(list('0123456789')): str_number = ''.join(digits) if is_pandigital_0to9(str_number): if has_substring_divisibility(str_number): print(f"Adding a new number = {int(str_number)}") candidates.append(int(str_number)) return candidates return [ int(str_number) for str_number in itertools.permutations('0123456789') if is_pandigital_0to9(str_number) and has_substring_divisibility(str_number) ] # In[80]: candidates = find_all_substring_divisible_0to9_pandigital_numbers() print(f"There is {len(candidates)} such numbers, their sum = {sum(candidates)}.") # --- # ## [Problem 44: pentagon numbers](https://projecteuler.net/problem=44) # In[92]: def nth_pentagonal_number(n: int) -> int: # if n is odd, 3*n-1 is even, so n * (3n-1) is even # if n is even, n * (3n-1) is even # so the division by two is always an integer division return n * (3*n - 1) // 2 # From the formula of $P_n = n(3n-1)/2$, it is clear that for any fixed $i$, $|P_{i+k}-P_i| = (i+k)(3(i+k)-1)/2 - i(3i-1)/2 = ( i(3i-1) + k(3i-1) + i(3k-1) + k(3k-1) - i(3i-1) )/2 = (k(3i-1) + i(3k-1) + k(3k-1))/2$ is an increasing function of $k$, and also for any fixed $k$, an increasing function of $i$. # Let's try naively to compute the first 1 million pentagonal numbers, put them in a set $\mathcal{P}$ then check if there are any $P_j,P_k$ (with $k > j$) that satisfy $P_j+P_k \in \mathcal{P}$ and $P_k - P_j \in \mathcal{P}$. # # My intuition is that the solution that minimizes $D=|P_k-P_j|$ will be found within the first 1 million pentagonal numbers, so we don't need to mathematically bound the value of the largest $k,j$. # In[93]: def compute_all_pentagonal_numbers(max_n: int) -> List[int]: return [nth_pentagonal_number(n) for n in range(1, int(max_n) + 1)] # In[98]: get_ipython().run_cell_magic('time', '', 'onemillion_pentagonal_numbers = compute_all_pentagonal_numbers(1e5)') # In[99]: def try_all_pairs_of_pentagonal_numbers(pentagonal_numbers: List[int]) -> List[int]: pairs_that_are_good = [] pentagonal_numbers_set = set(pentagonal_numbers) n = len(pentagonal_numbers) for j in range(n): Pj = pentagonal_numbers[j] for k in range(j + 1, n): Pk = pentagonal_numbers[k] their_sum = Pk + Pj if their_sum not in pentagonal_numbers_set: continue their_difference = Pk - Pj if their_difference not in pentagonal_numbers_set: continue print(f"Found another pair P{j} = {Pj}, P{k} = {Pk} whose sum = {their_sum} and difference = {their_difference} are pentagonal numbers too.") pairs_that_are_good.append((their_difference, Pj, Pk, their_sum)) return pairs_that_are_good # In[100]: get_ipython().run_cell_magic('time', '', 'try_all_pairs_of_pentagonal_numbers(onemillion_pentagonal_numbers)') # I got bored of waiting, and tadam the first solution is the good one! # --- # ## [Problem 45: Triangular, pentagonal, and hexagonal](https://projecteuler.net/problem=45) # # Ideally, this problem will be easy to solve if there is a $\mathcal{O}(1)$ method to check if a number is a triangle, pentagonal and hexagonal number. # I don't see an easy method, so let's just compute a lot of each of them, and try. # In[108]: def nth_triangle_number(n: int) -> int: return n * (n+1) // 2 def nth_pentagonal_number(n: int) -> int: return n * (3*n - 1) // 2 def nth_hexagonal_number(n: int) -> int: return n * (2*n - 1) # In[109]: assert nth_triangle_number(285) == 40755 assert nth_pentagonal_number(165) == 40755 assert nth_hexagonal_number(143) == 40755 # Now we can generate a set of a lot of triangle numbers, pentagonal numbers, hexagonal numbers, and find their intersection. # In[110]: {1, 3, 6, 10, 15} & {1, 5, 12, 22, 35} & {1, 6, 15, 28, 45} # In[111]: def intersection_of_sets_of_three_kinds_of_numbers(max_n: int) -> List[int]: set_triangles = { nth_triangle_number(n) for n in range(1, max_n + 1) } set_pentagonals = { nth_pentagonal_number(n) for n in range(1, max_n + 1) } set_hexagonals = { nth_hexagonal_number(n) for n in range(1, max_n + 1) } intersection = set_triangles & set_pentagonals & set_hexagonals return sorted(list(intersection)) # First, let's check that 40755=$T_{285}=P_{165}=H_{143}$ can indeed be found using this approach: # In[112]: intersection_of_sets_of_three_kinds_of_numbers(300) # How can we know when to stop? # I simply tried augmenting `max_n` until I found the next triangle number that is also pentagonal and hexagonal. # In[119]: intersection_of_sets_of_three_kinds_of_numbers(100000) # --- # ## [Problem 46: Goldbach's other conjecture](https://projecteuler.net/problem=46) # In[120]: get_ipython().run_line_magic('load_ext', 'Cython') # In[121]: get_ipython().run_cell_magic('cython', '', 'import math\n\ndef erathostene_sieve(int n):\n cdef list primes = [False, False] + [True] * (n - 1) #\xa0from 0 to n included\n cdef int max_divisor = math.floor(math.sqrt(n))\n cdef int i = 2\n for divisor in range(2, max_divisor + 1):\n if primes[divisor]:\n number = 2*divisor\n while number <= n:\n primes[number] = False\n number += divisor\n return primes') # For 10 million primes, this takes about 3 seconds on my machine, thanks to Cython being faster than naive Python code. # In[122]: sieve10million = erathostene_sieve(int(1e7)) primes_upto_1e7 = [p for p,b in enumerate(sieve10million) if b] # In[126]: def is_integer_square(n: int) -> bool: sqrt_n = math.sqrt(n) return int(sqrt_n)**2 == n # In[127]: assert not is_integer_square(24) assert is_integer_square(25) assert not is_integer_square(26) # In[144]: def is_goldbach(n: int, primes: List[int], verbose: bool=True) -> bool: # n has to be odd to be Goldbach if n % 2 == 0: return False # n has to not be prime if n in primes: return False # now we can test possible candidate primes for p in primes: if p >= n: return False # if n = p + 2*(a square) remainder = n - p if remainder % 2 != 0: continue candidate_square = remainder // 2 if is_integer_square(candidate_square): i = int(math.sqrt(candidate_square)) if verbose: print(f"n = {n} was found to satisfy Goldbach's property with prime p = {p} and i = {i} (n = p + 2 i²)") return True return False # In[145]: count = 0 max_n = 50 for n in range(1, max_n + 1): if is_goldbach(n, primes_upto_1e7, verbose=True): count += 1 print(f"For numbers n = 1 ... {max_n} we found {count} numbers satisfying Goldbach's conjecture.") # In[146]: def find_smallest_odd_composite_non_Goldbach_number(primes: List[int]) -> int: n = 7 max_n = max(primes) while n <= max_n: n += 2 # let's try all the odd numbers from 9 if n in primes: continue # we skip the primes if not is_goldbach(n, primes, verbose=False): return n raise ValueError(f"Not enough prime numbers in the list (max = {max_n}), no counter example was found...") # In[147]: find_smallest_odd_composite_non_Goldbach_number(primes_upto_1e7) # It was that fast, about 24 seconds, but it's enough. # There is surly plenty of solutions to solve this problem faster, but I don't care. # --- # ## [Problem 47: distinct primes factors](https://projecteuler.net/problem=47) # In[183]: def how_many_times_x_divides_y(x: int, y: int) -> int: c = 0 while x % y == 0: x = x // y c += 1 return c # In[187]: assert 2 == how_many_times_x_divides_y(100, 10) assert 2 == how_many_times_x_divides_y(644, 2) assert 1 == how_many_times_x_divides_y(644, 7) assert 1 == how_many_times_x_divides_y(644, 23) assert 0 == how_many_times_x_divides_y(644, 39) # In[ ]: from functools import lru_cache from typing import Tuple @lru_cache(maxsize=None) def number_of_distinct_prime_factors(n: int, primes: Tuple[int]) -> int: return len([p for p in primes if p <= n and n % p == 0]) # In[192]: def distinct_prime_factors(n: int, primes: Tuple[int]) -> List[int]: prime_factors = [p for p in primes if p <= n and n % p == 0] return prime_factors def str_prime_decomposition(n: int, primes: Tuple[int]) -> str: prime_factors = distinct_prime_factors(n, primes) prime_decomposition = { p: how_many_times_x_divides_y(n, p) for p in prime_factors } return ' * '.join( f"{p}**{power}" if power >= 2 else f"{p}" for (p, power) in prime_decomposition.items() ) # In[168]: primes_upto_1e7 = tuple(primes_upto_1e7) # In[169]: assert 2 == number_of_distinct_prime_factors(14, primes_upto_1e7) assert 2 == number_of_distinct_prime_factors(15, primes_upto_1e7) # In[195]: for n in [14, 15]: print(f"{n} = {str_prime_decomposition(n, primes_upto_1e7)}") # In[170]: assert 3 == number_of_distinct_prime_factors(644, primes_upto_1e7) assert 3 == number_of_distinct_prime_factors(645, primes_upto_1e7) assert 3 == number_of_distinct_prime_factors(646, primes_upto_1e7) # In[194]: for n in [644, 645, 646]: print(f"{n} = {str_prime_decomposition(n, primes_upto_1e7)}") # Now let's find the first number to have $x$ consecutive integers each having $y$ prime factors: # In[198]: def find_first_number_to_have_x_consecutive_int_with_y_prime_factors(max_n: int, x: int, y: int, primes: List[int]) -> int: n = 1 primes = tuple([p for p in primes if y*p <= max_n]) while n <= max_n: can_still_be_good = True number = n while number < n+x and can_still_be_good: if number_of_distinct_prime_factors(number, primes) != y: can_still_be_good = False # we can skip a few values for n n = number number += 1 if can_still_be_good: print(f"Proof that n = {n} has x = {x} consecutive integers with each y = {y} distinct prime factors:") for i in range(n, n + x): print(f" i = {i} has prime factors {distinct_prime_factors(i, primes)}") print(f" Proof: i = {i} == {str_prime_decomposition(i, primes)}") return n n += 1 raise ValueError(f"Not enough prime numbers in the list (max = {max_n}), no counter example was found...") # Now let's try with $x=y=2$ and $x=y=3$, to check the example given in the problem statement, before trying with $x=y=4$: # In[199]: find_first_number_to_have_x_consecutive_int_with_y_prime_factors(100, 2, 2, primes_upto_1e7) # In[200]: find_first_number_to_have_x_consecutive_int_with_y_prime_factors(1000, 3, 3, primes_upto_1e7) # Now to answer the problem, let's compute for $x=y=4$: # In[201]: get_ipython().run_cell_magic('time', '', 'find_first_number_to_have_x_consecutive_int_with_y_prime_factors(150000, 4, 4, primes_upto_1e7)') # --- # ## [Problem 49: Prime permutations](https://projecteuler.net/problem=49) # The arithmetic sequence, 1487, 4817, 8147, in which each of the terms increases by 3330, is unusual in two ways: (i) each of the three terms are prime, and, (ii) each of the 4-digit numbers are permutations of one another. # # There are no arithmetic sequences made up of three 1-, 2-, or 3-digit primes, exhibiting this property, but there is one other 4-digit increasing sequence. # # *What 12-digit number do you form by concatenating the three terms in this sequence?* # # I won't try to solve this efficiently: let's build a set of all the primes smaller than 9999 (to check in efficiently if a number is prime), then try all 4-digit prime numbers $p$ (ie. larger than 1000 and smaller than 9999) and keep the ones that satisfy this property. # In[202]: get_ipython().run_line_magic('load_ext', 'Cython') # In[203]: get_ipython().run_cell_magic('cython', '', 'import math\n\ndef erathostene_sieve(int n):\n cdef list primes = [False, False] + [True] * (n - 1) #\xa0from 0 to n included\n cdef int max_divisor = math.floor(math.sqrt(n))\n cdef int i = 2\n for divisor in range(2, max_divisor + 1):\n if primes[divisor]:\n number = 2*divisor\n while number <= n:\n primes[number] = False\n number += divisor\n return primes') # In[206]: sieve10000 = erathostene_sieve(10000+1) primes_upto_10000 = [p for p,b in enumerate(sieve10000) if b] print(f"There are {len(primes_upto_10000)} prime numbers smaller than 10000") # In[209]: def are_permutations(x: int, y: int, z: int) -> bool: sorted_str_x = sorted(list(str(x))) sorted_str_y = sorted(list(str(y))) if sorted_str_x != sorted_str_y: return False sorted_str_z = sorted(list(str(z))) if sorted_str_x != sorted_str_z: return False return True # In[211]: assert not are_permutations(1481, 4817, 8147) assert are_permutations(1487, 4817, 8147) assert not are_permutations(1487, 4317, 8147) # In[219]: from typing import List, Tuple def try_all_prime_numbers_from_1000_to_9999(primes: List[int]=primes_upto_10000) -> List[Tuple[int, int, int]]: primes = [p for p in primes if 1000 <= p <= 9999] primes_set = set(primes) solutions = [] for prime in primes: prime_1 = prime for delta in range(1, (9999-prime)//2): prime_2 = prime + delta if prime_2 not in primes_set: continue # next delta prime_3 = prime + 2*delta if prime_3 not in primes_set: continue # next delta if are_permutations(prime_1, prime_2, prime_3): solutions.append((prime_1, prime_2, prime_3)) if len(solutions) >= 2: return solutions return solutions # In[220]: get_ipython().run_cell_magic('time', '', 'try_all_prime_numbers_from_1000_to_9999()') # So the solution is the concatenation of 2969, 6299, 9629: 296962999629. # --- # ## [Problem 50: Consecutive prime sum](https://projecteuler.net/problem=50) # The prime 41, can be written as the sum of six consecutive primes: $41 = 2 + 3 + 5 + 7 + 11 + 13$ # # This is the longest sum of consecutive primes that adds to a prime below one-hundred. # # The longest sum of consecutive primes below one-thousand that adds to a prime, contains 21 terms, and is equal to 953. # # *Which prime, below one-million, can be written as the sum of the most consecutive primes?* # # I don't see an efficient way of checking this property. # In[202]: get_ipython().run_line_magic('load_ext', 'Cython') # In[203]: get_ipython().run_cell_magic('cython', '', 'import math\n\ndef erathostene_sieve(int n):\n cdef list primes = [False, False] + [True] * (n - 1) #\xa0from 0 to n included\n cdef int max_divisor = math.floor(math.sqrt(n))\n cdef int i = 2\n for divisor in range(2, max_divisor + 1):\n if primes[divisor]:\n number = 2*divisor\n while number <= n:\n primes[number] = False\n number += divisor\n return primes') # In[266]: sieve10million = erathostene_sieve(int(1e7)) primes_upto_10million = [p for p,b in enumerate(sieve10million) if b] print(f"There are {len(primes_upto_10million)} prime numbers smaller than 1 million") # In[267]: def sum_of_k_consecutive_primes(k: int, primes: List[int] = primes_upto_10million) -> List[int]: set_of_primes = set(primes) dict_of_sum_of_k_consecutive_primes = { # sum of k consecutive primes pi + pi+1 + pi+2 + ... + pi+k starting at pi sum(primes[j] for j in range(i, i + k)): tuple(primes[j] for j in range(i, i + k)) for i in range(len(primes) - k + 1) if sum(primes[j] for j in range(i, i + k)) in set_of_primes } set_of_sum_of_k_consecutive_primes = set(dict_of_sum_of_k_consecutive_primes.keys()) return set_of_primes & set_of_sum_of_k_consecutive_primes, dict_of_sum_of_k_consecutive_primes # In[268]: def find_largest_possible_k(primes: List[int] = primes_upto_10million) -> int: k = 2 max_prime = max(primes) while sum(primes[:k]) <= max_prime: k += 1 return k # In[274]: k_max = find_largest_possible_k(primes_upto_1million) print(f"Maximum k = {k_max}") # In[275]: def find_largest_k_that_gives_at_least_one_sum_of_k_consecutive_primes(kmax: int, kmin: int, primes: List[int] = primes_upto_10million) -> int: k = kmax while k >= kmin: print(f"Trying for {k}...") at_least_one_sum_of_k_consecutive_primes, dict_answer = sum_of_k_consecutive_primes(k, primes=primes) if at_least_one_sum_of_k_consecutive_primes: print(f"For k = {k}, intersection = {at_least_one_sum_of_k_consecutive_primes}...") print(f"For these primes, here are their sum: {dict_answer}") print(f"The largest of such prime is {max(at_least_one_sum_of_k_consecutive_primes)}") return k k -= 1 raise ValueError(f"Couldn't find a k in [{kmin},{kmax}] such that there is at least one prime that is the sum of k consecutive primes.") # In[271]: get_ipython().run_cell_magic('time', '', 'find_largest_k_that_gives_at_least_one_sum_of_k_consecutive_primes(6, 2, [p for p in primes_upto_10million if p <= 100])') # In[272]: get_ipython().run_cell_magic('time', '', 'find_largest_k_that_gives_at_least_one_sum_of_k_consecutive_primes(21, 2, [p for p in primes_upto_1million if p <= 1000])') # In[259]: k_max # In[276]: get_ipython().run_cell_magic('time', '', 'find_largest_k_that_gives_at_least_one_sum_of_k_consecutive_primes(k_max, 21, primes_upto_1million)') # Done. # I'm impressed to see that the naive majoration of $k$ (that gave $k_\max=547$) was very tight, as the answer was found to be $543$. # The approach I used to compute this majoration of $k$ was the most naive one: $k$ is the smallest value such that $p_1+p_2+...+p_k = 2 + 3 + 5 + ... + p_k$ the sum of the first $k$ prime numbers is larger than $1000000$ (one million). # --- # ## Continue # In[ ]: