Solve problem 214

python/e214.py

@@ -0,0 +1,51 @@
+from lib_prime import primes
+from functools import cache
+from math import gcd
+
+
+def get_phi_table(n: int, ps: list[int]) -> list[int]:
+    tbl = [-1] * (n + 2)
+    tbl[1] = 1
+    queue = [(1, 1, 0)]
+
+    while queue:
+        x, phi, pi = queue.pop()
+        for pin in range(pi, len(ps)):
+            p = ps[pin]
+            tbl[p] = p - 1
+            xn = x * p
+            if xn > n:
+                break
+            d = gcd(x, p)
+            phid = tbl[d]
+            assert phid != -1
+            phin = (phi * (p - 1) * d) // phid
+            tbl[xn] = phin
+            queue.append((xn, phin, pin))
+    return tbl
+
+
+def euler_214():
+    n = 40 * 10**6
+    ps = primes(n)
+    tbl = get_phi_table(n, ps)
+
+    @cache
+    def chain(n: int) -> int:
+        if n == 1:
+            return 1
+        n = tbl[n]
+        return 1 + chain(n)
+
+    s = 0
+    for i in range(len(ps)):
+        p = ps[i]
+        if chain(p) == 25:
+            s += p
+    return s
+
+
+if __name__ == "__main__":
+    solution = euler_214()
+    print(f"e214.py: {solution}")
+    assert solution == 1677366278943

python/e231.py

@@ -0,0 +1,25 @@
+from lib_prime import prime_factors_count
+from collections import defaultdict
+
+
+def sum_factors(n: int, k: int) -> int:
+    factors_num = defaultdict(int)
+    for i in range(n - k + 1, n + 1):
+        for factor, count in prime_factors_count(i).items():
+            factors_num[factor] += count
+    for i in range(2, k + 1):
+        for factor, count in prime_factors_count(i).items():
+            factors_num[factor] -= count
+    return sum(factor * count for factor, count in factors_num.items())
+
+
+def euler_231():
+    assert sum_factors(10, 3) == 14
+    return sum_factors(20 * 10**6, 15 * 10**6)
+
+
+if __name__ == "__main__":
+    solution = euler_231()
+    print(f"e231.py: {solution}")
+    # assert(solution == 0)
+

python/e265.py

@@ -0,0 +1,47 @@
+def is_valid(circle_with_seen, n):
+    circle, seen = circle_with_seen
+    if len(circle) < n:
+        return True
+
+    end = tuple(circle[-n:])
+    if end in seen:
+        return False
+    seen.add(end)
+    return True
+
+
+def overlap_valid(circle_with_seen, n):
+    circle, seen = circle_with_seen
+    for i in range(1, n):
+        end = tuple(circle[-i:] + circle[: n - i])
+        if end in seen:
+            return False
+        seen.add(end)
+    return True
+
+
+def s(n):
+    length = 2**n
+    start = [0 for _ in range(n)]
+    xs = [(start, set([tuple(start)]))]
+    for _ in range(length - n):
+        nxs = []
+        for circle, seen in xs:
+            for d in [0, 1]:
+                nx = (circle + [d], seen.copy())
+                if is_valid(nx, n):
+                    nxs.append(nx)
+        xs = nxs
+    circles = [x[0] for x in xs if overlap_valid(x, n)]
+    return sum(int("".join(map(str, circle)), 2) for circle in circles)
+
+
+def euler_265():
+    assert s(3) == 52
+    return s(5)
+
+
+if __name__ == "__main__":
+    solution = euler_265()
+    print(f"e265.py: {solution}")
+    assert solution == 209110240768

python/e274.py

