The Hexadecimal Digit Canon Challenge: Solution

from __future__ import annotations

# We only need very lightweight helpers here; this is intentionally a
# *naive* reference implementation that favors clarity over performance.


def _format_base(d: int, base: int) -> str:
    """Return the base-`base` representation of `d` as a string.

    We restrict ourselves to bases 8, 10, and 16:
      - base 8  -> octal digits 0–7
      - base 10 -> decimal digits 0–9
      - base 16 -> hexadecimal digits 0–9, a–f

    The returned string uses lowercase letters for hexadecimal digits.
    """
    if base == 10:
        # Decimal needs no special formatting
        return str(d)
    if base == 8:
        # Python's built-in octal formatter
        return format(d, "o")
    if base == 16:
        # Python's built-in hexadecimal formatter (lowercase)
        return format(d, "x")

    # Any other base is deliberately rejected to keep the code simple
    raise ValueError(f"Unsupported base={base}. Use one of: 8, 10, 16.")


def f_of_d(d: int, base: int) -> int:
    """Compute f(d) in the given base.

    This follows the problem definition *literally*:

      1) Write d in the chosen base.
      2) Sort the digits in ascending order.
      3) Remove all zero digits.
      4) Interpret the remaining digits again as a number in that base.

    Example (base 16):
      d = 0x3A04  -> "3a04"
      digits     -> ["3", "a", "0", "4"]
      sorted     -> ["0", "3", "4", "a"]
      no zeros   -> ["3", "4", "a"]
      result     -> int("34a", 16)
    """
    # Convert d to a string in the given base
    s = _format_base(d, base)

    # Remove all '0' digits, because f(d) discards them by definition.
    # The result is a list of digit *characters*.
    digits_no_zero = [ch for ch in s if ch != "0"]

    # Sort remaining digits in ascending order.
    #
    # Important subtlety:
    # - For base 8 and 10, lexicographic order equals numeric order.
    # - For base 16, Python uses characters '0'–'9','a'–'f', and their
    #   lexicographic order also matches numeric digit order.
    digits_no_zero.sort()

    # There is always at least one non-zero digit:
    # any representation of a positive integer contains at least one
    # non-zero digit in any base.
    canonical_str = "".join(digits_no_zero)

    # Convert back to an integer, interpreting the string in the same base
    return int(canonical_str, base)


def S_naive(
    n: int,
    base: int,
    debug_target: int | None = None,
) -> tuple[int, list[int], int, int]:
    """Naively compute S(n) by brute force.

    We enumerate *all* positive integers d with at most n digits in the
    given base, i.e.

        1 <= d < base**n

    For each such d, we compute f(d) and accumulate the sum.

    This function is intentionally slow and is meant only as:
      - a correctness reference
      - a debugging / sanity-check tool

    Returns:
      total:
          The value of S(n) as a Python integer.
      debug_hits:
          All d such that f(d) == debug_target (empty if debug_target is None).
      num_values:
          Total number of integers d that were enumerated.
      num_distinct:
          Number of distinct values taken by f(d).
    """
    if n <= 0:
        raise ValueError("n must be a positive integer.")
    if base not in (8, 10, 16):
        raise ValueError("base must be one of: 8, 10, 16.")

    # All positive integers with <= n digits in base `base` satisfy
    #   1 <= d < base**n
    upper_exclusive = base**n

    debug_hits: list[int] = []

    # Track which f(d) values we have seen so far
    seen: set[int] = set()

    total = 0           # running sum S(n)
    num_distinct = 0    # number of distinct f(d) values

    for d in range(1, upper_exclusive):
        # Compute f(d) directly from the definition
        f_d = f_of_d(d, base)

        # Add contribution to the total sum
        total += f_d

        # Count distinct f(d) values (purely diagnostic)
        if f_d not in seen:
            seen.add(f_d)
            num_distinct += 1

        # Optional debugging hook:
        # record which original d map to a particular f(d)
        if debug_target is not None and f_d == debug_target:
            debug_hits.append(d)

    # Number of values enumerated (purely diagnostic)
    num_values = upper_exclusive - 1

    return total, debug_hits, num_values, num_distinct


