euler/910.py

#!/usr/bin/env python

#  runtime: 43s on pypy3 11
'''
From https://github.com/lucky-bai/projecteuler-solutions/issues/102 license unknown
Runtime: 34s on pypy3 3.11.11, Ubuntu 25.04, Lenovo ThinkPad P14s.
'''

from __future__ import annotations

"""
Standalone Project Euler 910 solver.

All supporting logic (5-adic helpers, Phi iterators, and the orbit solver)
is consolidated here so this file can run by itself:
    python pe_910_solver.py
"""

from dataclasses import dataclass
from functools import lru_cache
from math import comb
from typing import Callable, Dict, List, Tuple

# Problem-specific constants
PE_B = 345_678
PE_C = 9_012_345
PE_D = 678
MOD_DEFAULT = 10**9


# 5-adic unit helpers -----------------------------------------------------

def modulus_to_k(mod: int) -> int:
    """Return k such that mod = 5^k (assumes mod is a power of 5)."""
    k = 0
    temp = mod
    while temp % 5 == 0 and temp > 1:
        temp //= 5
        k += 1
    return k


def teichmuller_lift(x: int, mod: int) -> int:
    """Compute the Teichmüller lift of x modulo mod (mod = 5^k)."""
    z = x % mod
    k = modulus_to_k(mod)
    for _ in range(k):
        z = pow(z, 5, mod)
    return z


def unit_decompose(x: int, mod: int) -> Tuple[int, int]:
    """
    Decompose a 5-adic unit x modulo mod into (z, u) with
        x = z * (1 + 5u),  z^4 = 1,  u in Z/5^{k-1}Z.
    Returns (z, u).
    """
    z = teichmuller_lift(x, mod)
    inv_z = pow(z, -1, mod)
    y = (x * inv_z) % mod
    k = modulus_to_k(mod)
    if k == 0:
        return z, 0
    u = (y - 1) // 5
    return z, u


@lru_cache(maxsize=None)
def _binom_prefix(exp: int, k: int) -> Tuple[int, ...]:
    """Precompute (C(exp, r) mod 5^k) for r = 0..k."""
    mod = 5**k if k else 1
    coeffs: List[int] = [1]
    for r in range(1, k + 1):
        coeffs.append(comb(exp, r) % mod)
    return tuple(coeffs)


def one_plus_5u_pow(u: int, exp: int, mod: int) -> int:
    """
    Compute (1 + 5u)^exp modulo mod (mod = 5^k) via truncated binomial series.
    """
    k = modulus_to_k(mod)
    if k == 0:
        return 1 % mod
    coeffs = _binom_prefix(exp, k)
    res = 1
    base = (5 * u) % mod
    pow_term = base
    for r in range(1, k + 1):
        res = (res + (coeffs[r] * pow_term) % mod) % mod
        pow_term = (pow_term * base) % mod
    return res


def unit_pow(x: int, exp: int, mod: int) -> int:
    """Raise a 5-adic unit x to exponent exp modulo mod using the decomposition."""
    k = modulus_to_k(mod)
    if k == 0 or x % 5 == 0:
        return pow(x, exp, mod)
    z, u = unit_decompose(x, mod)
    z_pow = pow(z, exp % 4, mod)
    tup = one_plus_5u_pow(u, exp, mod)
    return (z_pow * tup) % mod


# Inner g-iterator for modulo 5^k ----------------------------------------