@@ -0,0 +1,57 @@
+from math import gcd
+from lib_prime import primes
+from lib_misc import modinv
+
+
+def find_m(p: int) -> int:
+    rs = []
+    test_n = p * 1337
+    for m in range(1, p):
+        test_f = f(test_n, m)
+        if test_f % p == 0:
+            rs.append(m)
+    assert len(rs) == 1
+    return rs[0]
+
+
+def calc_m(p: int) -> int:
+    """
+    (n // 10 + (n % 10) * m) % p == 0
+    a + b * m == 0 (mod p)
+        b * m == -a (mod p)
+            m == -a * b**-1 (mod p)
+    """
+    n = p * 1337
+    a = n // 10
+    b = (n % 10)
+    m = -a * modinv(b, p)
+    m %= p
+    return m
+
+
+def f(n: int, m: int = 34) -> int:
+    left = n // 10
+    right = n % 10 * m
+    r = left + right
+    return r
+
+
+def euler_274():
+    ps = primes(10**7)
+    sum_m = 0
+    for p in ps:
+        if gcd(p, 10) != 1:
+            continue
+        m = calc_m(p)
+        if p < 1000:
+            m_naiv = find_m(p)
+            assert m == m_naiv
+        sum_m += m
+    return sum_m
+
+
+if __name__ == "__main__":
+    solution = euler_274()
+    print(f"e274.py: {solution}")
+    assert solution == 1601912348822
+

python/e297.py

@@ -0,0 +1,116 @@
+from functools import cache
+from typing import Optional
+
+
+def fib(max_n: int) -> list[int]:
+    a, b = 1, 2
+    xs = [a]
+    while b <= max_n:
+        xs.append(b)
+        a, b = b, a + b
+    return xs
+
+
+FIBS = fib(10**18)
+
+
+@cache
+def zgo(rem: int, i: int, count: int) -> Optional[int]:
+    if rem == 0:
+        return count
+    if rem < 0 or rem < FIBS[i]:
+        return None
+
+    r1 = zgo(rem - FIBS[i], i + 2, count + 1)
+    if r1:
+        return r1
+    return zgo(rem, i + 1, count)
+
+
+@cache
+def z(n: int) -> int:
+    return zgo(n, 0, 0) or 0
+
+
+def zsum(n: int) -> int:
+    return sum(z(i) for i in range(1, n))
+
+
+def zxsgo(rem: int, i: int, xs: list[int]) -> Optional[list]:
+    if rem == 0:
+        return xs
+    if rem < 0 or rem < FIBS[i]:
+        return None
+
+    r1 = zxsgo(rem - FIBS[i], i + 2, xs + [FIBS[i]])
+    if r1:
+        return r1
+    return zxsgo(rem, i + 1, xs)
+
+
+def zxs(n: int) -> Optional[list]:
+    return zxsgo(n, 0, [])
+
+
+@cache
+def zsumr(n: int) -> int:
+    if n == 1:
+        return 0
+    if n == 2:
+        return 1
+    if n in FIBS:
+        fib_i = FIBS.index(n)
+        a, b = FIBS[fib_i - 2], FIBS[fib_i - 1]
+        # print(f"zsumr({a})={zsumr(a)} zsumr({b})={zsumr(b)} {b - a}")
+        res = 0
+
+        # Zeckendorf sum up the previous fibonacci
+        res += zsumr(b)
+
+        # Zeckendorf sum up to the pervious previous fibonacci. This works because
+        # if we look at the Zeckendorf representations between two fibonaccis it is
+        # always a combination of the previous fibonacci plus the representations
+        # starting from zero. Since we cannot use the previous fibonacci because of
+        # the *nonconsecutive* requirement, it has to be the second previous.
+        res += zsumr(a)
+
+        # However, then for all the zsumr(a) terms we considered, that fibonacci
+        # a was not part of the count, so we have to add it:
+        res += n - b
+
+        return res
+
+    for fib_i, fib in enumerate(FIBS):
+        if fib > n:
+            break
+    else:
+        assert False
+    fib_i -= 1
+    b = FIBS[fib_i]  # fib before n
+
+    res = 0
+    res += zsumr(b)  # Zeckendorf up to previous fib
+    res += zsumr(n - b)  # Zeckendorf from b to n without last term
+    res += n - b  # account for last term
+    return res
+
+
+def euler_297():
+    assert zsumr(3) == 2
+    assert zsumr(5) == 5
+    assert zsumr(13) == 20
+    assert zsumr(89) == 235
+    assert zsumr(100) == 261
+    assert zsumr(10**6) == 7894453
+
+    assert z(1) == 1
+    assert z(5) == 1
+    assert z(14) == 2
+    assert z(100) == 3
+    return zsumr(10**17)
+
+
+if __name__ == "__main__":
+    solution = euler_297()
+    print(f"e297.py: {solution}")
+    assert solution == 2252639041804718029