def format_answer(x: int, base: int) -> str:
    """Format the final answer in the requested output base.

    This is useful because the *internal* computation uses Python integers,
    but the puzzle statement asks for answers written in base 8/10/16.
    """
    if base == 10:
        return str(x)
    if base == 8:
        return format(x, "o")
    if base == 16:
        # Uppercase hex looks nicer and is common in math puzzles
        return format(x, "X")
    raise ValueError(f"Unsupported base={base}.")


# ---------------------------------------------------------------------
# Example usage
# ---------------------------------------------------------------------

base = 16          # choose 8, 10, or 16
n = 0x3              # compute S(n)

# Set to an integer (in base-10) to see which d map to that f(d),
# or set to None to disable debugging.
debug_target = 0xDF

ans, debug_hits, num_values, num_distinct = S_naive(
    n=n,
    base=base,
    debug_target=debug_target,
)

print("num_values   =", num_values)
print("num_distinct =", num_distinct)
print("debug_hits   =", len(debug_hits))
print("S(n)         =", format_answer(ans, base))

def distinct_counts_by_n(max_n: int, base: int) -> list[tuple[int, int]]:
    """Return [(n, num_distinct_for_S(n))] for n = 1..max_n.

    num_distinct_for_S(n) counts the number of distinct f(d) values
    among all d with at most n digits in the given base, i.e.

        1 <= d < base**n
    """
    results: list[tuple[int, int]] = []

    for n in range(1, max_n + 1):
        upper_exclusive = base**n
        distinct: set[int] = set()

        for d in range(1, upper_exclusive):
            distinct.add(f_of_d(d, base))

        results.append((n, len(distinct)))

    return results


# --------------------------------------------------
# Example: observe growth for small n
# --------------------------------------------------

base = 16  # choose 8, 10, or 16

for n, cnt in distinct_counts_by_n(max_n=5, base=base):
    print(f"n={n}  num_distinct={cnt}")

from __future__ import annotations

from math import comb


def num_distinct_closed_form(n: int, base: int) -> int:
    """Number of distinct f(d) values among all positive integers with <= n digits in `base`.

    Key idea:
      For a fixed digit-length k, the value f(d) is determined solely by how many times
      each *non-zero* digit appears (zeros are discarded).

      Let b = base. The non-zero digits are {1, 2, ..., b-1}.
      For length k we count solutions to

          c_1 + c_2 + ... + c_{b-1} = k,  with c_i >= 0,

      i.e. the number of multisets of size k drawn from (b-1) symbols:

          C(k + (b-1) - 1, (b-1) - 1) = C(k + b - 2, b - 2).

      Summing over k = 1..n gives:

          sum_{k=1..n} C(k + b - 2, b - 2)
          = C(n + b - 1, b - 1) - 1

      (The "-1" removes the k=0 case corresponding to the empty multiset.)
    """
    if n <= 0:
        return 0
    if base not in (8, 10, 16):
        raise ValueError("base must be one of 8, 10, 16.")

    # Closed form: C(n + base - 1, base - 1) - 1
    return comb(n + base - 1, base - 1) - 1


def print_distinct_table(n_max: int, base: int) -> None:
    """Print a small table of num_distinct for n=1..n_max."""
    print(f"Base = {base}")
    for n in range(1, n_max + 1):
        print(f"n={n:2d}  num_distinct={num_distinct_closed_form(n, base)}")


# --------------------------------------------------
# Example usage
# --------------------------------------------------

for base in (16, ): # could also add other bases like 8 or 10
    print_distinct_table(n_max=20, base=base)
    print()

from __future__ import annotations

import math
from collections import Counter, deque
from collections.abc import Iterable


