In [7]:

s = ''
for n in range(81, 100):
    if n < 10:
        n = '0' + str(n)
    s += '''[**{n}**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p{n}), '''.format(n=n)
s

Out[7]:

'[**81**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p81), [**82**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p82), [**83**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p83), [**84**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p84), [**85**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p85), [**86**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p86), [**87**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p87), [**88**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p88), [**89**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p89), [**90**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p90), [**91**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p91), [**92**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p92), [**93**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p93), [**94**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p94), [**95**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p95), [**96**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p96), [**97**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p97), [**98**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p98), [**99**](http://nbviewer.jupyter.org/github/ymer/euler-python/blob/master/Python.ipynb#p99), '

General¶

In [1]:

%load_ext autotime
from itertools import takewhile, count, combinations_with_replacement, combinations, product, permutations
from operator import itemgetter, mul
from collections import defaultdict, Counter
from math import sqrt, factorial, log, floor
from functools import reduce, lru_cache

#prod([1, 2, 3]) = 6
def prod(a):
    return reduce(mul, a)

Problem 1 ¶

Add all the natural numbers below one thousand that are multiples of 3 or 5.

In [2]:

sum([n for n in range(1000) if n % 3 == 0 or n % 5 == 0])

Out[2]:

time: 30 ms

Problem 2 ¶

Find the sum of all the even-valued terms in the Fibonacci sequence which do not exceed four million.

In [3]:

fsum = 0
a, b = 0, 1

while a < 4*10**6:
    if a % 2 == 0:
        fsum += a
    a, b = b, a + b

fsum

Out[3]:

time: 7 ms

Problem 4 ¶

Find the largest palindrome made from the product of two 3-digit numbers.

In [5]:

def is_palindrome(n):
    return str(n) == str(n)[::-1]

max(n1*n2 for n1, n2 in combinations_with_replacement(range(100, 1000), 2) if is_palindrome(n1*n2))

Out[5]:

time: 308 ms

Problem 5 ¶

What is the smallest number divisible by each of the numbers 1 to 20?

In [6]:

factordict = defaultdict(lambda: 0, {})
for i in range(2, 21):
    for prime, power in factorint(i).items():
        factordict[prime] = max(factordict[prime], power)

prod([prime ** power for prime, power in factordict.items()])

Out[6]:

232792560

time: 5 ms

Problem 6 ¶

Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

In [7]:

sum(range(101))**2 - sum([n**2 for n in range(101)])

Out[7]:

25164150

time: 12 ms

Problem 7 ¶

Find the 10001st prime.

In [8]:

from sympy import prime

prime(10001)

Out[8]:

time: 74 ms

Problem 8 ¶

Discover the largest product of 13 consecutive digits in the 1000-digit number.

In [9]:

s = open('p/p008_number.txt').read()

max([prod(int(digit) for digit in number) for number in [s[i:i+13] for i in range(len(s) - 13)]])

Out[9]:

23514624000

time: 9 ms

Problem 9 ¶

Find the only Pythagorean triplet, {a, b, c}, for which a + b + c = 1000.

In [10]:

for a, b in combinations(range(1, 1000), 2):
    c = 1000 - a -b
    if a**2 + b**2 == c**2:
        break
        
a*b*c

Out[10]:

31875000

time: 242 ms

Problem 10 ¶

Calculate the sum of all the primes below two million.

In [11]:

from sympy import sieve

sum(sieve.primerange(2, 2000000))

Out[11]:

142913828922

time: 745 ms

Problem 11 ¶

What is the greatest product of four numbers on the same straight line in the 20 by 20 grid?

In [12]:

with open('p/p011_grid.txt') as inp:
    a = [[int(elem) for elem in line.split()] for line in inp]

grid = [(x,y) for x in range(20) for y in range(20)]

def up(pos, n):
    for i in range(n):
        yield pos[0] -i, pos[1]

def right(pos, n):
    for i in range(n):
        yield pos[0], pos[1] +i

def up_right(pos, n):
    for i in range(n):
        yield pos[0] -i, pos[1] +i
        
def down_right(pos, n):
    for i in range(n):
        yield pos[0] +i, pos[1] +i

directions = [up, right, up_right, down_right]

def value(pos):
    return a[pos[0]][pos[1]] if pos in grid else 0

def adjacent_product():
    for pos, direction in product(grid, directions):
        yield prod([value(pos) for pos in direction(pos, 4)])

max(adjacent_product())

Out[12]:

70600674

time: 81 ms

Problem 12 ¶

What is the value of the first triangle number to have over five hundred divisors?

In [13]:

from sympy import divisors

def triangle_numbers():
    t = 0
    for n in count(1):
        t += n
        yield t

for number in triangle_numbers():
    if len(divisors(number)) > 500:
        break

number

Out[13]:

76576500

time: 1.01 s

Problem 13 ¶

Find the first ten digits of the sum of one-hundred 50-digit numbers.

In [14]:

numbers = [int(line) for line in open('p/p013_numbers.txt')]

int(str(sum(numbers))[:10])

Out[14]:

5537376230

time: 3 ms

Problem 14 ¶

Which starting number, under one million, produces the longest collatz chain?

In [15]:

@lru_cache(maxsize=None)
def chain_length(n):
    if n == 1:
        return 1
    if n % 2 == 0:
        return(chain_length(n / 2)) + 1
    else:
        return chain_length(n * 3 + 1) + 1
    
max(range(1,10**6), key=chain_length)

Out[15]:

time: 2.31 s

Problem 15 ¶

Starting in the top left corner in a 20 by 20 grid, how many routes are there to the bottom right corner?

Analysis: To get to the corner, you need to go down exactly 20 times in a total of 40 moves, with each step being a choice between two options. This is given by the binomial distribution.

In [16]:

from sympy import binomial

binomial(40, 20)

Out[16]:

137846528820

time: 3 ms

Problem 16 ¶

What is the sum of the digits of the number 2^1000?

In [17]:

sum(map(int, str(2**1000)))

Out[17]:

time: 19 ms

Problem 17 ¶

How many letters would be needed to write all the numbers in words from 1 to 1000?

In [18]:

import inflect
p = inflect.engine()

sum(len(p.number_to_words(n).replace('-', '').replace(' ', '')) for n in range(1, 1001))

Out[18]:

time: 76 ms

Problem 18 ¶

Find the maximum sum traveling from the top of the triangle to the base.

Analysis: Start at the bottom and go up. Assign each tile the value of itself plus the largest of the two tiles that can lead to it. This will be the maximum sum leading to that tile.

In [19]:

with open('p/p018_triangle.txt') as inp:
    a = [[int(elem) for elem in line.split()] for line in inp]

for row in reversed(range(len(a) -1)):
    for col in range(len(a[row])):
        a[row][col] += max(a[row+1][col], a[row+1][col +1])

a[0][0]

Out[19]:

time: 7 ms