python/e313.py

@@ -0,0 +1,99 @@
+from lib_prime import primes, is_prime
+from dataclasses import dataclass
+from heapq import heappop, heappush
+from collections import namedtuple
+from collections import defaultdict
+from math import sqrt
+
+P = namedtuple("P", ["x", "y"])
+
+State = namedtuple("State", ["f", "g", "red", "space"])
+
+
+@dataclass(order=True)
+class State:
+    g: int
+    red: P
+    space: P
+
+
+def s(m: int, n: int) -> int:
+    """What am I doing? It was necessary to find the pattern."""
+
+    def h(s: State) -> int:
+        return abs(m - s.space.x) + abs(n - s.space.y)
+
+    def f(s: State) -> int:
+        return s.g + h(s)
+
+    def done(s: State) -> bool:
+        return s.red.x == m - 1 and s.red.y == n - 1
+
+    # assert m < n, then we could just take result times 2
+    init = State(0, P(0, 0), P(m - 1, n - 1))
+    states = [(0, init)]
+    g_costs = {(init.red, init.space): init.g}
+
+    while states:
+        _, s = heappop(states)
+        if done(s):
+            return s.g
+
+        for dx, dy in [(-1, 0), (0, 1), (1, 0), (0, -1)]:
+            nx, ny = s.space.x + dx, s.space.y + dy
+            if not (0 <= nx < m and 0 <= ny < n):
+                continue
+            nstate = State(s.g + 1, s.red, P(nx, ny))
+            if s.red.x == nx and s.red.y == ny:
+                nstate.red = P(s.space.x, s.space.y)
+            if (
+                not (nstate.red, nstate.space) in g_costs
+                or g_costs[(nstate.red, nstate.space)] > nstate.g
+            ):
+                g_costs[nstate.red, nstate.space] = nstate.g
+                heappush(states, (f(nstate), nstate))
+    return -1
+
+
+def solve_naiv(p_max: int) -> int:
+    t = 0
+    p2_max = p_max * p_max
+    ps = primes(p_max)
+    p2s = set([p * p for p in ps])
+    for m in range(2, p2_max):
+        for n in range(m, p2_max):
+            r = s(m, n)
+            if r in p2s:
+                t += 2
+            if r > p2_max:
+                break
+    return t
+
+
+def solve(p_max: int) -> int:
+    t = 2
+    r_base = 0
+    delta = 1
+    three = 3
+    for p in range(3, p_max + 1, 2):
+        if is_prime(p):
+            t += r_base * 2
+        r_base += delta
+        three -= 1
+        if three == 0:
+            delta += 1
+            three = 3
+    return t
+
+
+def euler_313():
+    assert solve_naiv(5) == 4
+    assert solve(20) == 86
+    assert solve(100) == 5482
+    return solve(10**6)
+
+
+if __name__ == "__main__":
+    solution = euler_313()
+    print(f"e313.py: {solution}")
+    # assert solution == 2057774861813004

python/e315.py