def S(n: int, base: int, fmt: str) -> int:
    """Compute S(n) by enumerating distinct f(d) values and counting multiplicities.

    Structure (kept the same):
      1) Enumerate all distinct values of f(d) among numbers with <= n digits.
         Each such value is a nondecreasing sequence of non-zero digits (1..base-1).
      2) For each distinct output f_d, count how many original numbers d (with <= n digits)
         map to it, by:
           - counting multiset permutations of the k non-zero digits
           - inserting z zeros for z = 0..(n-k), while forbidding a leading zero

    Args:
        n: Maximum number of digits of d (in base `base`).
        base: Base b (e.g. 8, 10, 16).
        fmt: Python `format` specifier matching the base:
            - base=10 -> "d"
            - base=16 -> "x" (or "X")
            - base=8  -> "o"

    Returns:
        S(n) as an exact integer.
    """
    # -------------------------------------------------------------------------
    # Phase 1: enumerate all distinct f(d) values for numbers with <= n digits.
    #
    # A value f(d) is just the sorted non-zero digits of d. Therefore it is
    # fully determined by a nondecreasing digit string over {1,2,...,base-1}.
    # We generate all such strings of length 1..n with a FIFO queue (BFS).
    # -------------------------------------------------------------------------
    max_value_exclusive = base**n  # any integer with >= n digits in base `base` is >= base**n

    fifo: deque[int] = deque(range(1, base))  # seed with all 1-digit non-zero values
    all_f_values: list[int] = []

    while fifo:
        f_d = fifo.popleft()
        all_f_values.append(f_d)

        # Try to extend by one digit (append in base `base`).
        new_prefix = base * f_d
        if new_prefix >= max_value_exclusive:
            continue  # already length n, cannot extend further

        last_digit = f_d % base
        for next_digit in range(last_digit, base):
            fifo.append(new_prefix + next_digit)

    # Optional sanity check (if available).
    assert len(all_f_values) == num_distinct_closed_form(n, base)

    # -------------------------------------------------------------------------
    # Helper: multiset permutations count = k! / prod_i c_i!
    # where c_i are multiplicities of equal digits in the k-digit string f_d.
    # -------------------------------------------------------------------------
    def multinomial_coefficient(k: int, counts: Iterable[int]) -> int:
        denom = 1
        for c in counts:
            denom *= math.factorial(c)
        return math.factorial(k) // denom

    # -------------------------------------------------------------------------
    # Phase 2: sum contributions.
    #
    # For a fixed distinct f_d:
    #   - let k be its length (number of non-zero digits)
    #   - perms_nonzero counts how many different permutations of those k digits exist
    #
    # Now fix z zeros to insert (0 <= z <= n-k). The final length is L = k+z.
    # Constraint: d must be an L-digit number, so its first digit cannot be 0.
    #
    # Therefore we choose the z zero-positions among the remaining L-1 slots:
    #
    #   interleave_zeros(k, z) = C(L-1, z) = C(k+z-1, z).
    # -------------------------------------------------------------------------
    total = 0

    for f_d in all_f_values:
        f_d_str = format(f_d, fmt)  # base-appropriate digit string
        k = len(f_d_str)

        # Multiplicities of digits within f_d (e.g. "11227" -> counts [2,1,1,1]).
        counts = Counter(f_d_str).values()

        # Number of distinct permutations of these k digits.
        perms_nonzero = multinomial_coefficient(k, counts)

        multiplicity = 0
        for z in range(0, n - k + 1):
            # Ways to insert z zeros into total length L=k+z without leading zero:
            # choose z positions among the last (L-1) positions.
            interleave_zeros = math.comb(k + z - 1, z) if z > 0 else 1

            multiplicity += perms_nonzero * interleave_zeros

        total += f_d * multiplicity

    return total


# Example usage: compute S(n) for hexadecimal numbers
for n_test in (3, 5, 0xA):
    value = S(n_test, base=16, fmt="x")
    print(f"S(0x{format_answer(n_test, base)}) in base 16 = 0x{format_answer(value, base)}")

from __future__ import annotations

import math
from fractions import Fraction