Problem 19 ¶

How many Sundays fell on the first of the month during the twentieth century?

In [20]:

from datetime import datetime

sum(1 for year in range(1901, 2001) for month in range(1, 13) if datetime(year, month, 1).isoweekday() == 7)

Out[20]:

time: 10 ms

Problem 20 ¶

Find the sum of digits in 100!

In [21]:

from math import factorial

sum(map(int, str(factorial(100))))

Out[21]:

time: 18 ms

Problem 21 ¶

Evaluate the sum of all the amicable numbers under 10000.

In [22]:

from sympy.ntheory import divisors

def sum_of_proper_divisors(n):
    return sum(divisors(n)) - n

def amicable_numbers():
    for a in range(1, 10000):
        b = sum_of_proper_divisors(a)
        if sum_of_proper_divisors(b) == a and a != b:
            yield a
        
sum(amicable_numbers())

Out[22]:

time: 350 ms

Problem 22 ¶

What is the total of all the name scores in the file of first names?

In [23]:

words = open('p/p022_names.txt').read().strip('"').split('","')

def value(word):
    return sum(ord(char) - ord('A') + 1 for char in word)

sum((i+1) * value(word) for i, word in enumerate(sorted(words)))

Out[23]:

871198282

time: 19 ms

Problem 23 ¶

Find the sum of all the positive integers which cannot be written as the sum of two abundant numbers.

In [24]:

from sympy import divisors

def sum_of_proper_divisors(n):
    return sum(divisors(n)) - n

ceiling = 28123

abundant_numbers = [n for n in range(1, ceiling) if sum_of_proper_divisors(n) > n]

sum(set(range(ceiling + 1)) - set(i+j for i, j in combinations_with_replacement(abundant_numbers, 2)))

Out[24]:

time: 3.95 s

Problem 24 ¶

What is the millionth lexicographic permutation of the digits 0, 1, 2, 3, 4, 5, 6, 7, 8 and 9?

In [25]:

int(''.join(list(permutations('0123456789'))[10**6-1]))

Out[25]:

2783915460

time: 870 ms

Problem 25 ¶

What is the first term in the Fibonacci sequence to contain 1000 digits?

In [26]:

from sympy import fibonacci

next(i for i in count(1) if len(str(fibonacci(i))) >= 1000)

Out[26]:

time: 757 ms

Problem 26 ¶

Find the value of d < 1000 for which 1/d contains the longest recurring cycle.

Analysis: https://en.wikipedia.org/wiki/Repeating_decimal#Fractions_with_prime_denominators. The length of the repetend (period of the repeating decimal) of 1/p is equal to the order of 10 modulo p.

In [27]:

def cycle_length(d):
    for k in range(1, d):
        if 10 ** k % d == 1:
            return k
    return 0

max(range(1000), key=cycle_length)

Out[27]:

time: 647 ms

Problem 27 ¶

Find a quadratic formula with terms below 1000 that produces the maximum number of primes for consecutive values of n.

In [28]:

from sympy import primerange, isprime

def quadratic(n, a, b):
    return n **2 + a*n + b

def quadratic_expressions():
    for a, b in product(range(-999, 1000), primerange(2, 1000)):
        n = next(n for n in count(1) if not isprime(quadratic(n, a, b)))
        yield n, a*b

max(quadratic_expressions(), key=itemgetter(0))[1]

Out[28]:

-59231

time: 1.65 s

Problem 28 ¶

What is the sum of both diagonals in a 1001 by 1001 spiral?

In [29]:

def diagonals(size):
    n = 1
    yield n
    for squaresize, corner in product(range(3, size + 1, 2), range(4)):
        n += squaresize - 1
        yield n
            
sum(diagonals(1001))

Out[29]:

669171001

time: 7 ms

Problem 29 ¶

How many distinct terms are in the sequence generated by a^b for 2 ≤ a ≤ 100 and 2 ≤ b ≤ 100?

In [30]:

len(set(a**b for a, b in product(range(2, 101), range(2, 101))))

Out[30]:

time: 24 ms

Problem 30 ¶

Find the sum of all the numbers that can be written as the sum of fifth powers of their digits.

Analysis: The highest sum of fifth powers of digits for a number with 7 digits is 9^5 * 7 = 413343, which has only 6 digits. So the number must have at most 6 digits.

In [31]:

ceiling = 10**6

sum(n for n in range(2, ceiling) if n == sum(int(digit)**5 for digit in str(n)))

Out[31]:

time: 3.99 s

Problem 31 ¶

How many different ways can 2 pounds be made using any number of coins?

Analysis: Dynamic programming. Approach described here: https://en.wikipedia.org/wiki/Change-making_problem

In [32]:

ways = [1] + [0] * 200
for coin in [1, 2, 5, 10, 20, 50, 100, 200]:
    for i in range(len(ways) - coin):
        ways[i + coin] += ways[i]

ways[-1]

Out[32]:

time: 4 ms

Problem 32 ¶

Find the sum of all products whose multiplicand/multiplier/product identity can be written as a 1 through 9 pandigital.

Analysis: For a x b to be a 9-digit number, either a is 1-digit and b is 4-digit, or a is 2-digit and b is 3-digit. So a can be at most 2 digit and b at most 4 digit. (This could easily be optimized further.)

In [33]:

def pandigital_products():
    for i, j in product(range(1, 100), range(100, 10000)):
        digits = ''.join((str(i*j), str(i), str(j)))
        if ''.join(sorted(digits)) == '123456789':
            yield i*j

sum(set(pandigital_products()))

Out[33]:

time: 2.33 s

Problem 33 ¶

Discover all the fractions with an curious cancelling method.

In [34]:

from fractions import Fraction

def curious_fractions():
    for num, den in product(range(10, 100), range(10, 100)):
        digits_num, digits_den = ([int(str(n)[i]) for i in range(2)] for n in (num, den))
        if digits_num[1] == digits_den[0] and digits_den[1] != 0 and digits_num[0] / digits_den[1] == num / den:
            yield Fraction(num, den)

prod(list(curious_fractions())).denominator

Out[34]:

time: 45 ms

Problem 34 ¶

Find the sum of all numbers which are equal to the sum of the factorial of their digits.

Analysis: 9! x 8 has only 7 digits, so an upper bound is 9! * 7.

In [35]:

ceiling = 10**5

sum(n for n in range(3, ceiling) if sum(factorial(int(digit)) for digit in str(n)) == n)

Out[35]:

time: 242 ms

Problem 35 ¶

How many circular primes are there below one million?

In [36]:

from sympy import sieve

def rotations(n):
    n = str(n)
    for i in range(1, len(n)):
        yield int(n[i:] + n[:i])

primes = set(sieve.primerange(2, 10**6))
        
sum(1 for n in primes if all(rotation in primes for rotation in rotations(n)))

Out[36]:

time: 206 ms

Problem 36 ¶