@@ -0,0 +1,111 @@
+from lib_prime import is_prime_rabin_miller
+from typing import Optional
+from functools import cache
+
+
+def compute_deltas() -> list[list[int]]:
+    S = [
+        set([0, 1, 2, 4, 5, 6]),  # 0
+        set([2, 5]),  # 1
+        set([0, 2, 3, 4, 6]),  # 2
+        set([0, 2, 3, 5, 6]),  # 3
+        set([1, 2, 3, 5]),  # 4
+        set([0, 1, 3, 5, 6]),  # 5
+        set([0, 1, 3, 4, 5, 6]),  # 6
+        set([0, 1, 2, 5]),  # 7
+        set([0, 1, 2, 3, 4, 5, 6]),  # 8
+        set([0, 1, 2, 3, 5, 6]),  # 9
+    ]
+
+    xs = [[0] * 10 for _ in range(10)]
+    for a in range(10):
+        for b in range(10):
+            xs[a][b] = len(S[a] - S[b]) + len(S[b] - S[a])
+    return xs
+
+
+NUM_TO_SEGMENTS = [6, 2, 5, 5, 4, 5, 6, 4, 7, 6]
+DELTA_SEGMENTS = compute_deltas()
+
+
+@cache
+def digital_root(n: int) -> int:
+    root = 0
+    while n != 0:
+        d = n % 10
+        n = n // 10
+        root += d
+    return root
+
+
+@cache
+def total_segments(n: int) -> int:
+    count = 0
+    n_next = 0
+    while n != 0:
+        d = n % 10
+        n = n // 10
+        count += NUM_TO_SEGMENTS[d]
+        n_next += d
+    return count
+
+
+@cache
+def transition_segments(a: int, b: int) -> int:
+    count = 0
+    while a != 0 or b != 0:
+        assert a >= b
+        da, db = a % 10, b % 10
+        if b == 0:
+            # We are just turning off the remaining digits
+            count += NUM_TO_SEGMENTS[da]
+        else:
+            count += DELTA_SEGMENTS[da][db]
+        a, b = a // 10, b // 10
+    return count
+
+
+@cache
+def sam_count(n: int) -> int:
+    count = total_segments(n)
+    count *= 2  # turn off n
+    n_next = digital_root(n)
+    if n == n_next:
+        return count
+    return count + sam_count(n_next)
+
+
+@cache
+def max_count(n_current: Optional[int], n: int) -> int:
+    if n_current is None:
+        # First number, only turn number on
+        count = total_segments(n)
+        n_next = digital_root(n)
+        assert n != n_next
+        return count + max_count(n, n_next)
+    elif n_current < 10:
+        # Last number, only turn number off
+        assert n_current == n
+        return total_segments(n)
+
+    # This is the case where we have to compute the actual delta
+    count = transition_segments(n_current, n)
+    n_next = digital_root(n)
+    return count + max_count(n, n_next)
+
+
+def euler_315():
+    assert sam_count(137) == 40
+    assert max_count(None, 137) == 30
+
+    t = 0
+    for n in range(10**7, 2 * 10**7):
+        if is_prime_rabin_miller(n):
+            t += sam_count(n) - max_count(None, n)
+    return t
+
+
+if __name__ == "__main__":
+    solution = euler_315()
+    print(f"e315.py: {solution}")
+    assert solution == 13625242

python/e323.py

@@ -0,0 +1,77 @@
+from random import randint
+from math import factorial
+from collections import defaultdict
+from fractions import Fraction
+
+
+def binomial(n: int, k: int) -> Fraction:
+    return Fraction(factorial(n), (factorial(k) * factorial(n - k)))
+
+
+def f(k: int, n: int) -> Fraction:
+    # Probability that k out of n trials succeed with rate p per trial.
+    p = Fraction(1, 2)
+    return binomial(n, k) * p**k * (1 - p) ** (n - k)
+
+
+def count(exp_n: Fraction, step: int, dist) -> Fraction:
+    step += 1
+
+    # For each number `n` of remaining "ones" calculate the probability that
+    # 0..`n` ones are remanining for that probability. We can then calculate a
+    # new probability mass by multiplying the current with the new probability
+    # for a given `n`.
+    ndist = defaultdict(Fraction)
+    for n, p in dist:
+        for i in range(n + 1):
+            np = f(i, n)
+            ndist[n - i] += p * np
+
+    # Probability that we reach `0` at this step counts to total expected
+    # number of steps.
+    exp_n += ndist[0] * step
+    del ndist[0]
+
+    # Drop probabilities small enough that they seem englible for the precision
+    # we need which is 10 digits after the decimal point.
+    dist = []
+    for n, p in ndist.items():
+        if p > Fraction(1, 10**15):
+            dist.append((n, p))
+    if dist:
+        return count(exp_n, step, dist)
+    return exp_n
+
+
+def stackoverflow():
+    # I found this expression when googling (geminiing) an expression for this
+    # type of problem. However, it did not yield the correct solution because I
+    # did not use Fractions. I then came up with `count` and was somehow clever
+    # enough to realize that I would have to use Fractions and after I came
+    # back I updated stackoverflow accordingly just to see if I had implemented
+    # the expression correctly.
+
+    # https://mathoverflow.net/questions/41604/what-is-the-expected-maximum-out-of-a-sample-size-n-from-a-geometric-distribut
+    t = Fraction()
+    p = Fraction(1, 2)
+    q = 1 - p
+    n = 32
+    for i in range(1, n + 1):
+        c = binomial(n, i)
+        c *= (-1) ** (i + 1)
+        c *= Fraction(1, 1 - q**i)
+        t += c
+    return float(t)
+
+
+def euler_323():
+    n_stackoverflow = stackoverflow()
+    n = count(Fraction(), 0, [(32, Fraction(1, 1))])
+    assert round(float(n_stackoverflow), 10) == round(float(n), 10)
+    return round(float(n), 10)
+
+
+if __name__ == "__main__":
+    solution = euler_323()
+    print(f"e323.py: {solution}")
+    assert solution == 6.3551758451

