Solve Euler 174

README.md

@@ -1,5 +1,7 @@
 # Project Euler
 
+![Euler Progress Banner](https://projecteuler.net/profile/failx.png)
+
 Project Euler discourages sharing solutions to the problem. Since there are
 plenty solutions out there on the internet I violate this rule for now. My goal
 for now is to finish one hundred problems. After that I might change the

python/e088.py

@@ -0,0 +1,90 @@
+
+def product(numbers):
+    r = 1
+    for n in numbers:
+        r = r * n
+    return r
+
+def _incrementer():
+    """ Generator function that returns natural numbers where the higher
+    digits are smaller or equal than the lower digits. """
+    n_digits = 2
+    r = [1] * n_digits
+    while True:
+        yield r
+        incremented = False
+        for i in range(n_digits - 1):
+            if r[i] < r[i + 1]:
+                r[i] += 1
+                incremented = True
+                for j in range(i):
+                    r[j] = 1
+                break
+        if incremented:
+            continue
+        elif r[-1] < 2000:
+            r[-1] += 1
+            for j in range(n_digits - 1):
+                r[j] = 1
+        else:
+            n_digits += 1
+            r = [1] * n_digits
+
+for k in range(40, 50):
+    s = product_sum_number(k)
+        print("k={}: {} = {}".format(k, product(s), s))
+
+
+def product_sum(k):
+    # Create initial list of canditates
+    canditates = []
+    for n in range(2, k + 1):
+        rest_sum = k - n
+        if rest_sum + n * 2 < 2 ** n:
+            break
+        canditates.append(tuple([2 for _ in range(n)]))
+    processed = set(canditates)
+    solutions = []
+
+    while canditates:
+        numbers = canditates.pop()
+        processed.add(numbers)
+
+        s = k - len(numbers) + sum(numbers)
+        p = product(numbers)
+        if s == p:
+            solutions.append(numbers)
+        elif s > p:
+            for i in range(len(numbers) - 1):
+                if numbers[i] < numbers[i + 1]:
+                    n = list(numbers)
+                    n[i] += 1
+                    t = tuple(n)
+                    if not t in processed:
+                        canditates.append(t)
+
+            n = list(numbers)
+            n[-1] += 1
+            t = tuple(n)
+            if not t in processed:
+                canditates.append(t)
+        else:
+            pass
+    result = min([product(ns) for ns in solutions])
+    return result
+
+
+def euler_088():
+    ns = set()
+    for k in range(2, 12001):
+        r = product_sum_2d(k)
+        print("k={}: {}".format(k, r))
+        ns.add(r)
+    return sum(ns)
+
+
+if __name__ == "__main__":
+    solution = euler_088()
+    print("e088.py: " + str(solution))
+    assert(solution == 7587457)
+

python/e090.py

@@ -1,4 +1,3 @@
-#
 def choose(xs, n):
     """
     Computes r choose n, in other words choose n from xs.

python/e094.py

@@ -0,0 +1,94 @@
+
+def odd(n):
+    return n % 2 == 1
+
+def is_square(n):
+    """ From the internet. I should implement my own version of this. """
+    ## Trivial checks
+    if type(n) != int:  ## integer
+        return False
+    if n < 0:      ## positivity
+        return False
+    if n == 0:      ## 0 pass
+        return True
+
+    ## Reduction by powers of 4 with bit-logic
+    while n&3 == 0:
+        n=n>>2
+
+    ## Simple bit-logic test. All perfect squares, in binary,
+    ## end in 001, when powers of 4 are factored out.
+    if n&7 != 1:
+        return False
+
+    if n==1:
+        return True  ## is power of 4, or even power of 2
+
+    ## Simple modulo equivalency test
+    c = n%10
+    if c in {3, 7}:
+        return False  ## Not 1,4,5,6,9 in mod 10
+    if n % 7 in {3, 5, 6}:
+        return False  ## Not 1,2,4 mod 7
+    if n % 9 in {2,3,5,6,8}:
+        return False
+    if n % 13 in {2,5,6,7,8,11}:
+        return False
+
+    ## Other patterns
+    if c == 5:  ## if it ends in a 5
+        if (n//10)%10 != 2:
+            return False    ## then it must end in 25
+        if (n//100)%10 not in {0,2,6}:
+            return False    ## and in 025, 225, or 625
+        if (n//100)%10 == 6:
+            if (n//1000)%10 not in {0,5}:
+                return False    ## that is, 0625 or 5625
+    else:
+        if (n//10)%4 != 0:
+            return False    ## (4k)*10 + (1,9)
+
+    ## Babylonian Algorithm. Finding the integer square root.
+    ## Root extraction.
+    s = (len(str(n))-1) // 2
+    x = (10**s) * 4
+    A = {x, n}
+    while x * x != n:
+        x = (x + (n // x)) >> 1
+        if x in A:
+            return False
+        A.add(x)
+    return True
+
+
+def euler_094():
+    n_max = 10 ** 9
+    n, total_peri = 2, 0
+    while n < n_max:
+        n += 2
+        for b2, c in [(n, n + 1), (n, n - 1)]:
+            b = b2 // 2
+            a_squared = c * c - b * b
+
+            if not is_square(a_squared):
+                continue
+
+            peri = 2 * c + b2
+            # print("b2={} c={} peri={}".format(b2, c, peri))
+
+            if peri > n_max:
+                return total_peri
+
+            total_peri += peri
+
+            if b2 > 10000:
+                n = int(n * 3.725) - 5
+                if odd(n):
+                    n -= 1
+
+
+if __name__ == "__main__":
+    solution = euler_094()
+    print("e094.py: " + str(solution))
+    assert(solution == 518408346)
+

python/e095.py

@@ -0,0 +1,48 @@
+from lib_misc import sum_proper_divisors
+
+
+def euler_095():
+    NO_CHAIN = set()
+    PROPER_DIVISOR_CACHE = {}
+
+    def sum_proper_divisors_cached(n):
+        if n in PROPER_DIVISOR_CACHE:
+            return PROPER_DIVISOR_CACHE[n]
+        r = sum_proper_divisors(n)
+        PROPER_DIVISOR_CACHE[n] = r
+        return r
+
+    def chain(n):
+        chain = []
+        next_number = n
+
+        while next_number not in chain:
+            if next_number >= 10**6 or next_number == 0 or next_number in NO_CHAIN:
+                for c in chain:
+                    NO_CHAIN.add(c)
+                return None
+            chain.append(next_number)
+            next_number = sum_proper_divisors_cached(next_number)
+        chain.append(next_number)
+        return chain
+
+    chains = [[]]
+    for i in range(1, 10**6 + 1):
+        c = chain(i)
+        if c and c[0] == c[-1]:
+            if len(c) == len(chains[0]):
+                # There might be multiple different chains with the same length
+                # so keep all of them to get the smalles element.
+                chains.append(c)
+            elif len(c) > len(chains[0]):
+                chains = [c]
+
+    elems = [e for c in chains for e in c]
+    return min(elems)
+
+
+if __name__ == "__main__":
+    solution = euler_095()
+    print("e095.py: " + str(solution))
+    assert(solution == 14316)
+

python/e098.py

@@ -0,0 +1,88 @@
+
+def get_words():
+    with open("../txt/e098.txt", "r") as f:
+        words = list(map(lambda w: w.strip().lower(), f.readlines()))
+    words.sort(key=len, reverse=True)
+    return words
+
+
+def get_pair_subsets(xs):
+    """ Returns all subsets of size two. """
+    r = [(xs[i], xs[j])
+            for i in range(len(xs))
+            for j in range(i + 1, len(xs))]
+    return r
+
+
+def find_anagrams(xs):
+    anagrams = {}
+    for x in xs:
+        x_tuple = tuple(sorted(x))
+        try:
+            anagrams[x_tuple].append(x)
+        except KeyError:
+            anagrams[x_tuple] = [x]
+    anagrams = [a for a in anagrams.values() if len(a) > 1]
+    return anagrams
+
+
+def find_anagram_pairs(xs):
+    anagrams = find_anagrams(xs)
+    r = [pair
+         for anagram in anagrams
+         for pair in get_pair_subsets(anagram)]
+    return r
+
+
+def words_match_numbers(words, numbers):
+    letters_a, letters_b = words
+    digits_a, digits_b = numbers
+
+    digit_to_letter = {}
+    for l, d in zip(letters_a, digits_a):
+        if d in digit_to_letter and digit_to_letter[d] != l:
+            # Each digit can only be assigned to one letter
+            return False
+        digit_to_letter[d] = l
+
+
+    letter_to_digit = {l: d for l, d in zip(letters_a, digits_a)}
+    for l, d in zip(letters_b, digits_b):
+        if letter_to_digit[l] != d:
+            return False
+    return True
+
+
+def find_pairs(words, squares):
+    pairs = []
+    anagram_pairs = find_anagram_pairs(words)
+    square_pairs = find_anagram_pairs(squares)
+    for word_pair in anagram_pairs:
+        for number_pair in square_pairs:
+            if words_match_numbers(word_pair, number_pair):
+                # pairs.append([word_pair, number_pair])
+                pairs.append(number_pair)
+    return pairs
+
+
+def euler_098():
+    squares = [str(i * i) for i in range(200000)]
+    words = get_words()
+
+    largest_number = 0
+    for i in range(3, 9):
+        words_i = list(filter(lambda w: len(w) == i, words))
+        squares_i = list(filter(lambda w: len(w) == i, squares))
+        pairs = find_pairs(words_i, squares_i)
+        for pair in pairs:
+            n = max(map(int, pair))
+            if n > largest_number:
+                largest_number = n
+    return largest_number
+
+
+if __name__ == "__main__":
+    solution = euler_098()
+    print("e098.py: " + str(solution))
+    assert(solution == 18769)
+

python/e100.py

@@ -0,0 +1,47 @@
+from fractions import Fraction
+
+
+def euler_100():
+    F_ONE_HALF = Fraction(1, 2)
+    prev_blue_discs, blue_discs = 10, 10
+    while True:
+        # Brute force search blue discs
+        blue_discs += 1
+
+        # Binary search red discs
+        red_upper = blue_discs // 2
+        red_lower = 1
+        while red_upper != red_lower:
+            red_discs = red_lower + (red_upper - red_lower) // 2
+            total_discs = blue_discs + red_discs
+            f1 = Fraction(blue_discs, total_discs)
+            f2 = Fraction(blue_discs - 1, total_discs - 1)
+            f = f1 * f2
+            if f == F_ONE_HALF:
+                # End criterion by Euler Problem
+                if total_discs > 10**12:
+                    return blue_discs
+
+                # The ration between two solutions seems to be constant so we
+                # can get the next larger candidate for blue_discs much
+                # quicker.
+                ratio = blue_discs / prev_blue_discs
+                # print("blue_discs={} ratio={}".format(blue_discs, ratio))
+                prev_blue_discs = blue_discs
+                blue_discs = int(blue_discs * ratio)
+                break
+            elif f < F_ONE_HALF:
+                red_upper = red_discs
+            else:
+                red_lower = red_discs
+
+            # No solution for current number of blue discs
+            if abs(red_upper - red_lower) <= 1:
+                break
+
+
+if __name__ == "__main__":
+    solution = euler_100()
+    print("e100.py: " + str(solution))
+    assert(solution == 756872327473)
+

python/e101.py

@@ -0,0 +1,37 @@
+import numpy as np
+
+
+def rev(xs):
+    return list(reversed(xs))
+
+
+def bop(degree, f):
+
+    if degree == 1:
+        return 1
+    A = [rev([n ** x for x in range(degree)])
+         for n in range(1, degree + 1)]
+    y = np.array([f(n) for n in range(1, degree + 1)])
+    x = np.linalg.solve(A, y)
+    ns = rev([(degree + 1) ** x for x in range(degree)])
+    return np.dot(x, ns)
+
+
+def euler_101():
+    # max_degree = 4
+    # def u(x):
+        # return x * x * x
+
+    max_degree = 10
+    def u(n):
+        return 1 - n + n**2 - n**3 + n**4 - n**5 + \
+               n**6 - n**7 + n**8 - n**9 + n**10
+
+    return int(sum([bop(d, u) for d in range(1, max_degree + 1)]))
+
+
+if __name__ == "__main__":
+    solution = euler_101()
+    print("e101.py: " + str(solution))
+    assert(solution == 37076114526)
+

python/e102.py

@@ -0,0 +1,51 @@
+from math import acos, sqrt, pi, atan
+import matplotlib.pyplot as plt
+
+
+def triangle(euler_text_line):
+    ps = list(map(int, euler_text_line.split(",")))
+    ps = [(ps[0], ps[1]), (ps[2], ps[3]), (ps[4], ps[5])]
+    return ps
+
+
+def angle_x(x, y):
+    if x > 0 and y == 0:
+        return 0
+    elif x < 0 and y == 0:
+        return pi
+    elif x == 0 and y > 0:
+        return pi / 2
+    elif x == 0 and y < 0:
+        return 3 / 2 * pi
+    elif x > 0 and y > 0:
+        return atan(y / x)
+    elif x < 0 and y > 0:
+        return pi - atan(y / abs(x))
+    elif x < 0 and y < 0:
+        return pi + atan(abs(y) / abs(x))
+    elif x > 0 and y < 0:
+        return 2 * pi - atan(abs(y) / x)
+    else:
+        raise Exception("Unhandled case")
+
+
+def contains_origin(triangle):
+    a, b, c = triangle
+    angles = sorted([angle_x(*a), angle_x(*b), angle_x(*c)])
+    a = angles[1] - angles[0]
+    b = angles[2] - angles[1]
+    c = angles[0] + 2 * pi - angles[2]
+    return not (a + b < pi or b + c < pi or c + a < pi)
+
+
+def euler_102():
+    with open("../txt/e102.txt", "r") as f:
+        triangles = list(map(triangle, f.readlines()))
+    return len([1 for t in triangles if contains_origin(t)])
+
+
+if __name__ == "__main__":
+    solution = euler_102()
+    print("e102.py: " + str(solution))
+    assert(solution == 228)
+

python/e103.py

@@ -0,0 +1,83 @@
+
+def build_sets(sets, numbers):
+    """ Sets is a list of two tuples. Each field represents the sum of the two
+    sets. For each new number, we can add the number to set one, add the
+    number to set two, or don't add the number to any sets. If there are no
+    more numbers we return the list so that we can check if there are tuples
+    that have the same sum. """
+    if not numbers:
+        return sets
+    current_number, numbers = numbers[0], numbers[1:]
+    new_sets = []
+    for t in sets:
+        new_sets.append(t)
+        new_sets.append((t[0] + current_number, t[1]))
+        new_sets.append((t[0], t[1] + current_number))
+    return build_sets(new_sets, numbers)
+
+
+def is_condition_one_met(s):
+    tuples = build_sets([(0, 0)], s)
+    for a, b in tuples:
+        if a != 0 and b != 0 and a == b:
+            return False
+    return True
+
+
+def is_condition_two_met(s):
+    """ The second condition says that if |B| > |C| then sum(B) > sum(C). Since
+    each line is sorted we take the two smallest integers and compare it to the
+    biggest. Then we take the three smallest integers and compare it to the two
+    biggest, and so on. Comparing the smallest with the biggest integers is the
+    worst case and if there is no violation of the condition it will be met for
+    all other subset pairs. """
+    half_len = (len(s) - 1) // 2
+    for i in range(1, half_len + 1):
+        if not sum(s[:i + 1]) > sum(s[-i:]):
+            # print(f"{s[:i + 1]=} < {s[-i:]=}")
+            return False
+    return True
+
+
+def is_special_sum_set(s):
+    """ Check for the two conditions and return True if they are both met. """
+    if not is_condition_two_met(s):
+        return False
+    if not is_condition_one_met(s):
+        return False
+    return True
+
+
+def next_set_heuristic(xs):
+    """ Calculate the estimated next optimum based on the formul given in the
+    question. """
+    b = xs[(len(xs) - 1) // 2]
+    next_set = [b]
+    for a in xs:
+        next_set.append(a + b)
+    return next_set
+
+
+def euler_103():
+    e_max = 50
+    min_set = None
+    for a in range(1, e_max):
+        for b in range(a + 1, e_max):
+            for c in range(b + 1, e_max):
+                for d in range(c + 1, e_max):
+                    for e in range(d + 1, e_max):
+                        for f in range(e + 1, e_max):
+                            for g in range(e + 1, e_max):
+                                s = [a, b, c, d, e, f, g,]
+                                if is_special_sum_set(s):
+                                    s_sum = sum(s)
+                                    if min_set is None or s_sum < sum(min_set):
+                                        min_set = s
+    return int("".join(map(str, min_set)))
+
+
+if __name__ == "__main__":
+    solution = euler_103()
+    print("e103.py: " + str(solution))
+    assert(solution == 20313839404245)
+

python/e104.py

@@ -0,0 +1,37 @@
+import math
+from itertools import permutations
+
+
+def first_n_digits(num, n):
+    # https://stackoverflow.com/questions/41271299/how-can-i-get-the-first-two-digits-of-a-number
+    p = int(math.log(num, 10)) - n + 1
+    if p < 0:
+        return num
+    return num // 10 ** (int(math.log(num, 10)) - n + 1)
+
+
+def is_end_pandigital(n, ps):
+    return (n % 10**9) in ps
+
+
+def is_start_pandigital(n, ps):
+    start = first_n_digits(n, 9)
+    return start in ps
+
+
+def euler_104():
+    a, b = 1, 1
+    ps = set([int("".join(map(str, p))) for p in permutations(range(1, 10))])
+    k = 3
+    while True:
+        a, b = b, a + b
+        if is_end_pandigital(b, ps) and is_start_pandigital(b, ps):
+            return k
+        k += 1
+
+
+if __name__ == "__main__":
+    solution = euler_104()
+    print("e104.py: " + str(solution))
+    assert(solution == 329468)
+

python/e105.py

@@ -0,0 +1,74 @@
+
+def build_sets(sets, numbers):
+    """ Sets is a list of two tuples. Each field represents the sum of the two
+    sets. For each new number, we can add the number to set one, add the
+    number to set two, or don't add the number to any sets. If there are no
+    more numbers we return the list so that we can check if there are tuples
+    that have the same sum. """
+    if not numbers:
+        return sets
+    current_number, numbers = numbers[0], numbers[1:]
+    new_sets = []
+    for t in sets:
+        new_sets.append(t)
+        new_sets.append((t[0] + current_number, t[1]))
+        new_sets.append((t[0], t[1] + current_number))
+    return build_sets(new_sets, numbers)
+
+
+def is_condition_one_met(s):
+    tuples = build_sets([(0, 0)], s)
+    for a, b in tuples:
+        if a != 0 and b != 0 and a == b:
+            return False
+    return True
+
+
+def is_condition_two_met(s):
+    """ The second condition says that if |B| > |C| then sum(B) > sum(C). Since
+    each line is sorted we take the two smallest integers and compare it to the
+    biggest. Then we take the three smallest integers and compare it to the two
+    biggest, and so on. Comparing the smallest with the biggest integers is the
+    worst case and if there is no violation of the condition it will be met for
+    all other subset pairs. """
+    half_len = (len(s) - 1) // 2
+    for i in range(1, half_len + 1):
+        if not sum(s[:i + 1]) > sum(s[-i:]):
+            # print(f"{s[:i + 1]=} < {s[-i:]=}")
+            return False
+    return True
+
+
+def is_special_sum_set(s):
+    """ Check for the two conditions and return True if they are both met. """
+    if not is_condition_two_met(s):
+        return False
+    if not is_condition_one_met(s):
+        return False
+    return True
+
+
+def load_sets():
+    """ Parse sets from file into a list of lists. Each line contains a list of
+    comma separated integers. We sort the integers in each line. """
+    def line_to_set(line):
+        return sorted(map(int, line.strip().split(",")))
+
+    with open("../txt/e105.txt") as f:
+        return list(map(line_to_set, f))
+
+
+def euler_105():
+    xs = load_sets()
+    s = 0
+    for i, x in enumerate(xs):
+        if is_special_sum_set(x):
+            s += sum(x)
+    return s
+
+
+if __name__ == "__main__":
+    solution = euler_105()
+    print("e105.py: " + str(solution))
+    assert(solution == 73702)
+

python/e106.py

@@ -0,0 +1,44 @@
+def repeat(n, xs):
+    return [xs for _ in range(n)]
+
+
+def sequence(xss, acc):
+    if not xss:
+        return acc
+    return sequence(xss[1:], [a + [x] for x in xss[0] for a in acc])
+
+
+def is_balanced(xs):
+    """ Check if subsets are balanced. That means for every lower number there
+    is a higher number that balances it out. """
+    ones_balanced, twos_balanced = True, True
+    ones, twos = 0, 0
+    for x in xs:
+        if x == 1:
+            ones += 1
+            twos -= 1
+        elif x == 2:
+            ones -= 1
+            twos += 1
+        if ones < 0:
+            ones_balanced = False
+        if twos < 0:
+            twos_balanced = False
+    return ones_balanced or twos_balanced
+
+
+def euler_106():
+    n = 12
+    # Generate all possible subsets where both subsets have the same number of
+    # elements and at least two elements per subset.
+    xs = repeat(n, [0, 1, 2])
+    xs = [x for x in sequence(xs, [[]])
+          if x.count(1) > 1
+          if x.count(1) == x.count(2)]
+    return len([1 for x in xs if not is_balanced(x)]) // 2
+
+
+if __name__ == "__main__":
+    solution = euler_106()
+    print("e106.py: " + str(solution))
+    assert(solution == 21384)

python/e107.py

@@ -0,0 +1,89 @@
+from dataclasses import dataclass
+from typing import List
+
+
+@dataclass
+class Edge:
+    weight: int
+    source_id: int
+    target_id: int
+
+
+@dataclass
+class Vertex:
+    id: int
+    edges: List[Edge]
+
+
+def load_matrix(filename):
+    def line_to_row(line):
+        r = []
+        for field in line.strip().split(","):
+            if field == "-":
+                r.append(0)
+            else:
+                r.append(int(field))
+        return r
+
+    with open(filename, "r") as f:
+        m = [line_to_row(l) for l in f]
+    return m
+
+
+def matrix_to_graph(matrix):
+    n_nodes = len(matrix)
+    graph = []
+    for row_index in range(n_nodes):
+
+        edges = []
+        for col_index in range(n_nodes):
+            weight = matrix[row_index][col_index]
+            if weight != 0:
+                edges.append(Edge(weight, row_index, col_index))
+        graph.append(Vertex(row_index, edges))
+    return graph
+
+
+def total_weight(graph):
+    total_weight = 0
+    for vertex in graph:
+        for edge in vertex.edges:
+            total_weight += edge.weight
+    total_weight //= 2
+    return total_weight
+
+
+def prims_algo(graph):
+    """ Takes a graph and returns the edges that for a minimum spanning tree
+    for that graph. """
+    visited_id = set([0])
+    min_edges = []
+    while len(visited_id) != len(graph):
+        min_edge = None
+        for vertex_id in visited_id:
+            vertex = graph[vertex_id]
+            for edge in vertex.edges:
+                if edge.target_id in visited_id:
+                    pass
+                elif min_edge is None:
+                    min_edge = edge
+                elif edge.weight < min_edge.weight:
+                    min_edge = edge
+        visited_id.add(min_edge.target_id)
+        min_edges.append(min_edge)
+    return min_edges
+
+
+def euler_107():
+    m = load_matrix("../txt/e107.txt")
+    g = matrix_to_graph(m)
+    weight = total_weight(g)
+    new_weight = sum([e.weight for e in prims_algo(g)])
+    return weight - new_weight
+
+
+if __name__ == "__main__":
+    solution = euler_107()
+    print("e107.py: " + str(solution))
+    assert(solution == 259679)
+

python/e108.py

@@ -0,0 +1,43 @@
+from fractions import Fraction
+from lib_misc import proper_divisors
+from math import ceil
+
+
+def get_distinct_solutions(n):
+    """ Get the number of distinct solutions for 1/x + 1/y = 1/n. The idea is
+    that 1/n*2 + 1/n*2 = 1/n is always a solution. So what we do is
+    brute-force all 1/x for x in [2, n * 2]. We have a solution if the
+    difference between 1/n and 1/x has a numerator of 1. """
+    n_inv = Fraction(1, n)
+    n_distinct = 0
+    for x in range(2, n * 2 + 1):
+        x_inv = Fraction(1, x)
+        y_inv = n_inv - x_inv
+        if y_inv.numerator == 1:
+            n_distinct += 1
+    return n_distinct
+
+
+def get_distinct_solutions2(n):
+    ds = proper_divisors(n * n)
+    return ceil((len(ds) + 1) / 2)
+
+
+def euler_108():
+    d_max, n_prev = 0, 0
+    # I arrived at the starting values empirically by observing the deltas.
+    for n in range(1260, 1000000, 420):
+        d = get_distinct_solutions2(n)
+        if d > d_max:
+            # print("n={} d={} delta={}".format(n, d, n - n_prev))
+            n_prev = n
+            d_max = d
+        if d > 1000:
+            return n
+
+
+if __name__ == "__main__":
+    solution = euler_108()
+    print("e108.py: " + str(solution))
+    assert(solution == 180180)
+

python/e109.py

@@ -0,0 +1,73 @@
+
+def all_singles():
+    return list(range(1, 21)) + [25]
+
+
+def all_doubles():
+    return [n * 2 for n in range(1, 21)] + [50]
+
+
+def all_triples():
+    return [n * 3 for n in range(1, 21)]
+
+
+def euler_109():
+    count = 0
+    target_score = 100
+    for first in all_doubles():
+        # D
+        if first < target_score:
+            count += 1
+
+        # D T
+        for second in all_triples():
+            if first + second < target_score:
+                count += 1
+        # D D
+        for second in all_doubles():
+            if first + second < target_score:
+                count += 1
+        # D S
+        for second in all_singles():
+            if first + second < target_score:
+                count += 1
+
+        # D T D
+        for second in all_triples():
+            for third in all_doubles():
+                if first + second + third < target_score:
+                    count += 1
+        # D T S
+        for second in all_triples():
+            for third in all_singles():
+                if first + second + third < target_score:
+                    count += 1
+        # D D S
+        for second in all_doubles():
+            for third in all_singles():
+                if first + second + third < target_score:
+                    count += 1
+
+        # D T T
+        for second in all_triples():
+            for third in all_triples():
+                if second >= third and first + second + third < target_score:
+                    count += 1
+
+        # D D
+        for second in all_doubles():
+            for third in all_doubles():
+                if second >= third and first + second + third < target_score:
+                    count += 1
+        # S S
+        for second in all_singles():
+            for third in all_singles():
+                if second >= third and first + second + third < target_score:
+                    count += 1
+    return count
+
+
+if __name__ == "__main__":
+    solution = euler_109()
+    print("e109.py: " + str(solution))
+    assert(solution == 38182)

python/e110.py

@@ -0,0 +1,75 @@
+from lib_misc import get_item_counts
+from lib_prime import primes
+
+
+def divisors(counts):
+    r = 1
+    for c in counts:
+        r *= (c + 1)
+    return r
+
+
+def tau(factors):
+    orig_factors = factors
+    factors = factors + factors
+    counts = get_item_counts(factors)
+    r = divisors(counts.values())
+    p = 1
+    for f in orig_factors:
+        p *= f
+    return r, p
+
+
+def counters(digits, max_digit):
+    def incrementing_counters(curr, left, max_digit, result):
+        if left == 0:
+            result.append(curr)
+            return
+        start = 1 if not curr else curr[-1]
+        for i in range(start, max_digit + 1):
+            incrementing_counters(curr + [i], left - 1, max_digit, result)
+
+    result = []
+    incrementing_counters([], digits, max_digit, result)
+    return result
+
+
+def euler_110():
+    target = 1000
+    target = 4 * 10**6
+    threshold = (target * 2) - 1
+    psupper = primes(1000)
+
+    lowest_distinct = 0
+    lowest_number = 0
+
+    # find upper bound
+    for i in range(len(psupper)):
+        distinct, number = tau(psupper[:i])
+        if distinct > threshold:
+            # print(lowest_distinct, number)
+            lowest_distinct = distinct
+            lowest_number = number
+            psupper = psupper[:i]
+
+    for j in range(1, len(psupper)):
+        ps = psupper[:-j]
+        for prime_counts in counters(len(ps), 5):
+            prime_counts.reverse()
+            nps = []
+            i = 0
+            for i in range(len(prime_counts)):
+                nps += [ps[i]] * prime_counts[i]
+            nps += ps[i + 1:]
+            distinct, number = tau(nps)
+            if distinct > threshold and distinct < lowest_distinct:
+                lowest_distinct = distinct
+                lowest_number = number
+                # print(lowest_distinct, lowest_number)
+    return lowest_number
+
+
+if __name__ == "__main__":
+    solution = euler_110()
+    print("e110.py: " + str(solution))
+    assert(solution == 9350130049860600)

python/e111.py

@@ -0,0 +1,46 @@
+from lib_prime import is_prime_rabin_miller as is_prime
+
+
+def get_permutations(number: str, number_permutations: int, start_index = 0):
+    """ Returns number permutated at number_permutations locations. """
+    if number_permutations == 0:
+        yield int(number)
+    else:
+        numbers = list("0123456789")
+        base = list(str(number))
+        for i in range(start_index, len(base)):
+            digit_stored = base[i]
+            for n in numbers:
+                if n == digit_stored:
+                    continue
+                base[i] = n
+                if base[0] == "0":
+                    continue
+                for p in get_permutations("".join(base), number_permutations - 1, i):
+                    yield p
+            base[i] = digit_stored
+
+
+def euler_111():
+    result = 0
+    n = 10
+    for d in range(10):
+        subresult = 0
+        base = "".join([str(d) for _ in range(n)])
+        for p in get_permutations(base, 1):
+            if is_prime(p):
+                subresult += p
+        if subresult == 0:
+            # d in [0, 2, 8] requires two permutations to yield at least one prime
+            for p in get_permutations(base, 2):
+                if is_prime(p):
+                    subresult += p
+        assert subresult > 0, "More than two permutations required to yield prime"
+        result += subresult
+    return result
+
+
+if __name__ == "__main__":
+    solution = euler_111()
+    print("e111.py: " + str(solution))
+    assert(solution == 612407567715)

python/e112.py

@@ -0,0 +1,37 @@
+
+def is_increasing(n_str):
+    if len(n_str) < 2:
+        return True
+    if n_str[0] <= n_str[1]:
+        return is_increasing(n_str[1:])
+    return False
+
+
+def is_decreasing(n_str):
+    if len(n_str) < 2:
+        return True
+    if n_str[0] >= n_str[1]:
+        return is_decreasing(n_str[1:])
+    return False
+
+
+def is_bouncy(n_str):
+    return not (is_increasing(n_str) or is_decreasing(n_str))
+
+
+def euler_112():
+    n, n_bouncy = 0, 0
+    while True:
+        n += 1
+        if is_bouncy(str(n)):
+            n_bouncy += 1
+            ratio = n_bouncy / n
+            if ratio >= 0.99:
+                return n
+
+
+if __name__ == "__main__":
+    solution = euler_112()
+    print("e112.py: " + str(solution))
+    assert(solution == 1587000)
+

python/e113.py

@@ -0,0 +1,52 @@
+def is_non_bouncy(n: int):
+    digits = list(map(int, list(str(n))))
+    decreasing = True
+    increasing = True
+    for i in range(len(digits) - 1):
+        increasing = increasing and digits[i] <= digits[i + 1]
+        decreasing = decreasing and digits[i] >= digits[i + 1]
+    return increasing or decreasing
+
+
+def non_bouncy_below_naiv(digits: int):
+    threshold = 10 ** digits
+    non_bouncy = sum([1 for i in range(1, threshold) if is_non_bouncy(i)])
+    return non_bouncy
+
+
+def non_bouncy_below(digits: int):
+    # increasing
+    s = [1 for _ in range(1, 10)]
+    result = sum(s)
+    for _ in range(digits - 1):
+        s = [sum(s[i:]) for i in range(9)]
+        result += sum(s) - 9
+    # The -9 is to avoid double counting 11, 22, 33 for increasing and
+    # decreasing.
+
+    # decreasing
+    decreasing_count = 0
+    s = [i + 1 for i in range(1, 10)]
+    for _ in range(digits - 1):
+        # print(s)
+        decreasing_count += sum(s)
+        s = [sum(s[:i]) + 1 for i in range(1, 10)]
+    result += decreasing_count
+    return result
+
+
+def euler_113():
+    assert non_bouncy_below_naiv(6) == 12951
+    assert non_bouncy_below(6) == 12951
+    assert non_bouncy_below(10) == 277032
+    assert is_non_bouncy(134468)
+    assert is_non_bouncy(66420)
+    assert (not is_non_bouncy(155349))
+    return non_bouncy_below(100)
+
+
+if __name__ == "__main__":
+    solution = euler_113()
+    print("e113.py: " + str(solution))
+    assert(solution == 51161058134250)
+

python/e114.py

@@ -0,0 +1,51 @@
+from functools import lru_cache
+
+
+@lru_cache()
+def block_combinations(row_length: int):
+
+    if row_length == 0:
+        return 1
+
+    min_block_length = 3
+    result = 0
+
+    # Moving from left to right. Let's say we have four spaces (____), then
+    # there are three options
+
+    # Use a blank:
+    #   -___
+    result = block_combinations(row_length - 1)
+         
+    # Or add 3:
+    #   ooo_ 
+    # For this case, the issue is that if we call block_combinations again with
+    # _, then it could happen that we will up with further o and double count.
+    # Therefore, for all cases where the spaces are not all filled, we will add
+    # a space.
+    for new_block in range(min_block_length, row_length):
+        result += block_combinations(row_length - new_block - 1)
+
+    # Or add 4:
+    #   oooo
+    if row_length >= min_block_length:
+        result += 1
+         
+    return result
+
+
+def euler_114():
+    assert block_combinations(0) == 1
+    assert block_combinations(1) == 1
+    assert block_combinations(2) == 1
+    assert block_combinations(3) == 2
+    assert block_combinations(4) == 4
+    assert block_combinations(7) == 17
+    return block_combinations(50)
+
+
+if __name__ == "__main__":
+    solution = euler_114()
+    print("e114.py: " + str(solution))
+    # assert(solution == 0)
+

python/e115.py

@@ -0,0 +1,28 @@
+from functools import lru_cache
+
+
+@lru_cache()
+def block_combinations(row_length: int, min_block_length: int = 3):
+    if row_length == 0:
+        return 1
+    result = 0
+    result = block_combinations(row_length - 1, min_block_length)
+    for new_block in range(min_block_length, row_length):
+        result += block_combinations(row_length - new_block - 1, min_block_length)
+    if row_length >= min_block_length:
+        result += 1
+    return result
+
+
+def euler_115():
+    for n in range(1000):
+        if block_combinations(n, 50) > 10**6:
+            return n
+    return -1
+
+
+if __name__ == "__main__":
+    solution = euler_115()
+    print("e115.py: " + str(solution))
+    assert(solution == 168)
+

python/e116.py

@@ -0,0 +1,31 @@
+from functools import lru_cache
+
+
+@lru_cache()
+def block_combinations(row_length: int, block_length: int = 3):
+    if row_length < 0:
+        return 0
+
+    if row_length == 0:
+        return 1
+
+    result = block_combinations(row_length - 1, block_length)
+    result += block_combinations(row_length - block_length, block_length)
+    return result
+
+
+def solve(row_length: int):
+    r = block_combinations(row_length, 2) - 1
+    r += block_combinations(row_length, 3) - 1
+    r += block_combinations(row_length, 4) - 1
+    return r
+
+
+def euler_116():
+    return solve(50)
+
+
+if __name__ == "__main__":
+    solution = euler_116()
+    print("e116.py: " + str(solution))
+    assert(solution == 20492570929)

python/e117.py

@@ -0,0 +1,30 @@
+from functools import lru_cache
+
+
+@lru_cache()
+def block_combinations(row_length: int):
+    block_sizes = [2, 3, 4]
+    if row_length < 0:
+        return 0
+
+    if row_length == 0:
+        return 1
+
+    # space
+    result = block_combinations(row_length - 1)
+
+    # one of the blocks
+    for block_length in block_sizes:
+        result += block_combinations(row_length - block_length)
+
+    return result
+
+
+def euler_117():
+    return block_combinations(50)
+
+
+if __name__ == "__main__":
+    solution = euler_117()
+    print("e117.py: " + str(solution))
+    assert(solution == 100808458960497)

python/e118.py

@@ -0,0 +1,45 @@
+from itertools import combinations, permutations
+from lib_prime import is_prime_rabin_miller as is_prime
+
+
+def as_number(xs) -> int:
+    return int("".join(map(str, xs)))
+
+
+def list_diff(xs, ys):
+    xs = list(xs)
+    for y in ys:
+        xs.remove(y)
+    return xs
+
+
+def combs(xs, size):
+    return sorted([p
+        for cs in combinations(xs, size)
+        for p in permutations(cs)])
+
+
+def count(min_set_size: int, min_num: int, remaining: list[int]) -> int:
+    if not remaining:
+        return 1
+
+    result = 0
+    # There are no pandigital primes with 9 digits, so we only check till 8.
+    for set_size in range(min_set_size, 9):
+        for subset in combs(remaining, set_size):
+            n = as_number(subset)
+            if n > min_num and is_prime(n):
+                new_remaining = list_diff(remaining, subset)
+                result += count(set_size, n, new_remaining)
+    return result
+
+
+def euler_118():
+    return count(1, 0, [1, 2, 3, 4, 5, 6, 7, 8, 9])
+
+
+if __name__ == "__main__":
+    solution = euler_118()
+    print("e118.py: " + str(solution))
+    assert(solution == 44680)
+

python/e119.py

@@ -0,0 +1,31 @@
+def digit_sum(n: int) -> int:
+    return sum(map(int, str(n)))
+
+
+def is_digital_power_sum(base: int, power: int):
+    if digit_sum(base ** power) == base:
+        return True
+    return False
+
+
+def get_nth_power_sum(n: int):
+    power_sums = []
+    for base in range(2, 600):
+        for power in range(2, 50):
+            if is_digital_power_sum(base, power):
+                power_sums.append(base ** power)
+    power_sums = sorted(power_sums)
+    return power_sums[n - 1]
+
+
+def euler_119():
+    assert(is_digital_power_sum(28, 4) == True)
+    assert(get_nth_power_sum(2) == 512)
+    assert(get_nth_power_sum(10) == 614656)
+    return get_nth_power_sum(30)
+
+
+if __name__ == "__main__":
+    solution = euler_119()
+    print("e119.py: " + str(solution))
+    assert(solution == 248155780267521)

python/e120.py

@@ -0,0 +1,28 @@
+def remainder(a: int, n: int) -> int:
+    d = a * a
+    return (pow(a - 1, n, d) + pow(a + 1, n, d)) % d
+
+
+def max_remainder(a: int) -> int:
+    n, r_max = 1, 2
+    while True:
+        r = remainder(a, n)
+        if r == r_max:
+            break
+        if r > r_max:
+            r_max = r
+        n += 1
+    assert(r_max > 2)
+    return r_max
+
+
+def euler_120():
+    assert(remainder(7, 3) == 42)
+    assert(max_remainder(7) == 42)
+    return sum([max_remainder(n) for n in range(3, 1001)])
+
+
+if __name__ == "__main__":
+    solution = euler_120()
+    print("e120.py: " + str(solution))
+    assert(solution == 333082500)

python/e121.py

@@ -0,0 +1,29 @@
+from math import factorial
+from itertools import combinations
+
+
+def product(xs):
+    r = 1
+    for x in xs:
+        r *= x
+    return r
+
+
+def sub_odds(n_blue, n_total):
+    odds = [n for n in range(1, n_total + 1)]
+    return sum(map(product, combinations(odds, n_total - n_blue)))
+
+
+def euler_121():
+    n_turns = 15
+    n_to_win = n_turns // 2 + 1
+    odds = sum([sub_odds(n_blue, n_turns)
+                for n_blue in range(n_to_win, n_turns + 1)])
+    return int(factorial(n_turns + 1) / odds)
+
+
+if __name__ == "__main__":
+    solution = euler_121()
+    print("e121.py: " + str(solution))
+    assert(solution == 2269)
+

python/e122.py

@@ -0,0 +1,27 @@
+def euler_122():
+    upper = 201
+    m_k = {k: 0 for k in range(2, upper)}
+    sets = [[1]]
+    for _ in range(11):
+        new_sets = []
+        for s in sets:
+            for i in range(len(s)):
+                new_elem = s[i] + s[-1]
+                if new_elem in m_k and m_k[new_elem] == 0:
+                    m_k[new_elem] = len(s)
+                new_sets.append(s + [new_elem])
+
+        # For better performance, we would have to prune here.
+        sets = new_sets
+
+    r = 0
+    for k in range(2, upper):
+        assert m_k[k] != 0
+        r += m_k[k]
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_122()
+    print("e122.py: " + str(solution))
+    assert solution == 1582

python/e123.py

@@ -0,0 +1,40 @@
+from lib_prime import prime_nth
+
+
+def varray(i):
+    return i * 2 + 1
+
+
+def rem(i):
+    n = varray(i)
+    p = prime_nth(n)
+    return (pow(p - 1, n) + pow(p + 1, n)) % (p * p)
+
+
+def bin_search(lo, hi, target):
+    if hi - lo == 1:
+        return varray(hi)
+    i = lo + ((hi - lo) // 2)
+    r = rem(i)
+
+    # print(f"{i=:<6} {r=}")
+    if r < target:
+        return bin_search(i, hi, target)
+    elif r > target:
+        return bin_search(lo, i, target)
+    assert False
+
+
+def euler_123():
+    """
+    I found out that when n is even, the remainder is always 2. When the
+    remainder is odd, the series grows monotonically. That means we can use a
+    binary search on odd values to find the solution.
+    """
+    return bin_search(1000, 30000, 10**10)
+
+
+if __name__ == "__main__":
+    solution = euler_123()
+    print("e123.py: " + str(solution))
+    assert(solution == 21035)

python/e124.py

@@ -0,0 +1,24 @@
+from lib_prime import prime_factors
+
+
+def radical(n: int) -> int:
+    fs = prime_factors(n)
+    r = 1
+    for f in set(fs):
+        r *= f
+    return r
+
+
+def euler_124():
+    assert radical(504) == 42
+    xs = []
+    for n in range(1, 100001):
+        xs.append((radical(n), n))
+    xs = sorted(xs)
+    return xs[10000 - 1][1]
+
+
+if __name__ == "__main__":
+    solution = euler_124()
+    print("e124.py: " + str(solution))
+    assert(solution == 21417)

python/e125.py

@@ -0,0 +1,23 @@
+from lib_misc import is_palindrome
+
+
+def euler_125():
+    upper = 10**8
+    ns = set()
+    for i in range(1, upper):
+        i2 = i * i
+        if i2 >= upper:
+            break
+        for j in range(i + 1, upper):
+            i2 += (j * j)
+            if i2 >= upper:
+                break
+            if is_palindrome(i2):
+                ns.add(i2)
+    return sum(ns)
+
+
+if __name__ == "__main__":
+    solution = euler_125()
+    print("e125.py: " + str(solution))
+    assert(solution == 2906969179)

python/e126.py

@@ -0,0 +1,141 @@
+from functools import lru_cache
+from collections import defaultdict
+
+
+CUBE_THRESHOLD = 10000
+
+def base_cubes(xlen, ylen, zlen) -> set:
+    r = set()
+    for x in range(xlen):
+        for y in range(ylen):
+            for z in range(zlen):
+                r.add((x, y, z))
+    return r
+
+
+@lru_cache
+def neighbors(x, y, z) -> list:
+    dirs = [
+        (1, 0, 0),
+        (-1, 0, 0),
+        (0, 1, 0),
+        (0, -1, 0),
+        (0, 0, 1),
+        (0, 0, -1),
+    ]
+    return [(x + dx, y + dy, z + dz) for dx, dy, dz in dirs]
+
+
+@lru_cache
+def count_cubes_per_layer(x, y, z) -> list:
+    """ First iteration. Naively add outer layers. """
+    r = []
+    clayer = base_cubes(x, y, z)
+    players = set(clayer)
+    while True:
+        nlayer = set()
+        for c in clayer:
+            for nb in neighbors(*c):
+                if nb not in players:
+                    nlayer.add(nb)
+        clayer = nlayer
+        players |= clayer
+        layer_count = len(nlayer)
+        r.append(layer_count)
+        if layer_count > CUBE_THRESHOLD:
+            break
+    return r
+
+
+def count_cubes_per_layer2(x, y, z) -> list:
+    """ Second iteration. Realized that the delta between the layers grows by
+    the same value. """
+    r = []
+    clayer = base_cubes(x, y, z)
+    players = set(clayer)
+    layer_0 = 0
+    layer_1 = 0
+    for i in range(2):
+        nlayer = set()
+        for c in clayer:
+            for nb in neighbors(*c):
+                if nb not in players:
+                    nlayer.add(nb)
+        clayer = nlayer
+        players |= clayer
+        layer_count = len(nlayer)
+        if i == 0:
+            layer_0 = layer_count
+        elif i == 1:
+            layer_1 = layer_count
+        r.append(layer_count)
+        if layer_count > CUBE_THRESHOLD:
+            break
+
+    layer = r[-1]
+    layer_delta = (layer_1 - layer_0) + 8
+    while layer <= CUBE_THRESHOLD:
+        layer += layer_delta
+        r.append(layer)
+        layer_delta += 8
+    return r
+
+
+def count_cubes_per_layer3(x, y, z) -> list:
+    """ Third iteration. Calculate first two layers directly, and then use
+    constant delta approach from second iteration. """
+    r = []
+    layer_0 = 2 * (x * y + x * z + y * z)
+    layer_1 = layer_0 + 4 * (x + y + z)
+    r = [layer_0, layer_1]
+    layer = r[-1]
+    layer_delta = (layer_1 - layer_0) + 8
+    while layer <= CUBE_THRESHOLD:
+        layer += layer_delta
+        r.append(layer)
+        layer_delta += 8
+    return r
+
+
+def euler_126():
+    # We assume that the least n for which C(n) is equal to `target`
+    # is under `layer_min_threshold`. We got this my trial and error.
+    layer_min_threshold = 20000
+    target = 1000
+
+    layers = defaultdict(int)
+    for x in range(1, layer_min_threshold):
+        for y in range(1, x + 1):
+
+            # If we exceed the threshold, cut the loop short.
+            if 2 * x * y  > layer_min_threshold:
+                break
+
+            for z in range(1, y + 1):
+                layer_0 = 2 * (x * y + x * z + y * z)
+                layer_1 = layer_0 + 4 * (x + y + z)
+
+                layers[layer_0] += 1
+                layers[layer_1] += 1
+
+                layer = layer_1
+                layer_delta = (layer_1 - layer_0) + 8
+
+                # If we exceed the threshold, cut the loop short.
+                while layer < layer_min_threshold:
+                    layer += layer_delta
+                    layer_delta += 8
+                    layers[layer] += 1
+
+    # Find the actual least value n for which C(n) == target.
+    layermin = 10**9
+    for n in layers.keys():
+        if layers[n] == target and n < layermin:
+            layermin = n
+    return layermin
+
+
+if __name__ == "__main__":
+    solution = euler_126()
+    print("e126.py: " + str(solution))
+    assert(solution == 18522)

python/e127.py

@@ -0,0 +1,67 @@
+from functools import lru_cache
+from lib_prime import prime_factors
+
+
+@lru_cache(maxsize=200000)
+def unique_factors(n: int) -> set[int]:
+    return set(prime_factors(n))
+
+
+def euler_127():
+    c_max = 120000
+    r = 0
+
+    # for a in [5]:
+    for a in range(1, c_max):
+        fsa = unique_factors(a)
+        # for b in [27]:
+        for b in range(a + 1, c_max):
+            if a + b > c_max:
+                break
+            c = a + b
+            rad = 1
+
+            do_continue = False
+            for fa in fsa:
+                rad *= fa
+                if b % fa == 0:
+                    do_continue = True
+                    break
+                if c % fa == 0:
+                    do_continue = True
+                    break
+            if do_continue:
+                continue
+            if rad > c:
+                continue
+
+            do_continue = False
+            fsb = unique_factors(b)
+            for fb in fsb:
+                rad *= fb
+                if c % fb == 0:
+                    do_continue = True
+                    break
+            if do_continue:
+                continue
+
+            if rad >= c:
+                continue
+
+            fsc = unique_factors(c)
+            for fc in fsc:
+                rad *= fc
+            if rad >= c:
+                continue
+
+            # print(fsa, fsb, fsc)
+            # print(a, b, c)
+            r += c
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_127()
+    print("e127.py: " + str(solution))
+    assert(solution == 18407904)
+

python/e128.py

@@ -0,0 +1,211 @@
+from lib_prime import primes, is_prime_rabin_miller
+from functools import lru_cache
+
+
+def first_approach():
+    """ This was my original naiv approach that worked, but was way to slow
+    because we computed the whole hex tile grid. """
+    ps = set(primes(1000000))
+    hextiles = {(0, 0): 1}
+    add_ring(1, hextiles)
+    add_ring(2, hextiles)
+
+    # Visualization of how the coord system for the hex tile grid works.
+    #                 
+    #     2 1 0 1 2  
+    #  4      8          
+    #  3    9  19        
+    #  2 10   2  18      
+    #  1    3   7        
+    #  0 11   1  17      
+    #  1    4   6        
+    #  2 12   5  16      
+    #  3   13  15        
+    #  4     14          
+    #                 
+    #                 
+
+    pds = [1, 2, 8]
+    for n in range(0, 100):
+        ring_coord = ring_start(n)
+        print("  ", ring_coord, hextiles[ring_coord])
+        if n % 10 == 0:
+            print(n)
+        add_ring(n + 1, hextiles)
+        for coord in ring_coords(n):
+            pdv = pd(coord, hextiles, ps)
+            if pdv == 3:
+                v = hextiles[coord]
+                assert v == first_number_ring(n) or v == first_number_ring(n + 1) - 1
+                pds.append(hextiles[coord])
+    print(len(pds))
+    print(sorted(pds))
+    target = 10
+    if len(pds) >= target:
+        return sorted(pds)[target - 1]
+    return None
+
+
+def ring_start(n):
+    return (-n * 2, 0)
+
+
+@lru_cache
+def first_number_ring(n):
+    if n == 0:
+        return 1
+    elif n == 1:
+        return 2
+    else:
+        return first_number_ring(n - 1) + (n - 1) * 6
+
+
+def get_nbvs_ring_start(n):
+    #    8        1
+    #  9   19   2   6   
+    #    2        0  
+    #  3   7    3   5
+    #    l        4 
+    v0 = first_number_ring(n)
+    v1 = first_number_ring(n + 1)
+    v2 = v1 + 1
+    v3 = v0 + 1
+    v4 = first_number_ring(n - 1)
+    v5 = v1 - 1
+    v6 = first_number_ring(n + 2) - 1
+    nbvs = [v1, v2, v3, v4, v5, v6]
+    # print(v0, nbvs)
+    return nbvs
+
+
+def pd_ring_start(n):
+    v0 = first_number_ring(n)
+    nbvs = get_nbvs_ring_start(n)
+    r = 0
+    for v in nbvs:
+        if is_prime_rabin_miller(abs(v0 - v)):
+            r += 1
+    return r
+
+
+def get_nbvs_ring_prev(n):
+    """ Get neighbors and values for tile before top tile. """
+    v0 = first_number_ring(n + 1) - 1
+    v1 = first_number_ring(n + 2) - 1
+    v2 = first_number_ring(n)
+    v3 = first_number_ring(n - 1)
+    v4 = v2 - 1
+    v5 = first_number_ring(n + 1) - 2
+    v6 = v1 - 1
+    nbvs = [v1, v2, v3, v4, v5, v6]
+    # print(v0, nbvs)
+    return nbvs
+
+
+def pd_ring_prev(n):
+    v0 = first_number_ring(n + 1) - 1
+    nbvs = get_nbvs_ring_prev(n)
+    r = 0
+    for v in nbvs:
+        if is_prime_rabin_miller(abs(v0 - v)):
+            r += 1
+    return r
+
+
+def add_ring(n, numbers):
+    """ Adds ring n coords and values to numbers. """
+    if n == 0:
+        return
+    first_coord = ring_start(n)
+    current_number = first_number_ring(n)
+    numbers[first_coord] = current_number
+    current_coord = tuple(first_coord)
+
+    for ro, co in [(1, -1), (2, 0), (1, 1), (-1, 1), (-2, 0), (-1, -1)]:
+        for _ in range(n):
+            current_coord = (current_coord[0] + ro, current_coord[1] + co)
+            current_number += 1
+            numbers[current_coord] = current_number
+
+    # Reset first coord which is overriden.
+    numbers[first_coord] = first_number_ring(n)
+
+
+def get_neighbor_values(coord, numbers):
+    neighbors = []
+
+    for ro, co in [(-2, 0), (-1, 1), (1, 1), (2, 0), (1, -1), (-1, -1)]:
+        nc = (coord[0] + ro, coord[1] + co)
+        neighbors.append(numbers[nc])
+    return neighbors
+
+
+def ring_coords(n):
+    """ Returns coords for ring n. """
+    if n == 0:
+        yield (0, 0)
+    current_coord = ring_start(n)
+    for ro, co in [(1, -1), (2, 0), (1, 1), (-1, 1), (-2, 0), (-1, -1)]:
+        for _ in range(n):
+            current_coord = (current_coord[0] + ro, current_coord[1] + co)
+            yield current_coord
+
+
+def pd(coord, numbers, primes):
+    prime_delta = 0
+    v = numbers[coord] 
+    for nbv in get_neighbor_values(coord, numbers):
+        if abs(v - nbv) in primes:
+            prime_delta += 1
+    return prime_delta
+
+
+def test_get_nbvs_functions():
+    """ Make sure that the function to compute the neighbor values works
+    correctly by comparing the values to the original grid based approach. """
+    hextiles = {(0, 0): 1}
+    add_ring(1, hextiles)
+    add_ring(2, hextiles)
+
+    for n in range(2, 30):
+        add_ring(n + 1, hextiles)
+        nbvs1 = get_nbvs_ring_start(n)
+        ring_coord = ring_start(n)
+        nbvs2 = get_neighbor_values(ring_coord, hextiles)
+        assert sorted(nbvs1) == sorted(nbvs2)
+
+        nbvs1 = get_nbvs_ring_prev(n)
+        ring_coord = (ring_start(n)[0] + 1, ring_start(n)[1] + 1)
+        nbvs2 = get_neighbor_values(ring_coord, hextiles)
+        assert sorted(nbvs1) == sorted(nbvs2)
+
+
+def euler_128():
+    # Conjecture: PD3s only occur at ring starts or at ring start minus one.
+    # Under this assumption, we can significantly reduce the search space.
+    # The only challenge is to find a way to directly compute the relevant
+    # coords and neighbor values.
+
+    test_get_nbvs_functions()
+
+    target = 2000
+    pds = [1, 2]
+    for n in range(2, 80000):
+        if pd_ring_start(n) == 3:
+            v0 = first_number_ring(n)
+            pds.append(v0)
+
+        if pd_ring_prev(n) == 3:
+            v0 = first_number_ring(n + 1) - 1
+            pds.append(v0)
+
+    assert len(pds) > target
+    return sorted(pds)[target - 1]
+
+
+if __name__ == "__main__":
+    solution = euler_128()
+    print("e128.py: " + str(solution))
+    assert(solution == 14516824220)
+
+

python/e129.py

@@ -0,0 +1,58 @@
+from math import gcd
+from functools import lru_cache
+
+
+@lru_cache
+def r(n):
+    assert n > 0
+    if n == 1:
+        return 1
+    else:
+        return 10**(n - 1) + r(n - 1)
+
+
+def r_closed_form(n):
+    assert n > 0
+    return (10**n - 1) // 9
+
+
+def r_modulo_closed_form(n, m):
+    assert n > 0 and m > 0
+    return ((pow(10, n, 9 * m) - 1) // 9) % m
+
+
+def a(n):
+    assert gcd(n, 10) == 1
+    k = 1
+    while True:
+        if r_modulo_closed_form(k, n) == 0:
+            return k
+        k += 1
+
+
+def euler_129():
+
+    # Comparing naiv to efficient implementations to find potential bugs.
+    for n in range(5, 1000):
+        assert r(n) == r_closed_form(n)
+        for m in range(2, 20):
+            assert (r_closed_form(n) % m)  == r_modulo_closed_form(n, m)
+
+    assert a(7) == 6
+    assert a(41) == 5
+    assert a(17) == 16
+
+    target = 10**6
+    delta = 200
+    for n in range(target, target + delta):
+        if gcd(n, 10) == 1:
+            if (av := a(n)):
+                if av > target:
+                    return n
+    return None
+
+
+if __name__ == "__main__":
+    solution = euler_129()
+    print("e129.py: " + str(solution))
+    assert(solution == 1000023)

python/e130.py

@@ -0,0 +1,32 @@
+from math import gcd
+from lib_prime import is_prime_rabin_miller
+
+
+def r_modulo_closed_form(n, m):
+    assert n > 0 and m > 0
+    return ((pow(10, n, 9 * m) - 1) // 9) % m
+
+
+def a_special(n):
+    k = n - 1
+    if r_modulo_closed_form(k, n) == 0:
+        return k
+    return None
+
+
+def euler_130():
+    c, r = 0, 0
+    for n in range(2, 10**9):
+        if gcd(n, 10) == 1 and not is_prime_rabin_miller(n):
+            if a_special(n) is not None:
+                r += n
+                c += 1
+            if c == 25:
+                break
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_130()
+    print("e130.py: " + str(solution))
+    assert(solution == 149253)

python/e131.py

@@ -0,0 +1,29 @@
+from lib_prime import primes
+
+
+def is_cube(n):
+    if n == 0:
+        return False
+    return round(n ** (1/3)) ** 3 == n
+
+
+def euler_131():
+    r = 0
+    ps = primes(10**6)
+    minm = 1
+    for p in ps:
+        for m in range(minm, minm + 20):
+            n = m * m * m
+            x = n * n * n + n * n * p
+            if is_cube(x):
+                # print(p, n)
+                r += 1
+                minm = m
+                break
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_131()
+    print("e131.py: " + str(solution))
+    assert(solution == 173)

python/e132.py

@@ -0,0 +1,32 @@
+from lib_prime import primes
+
+
+def r_modulo_closed_form(n, m):
+    assert n > 0 and m > 0
+    return ((pow(10, n, 9 * m) - 1) // 9) % m
+
+
+def is_factor(n, m):
+    return pow(10, n, 9 * m) == 1
+
+
+def repunit_factor_sum(n):
+    ps = []
+    for p in primes(10**6):
+        assert is_factor(n, p) == (r_modulo_closed_form(n, p) == 0)
+        if r_modulo_closed_form(n, p) == 0:
+            ps.append(p)
+        if len(ps) == 40:
+            break
+    return sum(ps)
+
+
+def euler_132():
+    assert repunit_factor_sum(10) == 9414
+    return repunit_factor_sum(10**9)
+
+if __name__ == "__main__":
+    solution = euler_132()
+    print("e132.py: " + str(solution))
+    assert solution == 843296
+

python/e133.py

@@ -0,0 +1,23 @@
+from lib_prime import primes
+
+
+def r_modulo_closed_form(n, m):
+    assert n > 0 and m > 0
+    return ((pow(10, n, 9 * m) - 1) // 9) % m
+
+
+def euler_132():
+    n = 10**20
+    r = 0
+    for p in primes(100_000):
+        if r_modulo_closed_form(n, p) == 0:
+            pass
+        else:
+            r += p
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_132()
+    print("e132.py: " + str(solution))
+    assert solution == 453647705

python/e134.py

@@ -0,0 +1,72 @@
+from lib_prime import primes
+from math import log, ceil
+
+
+def s(p1, p2):
+    p1l = ceil(log(p1, 10))
+    base = 10**p1l
+    for lhs in range(base, 10**12, base):
+        r = lhs + p1
+        if r % p2 == 0:
+            return r
+    assert False
+
+
+def s2(p1, p2):
+    # Invert, always invert... but still slow, haha
+    p1l = ceil(log(p1, 10))
+    base = 10**p1l
+    c = p2
+    while True:
+        if c % base == p1:
+            return c
+        c += p2
+
+
+def ext_gcd(a, b):
+    if a == 0:
+        return (b, 0, 1)
+    else:
+        gcd, x, y = ext_gcd(b % a, a)
+        return (gcd, y - (b // a) * x, x)
+
+
+def modinv(a, b):
+    gcd, x, _ = ext_gcd(a, b)
+    if gcd != 1:
+        raise Exception('Modular inverse does not exist')
+    else:
+        return x % b
+
+
+def s3(p1, p2):
+    # Okay with some math we are fast.
+    d = 10**ceil(log(p1, 10))
+    dinv = modinv(d, p2)
+    k = ((p2 - p1) * dinv) % p2
+    r = k * d + p1
+    return r
+
+
+def euler_134():
+    ps = primes(10000)
+    for i in range(2, len(ps) - 1):
+        p1, p2 = ps[i], ps[i + 1]
+        sc = s(p1, p2)
+        sc2 = s2(p1, p2)
+        sc3 = s3(p1, p2)
+        assert sc == sc2 and sc2 == sc3
+
+    r = 0
+    ps = primes(10**6 + 10) # p1 < 10**6 but p2 is the first p > 10**6 
+    for i in range(2, len(ps) - 1):
+        p1, p2 = ps[i], ps[i + 1]
+        sc = s3(p1, p2)
+        r += sc
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_134()
+    print("e134.py: " + str(solution))
+    assert solution == 18613426663617118

python/e135.py

@@ -0,0 +1,37 @@
+def euler_135():
+    n = 10**7
+    threshold = 10**6
+    lookup = {i: 0 for i in range(1, threshold)}
+    firstpositivestart = 1
+    for x in range(1, n):
+        firstpositive = False
+        for dxy in range(firstpositivestart, x // 2 + 1):
+            z = x - 2 * dxy
+            if z <= 0:
+                break
+            y = x - dxy
+            xyz = x*x - y*y - z*z
+            if xyz <= 0:
+                continue
+
+            # Start when xyz is already greater than 0.
+            if not firstpositive:
+                firstpositivestart = dxy
+                firstpositive = True
+
+            # If xyz is greater than threshold there are no more solutions.
+            if xyz >= threshold:
+                break
+            lookup[xyz] += 1
+
+    r = 0
+    for _, v in lookup.items():
+        if v == 10:
+            r += 1
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_135()
+    print("e135.py: " + str(solution))
+    assert solution == 4989

python/e136.py

@@ -0,0 +1,37 @@
+def euler_136():
+    n = 10**8
+    threshold = 50*10**6
+    lookup = [0 for _ in range(0, threshold)]
+    firstpositivestart = 1
+    for x in range(1, n):
+        firstpositive = False
+        for dxy in range(firstpositivestart, x // 2 + 1):
+            z = x - 2 * dxy
+            if z <= 0:
+                break
+            y = x - dxy
+            xyz = x*x - y*y - z*z
+            if xyz <= 0:
+                continue
+
+            # Start when xyz is already greater than 0.
+            if not firstpositive:
+                firstpositivestart = dxy
+                firstpositive = True
+
+            # If xyz is greater than threshold there are no more solutions.
+            if xyz >= threshold:
+                break
+            lookup[xyz] += 1
+
+    r = 0
+    for v in lookup:
+        if v == 1:
+            r += 1
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_136()
+    print("e136.py: " + str(solution))
+    assert solution == 2544559

python/e138.py

@@ -0,0 +1,32 @@
+from math import gcd
+
+
+def euler_138():
+    count = 12
+    r, m, n = 0, 2, 1
+    while True:
+        if count == 0:
+            break
+
+        if (m - n) % 2 == 1 and gcd(m, n) == 1:  # Ensure m and n are coprime and not both odd
+            a = m * m - n * n
+            b = 2 * m * n
+            c = m * m + n * n
+            if a > b:
+                a, b = b, a
+            assert a * a + b * b == c * c
+            if a * 2 - 1 == b or a * 2 + 1 == b:
+                # print(f"{m=} {n=} l={c}")
+                r += c # c is L
+                count -= 1
+                n = m
+                m *= 4
+        m += 1
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_138()
+    print("e138.py: " + str(solution))
+    assert solution == 1118049290473932
+

python/e139.py

@@ -0,0 +1,35 @@
+from math import gcd
+
+
+def euler_139():
+    limit = 100_000_000
+    r, m, n = 0, 2, 1
+    while True:
+        a = m * m - 1 * 1
+        b = 2 * m * 1
+        c = m * m + 1 * 1
+        if a + b + c > limit:
+            return r
+
+        for n in range(1, m):
+            if (m - n) % 2 == 1 and gcd(m, n) == 1:  # Ensure m and n are coprime and not both odd
+                a = m * m - n * n
+                b = 2 * m * n
+                c = m * m + n * n
+                if a > b:
+                    a, b = b, a
+
+                ka, kb, kc = a, b, c 
+                while ka + kb + kc < limit:
+                    assert ka * ka + kb * kb == kc * kc
+                    if kc % (kb - ka) == 0:
+                        r += 1
+                    ka, kb, kc = ka + a, kb + b, kc + c
+        m += 1
+
+
+if __name__ == "__main__":
+    solution = euler_139()
+    print("e139.py: " + str(solution))
+    assert solution == 10057761
+

python/e141.py

@@ -0,0 +1,80 @@
+import concurrent.futures
+from math import isqrt, gcd
+
+
+def issqrt(n):
+    isq = isqrt(n)
+    return isq * isq == n
+
+
+def have_same_ratio(a, b, c, d):
+    if b == 0 or d == 0:
+        raise ValueError("Denominators must be non-zero")
+    return a * d == b * c
+
+
+def check_range(start, end):
+    local_sum = 0
+    for i in range(start, end):
+        if i % 10000 == 0:
+            print(i)
+        n = i * i
+        for d in range(1, i):
+            q = n // d
+            r = n % d
+            if r == 0:
+                continue
+            assert r < d and d < q
+            if have_same_ratio(d, r, q, d):
+                print(f"{n=} {i=} {r=} {d=} {q=}")
+                local_sum += n
+                break
+    return local_sum
+
+
+def euler_141_brute_force():
+    s = 0
+    m = 1_000_000
+    range_size = 10_000
+
+    with concurrent.futures.ProcessPoolExecutor() as executor:
+        futures = []
+        for start in range(1, m, range_size):
+            end = min(start + range_size, m)
+            futures.append(executor.submit(check_range, start, end))
+        for future in concurrent.futures.as_completed(futures):
+            s += future.result()
+    return s
+
+
+def euler_141():
+    r = 0
+    NMAX = 10**12
+    a = 1
+    while True:
+        if a**3 > NMAX:
+            break
+        for b in range(1, a):
+            if a**3 * b**2 > NMAX:
+                break
+
+            if gcd(a, b) != 1:
+                continue
+
+            c = 1
+            while True:
+                n = a**3 * b * c**2 + c * b**2
+                if n > NMAX:
+                    break
+                if issqrt(n):
+                    # print(n)
+                    r += n
+                c += 1
+        a += 1
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_141()
+    print("e141.py: " + str(solution))
+    assert solution == 878454337159

python/e142.py

@@ -0,0 +1,38 @@
+from math import isqrt
+
+
+def iq(n):
+    isq = isqrt(n)
+    return isq * isq == n
+
+
+def euler_142():
+    # x + y x - y x + z x - z y + z y - z
+    NMAX = 1_000_000
+    squares = [n * n for n in range(1, NMAX)]
+    for x in range(1, NMAX):
+        for s1 in squares:
+            y = x - s1 # s1 = x - y
+            if y <= 0:
+                break
+            if not iq(x + y): # s2 = x + y
+                continue
+            for s3 in squares:
+                z = y - s3 # s3 = y - z
+                if x - z <= 0:
+                    continue
+                if z <= 0:
+                    break
+                if not iq(y + z): # s4 = y + z
+                    continue
+                if iq(x + z) and iq(x - z):
+                    print(x, y, z, x + y + z)
+                    return x + y + z
+    return None
+
+
+if __name__ == "__main__":
+    solution = euler_142()
+    print("e142.py: " + str(solution))
+    assert solution == 1006193
+

python/e143.py

@@ -0,0 +1,127 @@
+import math
+import concurrent.futures
+from collections import defaultdict
+
+
+def attempt_1(a, b, c):
+    """ For the triangle with the sites a, b, c, calculate
+        the coordinates of C where we define A to be at (0, 0),
+        and B at (c, 0). """
+    assert a + b > c
+    assert a + c > b
+    assert b + c > a
+    S2 = [s**2 for s in range(120000)]
+    S4 = [s**4 for s in range(120000)]
+
+    a2 = S2[a]
+    b2 = S2[b]
+    c2 = S2[c]
+
+    a4 = S4[a]
+    b4 = S4[b]
+    c4 = S4[c]
+
+    # I've derived this myself but unfortunately I was completely on the wrong
+    # track.
+    x = -3 * (a4 + b4 + c4) + 6 * (a2*b2 + a2*c2 + b2*c2)
+    y = math.isqrt(x)
+    # y = math.sqrt(-3 * (a4 + b4 + c4) + 6 * (a2*b2 + a2*c2 + b2*c2))
+    if not y * y == x:
+        return "no int"
+
+    d = math.isqrt((c2 + b2 + a2 + y) // 2)
+    if d * d * 2 == c2 + b2 + a2 + y:
+        return d
+    else:
+        return "no int"
+
+
+def is_square(n):
+    i = math.isqrt(n)
+    return i * i == n
+
+
+def check_range(start, end, limit):
+    # When you realize that the angles in the middle are all 120 degress, you
+    # can derive the following formula with the law of cosines:
+    # c**2 = p**2 + p * r + r**2
+    rs = set()
+    for q in range(start, end):
+        for r in range(1, q):
+            if q + r > limit:
+                break
+            a2 = q**2 + q * r + r**2
+            if not is_square(a2):
+                continue
+            for p in range(1, r):
+                if q + r + p > limit:
+                    break
+                b2 = p**2 + p * q + q**2
+                if not is_square(b2):
+                    continue
+
+                c2 = p**2 + p * r + r**2
+                if is_square(c2):
+                    d = p + q + r
+                    # assert 3*(a**4 + b**4 + c**4 + d**4) == (a**2 + b**2 + c**2 + d**2)**2
+                    print(math.sqrt(a2), math.sqrt(b2), math.sqrt(c2), d)
+                    if d < limit:
+                        rs.add(d)
+    return rs
+
+
+def euler_143_brute_foce():
+    m = 120_000
+    range_size = 1000
+    s = set()
+    with concurrent.futures.ProcessPoolExecutor() as executor:
+        futures = []
+        for start in range(1, m, range_size):
+            end = min(start + range_size, m)
+            futures.append(executor.submit(check_range, start, end, m))
+        for future in concurrent.futures.as_completed(futures):
+            s |= future.result()
+    return sum(s)
+
+
+def euler_143():
+    pairs = defaultdict(set)
+    limit = 120_000
+
+    # Instead of bruteforcing, we can calculate 120 degree triangles directly.
+    #
+    # https://en.wikipedia.org/wiki/Integer_triangle#Integer_triangles_with_a_120%C2%B0_angle
+    #
+    # We then map one 120 angle adjacent site to another one.
+    for m in range(1, math.isqrt(limit)):
+        for n in range(1, m):
+            a = m**2 + m * n + n**2
+            b = 2*m*n + n**2
+            c = m**2 - n**2
+            assert a**2 == b**2 + b * c + c**2
+
+            k = 1
+            while k * (b + c) < limit:
+                pairs[k * b].add(k * c)
+                pairs[k * c].add(k * b)
+                k += 1
+
+    # Which ultimately allows us to construct all the triangles by iterating
+    # over the sides and looking up the respective next side.
+    xs = set()
+    for q in pairs:
+        for r in pairs[q]:
+            for p in pairs[r]:
+                if q in pairs[p]:
+                    s = q + r + p
+                    if s < limit:
+                        xs.add(s)
+
+    return sum(xs)
+
+
+if __name__ == "__main__":
+    solution = euler_143()
+    print("e143.py: " + str(solution))
+    assert solution == 30758397
+

python/e144.py

@@ -0,0 +1,129 @@
+import math
+from dataclasses import dataclass
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+@dataclass
+class P:
+    x: float
+    y: float
+
+
+@dataclass
+class Line:
+    m: float
+    n: float
+
+    def plot(self, x_min: float, x_max: float, color='r'):
+        x = np.linspace(x_min, x_max, 10**5)
+        y = self.m * x + self.n
+        plt.plot(x, y, color, label='') 
+
+
+def line_from_two_points(p1: P, p2: P) -> Line:
+    m = (p2.y - p1.y) / (p2.x - p1.x)
+    n = ((p1.y * p2.x) - (p2.y * p1.x)) / (p2.x - p1.x)
+    return Line(m, n)
+
+
+def line_from_one_point(p: P, m: float) -> Line:
+    n = -m * p.x + p.y
+    return Line(m, n)
+
+
+def reflected_slope(m1, m2):
+    """ Math """
+    alpha1 = math.atan(m1)
+    alpha2 = -math.atan(1/m2)
+    reflected_angle = alpha1 - 2*alpha2
+    m3 = math.tan(reflected_angle)
+    return -m3
+
+
+def get_next_point(start_point: P, line: Line) -> P:
+    # With y = m*x + n into 4x**2 + y**2 = 10 get
+    # quadratic equation x**2 + 2mn * x + (n**2 - 10) = 0.
+
+    # Solve quadratic equation
+    a = 4 + line.m ** 2
+    b = 2 * line.m * line.n
+    c = line.n**2 - 100
+    s = math.sqrt(b**2 - 4*a*c)
+
+    # Pick correct x 
+    EPSILON = 10e-9
+    x1 = (-b + s) / (2*a)
+    x2 = (-b - s) / (2*a)
+    x = x2 if abs(x1 - start_point.x) < EPSILON else x1
+
+    # Calculate corresponding y and return point
+    y = x * line.m + line.n
+    return P(x, y)
+
+
+def calc_slope_ellipse(p: P) -> float:
+    """ Given by Project Euler. """
+    return -4 * p.x / p.y
+
+
+def plot_first_ten():
+    x = np.linspace(-10, 10, 10**5)
+    y_positive = np.sqrt(100 - 4*x**2)
+    y_negative = -np.sqrt(100 - 4*x**2)
+
+    plt.figure(figsize=(6,6)) 
+    plt.plot(x, y_positive, 'b', label='y=sqrt(10-4x^2)') 
+    plt.plot(x, y_negative, 'b', label='y=-sqrt(10-4x^2)') 
+    plt.xlabel('x') 
+    plt.ylabel('y') 
+    plt.title('Plot of the ellipse 4x^2 + y^2 = 10') 
+    plt.axis('equal') 
+
+    current_point = P(0.0, 10.1)
+    next_point = P(1.4, -9.6)
+    line = line_from_two_points(current_point, next_point)
+    next_point_computed = get_next_point(current_point, line)
+    line.plot(0, 1.4)
+    assert(next_point == next_point_computed)
+    current_point = next_point
+
+    for _ in range(10):
+        m1 = line.m
+        m2 = calc_slope_ellipse(current_point)
+        m3 = reflected_slope(m1, m2)
+        line = line_from_one_point(current_point, m3)
+        next_point = get_next_point(current_point, line)
+        x_min = min(current_point.x, next_point.x)
+        x_max = max(current_point.x, next_point.x)
+        line.plot(x_min, x_max, 'g')
+        current_point = next_point
+
+    plt.legend()
+    plt.show()
+
+
+def euler_144():
+    plot_first_ten()
+    current_point = P(0.0, 10.1)
+    next_point = P(1.4, -9.6)
+    line = line_from_two_points(current_point, next_point)
+    current_point = next_point
+    counter = 0
+    while True:
+        counter += 1
+        m1 = line.m
+        m2 = calc_slope_ellipse(current_point)
+        m3 = reflected_slope(m1, m2)
+        line = line_from_one_point(current_point, m3)
+        next_point = get_next_point(current_point, line)
+        current_point = next_point
+        if -0.01 <= current_point.x and current_point.x <= 0.01 and current_point.y > 0:
+            break
+    return counter
+
+
+if __name__ == "__main__":
+    solution = euler_144()
+    print("e144.py: " + str(solution))
+    assert(solution == 354)

python/e146.py

@@ -0,0 +1,49 @@
+from lib_prime import is_prime_rabin_miller
+import concurrent.futures
+
+
+def follows_prime_pattern(n):
+    n2 = n * n
+    prev_prime = None
+    for i in [1, 3, 7, 9, 13, 27]:
+        p = n2 + i
+        if not is_prime_rabin_miller(p):
+            return False
+        if prev_prime is not None:
+            for x in range(prev_prime + 2, p, 2):
+                if is_prime_rabin_miller(x):
+                    return False
+        prev_prime = p
+    return True
+
+
+def check_range(n_min, n_max):
+    s = 0
+    for n in range(n_min - (n_min % 10), n_max, 10):
+        if n % 3 == 0:
+            continue
+        if follows_prime_pattern(n):
+            print(n)
+            s += n
+    return s
+
+
+def euler_146():
+    s = 0
+    m = 150_000_000
+    range_size = 500_000
+
+    with concurrent.futures.ProcessPoolExecutor() as executor:
+        futures = []
+        for start in range(1, m, range_size):
+            end = min(start + range_size, m)
+            futures.append(executor.submit(check_range, start, end))
+        for future in concurrent.futures.as_completed(futures):
+            s += future.result()
+    return s
+
+
+if __name__ == "__main__":
+    solution = euler_146()
+    print("e146.py: " + str(solution))
+    assert solution == 676333270

python/e147.py

@@ -0,0 +1,35 @@
+from math import ceil
+
+
+def combs_orth(rows, cols):
+    count = 0
+    for r in range(1, rows + 1):
+        for c in range(1, cols + 1):
+            cr = ceil(rows / r)
+            cc = ceil(cols / c)
+            count += cr * cc
+    return count
+
+
+def combs_vert(rows, cols):
+    return 0
+
+
+def combs(rows, cols):
+    return combs_orth(rows, cols) + combs_vert(rows, cols)
+
+
+def euler_147():
+    # assert combs(1, 1) == 1
+    # assert combs(2, 1) == 4
+    # assert combs(3, 1) == 8
+    # assert combs(2, 2) == 18
+    assert combs(2, 3) == 37 - 19
+    return 0
+
+
+if __name__ == "__main__":
+    solution = euler_147()
+    print("e147.py: " + str(solution))
+    # assert(solution == 0)
+

python/e148.py

@@ -0,0 +1,20 @@
+def count(i, s, c, depth):
+    if depth == 0:
+        return s + i, c + 1
+    for i in range(i, i * 8, i):
+        s, c = count(i, s, c, depth - 1)
+        if c == 10**9:
+            break
+    return s, c
+
+
+def euler_148():
+    s, c = count(1, 0, 0, 11)
+    assert c == 10**9
+    return s
+
+
+if __name__ == "__main__":
+    solution = euler_148()
+    print("e148.py: " + str(solution))
+    assert solution == 2129970655314432

python/e149.py

@@ -0,0 +1,105 @@
+def gen_numbers():
+    numbers = []
+    for k in range(1, 56):
+        n = (100_003 - 200_003*k + 300_007*k**3) % 1_000_000 - 500_000
+        numbers.append(n)
+    for k in range(56, 2000 * 2000 + 1):
+        n = (numbers[k-24-1] + numbers[k-55-1] + 1_000_000) % 1_000_000 - 500_000
+        numbers.append(n)
+    return numbers
+
+
+def split(xs, n):
+    return [xs[i:i + n] for i in range(0, len(xs), n)]
+
+
+def diags(rows):
+    n = len(rows)
+    diags = []
+    for col_i in range(-n + 1, n):
+        d = []
+        for row_i in range(n):
+            if col_i >= 0 and col_i < n:
+                d.append(rows[row_i][col_i])
+            col_i += 1
+        diags.append(d)
+    return diags
+
+
+def anti_diags(rows):
+    n = len(rows)
+    diags = []
+    for col_i in range(-n + 1, n):
+        d = []
+        for row_i in range(n - 1, -1, -1):
+            if col_i >= 0 and col_i < n:
+                d.append(rows[row_i][col_i])
+            col_i += 1
+        diags.append(d)
+    return diags
+
+
+def shorten(xs):
+    r = []
+    i = 0
+    while i < len(xs):
+        r.append(0)
+        if xs[i] < 0:
+            while i < len(xs) and xs[i] < 0:
+                r[-1] += xs[i]
+                i += 1
+        else:
+            while i < len(xs) and xs[i] >= 0:
+                r[-1] += xs[i]
+                i += 1
+    if r[0] <= 0:
+        r = r[1:]
+    if r and r[-1] <= 0:
+        r = r[:-1]
+    return r
+
+
+def max_subarray(numbers):
+    """Find the largest sum of any contiguous subarray.
+
+    https://en.wikipedia.org/wiki/Maximum_subarray_problem
+
+    Emberassing that I didn't come up with that myself...
+    """
+    best_sum = -10**12
+    current_sum = 0
+    for x in numbers:
+        current_sum = max(x, current_sum + x)
+        best_sum = max(best_sum, current_sum)
+    return best_sum
+
+
+def max_subarray_naiv(xs):
+    """ Naiv implementation in O(n^3) or something. """
+    r = 0
+    xs = shorten(xs)
+    for width in range(1, len(xs) + 1, 2):
+        for i in range(len(xs)):
+            r = max(r, sum(xs[i:i+width]))
+    return r
+
+
+def euler_149():
+    ns = gen_numbers()
+    assert ns[10-1] == -393027
+    assert ns[100-1] == 86613
+
+    rows = split(ns, 2000)
+    cols = list(map(list, zip(*rows)))
+
+    r = 0
+    for elems in [rows, cols, diags(rows), anti_diags(rows)]:
+        for xs in elems:
+            r = max(max_subarray(xs), r)
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_149()
+    print("e149.py: " + str(solution))
+    assert solution == 52852124

python/e150.py

@@ -0,0 +1,61 @@
+def generate_triangle():
+    t = 0
+    ts = [[]]
+    i_in_row, max_in_row = 0, 1
+    for _ in range(1, 500_500 + 1):
+        t = (615_949 * t + 797_807) % 2**20
+        s = t - 2**19
+        ts[-1].append(s)
+        i_in_row += 1
+        if i_in_row == max_in_row:
+            i_in_row = 0
+            max_in_row += 1
+            ts.append([])
+    ts.pop() # remove empty list
+    return ts
+
+
+def array_to_prefix_sum_array(xs):
+    ys = []
+    current_sum = 0
+    for x in xs:
+        current_sum += x
+        ys.append(current_sum)
+    return ys
+
+
+def range_sum(prefix_sum, i, j):
+    if i == 0:
+        return prefix_sum[j]
+    else:
+        return prefix_sum[j] - prefix_sum[i - 1]
+
+
+def find_min_triangle(row_i, col_i, sum_triangle):
+    current_sum, min_sum = 0, float("inf")
+    left_i, right_i = col_i, col_i
+    while row_i < len(sum_triangle):
+        current_sum += range_sum(sum_triangle[row_i], left_i, right_i)
+        row_i += 1
+        right_i += 1
+        min_sum = min(current_sum, min_sum)
+    return min_sum
+
+
+def euler_150():
+    triangle = generate_triangle()
+    # triangle = [[15], [-14, -7], [20, -13, -5], [-3, 8, 23, -26], [1, -4 , -5, -18, 5], [-16, 31, 2, 9, 28, 3]]
+    sum_triangle = [array_to_prefix_sum_array(xs) for xs in triangle]
+
+    min_triangle = float("inf")
+    for row_i in range(len(triangle)):
+        for col_i in range(len(triangle[row_i])):
+            potential_min = find_min_triangle(row_i, col_i, sum_triangle)
+            min_triangle = min(min_triangle, potential_min)
+    return min_triangle
+
+
+if __name__ == "__main__":
+    solution = euler_150()
+    print("e150.py: " + str(solution))
+    assert solution == -271248680

python/e151.py

@@ -0,0 +1,31 @@
+from fractions import Fraction
+from functools import lru_cache
+
+
+@lru_cache
+def expected_singles(sheets):
+    count = 0
+    if len(sheets) == 1 and sheets[0] == 5:
+        return 0
+    elif len(sheets) == 1:
+        count += 1
+
+    weight = Fraction(1, len(sheets))
+    for sheet in sheets:
+        nsheets = list(sheets)
+        nsheets.remove(sheet)
+        for nsheet in range(sheet, 5):
+            nsheets.append(nsheet + 1)
+        count += weight * expected_singles(tuple(sorted(nsheets)))
+    return count
+
+
+def euler_151():
+    return round(float(expected_singles(tuple([2, 3, 4, 5]))), 6)
+
+
+if __name__ == "__main__":
+    solution = euler_151()
+    print("e151.py: " + str(solution))
+    assert(solution == 0.464399)
+

python/e155.py

@@ -0,0 +1,40 @@
+from fractions import Fraction
+from functools import lru_cache
+
+C = 60
+
+def parallel(a, b):
+    return Fraction(a * b, a + b)
+
+
+def series(a, b):
+    return a + b
+    
+
+@lru_cache()
+def caps(n):
+    if n == 1:
+        return set([C])
+
+    cs = set()
+    for size_a in range(1, (n // 2) + 1):
+        size_b = n - size_a
+        for a in caps(size_a):
+            for b in caps(size_b):
+                cs.add(a)
+                cs.add(parallel(a, b))
+                cs.add(series(a, b))
+    return cs
+
+
+def euler_155():
+    assert len(caps(1)) == 1
+    assert len(caps(2)) == 3
+    assert len(caps(3)) == 7
+    return len(caps(18))
+
+
+if __name__ == "__main__":
+    solution = euler_155()
+    print("e155.py: " + str(solution))
+    assert solution == 3857447

python/e161.py

@@ -0,0 +1,108 @@
+from copy import copy
+from typing import Optional, Tuple
+
+
+CACHE = {}
+EMPTY = '⬛'
+FIELDS = ['🔴', '🟠', '🟡', '🟢', '🔵', '🟣', '⚪', '🟤', '⭐']
+
+
+# We use [row, col] indexing.
+BLOCKS = [
+    ((0, 0), (0, 1), (0, 2)), # dark blue
+    ((0, 0), (1, 0), (2, 0)), # black
+    ((0, 0), (0, 1), (1, 0)), # red
+    ((0, 0), (0, 1), (1, 1)), # green
+    ((0, 0), (1, 0), (1, 1)), # blue
+    ((0, 0), (1, 0), (1, -1)), # orange
+]
+
+
+class Field:
+
+    def __init__(self, N_COLS, N_ROWS):
+        self.N_COLS = N_COLS
+        self.N_ROWS = N_ROWS
+        self.state = 0
+
+    def __hash__(self):
+        return hash(self.state)
+
+    def get_first_empty(self) -> Optional[Tuple[int, int]]:
+        for row in range(self.N_ROWS):
+            for col in range(self.N_COLS):
+                if self.is_empty(row, col):
+                    return (row, col)
+        return None
+
+    def is_empty(self, row: int, col: int) -> bool:
+        index = row * self.N_COLS + col
+        result = not bool(self.state & (1 << index))
+        return result
+
+    def set(self, row: int, col: int):
+        index = row * self.N_COLS + col
+        self.state |= (1 << index)
+
+    def __getitem__(self, row: int):
+        return [self.is_empty(row, col) for col in range(self.N_COLS)]
+
+    def print(self):
+        for row in range(self.N_ROWS):
+            row_str = ""
+            for col in range(self.N_COLS):
+                row_str += '⬛' if self.is_empty(row, col) else '⚪'
+            print(row_str)
+
+
+def fits(field, block, coord):
+    N_ROWS = field.N_ROWS
+    N_COLS = field.N_COLS
+
+    for rel_coord in block:
+        abs_row = coord[0] + rel_coord[0]
+        abs_col = coord[1] + rel_coord[1]
+
+        if abs_row < 0 or abs_col < 0 or abs_row >= N_ROWS or abs_col >= N_COLS:
+            return None
+
+        if not field.is_empty(abs_row, abs_col):
+            return None
+
+    new_field = copy(field)
+    for rel_coord in block:
+        abs_row = coord[0] + rel_coord[0]
+        abs_col = coord[1] + rel_coord[1]
+        new_field.set(abs_row, abs_col)
+
+    return new_field
+
+
+def count(field):
+    global CACHE
+
+    if hash(field) in CACHE:
+        return CACHE[hash(field)]
+
+    first_empty = field.get_first_empty()
+    if first_empty is None:
+        return 1
+    
+    result = 0
+    for block in BLOCKS:
+        if (new_field := fits(field, block, first_empty)) is not None:
+            result += count(new_field)
+
+    CACHE[hash(field)] = result
+    return result
+
+
+def euler_161():
+    return count(Field(9, 12))
+
+
+if __name__ == "__main__":
+    solution = euler_161()
+    print("e161.py: " + str(solution))
+    assert(solution == 20574308184277971)
+

python/e162.py

@@ -0,0 +1,39 @@
+from functools import lru_cache
+
+
+@lru_cache
+def combs(xs, first=False):
+    if xs == (0, 0, 0, 0):
+        return 1
+    
+    c = 0
+    z, o, a, r = xs
+
+    if (not first) and z > 0:
+        c += combs((z - 1, o, a, r))
+    if o > 0:
+        c += combs((z, o - 1, a, r))
+    if a > 0:
+        c += combs((z, o, a - 1, r))
+    if r > 0:
+        c += 13 * combs((z, o, a, r - 1))
+    return c
+
+
+def euler_162():
+    t = 0
+    for n in range(17):
+        for z in range(1, n + 1):
+            for o in range(1, n + 1):
+                for a in range(1, n + 1):
+                    r = n - z - o - a
+                    if not r >= 0:
+                        continue
+                    t += combs((z, o, a, r), True)
+    return hex(t)[2:].upper()
+
+
+if __name__ == "__main__":
+    solution = euler_162()
+    print("e162.py: " + str(solution))
+    assert(solution == "3D58725572C62302")

python/e164.py

@@ -0,0 +1,25 @@
+from functools import lru_cache
+
+
+@lru_cache
+def count(last, left):
+    if left == 0:
+        return 1
+    c = 0
+    for d in range(10):
+        if sum(last) + d <= 9:
+            c += count((last[-1], d), left - 1)
+    return c
+
+
+def euler_164():
+    c = 0
+    for d in range(1, 10):
+        c += count((d, ), 20-1)
+    return c
+
+
+if __name__ == "__main__":
+    solution = euler_164()
+    print("e164.py: " + str(solution))
+    assert solution == 378158756814587

python/e173.py

@@ -0,0 +1,57 @@
+
+
+def tiles(size):
+    """ Returns the number of tiles for a square of the given size. """
+    return size * 4 - 4
+
+
+def tiles_for_outer(size):
+    """ Returns the number of tiles for all possible laminae given the size of
+    the outermost ring. """
+
+    n_outer = tiles(size)
+    r = [n_outer]
+    for ni in range(size - 2, 2, -2):
+        n_outer += tiles(ni)
+        r.append(n_outer)
+    return r
+
+
+def count_tiles(size, n_max):
+
+    """ Count how many possible laminae with a size smaller the `n_max` are
+    possible for the given outermost ring `size`. """
+
+    n_outer = tiles(size)
+    if n_outer > n_max:
+        return 0
+
+    count = 1
+    for ni in range(size - 2, 2, -2):
+        n_outer += tiles(ni)
+        if n_outer > n_max:
+            break
+        count += 1
+    return count
+
+
+def euler_173():
+    assert tiles(3) == 8
+    assert tiles(4) == 12
+    assert tiles(6) == 20
+    assert tiles(9) == 32
+
+    count = 0
+    n_max = 10**6
+    edge_length = 3
+    for edge_length in range(3, n_max):
+        if tiles(edge_length) > n_max:
+            break
+        count += count_tiles(edge_length, n_max)
+    return count
+
+
+if __name__ == "__main__":
+    solution = euler_173()
+    print("e173.py: " + str(solution))
+    assert solution == 1572729

python/e174.py

@@ -0,0 +1,42 @@
+from collections import defaultdict
+
+
+def tiles(size):
+    """ Returns the number of tiles for a square of the given size. """
+    return size * 4 - 4
+
+
+def tiles_for_outer(size, ts, limit):
+    """ Adds t to ts for given outer size as long as t is
+    smaller than the limit. """
+    t = tiles(size)
+    if t > limit:
+        return
+    ts[t] += 1
+    for ni in range(size - 2, 2, -2):
+        t += tiles(ni)
+        if t > limit:
+            break
+        ts[t] += 1
+
+
+def euler_173():
+    ts = defaultdict(int)
+    limit = 1_000_000
+
+    for outer_size in range(3, limit // 4):
+        tiles_for_outer(outer_size, ts, limit)
+
+    c = 0
+    for n in range(1, 11):
+        for v in ts.values():
+            if v == n:
+                c += 1
+
+    return c
+
+
+if __name__ == "__main__":
+    solution = euler_173()
+    print("e173.py: " + str(solution))
+    assert solution == 209566

python/e179.py

@@ -0,0 +1,33 @@
+from lib_prime import get_divisors_count
+
+
+def euler_179_orig():
+    r = 0
+    for n in range(2, 10**7):
+        if n % 10**5 == 0:
+            print(n)
+        if get_divisors_count(n) == get_divisors_count(n + 1):
+            r += 1
+    return r
+
+
+def euler_179():
+    """ Much faster version. """
+    upper = 10**7
+    ndivs = [1 for _ in range(0, upper + 2)]
+    for d in range(2, upper // 2):
+        for i in range(d, upper + 2, d):
+            ndivs[i] += 1
+
+    r = 0
+    for i in range(2, upper + 1):
+        if ndivs[i] == ndivs[i + 1]:
+            r += 1
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_179()
+    print("e179.py: " + str(solution))
+    assert(solution == 986262)
+

python/e187.py

@@ -0,0 +1,19 @@
+from lib_prime import primes
+
+
+def euler_187():
+    r = 0
+    upper = 10**8
+    ps = primes(upper // 2)
+    for i in range(len(ps)):
+        for j in range(i, len(ps)):
+            if ps[i] * ps[j] >= upper:
+                break
+            r += 1
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_187()
+    print("e187.py: " + str(solution))
+    assert(solution == 17427258)

python/e188.py

@@ -0,0 +1,38 @@
+from lib_prime import prime_factors
+
+
+def phi(n):
+    # Calculate phi using Euler's totient function
+    c = n
+    fs = set(prime_factors(n))
+    for f in fs:
+        # same as c *= (1 - 1/f)
+        c *= (f - 1)
+        c //= f
+    return c
+
+
+def tetmod(a, k, m):
+    # https://stackoverflow.com/questions/30713648/how-to-compute-ab-mod-m
+    assert k > 0
+    if k == 1:
+        return a
+    if m == 1:
+        if a % 2 == 0:
+            return 0
+        else:
+            return 1
+    phi_m = phi(m)
+    return pow(a, phi_m + tetmod(a, k - 1, phi_m), m)
+
+
+def euler_188():
+    assert phi(100) == 40
+    assert tetmod(14, 2016, 100) == 36
+    return tetmod(1777, 1855, 10**8)
+
+
+if __name__ == "__main__":
+    solution = euler_188()
+    print("e188.py: " + str(solution))
+    assert solution == 95962097

python/e191.py

@@ -0,0 +1,38 @@
+from functools import lru_cache
+
+
+@lru_cache
+def count(days_left, absent_count):
+    if absent_count > 1:
+        return 0
+    if days_left == 0:
+        return 1
+    if days_left < 0:
+        return 0
+    c = 0
+    c += count(days_left - 1, absent_count + 1) # "a"
+    c += count(days_left - 2, absent_count + 1) # "la"
+    c += count(days_left - 3, absent_count + 1) # "lla"
+
+    c += count(days_left - 1, absent_count)     # "o"
+    c += count(days_left - 2, absent_count)     # "lo"
+    c += count(days_left - 3, absent_count)     # "llo"
+
+    if days_left == 2:
+        c += 1 # "ll"
+
+    if days_left == 1:
+        c += 1 # "l"
+
+    return c
+
+
+def euler_191():
+    return count(30, 0)
+
+
+if __name__ == "__main__":
+    solution = euler_191()
+    print("e191.py: " + str(solution))
+    assert(solution == 1918080160)
+

python/e195.py

@@ -0,0 +1,50 @@
+import math
+
+
+def inradius(a, b, c):
+    s = (a + b + c) / 2
+    area = (s * (s - a) * (s - b) * (s - c))
+    inradius = area / s**2
+    return inradius
+
+
+def euler_195():
+    count = 0
+    limit = 1053779
+    limit2 = limit * limit
+    for m in range(0, limit2):
+        for n in range(1, m // 2 + 1):
+            if math.gcd(m, n) == 1:
+                # Generator function from wiki for 60 degree eisenstein
+                # triangles.
+                a = m**2 - m*n + n**2
+                b = 2*m*n - n**2
+                c = m**2 - n**2
+                gcd = math.gcd(a, b, c)
+                a, b, c = a // gcd, b // gcd, c // gcd
+
+                # if gcd == 3 and inradius(a, b, c) > limit2:
+                #     break
+                if inradius(a / 3, b / 3, c / 3) > limit2:
+                    if n == 1:
+                        return count
+                    break
+
+                if a == b and b == c:
+                    continue
+
+                k = 1
+                while True:
+                    ir = inradius(a*k, b*k, c*k)
+                    if ir > limit2:
+                        break
+                    count += 1
+                    k += 1
+
+    return count
+
+
+if __name__ == "__main__":
+    solution = euler_195()
+    print("e195.py: " + str(solution))
+    assert solution == 75085391 

python/e203.py

@@ -0,0 +1,54 @@
+from lib_prime import primes
+from math import sqrt
+
+
+def pascal_numbers_till(row):
+    def pascal_next(xs):
+        r = [xs[0]]
+        for i in range(len(xs) - 1):
+            r.append(xs[i] + xs[i + 1])
+        r.append(xs[-1])
+        return r
+
+    all = set([1])
+    xs = [1]
+    for _ in range(row - 1):
+        xs = pascal_next(xs)
+        all |= set(xs)
+    return sorted(list(set(all)))
+
+
+def is_square_free(ps, n):
+    for p in ps:
+        if n % p == 0:
+            return False
+    return True
+
+
+def solve(n):
+    r = 0
+    pascals = pascal_numbers_till(n)
+    ps = [p * p for p in primes(int(sqrt(max(pascals))))]
+    for n in pascals:
+        if is_square_free(ps, n):
+            r += n
+    return r
+
+
+def test():
+    assert pascal_numbers_till(1) == [1]
+    assert pascal_numbers_till(2) == [1]
+    assert pascal_numbers_till(3) == [1, 2]
+    assert pascal_numbers_till(8) == [1, 2, 3, 4, 5, 6, 7, 10, 15, 20, 21, 35]
+    assert solve(8) == 105
+
+
+def euler_203():
+    test()
+    return solve(51)
+
+
+if __name__ == "__main__":
+    solution = euler_203()
+    print("e203.py: " + str(solution))
+    assert solution == 34029210557338

python/e204.py

@@ -0,0 +1,38 @@
+from lib_prime import primes
+from functools import lru_cache
+
+
+@lru_cache
+def next_prime(n):
+    if n == 1:
+        return 2
+    ps = primes(105)
+    assert n < ps[-1]
+    for i, p in enumerate(ps):
+        if p == n and i + 1 < len(ps):
+            return ps[i + 1]
+    assert False
+
+
+def count(current_product, max_product, current_prime, max_prime):
+    if current_product > max_product:
+        return 0
+    r = 1
+
+    if current_prime > 1:
+        r += count(current_product * current_prime, max_product, current_prime, max_prime)
+
+    while (current_prime := next_prime(current_prime)) <= max_prime:
+        r += count(current_product * current_prime, max_product, current_prime, max_prime)
+    return r
+
+
+def euler_204():
+    assert count(1, 10**8, 1, 5) == 1105
+    return count(1, 10**9, 1, 100)
+
+
+if __name__ == "__main__":
+    solution = euler_204()
+    print("e204.py: " + str(solution))
+    assert solution == 2944730

python/e205.py

@@ -0,0 +1,34 @@
+from lib_misc import get_item_counts
+
+
+def acc(scores, faces, rem):
+    if rem == 0:
+        return scores
+    nscores = []
+    for f in faces:
+        for s in scores:
+            nscores.append(f + s)
+    return acc(nscores, faces, rem - 1)
+
+
+def euler_205():
+    scores_four = acc([0], [1, 2, 3, 4], 9)
+    scores_six = acc([0], [1, 2, 3, 4, 5, 6], 6)
+    all_combinations = len(scores_four) * len(scores_six)
+
+    four_won_count = 0
+    for four_value, four_count in get_item_counts(scores_four).items():
+        current_six_count = 0
+        for six_value, six_count in get_item_counts(scores_six).items():
+            if four_value > six_value:
+                current_six_count += six_count
+        four_won_count += four_count * current_six_count
+
+    p = four_won_count / all_combinations
+    return f"{round(p, 7)}"
+
+
+if __name__ == "__main__":
+    solution = euler_205()
+    print("e205.py: " + str(solution))
+    # assert(solution == 0)

python/e206.py

@@ -2,12 +2,13 @@ import math
 
 
 def is_square(apositiveint):
-    """ From: https://stackoverflow.com/a/2489519 """
+    """From: https://stackoverflow.com/a/2489519"""
     x = apositiveint // 2
     seen = set([x])
     while x * x != apositiveint:
         x = (x + (apositiveint // x)) // 2
-         if x in seen: return False
+        if x in seen:
+            return False
         seen.add(x)
     return True
 
@@ -16,7 +17,7 @@ def canditates(current_value, depth):
     if depth <= 0:
         yield current_value
     else:
-        d = 10 ** depth
+        d = 10**depth
         for i in range(10):
             for new_value in canditates(current_value - i * d, depth - 2):
                 yield new_value
@@ -32,5 +33,4 @@ def euler_206():
 if __name__ == "__main__":
     solution = euler_206()
     print("e206.py: " + str(solution))
-    assert(solution == 1389019170)
+    assert solution == 1389019170
-

python/e243.py

@@ -0,0 +1,40 @@
+from math import gcd
+from fractions import Fraction
+
+
+def product(xs):
+    r = 1
+    for x in xs:
+        r *= x
+    return r
+
+
+def ratio(d):
+    r = 0
+    for i in range(1, d):
+        if gcd(i, d) == 1:
+            r += 1
+    return Fraction(r, (d - 1))
+
+
+def euler_243():
+    target = Fraction(15499, 94744)  # 0.1635881955585578
+
+    # The more prime factors, the lower the initial fraction gets. We figure
+    # out empirically that using primes up to 23 yields results close to the
+    # target fraction. We then iterate over the multiples to find the
+    # solution.
+    factors = [2, 3, 5, 7, 11, 13, 17, 19, 23]
+
+    for mul in range(1, 10):
+        d = mul * product(factors)
+        r = ratio(d)
+        if r < target:
+            return d
+    return None
+
+
+if __name__ == "__main__":
+    solution = euler_243()
+    print("e243.py: " + str(solution))
+    assert solution == 892371480

python/e293.py

@@ -0,0 +1,70 @@
+from lib_prime import primes, is_prime_rabin_miller
+from math import isqrt
+
+
+def next_larger_prime(n):
+    """ Returns the next prime larger or equal to n. """
+    if n == 1:
+        return 2
+
+    n += 1
+    if n % 2 == 0:
+        n += 1
+
+
+    for _ in range(1000):
+        if is_prime_rabin_miller(n):
+            return n
+        n += 2
+    assert False
+
+
+def m(n):
+    higher_prime = next_larger_prime(n + 1)
+    m = higher_prime - n
+    return m
+
+
+def admissible(xs, threshold):
+    found = set()
+    numbers = [(xs[0], 0)]
+
+    while numbers:
+        value, index = numbers.pop()
+        if value >= threshold:
+            continue
+        if index == len(xs) - 1:
+            found.add(value)
+        numbers.append((value * xs[index], index))
+        if index < len(xs) - 1:
+            numbers.append((value * xs[index + 1], index + 1))
+    return found
+
+
+def euler_293():
+    n = 10**9
+    assert next_larger_prime(1) == 2
+    assert next_larger_prime(7) == 11
+    assert m(630) == 11
+
+    m_set = set()
+    ps = primes(isqrt(n) + 10*8)
+
+    pow2 = 2
+    while pow2 < n:
+        mv = m(pow2)
+        m_set.add(mv)
+        pow2 *= 2
+
+    for i in range(len(ps)):
+        xs = ps[0:i + 1]
+        for x in admissible(xs, n):
+            m_set.add(m(x))
+
+    return sum(m_set)
+
+
+if __name__ == "__main__":
+    solution = euler_293()
+    print("e293.py: " + str(solution))
+    assert solution == 2209

python/e301.py

@@ -0,0 +1,70 @@
+from functools import lru_cache
+
+
+@lru_cache
+def x_naiv(s):
+    zcount = s.count(0)
+    if zcount == 3:
+        return 0
+    elif zcount == 2:
+        return 1
+    elif zcount == 1:
+        sl = list(s)
+        sl.remove(0)
+        if sl[0] == sl[1]:
+            return 0
+
+    s = list(s)
+    for i in range(len(s)):
+        for v in range(1, s[i] + 1):
+            ns = list(s)
+            ns[i] -= v
+            if x_naiv(tuple(sorted(ns))) == 0: 
+                # Other player loses.
+                return 1
+    return 0
+
+
+def x_nim_sum(s):
+    """ I wasn't able to deduce this myself quickly and looked it up instead. """
+    return 0 if (s[0] ^ s[1] ^ s[2]) == 0 else 1
+
+
+def test_nim_functions():
+    assert(x_naiv(tuple([0, 0, 0])) == 0)
+    assert(x_naiv(tuple([1, 0, 0])) == 1)
+    assert(x_naiv(tuple([1, 1, 0])) == 0)
+    assert(x_naiv(tuple([1, 1, 1])) == 1)
+    assert(x_naiv(tuple([1, 2, 3])) == 0)
+
+    for n1 in range(0, 8):
+        for n2 in range(1, n1 + 1):
+            for n3 in range(1, n2 + 1):
+                s = tuple([n1, n2, n3])
+                assert x_naiv(s) == x_nim_sum(s)
+
+
+def euler_301():
+    test_nim_functions()
+
+    # log(2**30, 10) ~= 9 -> 1 billion easy brute force
+    # Probably a more efficient computation is possible by counting the permutations
+    # where all binary digits are zero.
+    # 0 0 0 -> 0
+    # 0 0 1 -> 0
+    # 0 1 1 -> 0
+
+    r = 0
+    for n in range(1, 2**30 + 1):
+        s = tuple([n, 2 * n, 3 * n])
+        xr = x_nim_sum(s)
+        if xr == 0:
+            r += 1
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_301()
+    print("e301.py: " + str(solution))
+    # assert(solution == 0)
+

python/e345.py

@@ -0,0 +1,76 @@
+from itertools import permutations
+from heapq import heappush, heappop
+
+
+m1 = """  7  53 183 439 863
+        497 383 563  79 973
+        287  63 343 169 583
+        627 343 773 959 943
+        767 473 103 699 303"""
+
+
+m2 = """ 7  53 183 439 863 497 383 563  79 973 287  63 343 169 583
+627 343 773 959 943 767 473 103 699 303 957 703 583 639 913
+447 283 463  29  23 487 463 993 119 883 327 493 423 159 743
+217 623   3 399 853 407 103 983  89 463 290 516 212 462 350
+960 376 682 962 300 780 486 502 912 800 250 346 172 812 350
+870 456 192 162 593 473 915  45 989 873 823 965 425 329 803
+973 965 905 919 133 673 665 235 509 613 673 815 165 992 326
+322 148 972 962 286 255 941 541 265 323 925 281 601  95 973
+445 721  11 525 473  65 511 164 138 672  18 428 154 448 848
+414 456 310 312 798 104 566 520 302 248 694 976 430 392 198
+184 829 373 181 631 101 969 613 840 740 778 458 284 760 390
+821 461 843 513  17 901 711 993 293 157 274  94 192 156 574
+ 34 124   4 878 450 476 712 914 838 669 875 299 823 329 699
+815 559 813 459 522 788 168 586 966 232 308 833 251 631 107
+813 883 451 509 615  77 281 613 459 205 380 274 302  35 805
+"""
+
+
+def brute_force(m):
+    idxs = [i for i in range(len(m))]
+    max_sum = 0
+    for idx in permutations(idxs):
+        temp_sum = 0
+        for i, j in enumerate(idx):
+            temp_sum += m[j][i]
+        if temp_sum > max_sum:
+            max_sum = temp_sum
+    return max_sum
+
+
+def a_star(m):
+    # init heuristic function
+    h = {}
+    for i in range(len(m)):
+        h[i] = sum(map(min, m[i:]))
+        h[i + 1] = 0
+
+    states = []
+    for i, r in enumerate(m[0]):
+        state = (r + h[1], r, 1, [i])
+        heappush(states, state)
+
+    while states:
+        _, gscore, index, seen = heappop(states)
+        if index == len(m):
+            return -gscore
+
+        for row_index, row_score in enumerate(m[index]):
+            if row_index in seen:
+                continue
+            ngscore = gscore + row_score
+            nfscore = ngscore + h[index + 1]
+            nstate = (nfscore, ngscore, index + 1, seen + [row_index])
+            heappush(states, nstate)
+
+
+def euler_345():
+    m = list(zip(*map(lambda line: list(map(lambda n: -int(n), line.split())), m2.splitlines())))
+    return a_star(m)
+
+
+if __name__ == "__main__":
+    solution = euler_345()
+    print("e345.py: " + str(solution))
+    assert(solution == 13938)

python/e346.py

@@ -0,0 +1,23 @@
+
+def euler_346():
+    n_max = 10**12
+    repunits = set()
+    result = 1  # Our algorithm doesn't catch 1 so we add it here
+    for base in range(2, n_max):
+        n = base ** 0 + base ** 1
+        for exp in range(2, n_max):
+            n += base**exp
+            if exp == 2 and n > n_max:
+                # There are no more repunits below n_max.
+                return result
+            if n >= n_max:
+                break
+            if not n in repunits:
+                repunits.add(n)
+                result += n
+
+
+if __name__ == "__main__":
+    solution = euler_346()
+    print("e346.py: " + str(solution))
+    assert(solution == 336108797689259276)

python/e347.py

@@ -0,0 +1,39 @@
+from lib_prime import primes
+from typing import Optional, Tuple
+
+
+def s(n):
+    ps = primes(n)
+
+    sum = 0
+    for p1_index in range(len(ps)):
+        p1 = ps[p1_index]
+        for p2_index in range(p1_index + 1, len(ps)):
+            p2 = ps[p2_index]
+            if p1 * p2 > n:
+                break
+            largest: Optional[Tuple[int, int, int]] = None
+            p1e = 1
+            while p1**p1e * p2 <= n:
+                p2e = 1
+                m = p1**p1e * p2**p2e
+                while m <= n:
+                    if largest is None or m > largest[2]:
+                        largest = (p1, p2, m)
+                    p2e += 1
+                    m = p1**p1e * p2**p2e
+                p1e += 1
+            assert largest is not None
+            sum += largest[2]
+    return sum
+
+
+def euler_347():
+    assert s(100) == 2262
+    return s(10_000_000)
+
+
+if __name__ == "__main__":
+    solution = euler_347()
+    print("e347.py: " + str(solution))
+    assert solution == 11109800204052

python/e357.py

@@ -0,0 +1,41 @@
+def primes_sieve(nmax):
+    a = [0 for _ in range(nmax)]
+    for i in range(2, nmax):
+        if a[i] == 0:
+            for j in range(i + i, nmax, i):
+                a[j] = i
+    return a
+
+
+def is_number(n, primes):
+    f = n // 2
+    if primes[f] == 0:
+        for d in [2, f]:
+            if primes[d + n // d] != 0:
+                return False
+        return True
+
+    for d in range(1, n // 2 + 1):
+        if n % d == 0:
+            if primes[d + n // d] != 0:
+                return False
+        if d * d > n:
+            break
+    return True
+
+
+def euler_357():
+    r = 1 # n = 1 is a number
+    ps = primes_sieve(100_000_000)
+    for p in range(3, len(ps)):
+        if ps[p] == 0:
+            n = p - 1
+            if is_number(n, ps):
+                r += n
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_357()
+    print("e357.py: " + str(solution))
+    assert solution == 1739023853137

python/e493.py

@@ -0,0 +1,63 @@
+from functools import lru_cache
+
+
+BALLS_PER_COLOR = 10
+N_COLORS = 7
+TO_DRAW = 20
+
+
+@lru_cache(maxsize=10**9)
+def count(urn, to_draw):
+    if to_draw == 0:
+        distinct_colors = sum([1 if c != BALLS_PER_COLOR else 0 for c in urn])
+        return distinct_colors
+
+    remaining_balls = sum(urn)
+    expected = 0
+    for color_i, color_remaining in enumerate(urn):
+        if color_remaining == 0:
+            continue
+        new_urn = list(urn)
+        new_urn[color_i] -= 1
+        expected_for_color = count(tuple(sorted(new_urn)), to_draw - 1)
+        probability = color_remaining / remaining_balls
+        expected += (probability * expected_for_color)
+    return expected
+
+
+def binomial_coefficient(n, k):
+    if k < 0 or k > n:
+        return 0
+    if k == 0 or k == n:
+        return 1
+    k = min(k, n - k)
+    result = 1
+    for i in range(k):
+        result *= n - i
+        result //= i + 1
+    return result
+
+
+def math():
+    # Find probability that a given color is chose in TO_DRAW draws.
+    combs_without_color = binomial_coefficient((N_COLORS - 1) * BALLS_PER_COLOR, TO_DRAW)
+    combs_with_color = binomial_coefficient(N_COLORS * BALLS_PER_COLOR, TO_DRAW)
+    p_single_color = 1 - (combs_without_color / combs_with_color)
+    # Then, expected colors is that probability times the number of colors.
+    return p_single_color * N_COLORS
+
+
+def euler_493():
+    urn = [BALLS_PER_COLOR for _ in range(N_COLORS)]
+    c = count(tuple(urn), TO_DRAW)
+    c = round(c, 9)
+    m = round(math(), 9)
+    assert c == m
+    return c
+
+
+if __name__ == "__main__":
+    solution = euler_493()
+    print("e493.py: " + str(solution))
+    # assert(solution == 0)
+

python/e684.py

@@ -0,0 +1,54 @@
+def r_modulo_closed_form(n, m):
+    # From e132.py
+    assert n > 0 and m > 0
+    return ((pow(10, n, 9 * m) - 1) // 9) % m
+
+
+def s(n, mod):
+    a = n % 9
+    b = n // 9
+    r = (a + 1) * pow(10, b, mod) - 1
+    return r % mod
+
+
+def s_big(n, mod):
+    # s = (n % 9 + 1) * 10^(n/9) - 1
+    # S = sum(k=0, n/9) (1 * 10^k + 2 * 10^k + ... + 9 * 10^k) - n
+    #   = 45 * sum(k=0, n/9) 10^k - n
+    #   = 45 * r(n/9) - n
+    # Plus some extra math to make n % 9 == 8 so that the above works.
+
+    r = -1
+    while n % 9 != 8:
+        n += 1
+        r -= s(n, mod)
+
+    r += (r_modulo_closed_form(n // 9 + 1, mod) * 45) % mod
+    r -= n
+    return r % mod
+
+
+def s_big_naiv(n, mod):
+    r = 0
+    for i in range(1, n + 1):
+        r += s(i, mod)
+        if mod is not None:
+            r %= mod
+    return r
+
+
+def euler_684():
+    assert s_big_naiv(20, 11) == s_big(20, 11)
+    r = 0
+    mod = 1_000_000_007
+    a, b = 0, 1
+    for _ in range(2, 91):
+        a, b = b, a + b
+        r += s_big(b, mod)
+    return r % mod
+
+
+if __name__ == "__main__":
+    solution = euler_684()
+    print("e684.py: " + str(solution))
+    assert(solution == 922058210)

python/e686.py

@@ -0,0 +1,29 @@
+def p(leading_digits, count):
+    n_cutoff = 18
+    leading_digits_str = str(leading_digits)
+    n_digits = len(leading_digits_str)
+    j, current_leading, found = 1, 2, 0
+
+    while found < count:
+        current_leading *= 2
+        j += 1
+        current_leading_str = str(current_leading)
+        if  current_leading_str[:n_digits] == leading_digits_str:
+            found += 1
+        current_leading = int(current_leading_str[:n_cutoff])
+    return j
+
+
+def euler_686():
+    assert p(12, 1) == 7
+    assert p(12, 2) == 80
+    assert p(123, 45) == 12710
+    assert p(123, 100) == 28343
+    assert p(123, 1000) == 284168
+    return p(123, 678910)
+
+
+if __name__ == "__main__":
+    solution = euler_686()
+    print("e686.py: " + str(solution))
+    assert solution == 193060223

python/e700.py

@@ -0,0 +1,16 @@
+def euler_700():
+    mod = 4503599627370517
+    mul = 1504170715041707
+    result = mul
+    while mul > 1:
+        while mod > mul:
+            mod = mod % mul
+        mul -= mod
+        result += mul
+    return result
+
+
+if __name__ == "__main__":
+    solution = euler_700()
+    print("e700.py: " + str(solution))
+    assert solution == 1517926517777556

python/e719.py

@@ -0,0 +1,69 @@
+from typing import List, Optional
+from math import sqrt
+from functools import lru_cache
+
+
+def to_int(digits: List[int]) -> int:
+    return int("".join(map(str, digits)))
+
+
+def int_len(n: int) -> int:
+    return len(str(n))
+
+
+@lru_cache
+def digits_make_sum(digits, remainder: int) -> Optional[List]:
+    if len(digits) == 0:
+        if remainder == 0:
+            return []
+        return None
+
+    remainder_digits_count = int_len(remainder)
+    if remainder_digits_count > len(digits):
+        return None
+
+    if to_int(digits) < remainder:
+        return None
+
+    for group_size in range(1, remainder_digits_count + 1): 
+        digits_list = list(digits)
+        rest_value = remainder - to_int(digits_list[:group_size])
+        rest_digits = tuple(digits_list[group_size:])
+        r = digits_make_sum(rest_digits, rest_value)
+        if type(r) is list:
+            return [digits_list[:group_size]] + r
+    return None
+
+
+def digits_sum_to(digits, number: int):
+    r = digits_make_sum(digits, number)
+    if type(r) is list:
+        # print(digits, r, number)
+        return True
+    return False
+
+
+def t(limit: int) -> int:
+    r = 0
+    for base in range(4, int(sqrt(limit)) + 1):
+        n = base * base
+        n_digits = tuple(map(int, list(str(n))))
+        if digits_sum_to(n_digits, base):
+            r += n
+    return r
+
+
+def euler_719():
+    assert(t(10**4) == 41333)
+    assert(digits_sum_to((1, 0, 0), 10) == True)
+    assert(digits_sum_to((1, 8), 9) == True)
+    assert(digits_sum_to((2, 5), 5) == False)
+    assert(digits_sum_to((8, 2, 8, 1), 91) == True)
+    assert(t(10**6) == 10804656)
+    return t(10**12)
+
+
+if __name__ == "__main__":
+    solution = euler_719()
+    print("e719.py: " + str(solution))
+    assert(solution == 128088830547982)

python/e751.py

@@ -0,0 +1,35 @@
+from math import floor
+
+
+def next_b(b_prev):
+    b = floor(b_prev) * (b_prev - floor(b_prev) + 1)
+    return b
+
+
+def euler_751():
+
+    min_found = 2
+    theta_str_orig = "2."
+    while len(theta_str_orig) < 26:
+        for n in range(min_found, 1000):
+            theta_str_new = str(theta_str_orig) + str(n)
+            b = float(theta_str_new)
+            theta_str = theta_str_new.replace("2.", "")
+            while theta_str != "":
+                b = next_b(b)
+                a_str = str(floor(b))
+                # print(a_str, theta_str)
+                if theta_str.startswith(a_str):
+                    theta_str = theta_str.replace(a_str, "", 1)
+                else:
+                    break
+            if theta_str == "":
+                theta_str_orig = theta_str_new
+                break
+    return theta_str_orig
+
+
+if __name__ == "__main__":
+    solution = euler_751()
+    print("e751.py: " + str(solution))
+    assert solution == "2.223561019313554106173177"

python/e808.py

@@ -0,0 +1,25 @@
+from lib_prime import primes
+from lib_misc import is_palindrome_string
+
+
+def euler_808():
+    r = []
+    squared_primes = set()
+    ps = primes(100000000)
+    for p in ps:
+        p_squared = p * p
+        p_squared_str = str(p_squared)
+        if not is_palindrome_string(p_squared_str):
+            squared_primes.add(p_squared)
+            p_squared_reverse = int(p_squared_str[::-1])
+            if p_squared_reverse in squared_primes:
+                r += [p_squared, p_squared_reverse]
+                if len(r) == 50:
+                    break
+    return sum(r)
+
+
+if __name__ == "__main__":
+    solution = euler_808()
+    print("e808.py: " + str(solution))
+    assert(solution == 3807504276997394)

python/e816.py

@@ -0,0 +1,39 @@
+from math import sqrt
+
+
+def distance(p1, p2) -> float:
+    return (p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2
+
+
+def points(k: int):
+    s = 290797
+    mod = 50515093
+    # s = 5
+    # mod = 23
+    points = []
+    for _ in range(k):
+        sn = (s ** 2) % mod
+        p = (s, sn)
+        points.append(p)
+        s = (sn ** 2) % mod
+    return points
+
+
+def euler_816():
+    ps = sorted(points(2000000), key=lambda p: p[0])
+    d_min = float("inf")
+    sweep_dist = 5  # heuristic
+    for i in range(len(ps)):
+        for j in range(i, min(i + sweep_dist, len(ps))):
+            if i == j:
+                continue
+            d  = distance(ps[i], ps[j])
+            if d < d_min:
+                d_min = d
+    return round(sqrt(d_min), 9)
+
+
+if __name__ == "__main__":
+    solution = euler_816()
+    print("e816.py: " + str(solution))
+    assert(solution == 20.880613018)

python/e836.py

@@ -0,0 +1,27 @@
+import re
+
+
+TEXT = """
+<p>Let $A$ be an <b>affine plane</b> over a <b>radically integral local field</b> $F$ with residual characteristic $p$.</p>
+
+<p>We consider an <b>open oriented line section</b> $U$ of $A$ with normalized Haar measure $m$.</p>
+
+<p>Define $f(m, p)$ as the maximal possible discriminant of the <b>jacobian</b> associated to the <b>orthogonal kernel embedding</b> of $U$ <span style="white-space:nowrap;">into $A$.</span></p>
+
+<p>Find $f(20230401, 57)$. Give as your answer the concatenation of the first letters of each bolded word.</p>
+"""
+
+
+def euler_836():
+    r = ''
+    for m in re.findall(r'<b>(.*?)</b>', TEXT):
+        words = m.split(' ')
+        for word in words:
+            r += word[0]
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_836()
+    print("e836.py: " + str(solution))
+    assert(solution == 'aprilfoolsjoke')

python/e853.py

@@ -0,0 +1,40 @@
+def fib():
+    a, b = 1, 1
+    yield a
+    yield b
+    while True:
+        a, b = b, a + b
+        yield b
+
+
+def pisano_period(fibs, n):
+    for period in range(2, 130):
+        if fibs[0] % n == fibs[period] % n and fibs[1] % n == fibs[period + 1] % n:
+            return period
+    return None
+
+
+def euler_853():
+    r = 0
+    fs = fib()
+    fibs = [next(fs) for _ in range(400)]
+    upper = 10**9
+    pi_target = 120
+    for n in range(2, upper + 1):
+        # Check if pi_target is a period.
+        if fibs[0] % n == fibs[pi_target] % n and fibs[1] % n == fibs[pi_target + 1] % n:
+            # Make sure that no lower period than pi_target exits.
+            if pisano_period(fibs, n) == pi_target:
+                r += n
+                # print(n, prime_factors(n))
+                # I have noticed that the highest prime factor is always 61 or
+                # 2521, so we could probably also get the solution by
+                # enumerating the prime factors that yield a sum smaller than
+                # upper and have 61 or 2521 as their last factor.
+    return r
+
+
+if __name__ == "__main__":
+    solution = euler_853()
+    print("e853.py: " + str(solution))
+    assert solution == 44511058204

python/lib_bitarray.py

@@ -0,0 +1,39 @@
+class Bitarray:
+    def __init__(self, num_bits):
+        self.word_size = 32
+        self.word_mask = 2**self.word_size - 1
+        num_words = num_bits // self.word_size
+        if num_bits % self.word_size != 0:
+            num_words += 1
+        self.num_bits = num_bits
+        self.num_words = num_words
+        self.array = [0] * num_words
+
+    def setall(self, target: bool):
+        if target:
+            self.array = [self.word_mask] * self.num_words
+        else:
+            self.array = [0] * self.num_words
+
+    def set(self, bit, target: bool):
+        assert(bit < self.num_bits)
+        word_index = bit // self.word_size
+        bit_index = bit % self.word_size
+        value = self.array[word_index]
+        if target:
+            value = (value | (1 << bit_index))
+        else:
+            value = (value & (~(1 << bit_index)))
+        self.array[word_index] = value
+
+    def get(self, bit) -> bool:
+        assert(bit < self.num_bits)
+        word_index = bit // self.word_size
+        bit_index = bit % self.word_size
+        return bool(self.array[word_index] & (1 << bit_index))
+
+    def __getitem__(self, bit: int) -> bool:
+        return self.get(bit)
+
+    def __setitem__(self, bit: int, value: bool):
+        self.set(bit, value)

python/lib_create_templates.py

@@ -1,7 +1,6 @@
-FROM = 206
+import sys
-TILL = 206
 
-template = """
+TEMPLATE = """
 def euler_XXX():
     return 0
 
@@ -9,17 +8,22 @@ def euler_XXX():
 if __name__ == "__main__":
     solution = euler_XXX()
     print("eXXX.py: " + str(solution))
-    assert(solution == 0)
+    # assert(solution == 0)
+
 """
 
 
-for i in range(FROM, TILL + 1):
+def main():
-    if i < 100:
+    try:
-        e = "0" + str(i)
+        e = sys.argv[1]
-    else:
+    except IndexError:
-        e = str(i)
+        print("Provide Euler problem number for which to create template.")
+        return
     filename = "e" + e + ".py"
     with open(filename, "w") as f:
-        s = template.replace("XXX", e)
+        s = TEMPLATE.replace("XXX", e)
         f.write(s)
 
+
+if __name__ == "__main__":
+    main()

python/lib_misc.py

@@ -218,3 +218,15 @@ def is_permutation(n, p):
     for p_n in str(p):
         digit_counts_p[int(p_n)] += 1
     return digit_counts_n == digit_counts_p
+
+
+def cache(f):
+    cache = {}
+    def func_cached(*args):
+        if args in cache:
+            return cache[args]
+        r = f(*args)
+        cache[args] = r
+        return r
+    return func_cached
+

python/lib_prime.py

@@ -1,4 +1,5 @@
 from functools import lru_cache
+from lib_bitarray import Bitarray
 try:
     from lib_misc import get_item_counts
     from lib_misc import product
@@ -65,12 +66,42 @@ def is_prime(n):
     return True
 
 
+def is_prime_rabin_miller(number):
+    """ Rabin-Miller Primality Test """
+    witnesses = [2, 3, 5, 7, 11, 13, 17, 23, 29, 31, 37, 41, 43, 47, 53]
+
+    if number < 2:
+        return False
+
+    if number in witnesses:
+        return True
+
+    if number % 6 not in [1,5]:
+        return False
+
+    r, s = 1, number - 1
+    while s % 2 == 0:
+        s //= 2
+        r += 1
+
+    for witness in witnesses:
+        remainder = pow(witness, s, number)
+        if remainder == 1:
+            continue
+        for _ in range(1, r):
+            if remainder == number - 1:
+                break
+            remainder = pow(remainder, 2, number)
+        else:
+            return False
+    return True
+
+
 def prime_nth(n):
     """Returns the nth prime number. The first number
     is 1 indexed, i.e. n = 1 -> 2, n = 2 -> 3, etc.
 
     :param n:
-
     """
     if n == 1:
         return 2
@@ -97,12 +128,8 @@ def primes(n_max):
 
     :param n_max:
     """
-    try:
+    b = Bitarray(n_max)
-        import bitarray
-        b = bitarray.bitarray(n_max)
     b.setall(True)
-    except ModuleNotFoundError:
-        b = [True for _ in range(n_max)]
     n = 1
     b[n - 1] = False
     while n * n <= n_max:

txt/e105.txt

@@ -0,0 +1,100 @@
+81,88,75,42,87,84,86,65
+157,150,164,119,79,159,161,139,158
+673,465,569,603,629,592,584,300,601,599,600
+90,85,83,84,65,87,76,46
+165,168,169,190,162,85,176,167,127
+224,275,278,249,277,279,289,295,139
+354,370,362,384,359,324,360,180,350,270
+599,595,557,298,448,596,577,667,597,588,602
+175,199,137,88,187,173,168,171,174
+93,187,196,144,185,178,186,202,182
+157,155,81,158,119,176,152,167,159
+184,165,159,166,163,167,174,124,83
+1211,1212,1287,605,1208,1189,1060,1216,1243,1200,908,1210
+339,299,153,305,282,304,313,306,302,228
+94,104,63,112,80,84,93,96
+41,88,82,85,61,74,83,81
+90,67,84,83,82,97,86,41
+299,303,151,301,291,302,307,377,333,280
+55,40,48,44,25,42,41
+1038,1188,1255,1184,594,890,1173,1151,1186,1203,1187,1195
+76,132,133,144,135,99,128,154
+77,46,108,81,85,84,93,83
+624,596,391,605,529,610,607,568,604,603,453
+83,167,166,189,163,174,160,165,133
+308,281,389,292,346,303,302,304,300,173
+593,1151,1187,1184,890,1040,1173,1186,1195,1255,1188,1203
+68,46,64,33,60,58,65
+65,43,88,87,86,99,93,90
+83,78,107,48,84,87,96,85
+1188,1173,1256,1038,1187,1151,890,1186,1184,1203,594,1195
+302,324,280,296,294,160,367,298,264,299
+521,760,682,687,646,664,342,698,692,686,672
+56,95,86,97,96,89,108,120
+344,356,262,343,340,382,337,175,361,330
+47,44,42,27,41,40,37
+139,155,161,158,118,166,154,156,78
+118,157,164,158,161,79,139,150,159
+299,292,371,150,300,301,281,303,306,262
+85,77,86,84,44,88,91,67
+88,85,84,44,65,91,76,86
+138,141,127,96,136,154,135,76
+292,308,302,346,300,324,304,305,238,166
+354,342,341,257,348,343,345,321,170,301
+84,178,168,167,131,170,193,166,162
+686,701,706,673,694,687,652,343,683,606,518
+295,293,301,367,296,279,297,263,323,159
+1038,1184,593,890,1188,1173,1187,1186,1195,1150,1203,1255
+343,364,388,402,191,383,382,385,288,374
+1187,1036,1183,591,1184,1175,888,1197,1182,1219,1115,1167
+151,291,307,303,345,238,299,323,301,302
+140,151,143,138,99,69,131,137
+29,44,42,59,41,36,40
+348,329,343,344,338,315,169,359,375,271
+48,39,34,37,50,40,41
+593,445,595,558,662,602,591,297,610,580,594
+686,651,681,342,541,687,691,707,604,675,699
+180,99,189,166,194,188,144,187,199
+321,349,335,343,377,176,265,356,344,332
+1151,1255,1195,1173,1184,1186,1188,1187,1203,593,1038,891
+90,88,100,83,62,113,80,89
+308,303,238,300,151,304,324,293,346,302
+59,38,50,41,42,35,40
+352,366,174,355,344,265,343,310,338,331
+91,89,93,90,117,85,60,106
+146,186,166,175,202,92,184,183,189
+82,67,96,44,80,79,88,76
+54,50,58,66,31,61,64
+343,266,344,172,308,336,364,350,359,333
+88,49,87,82,90,98,86,115
+20,47,49,51,54,48,40
+159,79,177,158,157,152,155,167,118
+1219,1183,1182,1115,1035,1186,591,1197,1167,887,1184,1175
+611,518,693,343,704,667,686,682,677,687,725
+607,599,634,305,677,604,603,580,452,605,591
+682,686,635,675,692,730,687,342,517,658,695
+662,296,573,598,592,584,553,593,595,443,591
+180,185,186,199,187,210,93,177,149
+197,136,179,185,156,182,180,178,99
+271,298,218,279,285,282,280,238,140
+1187,1151,890,593,1194,1188,1184,1173,1038,1186,1255,1203
+169,161,177,192,130,165,84,167,168
+50,42,43,41,66,39,36
+590,669,604,579,448,599,560,299,601,597,598
+174,191,206,179,184,142,177,180,90
+298,299,297,306,164,285,374,269,329,295
+181,172,162,138,170,195,86,169,168
+1184,1197,591,1182,1186,889,1167,1219,1183,1033,1115,1175
+644,695,691,679,667,687,340,681,770,686,517
+606,524,592,576,628,593,591,584,296,444,595
+94,127,154,138,135,74,136,141
+179,168,172,178,177,89,198,186,137
+302,299,291,300,298,149,260,305,280,370
+678,517,670,686,682,768,687,648,342,692,702
+302,290,304,376,333,303,306,298,279,153
+95,102,109,54,96,75,85,97
+150,154,146,78,152,151,162,173,119
+150,143,157,152,184,112,154,151,132
+36,41,54,40,25,44,42
+37,48,34,59,39,41,40
+681,603,638,611,584,303,454,607,606,605,596

txt/e107.txt

@@ -0,0 +1,40 @@
+-,-,-,427,668,495,377,678,-,177,-,-,870,-,869,624,300,609,131,-,251,-,-,-,856,221,514,-,591,762,182,56,-,884,412,273,636,-,-,774
+-,-,262,-,-,508,472,799,-,956,578,363,940,143,-,162,122,910,-,729,802,941,922,573,531,539,667,607,-,920,-,-,315,649,937,-,185,102,636,289
+-,262,-,-,926,-,958,158,647,47,621,264,81,-,402,813,649,386,252,391,264,637,349,-,-,-,108,-,727,225,578,699,-,898,294,-,575,168,432,833
+427,-,-,-,366,-,-,635,-,32,962,468,893,854,718,427,448,916,258,-,760,909,529,311,404,-,-,588,680,875,-,615,-,409,758,221,-,-,76,257
+668,-,926,366,-,-,-,250,268,-,503,944,-,677,-,727,793,457,981,191,-,-,-,351,969,925,987,328,282,589,-,873,477,-,-,19,450,-,-,-
+495,508,-,-,-,-,-,765,711,819,305,302,926,-,-,582,-,861,-,683,293,-,-,66,-,27,-,-,290,-,786,-,554,817,33,-,54,506,386,381
+377,472,958,-,-,-,-,-,-,120,42,-,134,219,457,639,538,374,-,-,-,966,-,-,-,-,-,449,120,797,358,232,550,-,305,997,662,744,686,239
+678,799,158,635,250,765,-,-,-,35,-,106,385,652,160,-,890,812,605,953,-,-,-,79,-,712,613,312,452,-,978,900,-,901,-,-,225,533,770,722
+-,-,647,-,268,711,-,-,-,283,-,172,-,663,236,36,403,286,986,-,-,810,761,574,53,793,-,-,777,330,936,883,286,-,174,-,-,-,828,711
+177,956,47,32,-,819,120,35,283,-,50,-,565,36,767,684,344,489,565,-,-,103,810,463,733,665,494,644,863,25,385,-,342,470,-,-,-,730,582,468
+-,578,621,962,503,305,42,-,-,50,-,155,519,-,-,256,990,801,154,53,474,650,402,-,-,-,966,-,-,406,989,772,932,7,-,823,391,-,-,933
+-,363,264,468,944,302,-,106,172,-,155,-,-,-,380,438,-,41,266,-,-,104,867,609,-,270,861,-,-,165,-,675,250,686,995,366,191,-,433,-
+870,940,81,893,-,926,134,385,-,565,519,-,-,313,851,-,-,-,248,220,-,826,359,829,-,234,198,145,409,68,359,-,814,218,186,-,-,929,203,-
+-,143,-,854,677,-,219,652,663,36,-,-,313,-,132,-,433,598,-,-,168,870,-,-,-,128,437,-,383,364,966,227,-,-,807,993,-,-,526,17
+869,-,402,718,-,-,457,160,236,767,-,380,851,132,-,-,596,903,613,730,-,261,-,142,379,885,89,-,848,258,112,-,900,-,-,818,639,268,600,-
+624,162,813,427,727,582,639,-,36,684,256,438,-,-,-,-,539,379,664,561,542,-,999,585,-,-,321,398,-,-,950,68,193,-,697,-,390,588,848,-
+300,122,649,448,793,-,538,890,403,344,990,-,-,433,596,539,-,-,73,-,318,-,-,500,-,968,-,291,-,-,765,196,504,757,-,542,-,395,227,148
+609,910,386,916,457,861,374,812,286,489,801,41,-,598,903,379,-,-,-,946,136,399,-,941,707,156,757,258,251,-,807,-,-,-,461,501,-,-,616,-
+131,-,252,258,981,-,-,605,986,565,154,266,248,-,613,664,73,-,-,686,-,-,575,627,817,282,-,698,398,222,-,649,-,-,-,-,-,654,-,-
+-,729,391,-,191,683,-,953,-,-,53,-,220,-,730,561,-,946,686,-,-,389,729,553,304,703,455,857,260,-,991,182,351,477,867,-,-,889,217,853
+251,802,264,760,-,293,-,-,-,-,474,-,-,168,-,542,318,136,-,-,-,-,392,-,-,-,267,407,27,651,80,927,-,974,977,-,-,457,117,-
+-,941,637,909,-,-,966,-,810,103,650,104,826,870,261,-,-,399,-,389,-,-,-,202,-,-,-,-,867,140,403,962,785,-,511,-,1,-,707,-
+-,922,349,529,-,-,-,-,761,810,402,867,359,-,-,999,-,-,575,729,392,-,-,388,939,-,959,-,83,463,361,-,-,512,931,-,224,690,369,-
+-,573,-,311,351,66,-,79,574,463,-,609,829,-,142,585,500,941,627,553,-,202,388,-,164,829,-,620,523,639,936,-,-,490,-,695,-,505,109,-
+856,531,-,404,969,-,-,-,53,733,-,-,-,-,379,-,-,707,817,304,-,-,939,164,-,-,616,716,728,-,889,349,-,963,150,447,-,292,586,264
+221,539,-,-,925,27,-,712,793,665,-,270,234,128,885,-,968,156,282,703,-,-,-,829,-,-,-,822,-,-,-,736,576,-,697,946,443,-,205,194
+514,667,108,-,987,-,-,613,-,494,966,861,198,437,89,321,-,757,-,455,267,-,959,-,616,-,-,-,349,156,339,-,102,790,359,-,439,938,809,260
+-,607,-,588,328,-,449,312,-,644,-,-,145,-,-,398,291,258,698,857,407,-,-,620,716,822,-,-,293,486,943,-,779,-,6,880,116,775,-,947
+591,-,727,680,282,290,120,452,777,863,-,-,409,383,848,-,-,251,398,260,27,867,83,523,728,-,349,293,-,212,684,505,341,384,9,992,507,48,-,-
+762,920,225,875,589,-,797,-,330,25,406,165,68,364,258,-,-,-,222,-,651,140,463,639,-,-,156,486,212,-,-,349,723,-,-,186,-,36,240,752
+182,-,578,-,-,786,358,978,936,385,989,-,359,966,112,950,765,807,-,991,80,403,361,936,889,-,339,943,684,-,-,965,302,676,725,-,327,134,-,147
+56,-,699,615,873,-,232,900,883,-,772,675,-,227,-,68,196,-,649,182,927,962,-,-,349,736,-,-,505,349,965,-,474,178,833,-,-,555,853,-
+-,315,-,-,477,554,550,-,286,342,932,250,814,-,900,193,504,-,-,351,-,785,-,-,-,576,102,779,341,723,302,474,-,689,-,-,-,451,-,-
+884,649,898,409,-,817,-,901,-,470,7,686,218,-,-,-,757,-,-,477,974,-,512,490,963,-,790,-,384,-,676,178,689,-,245,596,445,-,-,343
+412,937,294,758,-,33,305,-,174,-,-,995,186,807,-,697,-,461,-,867,977,511,931,-,150,697,359,6,9,-,725,833,-,245,-,949,-,270,-,112
+273,-,-,221,19,-,997,-,-,-,823,366,-,993,818,-,542,501,-,-,-,-,-,695,447,946,-,880,992,186,-,-,-,596,949,-,91,-,768,273
+636,185,575,-,450,54,662,225,-,-,391,191,-,-,639,390,-,-,-,-,-,1,224,-,-,443,439,116,507,-,327,-,-,445,-,91,-,248,-,344
+-,102,168,-,-,506,744,533,-,730,-,-,929,-,268,588,395,-,654,889,457,-,690,505,292,-,938,775,48,36,134,555,451,-,270,-,248,-,371,680
+-,636,432,76,-,386,686,770,828,582,-,433,203,526,600,848,227,616,-,217,117,707,369,109,586,205,809,-,-,240,-,853,-,-,-,768,-,371,-,540
+774,289,833,257,-,381,239,722,711,468,933,-,-,17,-,-,148,-,-,853,-,-,-,-,264,194,260,947,-,752,147,-,-,343,112,273,344,680,540,-