Solution

Answer: 789184709

Let the decimal expansion of $k$ be

$$k=\sum_{i=0}^{n-1} a_i 10^i, \qquad 0\le a_i\le 9.$$

We want to count the numbers $k<10^n$ such that

$k\equiv 0 \pmod{23}$,
$\sum a_i = 23$.

Denote this count by $S(n)$.

1. Mathematical analysis

1.1 Encoding the conditions

The digit-sum condition is

$$\sum_{i=0}^{n-1} a_i = 23.$$

The divisibility condition is

$$\sum_{i=0}^{n-1} a_i 10^i \equiv 0 \pmod{23}.$$

So every digit contributes simultaneously:

$+a_i$ to the digit sum,
$+a_i 10^i$ modulo $23$.

This strongly suggests a generating-function approach.

1.2 A bivariate generating function

For a fixed position $i$, the digit $d\in{0,\dots,9}$ contributes

$$x^d y^{d10^i}.$$

Thus position $i$ contributes the polynomial

$$P_i(x,y)=\sum_{d=0}^9 x^d y^{d10^i}.$$

Multiplying over all positions:

$$F_n(x,y)=\prod_{i=0}^{n-1} P_i(x,y).$$

Then:

exponent of $x$ records digit sum,
exponent of $y$ records residue modulo $23$.

Therefore,

$$S(n)= [x^{23}y^0],F_n(x,y),$$

where the exponent of $y$ is interpreted modulo $23$.

1.3 Periodicity modulo 23

A crucial observation:

$$10^{22}\equiv 1 \pmod{23},$$

because $23$ is prime and $10$ is a primitive element modulo $23$.

Hence the sequence $10^i \pmod{23}$ has period $22$.

Therefore the factors repeat every $22$ positions.

Let

$$n = 22q + r,\qquad 0\le r<22.$$

Then each residue class position appears either:

$q$ times, or
$q+1$ times.

Specifically, for $j=0,\dots,21$,

$$P_j(x,y)=\sum_{d=0}^9 x^d y^{d10^j},$$

and

$$F_n(x,y) = \prod_{j=0}^{21} P_j(x,y)^{,q+\mathbf 1_{j<r}}.$$

1.4 Truncation miracle

We only care about coefficient $x^{23}$.

Thus during all multiplications we may discard every term of degree $>23$.

So the entire computation lives in a tiny finite algebra:

$x$-degree: $0\ldots 23$,
residue modulo $23$: $0\ldots 22$.

Total states:

$$24\times 23 = 552.$$

Exponentiation becomes extremely fast.

2. Python implementation

MOD = 10**9

# --------------------------------------------------------------------
# Polynomial representation:
# dictionary mapping (digit_sum, residue_mod_23) -> coefficient
# --------------------------------------------------------------------

def multiply(A, B):
    """
    Multiply two truncated polynomials.

    State:
        (s, r)
    means:
        x^s y^r

    We truncate all terms with s > 23.
    """
    C = {}

    for (s1, r1), v1 in A.items():
        for (s2, r2), v2 in B.items():

            s = s1 + s2

            # We only need digit sum <= 23
            if s <= 23:
                r = (r1 + r2) % 23

                key = (s, r)

                C[key] = (C.get(key, 0) + v1 * v2) % MOD

    return C

def power(P, e):
    """
    Fast exponentiation of a polynomial.
    """
    R = {(0, 0): 1}

    while e > 0:

        if e & 1:
            R = multiply(R, P)

        P = multiply(P, P)

        e >>= 1

    return R

def S(n):
    """
    Compute S(n) mod 1e9.
    """

    q, r = divmod(n, 22)

    # Start with multiplicative identity
    ans = {(0, 0): 1}

    # Build the 22 periodic factors
    for j in range(22):

        w = pow(10, j, 23)

        # P_j(x,y)
        P = {}

        for d in range(10):
            key = (d, (d * w) % 23)
            P[key] = P.get(key, 0) + 1

        exponent = q + (1 if j < r else 0)

        ans = multiply(ans, power(P, exponent))

    return ans.get((23, 0), 0)

# --------------------------------------------------------------------
# Verification on given examples
# --------------------------------------------------------------------

print(S(9))      # 263626
print(S(42))     # 878570056  (matches the last 9 digits)

# Final computation
N = 11**12

print(S(N))

3. Code walkthrough

`multiply(A, B)`

Each polynomial is stored sparsely.

A term

(s, r): value

means:

$$\text{value}\cdot x^s y^r.$$

When multiplying:

$$x^{s_1}y^{r_1}\cdot x^{s_2}y^{r_2} = x^{s_1+s_2}y^{r_1+r_2}.$$

The residue is reduced modulo $23$.

We discard all terms with digit sum $>23$.

`power(P, e)`

Standard binary exponentiation.

This computes:

$$P^e$$

in $O(\log e)$ multiplications.

Since $11^{12}\approx 3.1\times 10^{12}$, only about $42$ squaring steps are needed.

`S(n)`

We use the periodicity:

$$10^{22}\equiv 1 \pmod{23}.$$

Thus there are only $22$ distinct position-types.

$$n=22q+r,$$

then factor $P_j$ appears:

$$q+\mathbf 1_{j<r}$$

times.

Finally we extract the coefficient:

ans[(23, 0)]

which corresponds exactly to:

digit sum $23$,
residue $0\pmod{23}$.

4. Verification

The program reproduces the supplied value:

$$S(9)=263626.$$

For $S(42)$, the problem statement gives

$$6377168878570056.$$

Modulo $10^9$, this is:

$$878570056,$$

which exactly matches the program output.

So the method is consistent with the provided checks.

Therefore,

$$S(11^{12}) \bmod 10^9 = 789184709.$$

Answer: 789184709

Project Euler Problem 294