def S(n: int, base: int) -> int:
    """Compute S(n) exactly using the A(k) / B(k) decomposition (rational DP).

    We use the identity

        S(n) = sum_{k=1..n}  k * A(k) * B(k),

    where

        A(k) = sum_{m=k..n} (m-1)! / (m-k)!,

    and B(k) is the weighted sum over all digit-multiplicity patterns
    (c_1, ..., c_{base-1}) with sum c_i = k:

        B(k) = sum_{c_1+...+c_{base-1}=k}  f(c_1,...,c_{base-1}) / prod_i c_i!.

    Instead of enumerating all (c_1,...,c_{base-1}), we compute B(k) via a DP
    over digits d = 1..base-1, maintaining two arrays:

      - W[t] = sum_{patterns of length t}  1 / prod c_i!
      - B[t] = sum_{patterns of length t}  f(pattern) / prod c_i!

    Both W and B are Fractions to keep the derivation exact.
    """
    # -------------------------------------------------------------------------
    # 1) Compute A[k] for k=1..n:
    #
    #    A(k) = sum_{m=k..n} (m-1)! / (m-k)!
    #
    # This quantity depends only on lengths, not on digit composition.
    # We store A[0]=0 unused for convenience.
    # -------------------------------------------------------------------------
    A: list[Fraction] = [Fraction(0, 1)] * (n + 1)

    for k in range(1, n + 1):
        # Sum over total lengths m (total digits, including zeros).
        # Note: (m-1)!/(m-k)! is an integer, but we keep Fraction for uniformity.
        for m in range(k, n + 1):
            A[k] += Fraction(math.factorial(m - 1), math.factorial(m - k))

    # -------------------------------------------------------------------------
    # 2) Compute B[k] for k=0..n via DP over digits d=1..base-1.
    #
    # DP state after processing digits 1..(d-1):
    #   W[t] = sum_{c_1+...+c_{d-1}=t}  1 / (c_1! ... c_{d-1}!)
    #   B[t] = sum_{c_1+...+c_{d-1}=t}  f(c_1,...,c_{d-1}) / (c_1! ... c_{d-1}!)
    #
    # Initialization (no digits chosen yet):
    #   - empty pattern of length 0 has weight 1
    #   - its value-sum is 0 (there is no number yet)
    # -------------------------------------------------------------------------
    W: list[Fraction] = [Fraction(0, 1)] * (n + 1)
    B: list[Fraction] = [Fraction(0, 1)] * (n + 1)
    W[0] = Fraction(1, 1)

    # Process each possible non-zero digit exactly once.
    for d in range(1, base):
        newW: list[Fraction] = [Fraction(0, 1)] * (n + 1)
        newB: list[Fraction] = [Fraction(0, 1)] * (n + 1)

        # We may append digit d exactly c times, where c >= 0.
        # If current non-zero length is t, then t+c must be <= n.
        for t in range(0, n + 1):
            if W[t] == 0 and B[t] == 0:
                # Small skip: no patterns of this length exist in the current DP state.
                continue

            for c in range(0, n - t + 1):
                # Each digit multiplicity contributes the factor 1/c!
                coef = Fraction(1, math.factorial(c))

                # Appending c digits to the right in base `base` multiplies by base^c.
                p_base = base**c

                # The suffix consisting of c copies of digit d has value:
                #   d * (1 + base + base^2 + ... + base^(c-1))
                # The geometric sum factor is:
                #   rep = (base^c - 1) / (base - 1)
                # Example (base 16): c=2 => rep = 1 + 16 = 0x11
                rep = (p_base - 1) // (base - 1)

                # Weight update:
                #   W_new[t+c] += W[t] * (1/c!)
                newW[t + c] += W[t] * coef

                # Value-sum update:
                #
                # For each old pattern/value x of length t:
                #   new_value = x * base^c + d * rep
                #
                # Aggregating over all patterns of length t:
                #   B[t] aggregates sum(x / weight_den)
                #   W[t] aggregates sum(1 / weight_den)
                #
                # Thus:
                #   B_new[t+c] += (B[t] * base^c + W[t] * d * rep) * (1/c!)
                newB[t + c] += (B[t] * p_base + W[t] * d * rep) * coef

        W, B = newW, newB

    # Now B[k] equals the desired B(k) after processing digits 1..base-1.

    # -------------------------------------------------------------------------
    # 3) Assemble S(n) = sum_{k=1..n} k * A[k] * B[k]
    # -------------------------------------------------------------------------
    total = Fraction(0, 1)
    for k in range(1, n + 1):
        total += k * A[k] * B[k]

    # The math guarantees the result is an integer.
    return int(total)


# Example usage: compute S(n) for hexadecimal numbers
# Example usage: compute S(n) for hexadecimal numbers
for n_test in (3, 5, 0xA, 0xAA):
    value = S(n_test, base=16)
    print(f"S(0x{format_answer(n_test, base)}) in base 16 = 0x{format_answer(value, base)}")

