Project Euler Problem 416

Solution

Answer: 898082747

Let $L=n-1$.

A complete outward or homeward trip is just a composition of $L$ into parts $1,2,3$.

For $m$ round-trips there are

$$p=2m$$

independent monotone paths from $0$ to $L$.

The key observation is that we can scan the board from left to right and only need to know, for each path, how far ahead its next landing point is.

For a single path at position $t$:

state $0$: the path lands on $t$,
state $1$: next landing is at $t+1$,
state $2$: next landing is at $t+2$.

If a path is currently in state $0$, then after landing it chooses a jump of length $1,2,3$, producing new states $0,1,2$ respectively.

For all $p$ paths together, the entire system is described only by the counts

$$(a,b,c),$$

where

$a$ paths are in state $0$,
$b$ paths are in state $1$,
$c$ paths are in state $2$,

with

$$a+b+c=p.$$

Thus the number of automaton states is

$$\binom{p+2}{2}.$$

For $p=20$, this is only $231$.

At an interior square $t$, the square is visited iff $a>0$.

We therefore augment the state with a bit telling whether we have already seen one unvisited square.

So the total number of states is only

$$231 \times 2 = 462.$$

This gives a finite linear transition system, so $F(m,n)$ can be computed by fast matrix exponentiation in $O(\log n)$.

Because the required modulus is

$$10^9 = 2^9\cdot 5^9,$$

we compute separately modulo $512$ and modulo $1953125$, where 64-bit matrix multiplication is safe, then combine with CRT.

The transition rule is:

If the current state is $(a,b,c)$, then among the $a$ landing paths choose

$x$ to jump length $1$,
$y$ to jump length $2$,
$z$ to jump length $3$,

with $x+y+z=a$.

The next state is

$$(b+x,; c+y,; z),$$

with multiplicity

$$\binom{a}{x,y,z}.$$

Python implementation:

import numpy as np
from math import comb

MOD1 = 512
MOD2 = 1953125

def solve_mod(mod):
    m = 10
    p = 2 * m
    L = 10**12 - 1

    # Enumerate all (a,b,c), a+b+c=p
    states = []
    idx = {}

    for a in range(p + 1):
        for b in range(p + 1 - a):
            c = p - a - b
            idx[(a, b, c)] = len(states)
            states.append((a, b, c))

    S = len(states) * 2   # extra bit = already used one unvisited square

    A = np.zeros((S, S), dtype=np.int64)
    B = np.zeros((S, S), dtype=np.int64)

    # Build transition matrices
    for si, (a, b, c) in enumerate(states):

        trans = []

        for x in range(a + 1):
            for y in range(a - x + 1):
                z = a - x - y

                mult = comb(a, x) * comb(a - x, y)
                mult %= mod

                nxt = (b + x, c + y, z)
                ni = idx[nxt]

                trans.append((ni, mult))

        for used in range(2):

            row = si * 2 + used

            for ni, mult in trans:

                # First step: no interior-square check
                A[row, ni * 2 + used] += mult
                A[row, ni * 2 + used] %= mod

                # Interior steps
                nused = used + (1 if a == 0 else 0)

                if nused <= 1:
                    B[row, ni * 2 + nused] += mult
                    B[row, ni * 2 + nused] %= mod

    start = idx[(p, 0, 0)] * 2

    vec = np.zeros(S, dtype=np.int64)
    vec[start] = 1

    # First step (position 0 is not interior)
    vec = (vec @ A) % mod

    # Remaining steps
    e = L - 1
    P = B.copy()

    while e:
        if e & 1:
            vec = (vec @ P) % mod

        e >>= 1

        if e:
            P = (P @ P) % mod

    ans = (vec[start] + vec[start + 1]) % mod
    return int(ans)

a1 = solve_mod(MOD1)
a2 = solve_mod(MOD2)

# CRT
inv = pow(MOD1, -1, MOD2)

answer = (
    a1
    + MOD1 * (((a2 - a1) * inv) % MOD2)
) % (10**9)

print(answer)

Verification on the examples:

$F(1,3)=4$
$F(1,4)=15$
$F(1,5)=46$
$F(2,3)=16$
$F(2,100)\bmod 10^9 = 429619151$

all match the statement.

Therefore,

Answer: 898082747