python/e336.py

@@ -0,0 +1,65 @@
+from itertools import permutations
+from string import ascii_uppercase
+from functools import cache
+
+
+@cache
+def count(xst: tuple) -> int:
+    expected = 0
+    for i in range(len(xst)):
+        if xst[i] == expected:
+            expected += 1
+        else:
+            # Carriage i is not the correct one
+            break
+    else:
+        # Train is fully sorted
+        return 0
+
+    # Find out where Simple Simon would split next; there are two cases.
+    xs = list(xst)
+
+    # If the next expected carriage is at the end of the list, Simon can just
+    # reverse from where we currently are.
+    if xs[-1] == expected:
+        xs[i:] = xs[i:][::-1]
+        return 1 + count(tuple(xs))
+
+    for j in range(i + 1, len(xst)):
+        if xst[j] == expected:
+            xs[j:] = xs[j:][::-1]
+            return 1 + count(tuple(xs))
+
+    assert False, "should not get here"
+
+
+def get_maximix_arrangements(n: int) -> list:
+    x = [i for i in range(n)]
+    xs = permutations(x)
+    max_xs = []
+    max_xs_count = 0
+    for x in xs:
+        c = count(x)
+        if c == max_xs_count:
+            max_xs.append(x)
+        elif c > max_xs_count:
+            max_xs = [x]
+            max_xs_count = c
+    return sorted(map(to_ascii, max_xs))
+
+
+def to_ascii(x) -> str:
+    return "".join([ascii_uppercase[d] for d in x])
+
+
+def euler_336():
+    xs = get_maximix_arrangements(6)
+    assert len(xs) == 24
+    assert xs[10 - 1] == "DFAECB"
+    return get_maximix_arrangements(11)[2011 - 1]
+
+
+if __name__ == "__main__":
+    solution = euler_336()
+    print(f"e336.py: {solution}")
+    assert solution == "CAGBIHEFJDK"

python/e348.py

@@ -0,0 +1,31 @@
+from functools import cache
+from collections import defaultdict
+
+
+@cache
+def is_palindrome(n: int) -> bool:
+    ns = str(n)
+    return ns == ns[::-1]
+
+
+def euler_348():
+    counts = defaultdict(int)
+
+    for a in range(30000):
+        for b in range(3000):
+            x = a**2 + b**3
+            if is_palindrome(x):
+                counts[x] += 1
+    rs = []
+    for k, v in counts.items():
+        if v == 4:
+            rs.append(k)
+    assert len(rs) == 5
+    return sum(rs)
+
+
+if __name__ == "__main__":
+    solution = euler_348()
+    print(f"e348.py: {solution}")
+    assert(solution == 1004195061)
+

python/e351.py