from __future__ import annotations

import math


def S(n: int, base: int) -> int:
    """Compute S(n) exactly using the integer-scaled DP (no Fractions).

    We use the decomposition
        S(n) = sum_{k=1..n} k * A(k) * B(k),

    but compute B(k) via scaled integers:
        B_scaled[k] = k! * B(k).

    DP meaning (after processing digits 1..d):
      - W_scaled[t] = t! * sum_{patterns length t} 1 / prod c_i!
      - B_scaled[t] = t! * sum_{patterns length t} f(pattern) / prod c_i!

    The digit-extension update (append digit d exactly c times) becomes integer-only:
        coef = C(t+c, c)
        W_scaled_new[t+c] += W_scaled[t] * coef
        B_scaled_new[t+c] += (B_scaled[t] * base^c + W_scaled[t] * d * rep(c)) * coef

    where rep(c) = (base^c - 1)/(base - 1) = 1 + base + ... + base^(c-1).

    Finally, since B_scaled[k] = k! * B(k), we have:
        S(n) = sum_{k=1..n} k * A(k) * B_scaled[k] / k!

    The division is exact by construction.
    """
    # -------------------------------------------------------------------------
    # 1) Precompute A[k] = sum_{m=k..n} (m-1)!/(m-k)!   (pure integers)
    # -------------------------------------------------------------------------
    A = [0] * (n + 1)
    for k in range(1, n + 1):
        # (m-1)!/(m-k)! = (m-k)(m-k+1)...(m-1) is an integer rising product of length (k-1).
        for m in range(k, n + 1):
            A[k] += math.factorial(m - 1) // math.factorial(m - k)

    # -------------------------------------------------------------------------
    # 2) Compute scaled B via DP over digits d=1..base-1.
    #
    # W_scaled[0]=1 represents the empty pattern (length 0) with weight 1.
    # B_scaled[0]=0 because the empty pattern has value 0.
    # -------------------------------------------------------------------------
    W_scaled = [0] * (n + 1)
    B_scaled = [0] * (n + 1)
    W_scaled[0] = 1

    for d in range(1, base):  # digits 1..9 (base10), 1..7 (base8), 1..15 (base16), ...
        newW = [0] * (n + 1)
        newB = [0] * (n + 1)

        # For each existing length t, we may add c copies of digit d, keeping total <= n.
        for t in range(0, n + 1):
            if W_scaled[t] == 0 and B_scaled[t] == 0:
                continue

            for c in range(0, n - t + 1):
                # coef = (t+c)!/(t! c!) = C(t+c, c)
                # This is exactly the scaling factor that removes the 1/c! fractions
                # after we multiply everything by (t+c)!.
                coef = math.comb(t + c, c)

                # Appending c digits shifts old values by base^c.
                p_base = base**c

                # Value of suffix "ddd...d" (c times) in base `base`:
                # d * (1 + base + ... + base^(c-1)) = d * ((base^c - 1)/(base - 1))
                rep = (p_base - 1) // (base - 1)

                # Update weights (scaled): only depends on how many patterns we had at length t.
                newW[t + c] += W_scaled[t] * coef

                # Update value-sums (scaled):
                # - old aggregated values shift by base^c
                # - plus the suffix contribution, scaled by the total weight mass W_scaled[t]
                newB[t + c] += (B_scaled[t] * p_base + W_scaled[t] * d * rep) * coef

        W_scaled, B_scaled = newW, newB

    # -------------------------------------------------------------------------
    # 3) Assemble S(n) using B_scaled[k] = k! * B(k)
    # -------------------------------------------------------------------------
    ans = 0
    for k in range(1, n + 1):
        ans += k * A[k] * B_scaled[k] // math.factorial(k)

    return ans


# Example usage: compute S(n) for hexadecimal numbers
for n_test in (3, 5, 0xA, 0xAA):
    value = S(n_test, base=16)
    print(f"S(0x{format_answer(n_test, base)}) in base 16 = 0x{format_answer(value, base)}")

mod = 0x1FFFFFFFFFFFFFFF
for n_test in (3, 5, 0xA, 0xAA):
    value = S(n_test, base=16) % mod
    print(f"S(0x{format_answer(n_test, base)}) in base 16 = 0x{format_answer(value, base)}")

