Solve Euler 174

python/e090.py

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

python/e108.py

@@ -1,4 +1,6 @@
 from fractions import Fraction
+from lib_misc import proper_divisors
+from math import ceil
 
 
 def get_distinct_solutions(n):
@@ -16,11 +18,16 @@ def get_distinct_solutions(n):
     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_solutions(n)
+        d = get_distinct_solutions2(n)
         if d > d_max:
             # print("n={} d={} delta={}".format(n, d, n - n_prev))
             n_prev = n

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

@@ -35,7 +35,7 @@ def euler_111():
             for p in get_permutations(base, 2):
                 if is_prime(p):
                     subresult += p
-        assert subresult > 0, "More than to permutations required to yield prime"
+        assert subresult > 0, "More than two permutations required to yield prime"
         result += subresult
     return result
 

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/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,13 +2,14 @@ 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
+        x = (x + (apositiveint // x)) // 2
-         if x in seen: return False
+        if x in seen:
-         seen.add(x)
+            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/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/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/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_prime.py

@@ -70,6 +70,9 @@ 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
 
@@ -99,7 +102,6 @@ def prime_nth(n):
     is 1 indexed, i.e. n = 1 -> 2, n = 2 -> 3, etc.
 
     :param n:
-
     """
     if n == 1:
         return 2