Find the sum of all numbers, less than one million, which are palindromic in base 10 and base 2.

In [37]:

def is_palindromic(n):
    return str(n) == str(n)[::-1]

def binary(n):
    return str(bin(n))[2:]

sum(n for n in range(10**6) if is_palindromic(n) and is_palindromic(binary(n)))

Out[37]:

time: 730 ms

Problem 37 ¶

Find the sum of the only eleven primes that are both truncatable from left to right and right to left.

In [38]:

from sympy import sieve, isprime

def is_truncatable(prime):
    s = str(prime) 
    return all(isprime(int(s[:i])) and int(isprime(s[i:])) for i in range(1, len(s)))

truncatable_primes = []
for prime in sieve:
    if prime > 7 and is_truncatable(prime):
        truncatable_primes.append(prime)
        
    if len(truncatable_primes) == 11:
        break

sum(truncatable_primes)

Out[38]:

time: 667 ms

Problem 38 ¶

What is the largest 1 to 9 pandigital 9-digit number that can be formed as the concatenated product of an integer with (1,2, ... , n)?

Analysis: Each increment of n adds at least as many digits as the number has. So n must be less than 10, and the number of digits in the base number multiplied by n must be at most 9. (This could easily be optimized further.)

In [39]:

def pandigital_concats():
    for n in range(2, 10):
        for i in range(1, 10**(9 // n)):
            s = ''.join(str(i * j) for j in range(1, n + 1))
            if ''.join(sorted(s)) == '123456789':
                yield int(s)

max(pandigital_concats())

Out[39]:

932718654

time: 50 ms

Problem 39 ¶

Find the perimeter <= 1000 for which there is the highest number of right-angled triangles with integer side lengths.

In [40]:

def perimeters():
    for a, b in product(range(500), range(500)):
        c = sqrt(a*a + b*b)
        if c == int(c) and a+b+c <= 1000:
            yield a+b+int(c)
            
Counter(perimeters()).most_common()[0][0]

Out[40]:

time: 181 ms

Problem 40 ¶

An irrational decimal fraction is created by concatenating the positive integers: If dn represents the nth digit, find the value of the following expression: d1 x d10 x d100 x d1000 x d10000 x d100000 x d1000000

In [41]:

fraction = '0.' + ''.join([str(n) for n in range(1, 200000)])

prod([int(fraction[10**power+1]) for power in range(7)])

Out[41]:

time: 71 ms

Problem 41 ¶

What is the largest n-digit pandigital prime that exists?

Analysis: 9 and 8 digit pandigital numbers are divisible by 3 and thus cannot be prime. We find the largest prime by starting from the highest pandigital number and going down until we find a prime.

In [42]:

from sympy import isprime

def number(s):
    return int(''.join(str(digit) for digit in s))

pandigital_numbers = [number(s) for ndigits in range(2, 8) for s in permutations(range(1, ndigits+1))]

max(n for n in pandigital_numbers if isprime(n))

Out[42]:

time: 82 ms

Problem 42 ¶

Find the number of triangle words in a list.

In [43]:

words = open('p/p042_words.txt').read().strip('"').split('","')

def word_value(word):
    return sum(ord(char) - ord('A') + 1 for char in word)

triangle_numbers = {0.5 * n * (n + 1) for n in range(30)}

sum(1 for word in words if word_value(word) in triangle_numbers)

Out[43]:

time: 256 ms

Problem 43 ¶

Find the sum of all pandigital numbers with a sub-string divisibility property.

In [44]:

from sympy import primerange

divisors = list(primerange(1, 18))

def pandigital_with_property():
    for perm in permutations('0123456789'):
        s = ''.join(perm)
        if not any(int(s[start+1: start+4]) % divisor for start, divisor in enumerate(divisors)) and s[0] != "0":
            yield int(s)
        
sum(pandigital_with_property())

Out[44]:

16695334890

time: 6.79 s

Problem 44 ¶

Find the pair of pentagonal numbers for which their sum and difference are pentagonal and the difference is minimized.

In [45]:

ceiling = 5000

pentagonals = set(n*(3*n-1)//2 for n in range(1, ceiling))

def pentagonal_pairs():
    for p1, p2 in combinations_with_replacement(pentagonals, 2):
        if abs(p1 - p2) in pentagonals and p1 + p2 in pentagonals:
            yield abs(p1 - p2)

min(pentagonal_pairs())

Out[45]:

time: 1.94 s

Problem 45 ¶

Find the next triangle number that is also pentagonal and hexagonal.

In [46]:

ceiling = 100000

tri = set([n*(n+1)//2 for n in range(ceiling)])
pent = set([n*(3*n-1)//2 for n in range(ceiling)])
hexa = set([n*(2*n-1) for n in range(ceiling)])

min(elem for elem in tri & pent & hexa if elem > 40755)

Out[46]:

1533776805

time: 93 ms

Problem 46 ¶

What is the smallest odd composite that is not a goldbach number (can be written as the sum of a prime and twice a square)?

In [47]:

from sympy import sieve

ceiling = 10000
primes = set(sieve.primerange(2, ceiling))
odd_composites = set(range(3, ceiling, 2)) - primes
squares = {i**2 for i in range(1, int(sqrt(ceiling)))}
goldbachs = {prime + 2*square for prime, square in product(primes, squares)}

min(list(odd_composites - goldbachs))

Out[47]:

time: 28 ms

Problem 47 ¶

Find the first four consecutive integers to have four distinct primes factors. What is the first of these numbers?

In [48]:

from sympy import factorint

nconsecutive = nfactors = 4
primefactorlist = []

for n in count(1):
    primefactorlist.append(factorint(n))
    if all(len(factors) == nfactors for factors in primefactorlist[n - nconsecutive: n+1]):
        break
        
n - nconsecutive + 1

Out[48]:

time: 1.65 s

Problem 48 ¶

Find the last ten digits of the series, 1^1 + 2^2 + 3^3 + ... + 1000^1000

In [49]:

int(str(sum([n**n for n in range(1, 1001)]))[-10:])

Out[49]:

9110846700

time: 11 ms

Problem 49 ¶

Find the members of a 4-digits sequence.

In [50]:

from sympy import isprime

for n1 in range(1000, 3340):
    n2, n3 = n1 + 3330, n1 + 6660
    if all(isprime(n) for n in [n1, n2, n3]) and sorted(str(n1))==sorted(str(n2))==sorted(str(n3)) and n1 != 1487:
        break
        
int(''.join([str(n1), str(n2), str(n3)]))

Out[50]:

296962999629

time: 38 ms

Problem 50 ¶

Which prime, below one-million, can be written as the sum of the most consecutive primes?

In [51]:

from sympy import sieve

ceiling = 10**6

primes = list(sieve.primerange(2, ceiling))
primeset = set(primes)

max_consecutive = 0
max_prime = 2
for i in range(len(primes) -1):
    consecutive = 0
    primesum = primes[i]
    while primesum < ceiling:
        if consecutive > max_consecutive and primesum in primeset:
            max_consecutive = consecutive
            max_prime = primesum
        consecutive += 1
        primesum += primes[i+consecutive]

max_prime

Out[51]:

time: 277 ms

Problem 51 ¶

Find the smallest prime which, by replacing part of the number (not necessarily adjacent digits) with the same digit, is part of an eight prime value family.

In [52]:

from sympy import sieve
from itertools import chain

ndigits = 6
family_size = 8
digits = [str(i) for i in range(10)]
primes = sieve.primerange(10**(ndigits-1), 10**ndigits)


def replace_xs(m):
    for i, prime in enumerate(m):
        if prime[0] == 'x':
            m[i] = prime.replace('x', '1')
        else:
            m[i] = prime.replace('x', '0') 
    return m


def subsets(s):
    return chain.from_iterable(combinations(s,n) for n in range(len(s)+1))


def replace_with_xs(prime, x_indices):
    p = list(str(prime))
    for i in x_indices:
        p[i] = 'x'
    return ''.join(p)


def primes_with_xs():
    """For each prime, find all possible x-containing numbers that could make that prime"""
    for prime, digit in product(primes, digits):
        for x_indices in subsets([i for i, elem in enumerate(str(prime)) if elem == digit]):
            yield replace_with_xs(prime, x_indices)
            
m = [prime for prime, copies in Counter(primes_with_xs()).items() if copies == family_size]

int(min(replace_xs(m)))

Out[52]:

time: 3.55 s

Problem 52 ¶

Find the smallest positive integer, x, such that 2x, 3x, 4x, 5x, and 6x, contain the same digits.

In [53]:

next(n for n in count(1) if all(sorted(str(n)) == sorted(str(n*power)) for power in range(2,7)))

Out[53]:

time: 405 ms

Problem 53 ¶

How many, not necessarily distinct, values of nCr, for 1 ≤ n ≤ 100, are greater than one-million?

In [54]:

sum(1 for n in range(101) for r in range(n) if factorial(n) / (factorial(r) * factorial(n-r)) > 1000000)

Out[54]:

time: 15 ms

Problem 54 ¶

How many poker hands does Player 1 win?

In [55]:

from itertools import groupby

def rank(hand):
    values = ['xx23456789TJQKA'.find(card[0]) for card in hand]
    suits = [card[1] for card in hand]
    flush = len(set(suits)) == 1
    straight = any(i for i in range(1, 11) if all(i+j in values for j in range(5)))
    counts = Counter(values).most_common()
    order = []
    for key, item in groupby(counts, itemgetter(1)):
        order += sorted([elem[0] for elem in item], reverse=True)
    if straight and flush:
        return [8] + order
    if counts[0][1] == 4:
        return [7] + order
    if counts[0][1] == 3 and counts[1][1] == 2:
        return [6] + order
    if flush:
        return [5] + order
    if straight:
        return [4] + order
    if counts[0][1] == 3:
        return [3] + order
    if counts[0][1] == 2 and counts[1][1] == 2:
        return [2] + order
    if counts[0][1] == 2:
        return [1] + order
    return [0] + order


def hands():
    with open('p/p054_poker.txt') as inp:
        for line in inp.readlines():
            yield line.split()[:5], line.split()[5:]

sum(rank(hand1) > rank(hand2) for hand1, hand2 in hands())

Out[55]:

time: 103 ms

Problem 55 ¶

How many Lychrel numbers are there below ten-thousand?

In [56]:

max_iterations = 50

def is_lychrel(n):
    for _ in range(max_iterations):
        n = str(int(n) + int(str(n)[::-1]))
        if n == n[::-1]:
            return False
    return True

sum(is_lychrel(n) for n in range(10000))

Out[56]:

time: 76 ms

Problem 56 ¶

Considering natural numbers of the form, a^b, where a, b < 100, what is the maximum digital sum?

In [57]:

max(sum(int(digit) for digit in str(a**b)) for a, b in product(range(100), range(100)))

Out[57]:

time: 248 ms

Problem 57 ¶

In the first one-thousand expansions of the square root of 2, how many fractions contain a numerator with more digits than denominator?

In [58]:

from fractions import Fraction

co = 0
expander = 2

for _ in range(1000):
    expander = Fraction(2 + 1/expander)
    fraction = Fraction(1 + 1/expander)

    if len(str(fraction.numerator)) > len(str(fraction.denominator)):
        co += 1
        
co

Out[58]:

time: 83 ms

Problem 58 ¶

What is the side length of the square spiral for which the ratio of primes along both diagonals first falls below 10%?

In [59]:

from sympy import isprime

diagonal_number = 1
primes = 0
for sidelength in count(3, 2):
    for corner in range(4):
        diagonal_number += sidelength - 1
        if isprime(diagonal_number):
            primes += 1

    if primes / (sidelength * 2 - 1) < 0.1:
        break

sidelength

Out[59]:

time: 685 ms

Problem 59 ¶

Find the key that decrypts a file.

Analysis: We check all possible keys, and identify the key that finds the most common words.

In [60]:

common_words = {'the', 'and', 'be', 'of', 'that', 'have', 'for', 'it', 'not', 'on', 'with', 'you'}
ciphertext = [int(s) for line in open("p/p059_cipher.txt") for s in line.split(',')]

max_words_found = 0
for key in product(range(97, 123), repeat=3):
    decrypted = [chr(ciphertext[i] ^ key[i % len(key)]) for i in range(len(ciphertext))]
    text = ''.join(decrypted).split()
    words_found = len([word for word in text if word in common_words])
    if words_found > max_words_found:
        max_words_found = words_found
        best_sum = sum([ciphertext[i] ^ key[i % len(key)] for i in range(len(ciphertext))])

best_sum

Out[60]:

time: 7.15 s

Problem 60 ¶

Find the lowest sum for a set of five primes for which any two primes concatenate to produce another prime.

Analysis: We add the primes recursively one at a time, so that all primes in the set concatenate with each other.

In [61]:

from sympy import sieve, isprime

primes = list(sieve.primerange(2, 10000))

def all_concats_are_prime(primelist, prime_index):
    primelist = [str(primes[i]) for i in primelist]
    p1 = str(primes[prime_index])
    return all(isprime(int(p1 + p2)) and isprime(int(p2 + p1)) for p2 in primelist)

def find_primeset(primelist):
    if len(primelist) == 5:
        return primelist

    start = primelist[-1] if primelist else 0

    for i in range(start, len(primes)):
        if all_concats_are_prime(primelist, i):
            result = find_primeset(primelist + [i])
            if result:
                return result

sum([primes[p] for p in find_primeset([])])

Out[61]:

time: 12.3 s

Problem 61 ¶

Find the sum of the only ordered set of six cyclic 4-digit numbers for which each polygonal type: triangle, square, pentagonal, hexagonal, heptagonal, and octagonal, is represented by a different number in the set.

Analysis: Build the sets recursively. Each new addition to the set must be cyclic with the previous number, and be of a new polygonal type.

In [62]:

low, high = 1000, 10000

def gen(func, low, high):
    for n in count(1):
        if func(n) > high:
            break
        if func(n) >= low:
            yield func(n)

polygons = {
    'triangle': set(gen(lambda n: n * (n + 1) // 2, low, high)),
    'square': set(gen(lambda n: n**2, low, high)),
    'pentagonal': set(gen(lambda n: n * (3 * n -1) // 2, low, high)),
    'hexagonal': set(gen(lambda n: n * (2*n -1), low, high)),
    'heptagonal': set(gen(lambda n: n * (5*n -3) // 2, low, high)),
    'octogonal': set(gen(lambda n: n * (3*n -2), low, high)),
}

def cycle(number):
    for i in range(10, 100):
        yield int(str(number)[-2:] + str(i))
    

def get_chain(starting_number, remaining_polygons, number):
    if not remaining_polygons:
        if starting_number in set(cycle(number)):
            return starting_number
        return

    for next_number, polygon in product(cycle(number), remaining_polygons):
        if next_number in polygons[polygon]:
            result = get_chain(starting_number, remaining_polygons.copy() - {polygon}, next_number)
            if result:
                return next_number + result

for number in polygons['pentagonal']:
    result = get_chain(number, set(polygons) - {'pentagonal'}, number)
    if result:
        break
        
result

Out[62]:

time: 140 ms

Problem 62 ¶

Find the smallest cube for which exactly five permutations of its digits are cube.

In [63]:

d = defaultdict(list)
for n in range(10000):
    d[str(sorted(str(n**3)))].append(n)

min(min(basenumbers)**3 for basenumbers in d.values() if len(basenumbers) == 5)

Out[63]:

127035954683

time: 54 ms

Problem 63 ¶

How many n-digit positive integers exist which are also an nth power?

In [64]:

sum(len(str(a**b)) == b for a in range(1,100) for b in range(1, 100))

Out[64]:

time: 17 ms

Problem 64 ¶

How many continued fractions for N ≤ 10000 have an odd period?

Analysis: Algorithm from https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Algorithm

In [65]:

from math import floor

def continued_fraction(n):
    m = 0
    d = 1
    a0 = a = floor(n ** 0.5)
    while 2 * a0 != a:
        m =  d * a - m
        d = (n - (m ** 2)) / d
        a = floor((a0 + m) / d)
        yield a
   
def non_squares_up_to(hi):
    for n in range(1, hi):
        if sqrt(n) != floor(sqrt(n)):
            yield n

sum(len(list(continued_fraction(n))) % 2 != 0 for n in non_squares_up_to(10000))

Out[65]:

time: 310 ms

Problem 65 ¶

Find the sum of digits in the numerator of the 100th convergent of the continued fraction for e.

In [66]:

def continued_fraction(i):
    if i % 3 == 2:
        return i // 3 * 2 + 2
    else:
        return 1

numerator, denominator = 2, 1
for i in range(1, 100):
    numerator, denominator = continued_fraction(i) * numerator + denominator, numerator

sum(int(digit) for digit in str(numerator))

Out[66]:

time: 7 ms

Problem 66 ¶

Solve diophantine equations.

Analysis: Method from https://en.wikipedia.org/wiki/Chakravala_method

In [67]:

def chakravala(N):
    sq = int(floor(sqrt(N)))
    a, b, k, m = sq, 1, sq**2 - N, sq
    while k != 1:
        m = ((m + sq) // k) * k - m 
        a, b, k = (a*m + N*b) // abs(k), (a + b*m) // abs(k), (m**2 - N) // k
    return a
 
max_x, max_D = 0, 0
for D in range(1, 1001):
    if floor(sqrt(D)) != sqrt(D):
        x = chakravala(D)
        if x > max_x:
            max_x, max_D = x, D
max_D

Out[67]:

time: 46 ms

Problem 67 ¶

Find the maximum sum traveling from the top of the triangle to the base.

Analysis: Reuse of code from problem 18.

In [68]:

with open('p/p067_triangle.txt') as inp:
    a = [[int(elem) for elem in line.split()] for line in inp]

for row in reversed(range(len(a) -1)):
    for col in range(len(a[row])):
        a[row][col] += max(a[row+1][col], a[row+1][col +1])

a[0][0]

Out[68]:

time: 22 ms

Problem 68 ¶

What is the maximum 16-digit string for a "magic" 5-gon ring?

Analysis: For each permutation, use only one of the five symmetrical configurations, and look only at permutations where all arms have equal sum.

In [69]:

def magic_strings():
    arms = [[0, 5, 6], [1, 6, 7], [2, 7, 8], [3, 8, 9], [4, 9, 5]]
    for ring in permutations(range(1, 11)):
        if all(ring[0] < ring[i] for i in range(1, 5)):
            if len(set(sum(ring[i] for i in arms[j]) for j in range(5))) == 1:
                s = "".join((str(ring[arms[i][j]]) for i in range(5) for j in range(3)))
                if len(s) == 16:
                    yield s

max(magic_strings())

Out[69]:

'6531031914842725'

time: 7.39 s

Problem 69 ¶

Find the value of n ≤ 1,000,000 for which n/φ(n) is a maximum.

Analysis: n/φ(n) is lowest when n has the most prime factors compared to its size. This is the case when n is a primorial

In [70]:

from sympy.ntheory.generate import primorial

primorial(max(takewhile(lambda n: primorial(n) <= 10**6, count(1))))

Out[70]:

time: 3 ms

Problem 70 ¶

Find the value of n, 1 < n < 10^7, for which φ(n) is a permutation of n and the ratio n/φ(n) produces a minimum

Analysis: φ is lowest when its the product of exactly two primes. The algorithm searches through pairs of primes and saves the best pair product.

In [71]:

from sympy import sieve

ceiling = 10**7

primes = sieve.primerange(2, int(sqrt(1.3 * int(ceiling))))

def phi_is_permutation():
    for p1, p2 in product(primes, repeat=2):
        n = p1 * p2
        phi = (1 - p1) * (1 - p2)
        if n <= ceiling and sorted(str(n)) == sorted(str(phi)):
            yield [n / phi, n]

min(phi_is_permutation(), key=itemgetter(0))[1]

Out[71]:

time: 598 ms

Problem 71 ¶

By listing the set of reduced proper fractions for d ≤ 1,000,000 in ascending order of size, find the numerator of the fraction immediately to the left of 3/7.

In [72]:

from fractions import Fraction

def find_nearest_fraction(fraction, ceiling):
    while 1:
        new_fraction = (fraction - Fraction(1, ceiling*10)).limit_denominator(ceiling)
        if new_fraction < fraction:
            return new_fraction

find_nearest_fraction(Fraction(3, 7), 10**6).numerator

Out[72]:

time: 4 ms

Problem 72 ¶

How many elements would be contained in the set of reduced proper fractions for d ≤ 1,000,000?

Analysis: Go through the primes in order. For each prime, update all the values of d and discount the denominators for which the fraction could be reduced by that prime.

In [73]:

ceiling = 10**7

d = list(range(ceiling+1))
for n in range(2, ceiling+1):
    if d[n] == n:
        for k in range(n, ceiling+1, n):
            d[k] -= d[k] // n
sum(d) - 1

Out[73]:

30396356427241

time: 15.1 s

Problem 73 ¶

How many fractions lie between 1/3 and 1/2 in the sorted set of reduced proper fractions for d ≤ 12,000?

In [74]:

from math import gcd

ceiling = 12000

sum(1 for denom in range(4, ceiling + 1) for num in range(denom//3 + 1, denom//2 + 1) if gcd(num, denom) == 1)

Out[74]:

time: 4.01 s

Problem 74 ¶

How many chains, with a starting number below one million, contain exactly sixty non-repeating terms?

In [75]:

chain_length = {
    169: 3,
    363601: 3,
    1454: 3,
    871: 2,
    45361: 2,
    872: 2,
    45362: 2
}

@lru_cache(2**20)
def get_chain(n):
    for steps in count(0):
        if n in chain_length:
            return steps + chain_length[n]
        next_n = sum(factorial(int(digit)) for digit in str(n))
        if next_n == n:
            return steps + 1
        n = next_n

for number in range(1, 10**6):
    chain_length[number] = get_chain(number)

sum(length == 60 for length in chain_length.values())

Out[75]:

time: 9.02 s

Problem 75 ¶

Given that L is the length of the wire, for how many values of L ≤ 1,500,000 can exactly one integer sided right angle triangle be formed?

Analysis: Formula from https://en.wikipedia.org/wiki/Pythagorean_triple#Generating_a_triple

In [76]:

from math import gcd

ceiling = 1500000

def perimeters():
    for n in range(1, int(sqrt(ceiling))):
        for m in range(n+1, int(sqrt(ceiling)), 2):
            if gcd(m, n) == 1:
                for k in count(1):
                    perimeter = k*(m**2 - n**2 + 2*m*n + m**2 + n**2)
                    if perimeter > ceiling:
                        break
                    yield perimeter

sum(value == 1 for value in Counter(perimeters()).values())

Out[76]:

time: 2.82 s

Problem 76 ¶

How many different ways can one hundred be written as a sum of at least two positive integers?

In [77]:

from sympy.functions.combinatorial.numbers import nT

nT(100) -1

Out[77]:

190569291

time: 47 ms

Problem 77 ¶

What is the first value which can be written as the sum of primes in over five thousand different ways?

Analysis: Reuse of code from problem 31.

In [78]:

from sympy import sieve

ceiling = 100

ways = [1] + [0] * ceiling
for prime in sieve.primerange(2, ceiling+1):
    for i in range(len(ways) - prime):
        ways[i + prime] += ways[i]
     
min(i for i, number_of_ways in enumerate(ways) if number_of_ways > 5000)

Out[78]:

time: 5 ms

Problem 78 ¶

Let p(n) represent the number of different ways in which n coins can be separated into piles. Find the least value of n for which p(n) is divisible by one million.

Analysis: Solution from https://en.wikipedia.org/wiki/Partition_(number_theory):
p(k) = p(k − 1) + p(k − 2) − p(k − 5) − p(k − 7) + p(k − 12) + p(k − 15) − p(k − 22) − ...

In [79]:

from itertools import cycle

partitions = [1]
pentagonals = sorted([n*(3*n - 1)//2 for n in range(-250, 250)])

for n in count(1):
    signs = cycle([1, 1, -1, -1])
    use_pentagonals = [p for p in pentagonals if 0 < p <= n]

    number_of_partitions = sum(sign * partitions[n - p] for sign, p in zip(signs, use_pentagonals))
    if number_of_partitions % 10**6 == 0:
        break
    
    partitions.append(number_of_partitions)

n

Out[79]:

time: 7.56 s

Problem 79 ¶

Analyse the file so as to determine the shortest possible secret passcode of unknown length.

Analysis: Sort the numbers in order of how many different numbers appear after it.

In [80]:

logins = set([line.strip() for line in open('p/p079_keylog.txt')])

numbers_pressed_after = defaultdict(set)
for login in logins:
    numbers_pressed_after[login[0]].add(login[1])
    numbers_pressed_after[login[0]].add(login[2])
    numbers_pressed_after[login[1]].add(login[2])
    numbers_pressed_after[login[2]].add('')
 
''.join([t[1] for t in sorted([(len(v), k) for k, v in numbers_pressed_after.items()], reverse=True)])

Out[80]:

'73162890'

time: 10 ms

Problem 80 ¶

For the first one hundred natural numbers, find the total of the digital sums of the first one hundred decimal digits for all the irrational square roots.

In [81]:

from decimal import Decimal, Context

sum(int(digit) for n in range(2, 101) for digit in str(Decimal(n).sqrt(Context(prec=102))).replace('.', '')[:-2])

Out[81]:

time: 18 ms

Problem 81 ¶

Find the minimal path sum in a matrix from the top left to the bottom right by only moving right and down.

In [82]:

import networkx as nx

matrix = {(row, col): int(n) for row, line in enumerate(open('p/p081_matrix.txt')) for col, n in enumerate(line.split(','))}

start, end = (0,-1), max(matrix)

matrix[start] = 0

directions = {(1,0), (0,1)}

def moves(tile):
    return {(tile[0] + move[0], tile[1] + move[1]) for move in directions} & set(matrix)

G = nx.DiGraph()
for tile in matrix:
    for adjacent_tile in moves(tile):
        G.add_edge(tile, adjacent_tile, weight = matrix[adjacent_tile])

nx.dijkstra_path_length(G, source=start, target=end)

Out[82]:

time: 2.51 s

Problem 82 ¶

Find the minimal path sum in a matrix from the left column to the right column, only moving up, down and right.

In [83]:

import networkx as nx

matrix = {(row, col): int(n) for row, line in enumerate(open('p/p081_matrix.txt')) for col, n in enumerate(line.split(','))}

nrows, ncols = max(matrix)

starts = {(row, -1) for row in range(nrows)}
ends = {(row, ncols) for row in range(nrows)}

for start in starts:
    matrix[start] = 0

directions = {(0, -1), (0, 1), (-1,0), (1,0)}

def moves(tile):
    return {(tile[0] + move[0], tile[1] + move[1]) for move in directions} & set(matrix)

G = nx.DiGraph()
for tile in matrix:
    for adjacent_tile in moves(tile) & set(matrix):
        G.add_edge(tile, adjacent_tile, weight = matrix[adjacent_tile])
        
paths = nx.single_source_dijkstra_path_length(G, source=(0, -1))

min(val for tile, val in paths.items() if tile in ends)

Out[83]:

time: 2.91 s

Problem 83 ¶

Find the minimal path sum ain matrix from the top left to the bottom right by moving left, right, up, and down.

In [84]:

import networkx as nx

matrix = {(row, col): int(n) for row, line in enumerate(open('p/p081_matrix.txt')) for col, n in enumerate(line.split(','))}

start, end = (0,-1), max(matrix)

matrix[start] = 0

directions = {(-1,0), (0,-1), (1,0), (0,1)}

def moves(tile):
    return {(tile[0] + move[0], tile[1] + move[1]) for move in directions} & set(matrix)

G = nx.DiGraph()
for tile in matrix:
    for adjacent_tile in moves(tile):
        G.add_edge(tile, adjacent_tile, weight = matrix[adjacent_tile])

nx.dijkstra_path_length(G, source=start, target=end)

Out[84]:

time: 2.04 s

Problem 84 ¶

Find the squares that you land most often on in a game of monopoly.

Analysis: We ignore the order of cards, and instead treat it as a random pick each turn. Likewise three double rolls in a row is treated as an independent event with likelihood 1/n^3.

In [85]:

from random import choice, randint

class Tiles():
    GO = 0
    CC1 = 2
    R1 = 5
    CH1 = 7
    JAIL = 10
    C1 = 11
    U1 = 12
    R2 = 15
    CC2 = 17
    CH2 = 22
    E3 = 24
    R3 = 25
    U2 = 28
    G2J = 30
    CC3 = 33
    R4 = 35
    CH3 = 36
    H2 = 39

dice_sides = 4
dice = range(1, dice_sides + 1)
community_chest = [Tiles.GO, Tiles.JAIL] + ['']*14

next_R = {
    Tiles.CH1: Tiles.R2,
    Tiles.CH2: Tiles.R3,
    Tiles.CH3: Tiles.R1
}

next_U = {
    Tiles.CH1: Tiles.U1,
    Tiles.CH2: Tiles.U2,
    Tiles.CH3: Tiles.U1
}

chance = [Tiles.GO, Tiles.JAIL, Tiles.C1, Tiles.E3, Tiles.H2, Tiles.R1, next_R, 
          next_R, next_U, 'Go back 3 squares'] + ['']*6

pos = Tiles.GO


def get_next(pos):
    
    roll = (choice(dice), choice(dice))

    if roll[0] == roll[1] and randint(1, dice_sides**2) == 1:
        return Tiles.JAIL
    
    pos += sum(roll)
    if pos >= 40:
        pos -= 40
    
    if pos in {Tiles.CC1, Tiles.CC2, Tiles.CC3}:
        card = choice(community_chest)
        if type(card) is int:
            return card
        
    if pos in {Tiles.CH1, Tiles.CH2, Tiles.CH3}:
        card = choice(chance)
        if type(card) is int:
            return card
        elif type(card) is dict:
            return card[pos]
        elif card == 'Go back 3 squares':
            pos -= 3 
            if pos < 0:
                pos += 40
    
    if pos == Tiles.G2J:
        return Tiles.JAIL
    
    return pos
        
c = Counter()
pos = Tiles.GO
for _ in range(10**5):
    pos = get_next(pos)
    c[pos] += 1

result = c.most_common(3)
''.join(str(item[0]) for item in result)

Out[85]:

'101516'

time: 499 ms

Problem 85 ¶

Although there exists no rectangular grid that contains exactly two million rectangles, find the area of the grid with the nearest solution.

In [86]:

def rectangles(x, y):
    return sum((x + 1 - a) * (y + 1 - b) for a, b in product(range(1, x+1), range(1, y+1)))

prod(min(product(range(1, 100), repeat=2), key=lambda sides: abs(rectangles(*sides) - 2*10**6)))

Out[86]:

time: 4.48 s

Problem 86 ¶

Investigate routes in a cuboid.

In [87]:

def solutions_with_largest_side_equal_to(a):
    for bc in range(1, 2 * a + 1):
        route = sqrt(a**2 + bc**2)
        if route == round(route):
            for c in range(1, bc//2 + 1):
                b = bc - c
                if a >= b:
                    yield a, b, c

solutions = []
for M in count(1):
    solutions += solutions_with_largest_side_equal_to(M)
    if len(solutions) > 1000000:
        break

M

Out[87]:

time: 4.81 s

Problem 87 ¶

How many numbers below fifty million can be expressed as the sum of a prime square, prime cube, and prime fourth power?

In [88]:

from sympy import sieve
from fractions import Fraction

ceiling = 5*10**7

def powers(power):
    max_root = round(pow(ceiling, Fraction(1, power))) + 1
    for n in sieve.primerange(2, max_root):
        yield n**power

len(set(sum(triplet) for triplet in product(powers(2), powers(3), powers(4)) if sum(triplet) < ceiling))

Out[88]:

time: 819 ms

Problem 88 ¶

What is the sum of all the minimal product-sum numbers for 2≤k≤12000?

In [89]:

min_solutions = defaultdict(lambda: 10**9)
ceiling = 12000

def get_solutions(nfactors, nproduct, nsum):
    k = nfactors + nproduct - nsum
    if k > ceiling:
        return
    min_solutions[k] = min(min_solutions[k], nproduct)

    for n in range(2, int(2*ceiling / nproduct) + 1):
        get_solutions(nfactors+1, n*nproduct, nsum+n)

get_solutions(1, 1, 1)

sum(set(val for k, val in min_solutions.items() if k >= 2))

Out[89]:

time: 6.94 s

Problem 89 ¶

Find the minimal form of roman numerals.

In [90]:

romans = [line.strip() for line in open('p/p089_roman.txt')]

savedict = {'VIIII': 1, 'IIII': 2, 'LXXXX': 1, 'XXXX': 2, 'DCCCC': 1, 'CCCC': 2}

sum(val for roman in romans for key, val in savedict.items() if key in roman)

Out[90]:

time: 6 ms

Problem 90 ¶

How many distinct arrangements of the two cubes allow for all of the square numbers to be displayed?

In [91]:

from itertools import chain

squares = {'01', '04', '06', '16', '25', '36', '46', '64', '81'}
numbers = [str(i) for i in range(9)] + ['6']


def configuration_works(cube1, cube2):
    possible_displays = {''.join(faces) for faces in chain(product(cube1, cube2), product(cube2, cube1))}
    return squares.issubset(possible_displays)


def get_solutions():
    for cube1, cube2 in combinations_with_replacement(combinations(numbers, 6), 2):      
        if configuration_works(cube1, cube2):
            yield (cube1, cube2)

len(list(get_solutions()))

Out[91]:

time: 448 ms

Problem 91 ¶

How many right triangles can be formed in a 50*50 grid?

In [92]:

from collections import namedtuple

Point = namedtuple('Point', 'x y')
origin = Point(0,0)


def length_of_side_squared(p1, p2):
    return (p2.x - p1.x)**2 + (p2.y - p1.y)**2


def right_angled_triangles(n):
    points = {Point(x,y) for x, y in product(range(0, n+1), repeat=2)} - {origin}
    
    for p1, p2 in combinations(points, 2):
        a2 = length_of_side_squared(p1, origin)
        b2 = length_of_side_squared(p2, origin)
        c2 = length_of_side_squared(p1, p2)
        if a2+b2==c2 or a2+c2==b2 or c2+b2==a2:
            yield([p1, p2, origin])

len(list(right_angled_triangles(50)))

Out[92]:

time: 13.3 s

Problem 92 ¶

A number chain is created by continuously adding the square of the digits in a number to form a new number until it has been seen before. How many starting numbers below ten million will arrive at 89?

In [93]:

@lru_cache(2**20)
def is89chain(n):
    if n == 1: return False
    if n == 89: return True
    return is89chain(sum((int(digit)**2 for digit in str(n))))

sum(is89chain(i) for i in range(2, 10**7))

Out[93]:

time: 57.1 s

Problem 93 ¶

Find the set of four distinct digits, a < b < c < d, for which the longest set of consecutive positive integers, 1 to n, can be obtained, giving your answer as a string: abcd.

In [94]:

all_operations = [lambda x, y: x + y, 
                  lambda x, y: x*y, 
                  lambda x, y: x - y, 
                  lambda x,y: x/y]


def targets(orig_digits):
    for a,b,c,d in permutations(orig_digits):
        for operations in product(all_operations, repeat=3):  
            target = operations[0](operations[1](operations[2](a, b), c), d)
            if target > 0 and round(target) == target:
                yield target

            target = operations[0](operations[1](a, b), operations[2](c, d))
            if target > 0 and round(target) == target:
                yield target

                
def consecutive_integers(integers):
    for i, elem in enumerate(set(integers)):
        if i+1 != elem:
            return i

result = max(combinations(range(1, 10), 4), key=lambda digits: consecutive_integers(targets(digits)))

''.join(str(digit) for digit in result)

Out[94]:

'1258'

time: 364 ms

Problem 94 ¶

Find the sum of the perimeters of all almost equilateral triangles with integral side lengths and area and whose perimeters do not exceed one billion.

In [95]:

from sympy.geometry import Triangle
from sympy.solvers.solvers import check_assumptions

def seq(n):
    if n == 0 or n == 1:
        return 1
    if n == 2:
        return 5

    return 3 * seq(n-1) + 3 * seq(n -2) - seq(n -3)


def get_perimeters(ceiling):
    for n in count(2):
        a = seq(n)
        if a*3 > ceiling:
            break
        triangle = Triangle(sss=(a, a, a+1))
        if check_assumptions(triangle.area, integer=True):
            yield triangle.perimeter
        triangle = Triangle(sss=(a, a, a-1))
        if check_assumptions(triangle.area, integer=True):
            yield triangle.perimeter

sum(get_perimeters(10**9))

Out[95]:

518408346

time: 929 ms

Problem 95 ¶

Find the smallest member of the longest amicable chain with no element exceeding one million.

In [96]:

from sympy import divisors
ceiling = 10**6

sumofdivisors = [1 for i in range(1, ceiling + 1)]
for n in range(2, ceiling + 1):
    for i in range(2*n, ceiling, n):
        sumofdivisors[i] += n

zeros = {0}

longest_chain = []
for start_n in range(10, ceiling):
    n = start_n
    chain = [n]
    while 1:
        n = sumofdivisors[n]

        if n in zeros:
            zeros |= set(chain)
            break

        if n == start_n:
            if len(chain) > len(longest_chain):
                longest_chain = chain
            break

        if n in chain or n > ceiling:
            break

        chain.append(n)

min(longest_chain)

Out[96]:

time: 11.3 s

Problem 96 ¶

Solve sudoku puzzles.

In [98]:

from dlxsudoku.sudoku import Sudoku

s = [l for l in open('p/p096_sudoku.txt').readlines() if l[0].isdigit()]
sudokus = [Sudoku(''.join(s[i*9:(i+1)*9])) for i in range(50)]
for sudoku in sudokus:
    sudoku.solve()
    
sum(int(sudoku.to_oneliner()[:3]) for sudoku in sudokus)

Out[98]:

time: 583 ms

Problem 97 ¶

Find the last ten digits of a prime number.

In [99]:

(28433 * pow(2, 7830457, 10**10) + 1) % 10**10

Out[99]:

8739992577

time: 8 ms

Problem 98 ¶

Find the largest square anagram word pairs in a list of words

In [100]:

words = [list(word.replace('"', '')) for word in open('p/p098_words.txt').readlines()[0].split(',')]

def is_anagram(w1, w2):
    return sorted(w1) == sorted(w2) and w1 != w2

def formed_number(w, d):
    return int(''.join([str(d[l]) for l in w]))

def is_square(n):
    return sqrt(n) == int(sqrt(n))


def square_anagrams():
    for w1, w2 in combinations(words, r=2):
        if is_anagram(w1, w2):
            for order in [order for numbers in combinations(range(10), len(w1)) for order in permutations(numbers)]:
                d = {w1[i]: order[i] for i in range(len(w1))}
                if d[w1[0]] != 0 and d[w2[0]] != 0:
                    nw1, nw2 = formed_number(w1, d), formed_number(w2, d)
                    if nw1 != nw2 and is_square(nw1) and is_square(nw2):
                        yield(nw1)
                        yield(nw2)

max(square_anagrams())

Out[100]:

time: 59.5 s

Problem 99 ¶

Determine which base and exponential pair has the greatest numerical value.

In [101]:

def get_log_values():
    for i, line in enumerate(open('p/p099_base_exp.txt')):
        base, exponent = [int(n) for n in line.split(',')]
        yield i+1, exponent * log(base)

max(get_log_values(), key=itemgetter(1))[0]

Out[101]:

time: 20 ms

Problem 100 ¶

By finding the first arrangement to contain over 10^12 = 1,000,000,000,000 discs in total, determine the number of blue discs that the box would contain.

Analysis: Sequence formula from https://oeis.org/A011900

In [102]:

from sympy import solve, Symbol

def seq(n):
    if n == 0 or n == -1:
        return 1
    return int(seq(n-1) * (seq(n-1) + 2) / seq(n-2))

total = Symbol('total')

for n in count(1):
    blue = seq(n)
    discs = max(solve(blue/total * (blue-1)/(total-1) -1/2, total))
    if discs > 10**12:
        break

blue

Out[102]:

756872327473

time: 1.83 s