from __future__ import annotations

import math
from sympy import totient


def S(n: int, base: int, mod: int) -> int:
    """Compute S(n) modulo `mod` using the modular A(k)/B(k) decomposition.

    This is the modular analogue of the integer-scaled DP:

        S(n) = sum_{k=1..n} k * A(k) * B(k)  (mod mod)

    where
        A(k) = sum_{m=k..n} (m-1)! / (m-k)!,
    and B(k) is obtained via a digit-DP over digits 1..(base-1).

    Modular inverses are computed using Euler's theorem:
        a^{-1} ≡ a^{φ(mod)-1} (mod mod),
    which requires gcd(a, mod) = 1 for all inverted values.
    """

    # -------------------------------------------------------------------------
    # Helpers
    # -------------------------------------------------------------------------
    phi_mod = int(totient(mod))

    def inv_mod(a: int) -> int:
        """Modular inverse using Euler's theorem (assumes gcd(a, mod) = 1)."""
        return pow(a % mod, phi_mod - 1, mod)

    def factorial_mod(m: int) -> int:
        """Compute m! modulo mod (naive O(m))."""
        acc = 1
        for i in range(1, m + 1):
            acc = (acc * i) % mod
        return acc

    # -------------------------------------------------------------------------
    # 1) Compute A(k) = sum_{m=k..n} (m-1)! / (m-k)!   (mod mod)
    # -------------------------------------------------------------------------
    A = [0] * (n + 1)

    for k in range(1, n + 1):
        for m in range(k, n + 1):
            numerator = factorial_mod(m - 1)
            denominator = factorial_mod(m - k)
            A[k] = (A[k] + numerator * inv_mod(denominator)) % mod

    # -------------------------------------------------------------------------
    # 2) Compute B(k) via digit DP (fraction-style, but done modulo mod)
    #
    # W[t] = sum of weights for patterns of length t
    # B[t] = weighted sum of values for patterns of length t
    # -------------------------------------------------------------------------
    W = [0] * (n + 1)
    B = [0] * (n + 1)

    W[0] = 1  # empty pattern has weight 1
    inv_base_minus_1 = inv_mod(base - 1)

    for d in range(1, base):  # digits 1..(base-1)
        newW = [0] * (n + 1)
        newB = [0] * (n + 1)

        for t in range(0, n + 1):
            if W[t] == 0 and B[t] == 0:
                continue

            for c in range(0, n - t + 1):
                # coef = 1 / c!  (mod mod)
                coef = inv_mod(factorial_mod(c))

                # base^c mod mod (shift old values)
                p_base = pow(base, c, mod)

                # rep(c) = (base^c - 1) / (base - 1)
                rep = (p_base - 1) * inv_base_minus_1 % mod

                # Update weights
                newW[t + c] = (newW[t + c] + W[t] * coef) % mod

                # Update value sums
                shifted_old = (B[t] * p_base) % mod
                suffix_part = (W[t] * d * rep) % mod
                newB[t + c] = (newB[t + c] + (shifted_old + suffix_part) * coef) % mod

        W, B = newW, newB

    # -------------------------------------------------------------------------
    # 3) Assemble final answer
    # -------------------------------------------------------------------------
    ans = 0
    for k in range(1, n + 1):
        ans = (ans + (k * A[k] % mod) * B[k]) % mod

    return ans


# Example usage: compute S(n) for hexadecimal numbers
base = 16
mod = 0x1FFFFFFFFFFFFFFF
for n_test in (3, 5, 0xA, 0xAA):
    value = S(n_test, base=base, mod=mod)
    print(f"S(0x{format_answer(n_test, base)}) in base 16 = 0x{format_answer(value, base)}")

from __future__ import annotations

from sympy import totient


