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/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/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/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/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

@@ -22,7 +22,7 @@ def euler_205():
         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)
+        four_won_count += four_count * current_six_count
 
     p = four_won_count / all_combinations
     return f"{round(p, 7)}"
@@ -32,4 +32,3 @@ 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/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/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/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

@@ -102,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