@@ -0,0 +1,108 @@
+from fractions import Fraction
+
+
+def plot(n):
+    """
+    First iteration with plotting so that I could understand
+    what's going on. Worked well to accomplish this plan.
+
+    1. Visualize till H(10) naively
+    2. See pattern
+    3. Generalize and compute solution
+    """
+    import numpy as np
+    import matplotlib
+    import matplotlib.pyplot as plt
+    from matplotlib.ticker import MultipleLocator
+
+    matplotlib.use("QtAgg")
+
+    _, ax = plt.subplots()
+    ax.set_aspect("equal")
+    dx = np.cos(np.pi / 3)
+    dy = np.sin(np.pi / 3)
+
+    c = 0
+    seen_set = set()
+    seen = []
+    for i in range(n + 1):
+        x, y = 2 * i, 0
+        for _ in range(i):
+            for xs, ys in seen:
+                if x % xs == 0 and x * ys == y * xs:
+                    assert Fraction(y, x) in seen_set
+                    ax.plot(x * dx, y * dy, "x", color="green")
+                    c += 1
+                    break
+            else:
+                ax.plot(x * dx, y * dy, "o", color="blue")
+                if x != 0 or y != 0:
+                    seen_set.add(Fraction(y, x))
+                    seen.append((x, y))
+            x -= 1
+            y += 1
+
+    ax.xaxis.set_major_locator(MultipleLocator(1))
+    ax.yaxis.set_major_locator(MultipleLocator(1))
+    ax.grid(True, alpha=0.2)
+    print(f"h({n}) = {c * 6}")
+    plt.show()
+
+
+def quick(n):
+    """
+    Improved version that uses slopes to identify
+    if a point is visible. Still way too slow for 10**9, though.
+    """
+    c = 0
+    seen_set = set()
+    for i in range(n + 1):
+        x, y = 2 * i, 0
+        for _ in range(i):
+            if Fraction(y, x) in seen_set:
+                c += 1
+            else:
+                if x != 0 or y != 0:
+                    seen_set.add(Fraction(y, x))
+            x -= 1
+            y += 1
+    return c * 6
+
+
+def h(n: int) -> int:
+    """Optimized version based on realization that the
+    totient, i.e., how many numbers don't share a non-one
+    divisor with the respective row can be used to calculate
+    how many points in that row are visible."""
+    totients = [i for i in range(n + 1)]
+
+    r = (n * (n + 1)) // 2  # All points in 1/6 of the hexagon
+    r -= 1  # totient of 1 is 1
+    for i in range(2, n + 1):
+        # Number had no factor so far and is therefore a prime.
+        if totients[i] == i:
+            # That phi(p) = p - 1:
+            totients[i] -= 1
+            # Update totients of other numbers of which i is a factor.
+            # This is the same as * (1 - 1 / i); since i is a factor
+            # //= i will always be accurate.
+            for j in range(i + i, n + 1, i):
+                totients[j] *= i - 1
+                totients[j] //= i
+        r -= totients[i]
+    # We only calculate the points for 1/6 of the grid.
+    return r * 6
+
+
+def euler_351():
+    assert h(5) == 30
+    assert h(10) == 138
+    assert h(1000) == 1177848
+    assert quick(1000) == 1177848
+    return h(100000000)
+
+
+if __name__ == "__main__":
+    solution = euler_351()
+    print(f"e351.py: {solution}")
+    assert solution == 11762187201804552

python/e401.py