class GCIterator:
    """
    Lazy binary-lifting iterator for the inner map g(c, x) = x^c (x + 1) mod 5^k.
    Provides fast evaluation of Phi_0(x) = g_c^{∘ b}(g_{c+1}(x)).
    """

    def __init__(self, mod: int):
        self.mod = mod
        self.jump_cache: List[Dict[int, int]] = [{} for _ in range(PE_B.bit_length() + 1)]
        self.k = modulus_to_k(mod)
        self.unit_cache: Dict[int, Tuple[int, int]] = {}
        self.cycle_cache: Dict[int, Tuple[int, int]] = {}

    def _pow_unit(self, x: int, exp: int) -> int:
        m = self.mod
        if m == 1:
            return 0
        if x % 5 == 0:
            return pow(x, exp, m)
        cache = self.unit_cache
        if x not in cache:
            cache[x] = unit_decompose(x, m)
        z, u = cache[x]
        z_pow = pow(z, exp % 4 or 4, m)
        tup = one_plus_5u_pow(u, exp, m)
        return (z_pow * tup) % m

    def g_c(self, x: int) -> int:
        m = self.mod
        pow_part = self._pow_unit(x, PE_C)
        return (pow_part * ((x + 1) % m)) % m

    def g_cplus1(self, x: int) -> int:
        m = self.mod
        pow_part = self._pow_unit(x, PE_C + 1)
        return (pow_part * ((x + 1) % m)) % m

    def jump(self, level: int, x: int) -> int:
        """Compute g_c^{∘ 2^level}(x) modulo mod, lazily cached."""
        x %= self.mod
        cache = self.jump_cache[level]
        if x in cache:
            return cache[x]
        if level == 0:
            val = self.g_c(x)
        else:
            mid = self.jump(level - 1, x)
            val = self.jump(level - 1, mid)
        cache[x] = val
        return val

    def iterate_steps(self, steps: int, x: int) -> int:
        """Apply g_c exactly `steps` times starting from x using binary lifting."""
        res = x % self.mod
        bit = 0
        s = steps
        while s:
            if s & 1:
                res = self.jump(bit, res)
            s >>= 1
            bit += 1
        return res

    def Phi0(self, x: int) -> int:
        """Phi_0(x) = h(b+1, c, x) modulo mod."""
        y0 = self.g_cplus1(x % self.mod)
        return self.iterate_steps(PE_B + 1, y0)


# Cycle-accelerated generic helpers --------------------------------------

def g_mod(c: int, x: int, mod: int) -> int:
    """g(c, x) = x^c * (x + 1) modulo mod."""
    return (pow(x, c, mod) * ((x + 1) % mod)) % mod


def h_mod(b: int, c: int, x: int, mod: int) -> int:
    """
    h(b, c, x):
      t = g(c+1, x);
      then apply g(c, ·) exactly b times.
    """
    t = g_mod(c + 1, x, mod)
    for _ in range(b):
        t = g_mod(c, t, mod)
    return t


def iterate_with_cycle(
    f: Callable[[int], int],
    x0: int,
    n: int,
    mod: int,
    max_steps: int | None = None,
) -> int:
    """
    Compute f^n(x0) modulo mod using Floyd's cycle-finding algorithm.

    Intended for small mod (e.g., 2^9 or 5^9) and huge n.
    """
    if n == 0:
        return x0 % mod

    tortoise = f(x0) % mod
    hare = f(f(x0) % mod) % mod
    steps = 0
    if max_steps is None:
        max_steps = mod * 4  # generous upper bound

    while tortoise != hare and steps < max_steps:
        tortoise = f(tortoise) % mod
        hare = f(f(hare) % mod) % mod
        steps += 1

    if steps >= max_steps:
        # Fallback: simple iteration (used only in tiny test cases)
        x = x0 % mod
        for _ in range(n):
            x = f(x) % mod
        return x

    # Find mu (cycle start)
    mu = 0
    tortoise = x0 % mod
    while tortoise != hare and mu < max_steps:
        tortoise = f(tortoise) % mod
        hare = f(hare) % mod
        mu += 1

    # Find lambda (cycle length)
    lam = 1
    hare = f(tortoise) % mod
    while tortoise != hare and lam < max_steps:
        hare = f(hare) % mod
        lam += 1

    if n < mu:
        x = x0 % mod
        for _ in range(n):
            x = f(x) % mod
        return x

    # Move to position mu
    x = x0 % mod
    for _ in range(mu):
        x = f(x) % mod

    rem = (n - mu) % lam
    for _ in range(rem):
        x = f(x) % mod
    return x