def S(n: int, base: int, mod: int) -> int:
    """Compute S(n) modulo `mod` using the modular A(k)/B(k) decomposition.

    Same logic as above version, but faster by caching:
      - factorials (0!..n!) modulo mod
      - inverse factorials via pow(fact[n], phi-1) and downward recurrence
      - powers base^c modulo mod for c=0..n
      - rep(c) = (base^c - 1)/(base-1) modulo mod for c=0..n

    NOTE: This assumes all inverses exist (e.g. mod is prime and mod > n, or at
    least gcd(denominator, mod)=1 for all denominators we invert).
    """
    phi_mod = int(totient(mod))

    # ---------------------------------------------------------------------
    # Precompute factorials and inverse factorials up to n
    # ---------------------------------------------------------------------
    fact = [1] * (n + 1)      # fact[i] = i! mod mod
    for i in range(1, n + 1):
        fact[i] = (fact[i - 1] * i) % mod

    # inv_fact[i] = (i!)^{-1} mod mod, computed in O(n) using one pow() + recurrence
    inv_fact = [1] * (n + 1)
    inv_fact[n] = pow(fact[n], phi_mod - 1, mod)  # (n!)^{-1}
    for i in range(n, 0, -1):
        # (i-1)!^{-1} = i!^{-1} * i
        inv_fact[i - 1] = (inv_fact[i] * i) % mod

    # Also cache (base-1)^{-1} mod mod, used in rep(c)
    inv_base_minus_1 = pow(base - 1, phi_mod - 1, mod)

    # ---------------------------------------------------------------------
    # Precompute powers and rep(c) for c=0..n
    # ---------------------------------------------------------------------
    pow_base = [1] * (n + 1)  # pow_base[c] = base^c mod mod
    for c in range(1, n + 1):
        pow_base[c] = (pow_base[c - 1] * base) % mod

    # rep[c] = (base^c - 1)/(base-1) mod mod
    rep = [0] * (n + 1)
    for c in range(0, n + 1):
        rep[c] = (pow_base[c] - 1) * inv_base_minus_1 % mod

    # coef[c] = 1/c! mod mod  (used very frequently in the DP)
    inv_c_fact = inv_fact  # alias to make intent clearer

    # ---------------------------------------------------------------------
    # 1) Compute A(k) = sum_{m=k..n} (m-1)! / (m-k)!   (mod mod)
    #    => (m-1)! * ((m-k)!)^{-1}
    # ---------------------------------------------------------------------
    A = [0] * (n + 1)
    for k in range(1, n + 1):
        acc = 0
        for m in range(k, n + 1):
            acc = (acc + fact[m - 1] * inv_fact[m - k]) % mod
        A[k] = acc

    # ---------------------------------------------------------------------
    # 2) Compute B(k) via digit DP (fraction-style, modulo arithmetic)
    # ---------------------------------------------------------------------
    W = [0] * (n + 1)
    B = [0] * (n + 1)
    W[0] = 1

    for d in range(1, base):
        newW = [0] * (n + 1)
        newB = [0] * (n + 1)

        for t in range(0, n + 1):
            wt = W[t]
            bt = B[t]
            if wt == 0 and bt == 0:
                continue

            # Append c copies of digit d, with 0 <= c <= n-t
            for c in range(0, n - t + 1):
                coef = inv_c_fact[c]         # 1/c!
                p = pow_base[c]              # base^c
                r = rep[c]                   # (base^c - 1)/(base-1)

                idx = t + c

                # newW[idx] += W[t] * coef
                newW[idx] = (newW[idx] + wt * coef) % mod

                # newB[idx] += (B[t] * base^c + W[t] * d * rep(c)) * coef
                shifted_old = (bt * p) % mod
                suffix_part = (wt * d) % mod
                suffix_part = (suffix_part * r) % mod
                newB[idx] = (newB[idx] + (shifted_old + suffix_part) * coef) % mod

        W, B = newW, newB

    # ---------------------------------------------------------------------
    # 3) Assemble final answer: sum_{k=1..n} k * A[k] * B[k]
    # ---------------------------------------------------------------------
    ans = 0
    for k in range(1, n + 1):
        ans = (ans + (k * A[k] % mod) * B[k]) % mod

    return ans


# Example usage (kept close to yours)
base = 16
mod  = 0x1FFFFFFFFFFFFFFF
for n_test in (3, 5, 0xA, 0xAA):
     value = S(n_test, base=base, mod=mod)
     print(f"S(0x{format_answer(n_test, base)}) in base 16 = 0x{format_answer(value, base)}")

value = S(0xAAA, base=base, mod=mod)
print(f"S(0xAAA) in base 16 = 0x{format_answer(value, base)}")