@@ -0,0 +1,59 @@
+def divisors(n: int) -> list[int]:
+    xs = set()
+    d = 1
+    while d * d <= n:
+        if n % d == 0:
+            xs.add(d)
+            xs.add(n // d)
+        d += 1
+    return sorted(list(xs))
+
+
+def sigma2(n: int) -> int:
+    return sum(d**2 for d in divisors(n))
+
+
+def sum2(a: int, b: int) -> int:
+    return sum(x**2 for x in range(a, b + 1))
+
+
+def sigma2caps_naive(n: int) -> int:
+    return sum(sigma2(i) for i in range(1, n + 1))
+
+
+def sum2(a: int, b: int) -> int:
+    # Calculate sum(x**2 for x in range(a, b + 1)) directly
+    # https://stackoverflow.com/questions/70902823/
+    r = b * (b + 1) * (2 * b + 1)
+    r -= (a - 1) * a * (2 * (a - 1) + 1)
+    r //= 6
+    return r
+
+
+def sigma2caps(n: int) -> int:
+    res = 0
+    mod = 10**9
+    a = 1
+    while a <= n:
+        d = n // a
+        b = n // d
+        s2 = sum2(a, b) % mod
+        res += (s2 * d) % mod
+        res %= mod
+        a = b + 1
+    return res
+
+
+def euler_401():
+    assert sum2(23, 1337) == sum(x**2 for x in range(23, 1337 + 1))
+    assert sigma2caps_naive(6) == 113
+    assert sigma2caps(6) == 113
+    assert sigma2caps(10**6) == 385965077
+    assert sigma2caps(10**9) == 73475174
+    return sigma2caps(10**15)
+
+
+if __name__ == "__main__":
+    solution = euler_401()
+    print(f"e401.py: {solution}")
+    assert solution == 281632621

python/e587.py

@@ -0,0 +1,61 @@
+from math import sqrt, pi, atan
+
+
+def get_intersection(n: int) -> tuple[int, int]:
+    """Compute the point where the line intersects the circle for a given
+    number of circles `n`.
+
+    Equation for circle: y = 1 - sqrt(2*x - x**2)
+    Equation for line: y = (1/n) * x
+    """
+
+    # Computed via Wolfram alpha; could be derived by setting the curves equal,
+    # isolating the root, and then applying the quadratic formula.
+    x_int = (-sqrt(2) * n ** (3 / 2) + n**2 + n) / (n**2 + 1)
+    y_int = x_int / n
+    return (x_int, y_int)
+
+
+def get_ratio_orange(n: int) -> float:
+    """Compute the ratio between the orange and the blue section. The blue
+    section is the L-section meaning the area outside of the quarter circle.
+    The orange area is the concave triangle created by the diagonal line."""
+
+    x, y = get_intersection(n)
+    area_blue = 1 - (pi / 4)
+
+    theta = atan((1 - x) / (1 - y))  # angle of concave triangle circle part
+    area_sector = theta / 2  # area of sector of concave triangle circle part
+
+    # area_orange = area triangle + area_rectangle + area_triangle - area_sector
+    area_orange = (0.5 * x * y) + ((1 - x) * y) + (0.5 * (1 - x) * (1 - y))
+    area_orange -= area_sector
+    return area_orange / area_blue
+
+
+def get_less_than(percentage: float) -> int:
+    ratio = percentage / 100
+    n_min, n_max = 1, 10000
+
+    while n_max - n_min > 1:
+        n_mid = (n_max + n_min) // 2
+        ratio_n = get_ratio_orange(n_mid)
+        if ratio_n < ratio:
+            n_max = n_mid
+        else:
+            n_min = n_mid
+
+    return n_max
+
+
+def euler_587():
+    assert get_ratio_orange(1) == 0.5
+    assert abs(get_ratio_orange(2) - 0.3646) <= 0.0001
+    assert get_less_than(10) == 15
+    return get_less_than(0.1)
+
+
+if __name__ == "__main__":
+    solution = euler_587()
+    print(f"e587.py: {solution}")
+    assert solution == 2240

python/e710.py

@@ -0,0 +1,95 @@
+MOD = 10**7
+SPACE = 10**7
+CACHE = [None for _ in range(SPACE * 2)]
+CACHE_TILL = [None for _ in range(SPACE)]
+CACHE_TILL_IS_TWO_EVEN = [None for _ in range(SPACE)]
+CACHE_TILL_IS_TWO_ODD = [None for _ in range(SPACE)]
+CACHE_T_FROM_ACC = [None for _ in range(SPACE * 2)]
+
+
+def t_from_acc(rem: int, is_two: bool) -> int:
+    offset = rem + SPACE * is_two
+    if (value := CACHE_T_FROM_ACC[offset]) is not None:
+        return value
+    if rem < 0:
+        return 0
+    s = t_from(rem, is_two) + t_from_acc(rem - 2, is_two)
+    s %= MOD
+    CACHE_T_FROM_ACC[offset] = s
+    return s
+
+
+def t_from(rem: int, is_two: bool) -> int:
+    offset = rem + SPACE * is_two
+    if (value := CACHE[offset]) is not None:
+        return value
+    s = 0
+    if rem < 0:
+        s = 0
+    elif rem == 1:
+        s = 0
+    elif rem == 0 and is_two:
+        s = 1
+    else:
+        s += t_from(rem - 2, is_two)  # d = 1
+        s %= MOD
+        s += t_from(rem - 4, True)  # d = 2
+        s %= MOD
+        s += t_from_acc(rem - 6, is_two)
+        s %= MOD
+        # Naive (too slow):
+        # for d in range(1, rem // 2 + 1):
+        #     s += t_from(rem - 2 * d, is_two or d == 2)
+        #     s %= MOD
+    CACHE[offset] = s % MOD
+    return s % MOD
+
+
+def t_old(n: int) -> int:
+    # print(f"t({n}) = ", end="")
+    s = 0
+    for i in range(0, n + 1):
+        rem = n - i
+        s += t_from(rem, i == 2)
+    return s % MOD
+
+
+def t(n: int) -> int:
+    s = CACHE_TILL[n - 4]
+    assert s is not None
+    s += t_from(n - 3, False)
+    CACHE_TILL[n - 3] = s % MOD
+    s += t_from(n - 2, True)
+    s %= MOD
+    s += t_from(n - 1, False)
+    s %= MOD
+    s += t_from(n, False)
+    s %= MOD
+    return s
+
+
+def euler_710():
+    for i in range(4):
+        CACHE_TILL[i] = 0
+
+    for n in range(4, 50):
+        t(n)
+    assert t(5) == 1
+    assert t(6) == 4
+    assert t(20) == 824
+    assert t(42) == 1999923
+
+    for n in range(10, SPACE):
+        tn = t(n)
+        if tn % 10**6 == 0:
+            # print(f"t({n}) = {tn} <<<")
+            return n
+        # if tn % 10000 == 0:
+        #     print(f"t({n}) = {tn}")
+    assert False, "Solution not found"
+
+
+if __name__ == "__main__":
+    solution = euler_710()
+    print(f"e710.py: {solution}")
+    assert solution == 1275000

python/e800.py

@@ -0,0 +1,68 @@
+from lib_prime import primes
+from math import log, ceil
+
+
+def c(n_base: int) -> int:
+    """Compute the number of prime pairs (p, q) where
+    p**q * q**p <= n_base**n_base and p < q"""
+    digits = ceil(log(n_base, 10)) + 1
+    n = n_base**n_base
+    n_log = log(n)
+    ps = primes(2 * 10**digits)
+    p, q = ps[0], ps[-1]
+    assert (p**q * q**p) > n, "largest prime too low"
+    j_max = len(ps) - 1
+
+    t = 0
+    for i in range(len(ps)):
+        p, q = ps[i], ps[i + 1]
+
+        # Crux:
+        #             (p**q * q**p) > n
+        # <=>       ln(p**q * q**p) > ln(n)
+        # <=>   ln(p**q) + ln(q**p) > ln(n)
+        # <=> q * ln(p) + p * ln(q) > ln(n)
+        # Much lower numbers and therefore faster to compute.
+        if q * log(p) + p * log(q) > n_log:
+            break
+
+        # Binary search to find the largest j for which p**q * q**p <= n:
+        j_lo, j_hi = i + 1, j_max
+        while (j_hi - j_lo) > 1:
+            j_mid = (j_hi + j_lo) // 2
+            p, q = ps[i], ps[j_mid]
+            if q * log(p) + p * log(q) > n_log:
+                j_hi = j_mid
+            else:
+                j_lo = j_mid
+        j_max = j_lo + 1
+
+        # All primes lower than ps[j_lo] have to be counted.
+        t += j_lo - i
+
+        # p_ref, q_ref = ps[i], ps[j_lo]
+        # print(p_ref, q_ref, "<=")
+
+        # Turns out the binary search is not necessary with the optimized
+        # version of the comparison and the following works just as fine.
+        # for j in range(j_max, i, -1):
+        #     p, q = ps[i], ps[j]
+        #     if q * log(p) + p * log(q) > n_log:
+        #         j_max = j
+        #     else:
+        #         assert p == p_ref
+        #         assert q == q_ref
+        #         t += (j - i)
+        #         break
+    return t
+
+
+def euler_800():
+    assert c(800) == 10790
+    return c(800800)
+
+
+if __name__ == "__main__":
+    solution = euler_800()
+    print(f"e800.py: {solution}")
+    assert(solution == 1412403576)

python/lib_create_templates.py

@@ -8,7 +8,7 @@ def euler_XXX():
 if __name__ == "__main__":
     solution = euler_XXX()
     print(f"eXXX.py: {solution}")
-    # assert(solution == 0)
+    # assert solution == 0
 
 """