def h_mod_cycle(b: int, c: int, x: int, mod: int) -> int:
    """Cycle-accelerated version of h_mod for large b and small mod."""
    x1 = g_mod(c + 1, x, mod)
    if b == 0:
        return x1
    f = lambda y: g_mod(c, y, mod)
    return iterate_with_cycle(f, x1, b, mod)


def Phi0_mod(E: int, b: int, c: int, mod: int = MOD_DEFAULT) -> int:
    """Phi_0(E) = h(b+1, c, E) modulo mod."""
    return h_mod(b + 1, c, E, mod)


def Phi0_mod_cycle(E: int, b: int, c: int, mod: int = MOD_DEFAULT) -> int:
    """
    Phi_0(E) using cycle detection on h_mod; suited for large b and small mod.
    """
    return h_mod_cycle(b + 1, c, E, mod)


def Phi_a_mod(a: int, E: int, b: int, c: int, mod: int = MOD_DEFAULT) -> int:
    """
    Compute Phi_a(E) via the outer update (naive for small parameters).
    """

    def phi0(x: int) -> int:
        return Phi0_mod(x, b, c, mod)

    if a == 0:
        return phi0(E)

    Phi_funcs: List[Callable[[int], int]] = [phi0]

    for _k in range(a):
        prev = Phi_funcs[-1]

        def make_next(prev_func: Callable[[int], int]) -> Callable[[int], int]:
            def next_func(x: int) -> int:
                s = (x * prev_func(x)) % mod
                for _ in range(b):
                    s = prev_func(s) % mod
                return s

            return next_func

        Phi_funcs.append(make_next(prev))

    return Phi_funcs[a](E)


def Phi_a_mod_cycle(a: int, E: int, b: int, c: int, mod: int = MOD_DEFAULT) -> int:
    """
    Cycle-accelerated Phi_a(E) for large b and small mod.
    """

    def phi0(x: int) -> int:
        return Phi0_mod_cycle(x, b, c, mod)

    if a == 0:
        return phi0(E)

    Phi_funcs: List[Callable[[int], int]] = [phi0]

    for _k in range(a):
        prev = Phi_funcs[-1]

        def make_next(prev_func: Callable[[int], int]) -> Callable[[int], int]:
            def next_func(x: int) -> int:
                s0 = (x * prev_func(x)) % mod

                def f(y: int) -> int:
                    return prev_func(y) % mod

                return iterate_with_cycle(f, s0, b, mod)

            return next_func

        Phi_funcs.append(make_next(prev))

    return Phi_funcs[a](E)


def T_mod(a: int, b: int, c: int, d: int, mod: int = MOD_DEFAULT) -> int:
    """T(a, b, c, d) modulo mod using the naive Phi_a recurrence."""
    return Phi_a_mod(a, d, b, c, mod)


def T_mod_cycle(a: int, b: int, c: int, d: int, mod: int = MOD_DEFAULT) -> int:
    """T(a, b, c, d) modulo mod using cycle-accelerated Phi_a."""
    return Phi_a_mod_cycle(a, d, b, c, mod)


# Orbit-focused Phi solver for 5^k ---------------------------------------

class JumpIter:
    """Binary-lifting iterator for an arbitrary function f: X -> X mod m."""

    def __init__(self, f: Callable[[int], int], mod: int, max_bits: int):
        self.f = f
        self.mod = mod
        self.jump: List[Dict[int, int]] = [dict() for _ in range(max_bits)]

    def _step(self, x: int) -> int:
        return self.f(x) % self.mod

    def _jump(self, bit: int, x: int) -> int:
        cache = self.jump[bit]
        xm = x % self.mod
        if xm in cache:
            return cache[xm]
        if bit == 0:
            val = self._step(xm)
        else:
            mid = self._jump(bit - 1, xm)
            val = self._jump(bit - 1, mid)
        cache[xm] = val
        return val

    def pow(self, x: int, n: int) -> int:
        res = x % self.mod
        bit = 0
        nn = n
        while nn:
            if nn & 1:
                res = self._jump(bit, res)
            nn >>= 1
            bit += 1
        return res


@dataclass
class OrbitPhiSolver:
    """
    Orbit-focused Phi_a solver for the PE 910 parameters modulo 5^k (k ≤ 9).
    Uses memoization and cycle detection along the orbit starting at d = 678.
    """

    mod: int
    cache: List[Dict[int, int]]
    g_iter: GCIterator

    def __init__(self, mod: int):
        self.mod = mod
        self.cache = []  # cache[a][x] = Phi_a(x)
        self.jump_iters: List[JumpIter | None] = []
        self.g_iter = GCIterator(mod)
        self.phi0_table = None
        if mod <= 2_000_000:  # up to 5^9 = 1,953,125
            self._build_phi0_table()

    def phi(self, a: int, x: int) -> int:
        """Phi_a(x) modulo self.mod (specialized to PE b,c)."""
        xm = x % self.mod
        self._ensure_level(a)
        cache_a = self.cache[a]
        if xm in cache_a:
            return cache_a[xm]

        if a == 0:
            if self.phi0_table is not None:
                val = self.phi0_table[xm]
            else:
                val = self.g_iter.Phi0(xm)
        else:
            prev = self.phi(a - 1, xm)
            s0 = (xm * prev) % self.mod
            it = self._ensure_jump_iter(a - 1)
            val = it.pow(s0, PE_B)

        cache_a[xm] = val
        return val

    def _build_phi0_table(self) -> None:
        """Precompute Phi0(x) for all residues modulo mod."""
        mod = self.mod
        table = [0] * mod
        for x in range(mod):
            table[x] = self.g_iter.Phi0(x)
        self.phi0_table = table

    def _ensure_level(self, a: int) -> None:
        while len(self.cache) <= a:
            self.cache.append({})
            self.jump_iters.append(None)

    def _ensure_jump_iter(self, level: int) -> JumpIter:
        self._ensure_level(level)
        it = self.jump_iters[level]
        if it is None:
            max_bits = PE_B.bit_length() + 1
            it = JumpIter(lambda y, lvl=level: self.phi(lvl, y), self.mod, max_bits)
            self.jump_iters[level] = it
        return it

    def U_sequence(self, a_max: int) -> List[int]:
        """Compute [Phi_a(PE_D) mod mod for a=0..a_max]."""
        return [self.phi(a, PE_D) for a in range(a_max + 1)]


# Top-level T(a,b,c,d) modulo specific powers -----------------------------

def T_mod_2pow9(a: int, b: int, c: int, d: int, e: int = 0) -> int:
    mod = 2**9
    m = mod
    Phi = [[0] * m for _ in range(a + 1)]

    for x in range(m):
        Phi[0][x] = Phi0_mod_cycle(x, b, c, mod)

    for k in range(0, a):
        prev = Phi[k]
        curr = Phi[k + 1]

        def f(y: int, table=prev) -> int:
            return table[y]

        for x in range(m):
            s0 = (x * prev[x]) % mod
            curr[x] = iterate_with_cycle(f, s0, b, mod)

    return (Phi[a][d % mod] + e) % mod


def T_mod_5pow9(a: int, b: int, c: int, d: int, e: int = 0) -> int:
    mod = 5**9
    if b == PE_B and c == PE_C:
        solver = OrbitPhiSolver(mod)
        return (solver.phi(a, d) + e) % mod
    return (T_mod_cycle(a, b, c, d, mod) + e) % mod


def T_mod_1e9(a: int, b: int, c: int, d: int, e: int = 0) -> int:
    mod2 = 2**9
    mod5 = 5**9
    t2 = T_mod_2pow9(a, b, c, d, e)
    t5 = T_mod_5pow9(a, b, c, d, e)
    inv_mod2_mod5 = pow(mod2, -1, mod5)
    k = ((t5 - t2) * inv_mod2_mod5) % mod5
    x = t2 + k * mod2
    return x % (mod2 * mod5)


def solve() -> int:
    a = 12
    b = PE_B
    c = PE_C
    d = PE_D
    e = 90
    return T_mod_1e9(a, b, c, d, e)


if __name__ == "__main__":
    print(solve())