Problem 324: Building a Tower

Project Euler Problem 324 Solution

EulerSolve provides an optimized solution for Project Euler Problem 324, Building a Tower, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We want to tile a $3\times 3\times n$ tower with $2\times 1\times 1$ blocks, allowing all axis-aligned orientations, and count the number of valid towers modulo $$M=100000007.$$ Let $f(n)$ be that count. The challenge is that the requested index is not merely large: it is $$N=10^{10000}.$$ So even a fast layer-by-layer dynamic program is still hopeless unless we compress the sequence very aggressively. Mathematical Approach 1) A parity observation. The tower volume is $9n$. Each block covers 2 unit cubes. Therefore $f(n)=0$ whenever $9n$ is odd, i.e. whenever $n$ is odd. This is not enough to solve the problem, but it is a good first sanity check on the sequence. 2) Slice the tower by layers. Process the tower one $3\times 3$ layer at a time. A single layer has 9 cells, so we encode a layer state by a 9-bit mask $$s\in[0,2^9).$$ A bit of $s$ is 1 if that cell is already occupied by a domino that started in the previous layer and sticks into the current one. Since there are 9 cells, the total number of masks is $$2^9=512.$$ 3) What one transition means. Fix an incoming mask $s$. We must finish filling the current $3\times 3$ layer completely....

Detailed mathematical approach

Problem Summary

Let $f(n)$ be the number of ways to fill a $3\times3\times n$ tower completely with blocks of size $2\times1\times1$. A block may be turned in any direction, so it covers two cells sharing a face, and rotations or reflections of the tower count as distinct. With $q=100000007$ the statement gives $f(2)=229$ and $f(4)=117805$, together with $f(10)\bmod q=96149360$, $f(10^3)\bmod q=24806056$ and $f(10^6)\bmod q=30808124$. We must find $f(10^{10000})\bmod q$. ^[1]

The whole solution in one sentence: count the ways one layer can meet the next, use those counts to work out the first $1301$ heights, read off the fixed rule every later height must obey, and apply it once to leap to the height asked for.

Mathematical Approach

1. Two tiny cases, and how large the answer is

The empty tower has one filling, the one using no block, so $f(0)=1$. A tower of height $n$ has $9n$ cells and every block covers $2$, so $9n$ must be even; as $9$ is odd, $f(n)=0$ for every odd $n$, which is why the statement quotes $f(2)$ and $f(4)$. Section 6 turns that parity into nineteen vanishing recurrence coefficients.

The modulus is no decoration either. As Section 7 shows, $f(n+2)/f(n)$ tends to $\lambda=552.0631464\ldots$ through even heights, so $f(n)$ has asymptotically $\tfrac12n\log_{10}\lambda\approx1.3709944\,n$ digits and the value asked for would need some $1.371\times10^{10000}$ of them. Only a residue can be asked for, and only an algebraic shortcut reaches that index.

2. Layers and states

Cut the tower into its $n$ layers of size $3\times3$ and fill them from the bottom up. Inside a layer a block can lie flat in $12$ places, the $6$ horizontal and $6$ vertical neighbouring pairs; otherwise it stands upright, its lower cell in this layer and its upper cell directly above.

So the only thing a layer passes upwards is which cells of the next layer are already occupied by upright blocks. Encode that set as a $9$-bit mask $s$, giving $2^9=512$ states, and let $T_{s,t}$ count the ways to fill one layer when the cells of $s$ arrive filled from below and exactly the cells of $t$ protrude upwards. A filling of the tower is a walk $\emptyset\to s_1\to\cdots\to s_{n-1}\to\emptyset$: nothing protrudes into the bottom layer and nothing sticks out of the top. Since the $(i,j)$ entry of the $n$-th power of a weighted adjacency matrix is the total weight of the walks of length $n$ from $i$ to $j$, ^[3]

$$f(n)=\bigl(T^{n}\bigr)_{\emptyset,\emptyset}.$$

3. What one entry of the matrix counts

Fix $s$ and $t$. Each cell is used in exactly one of three ways: it lies in $s$ and came filled from below, it lies in $t$ and carries an upright block, or a flat block covers it. The first two must not clash, so $T_{s,t}=0$ unless $s\cap t=\emptyset$, and the flat blocks must cover the rest exactly. A matching is a set of neighbouring cell pairs no two of which share a cell; it is perfect on a cell set $X$ when every cell of $X$ lies in a pair, the empty set having one such matching. Writing $\mu(X)$ for the number of perfect matchings on $X$,

$$T_{s,t}=\begin{cases}0,&s\cap t\ne\emptyset,\\\mu(\overline{s\cup t}),&s\cap t=\emptyset,\end{cases}$$

with $\overline{X}$ the complement of $X$ among the nine cells. Swapping $s$ and $t$ leaves both the test $s\cap t=\emptyset$ and the set $s\cup t$ untouched, so either case returns the same value: $T$ is symmetric. (The union alone does not fix the entry — for $u$ a single corner cell, $T_{\emptyset,u}=4$ while $T_{u,u}=0$ — but the disjointness test is symmetric too, which is all the argument needs.) Symmetry turns the second power into a sum of squares:

$$f(2)=\sum_{t}T_{\emptyset,t}T_{t,\emptyset}=\sum_{t}T_{\emptyset,t}^{2}.$$

The row $T_{\emptyset,\cdot}$ has $98$ non-zero entries: $77$ equal to $1$, thirteen to $2$, four to $3$ and four to $4$. Their sum $131$ counts the fillings of the bottom layer and, by the formula, the matchings of the $3\times3$ grid graph including the empty one, as a pass over all $2^{12}$ edge subsets confirms. Their squares sum to $77+4\cdot13+9\cdot4+16\cdot4=229$, the checkpoint $f(2)=229$. Also $T_{\emptyset,\emptyset}=0$: nine cells cannot be covered by flat blocks alone.

4. The same method on a tower we can check by hand

Shrink the cross-section to $2\times2$ and let $h(n)$ count fillings of a $2\times2\times n$ tower. A layer has four cells, so no parity obstruction, and $2^4=16$ states. Reading Section 3 for $s=\emptyset$: the state $\emptyset$ is reached in two ways, the two perfect matchings of the square; each of the four states made of one neighbouring pair in one way; and the full state in one way, using no flat block. The total $2+4+1=7$ counts the matchings of the $2\times2$ grid, whose edges form a cycle, and $h(1)=T_{\emptyset,\emptyset}=2$ — correct, since a $2\times2\times1$ slab takes two blocks either side by side or one behind the other. The sum of squares gives

$$h(2)=\sum_tT_{\emptyset,t}^2=2^2+4\cdot1^2+1^2=9,$$

and continuing the walk count yields $1,2,9,32,121,450,1681$ for $n=0,\dots,6$. Although the matrix has size $16$, these obey a rule of order only three,

$$h(n)=3h(n-1)+3h(n-2)-h(n-3),$$

as one checks by hand: $3\cdot9+3\cdot2-1=32$ and $3\cdot32+3\cdot9-2=121$. It holds for ever, as three counts show. Let $a_n$, $b_n$, $c_n$ count partial towers of height $n$ whose outgoing state is empty, one fixed neighbouring pair, and the full square; the four pairs have equal counts by symmetry, and $h(n)=a_n$. The transitions above give

$$\begin{aligned}a_{n+1}&=2a_n+4b_n+c_n,\\ b_{n+1}&=a_n+b_n,\\ c_{n+1}&=a_n,\end{aligned}$$

from $a_0=1$, $b_0=c_0=0$, and eliminating $b_n$ and $c_n$ leaves exactly that rule. Nothing shorter works: an order-two rule fitted to $1,2,9,32$ is forced to the coefficients $14/5$ and $17/5$, and then predicts $601/5$ where $h(4)=121$; an order-one rule gives $h(2)=4$. A matrix of size $16$ can thus generate a sequence of order $3$, and only computation reveals which order to expect.

5. Why a recurrence must exist, and how long it can be

By the Cayley–Hamilton theorem $T$ satisfies its own characteristic polynomial, of degree $512$; ^[4] multiplying that identity by $T^{m}$ and reading the $(\emptyset,\emptyset)$ entry turns its coefficients into a linear recurrence for $f$. Any sequence $u^{\top}A^nv$ with $A$ of size $d$ behaves so, and its order is the length of its shortest recurrence. ^[5]

The bound sharpens before any computation. The eight symmetries of the square permute the nine cells, hence the masks; $\mu$ is invariant under them, so each commutes with $T$, and the starting vector, concentrated on $\emptyset$, is fixed by all eight. Every vector the recursion produces is therefore constant on orbits — confirmed over $24$ layers — and the masks form $102$ orbits, so the order is at most $102$.

The exact order needs a certificate, available once Section 7 has produced the integer coefficients $a_j$. Put

$$Q(x)=x^{38}-\sum_{j=1}^{19}a_jx^{38-2j}.$$

Exact integer arithmetic gives $Q(T)e_{\emptyset}=0$ in all $512$ coordinates; multiplying on the left by $e_{\emptyset}^{\top}T^{n}$ proves the order-$38$ recurrence for every $n\ge0$, over the integers and so modulo $q$ too. For the lower bound let $H$ be the $38\times38$ matrix $H_{i,j}=f(i+j)$: its determinant $6286188$ modulo $q$ is non-zero, whereas a recurrence of order $L\lt38$ would make column $L$ of $H$ a combination of earlier columns. The order is therefore exactly $38$; neither certificate is a step the programs perform.

6. Only even lags survive

If $T_{s,t}\ne0$ the flat blocks cover $9-|s|-|t|$ cells, which must be even, so $|s|+|t|$ is odd at every non-zero entry. With $D$ diagonal, $D_{s,s}=(-1)^{|s|}$, we get $(DTD)_{s,t}=(-1)^{|s|+|t|}T_{s,t}$, that is

$$DTD=-T,\qquad D^2=I.$$

Two consequences. First $f(n)=(DT^{n}D)_{\emptyset,\emptyset}=(-1)^nf(n)$, re-proving $f(n)=0$ for odd $n$. Second, if a polynomial annihilates the sequence then so does the same polynomial in $-x$, up to sign; the minimal one $Q$ is monic of even degree $38$, so $Q(x)=Q(-x)$ and $Q$ is a polynomial in $x^2$. Writing $Q(x)=P(x^2)$ with $\deg P=19$, only even lags appear:

$$f(n)=\sum_{j=1}^{19}a_j\,f(n-2j).$$

Equivalently $g(m)=f(2m)$ has order $19$, and half the $38$ coefficients vanish before anything is computed.

7. Berlekamp–Massey finds the coefficients

A recurrence predicts each term from earlier ones, and Berlekamp–Massey repairs such a rule wherever it mispredicts. Holding coefficients $C_1,\dots,C_L$ it measures the discrepancy

$$d_n=f(n)+\sum_{i=1}^{L}C_if(n-i),$$

which is $0$ exactly when the rule gets term $n$ right. Otherwise a saved earlier rule is shifted, scaled by the ratio of $d_n$ to a stored non-zero discrepancy and added in, the order rising only when necessary; the result is the shortest rule fitting the prefix, at $O(K^2)$ field operations for $K$ terms. If the sequence has order at most $D$, its first $2D$ terms already determine the minimal rule. ^[6] The field is $\mathbb{F}_q$, the residues modulo the prime $q=100000007$; as $q-1=2\cdot491\cdot101833$, Fermat's little theorem gives the inverse of $a\ne0$ as $a^{q-2}$. Cayley–Hamilton supplies $D=512$, so $1024$ terms would suffice; the program feeds in $1301$.

The output has length $38$ with all nineteen odd-lag coefficients zero, as Section 6 predicts. The even-lag ones are residues, but they reduce an integer recurrence: solving the $19\times19$ system over the rationals with exact values of $g$ yields integers whose residues match the output,

$a_1=675$, $a_2=-73471$, $a_3=3221189$;
$a_4=-72583272$, $a_5=925908264$;
$a_6=-6971103216$, $a_7=31523058272$;
$a_8=-86171526770$, $a_9=142604534086$;
$a_j=-a_{19-j}$ for $10\le j\le18$, and $a_{19}=1$.

That mirror symmetry gives free checks. Summed over $j\le18$ it yields $0$, so $\sum_ja_j=a_{19}=1$ and $P(1)=0$; it says exactly $x^{19}P(1/x)=-P(x)$, equivalently $P(x)=(x-1)S(x)$ with $S$ palindromic of degree $18$. That pairs the roots as $r,1/r$ with $1$ paired to itself, but reality comes from the matrix: $g(m)=e_{\emptyset}^{\top}(T^2)^me_{\emptyset}$ with $T$ real symmetric and invertible, so $T^2$ is positive definite and an orthonormal eigenbasis writes $g(m)=\sum_jc_jr_j^{\,m}$ with $c_j\gt0$ for each contributing eigenvalue. The roots of $P$ are exactly those eigenvalues, hence real, positive and distinct, and the largest, $\lambda\approx552.063146417369$, is the limit of $g(m+1)/g(m)$ used in Section 1 — the ratio $g(12)/g(11)\approx552.063146414713$ lies within $2.7\times10^{-9}$ of it.

Two more checks cost nothing. As $a_{19}=1$ the rule runs backwards in integers, and $x^{19}P(1/x)=-P(x)$ makes the reversed index obey the same rule, so $g(m)$ and $g(-1-m)$ are both solutions; they agree for $0\le m\le20$ by exact substitution, and $19$ consecutive agreements fix a solution, so $g(-1-m)=g(m)$ always — the three-dimensional counterpart of Propp's reciprocity $T(m,-2-n)=\pm T(m,n)$ for $m\times n$ rectangles. ^[9] And OEIS A028452 publishes a recurrence of order $23$ whose characteristic polynomial is exactly $P(x)\,(x-1)^4$, a redundant factor its own degree-$19$ generating-function denominator betrays. ^[2]

8. Jumping to an astronomical index

Let $S$ be the shift operator, $(Sf)(n)=f(n+1)$. Section 6 says precisely $Q(S)f=0$ with $\deg Q=38$. Divide $x^N$ by $Q$ with remainder,

$$x^N\equiv\sum_{i=0}^{37}b_ix^i\pmod{Q(x)},$$

and substitute $S$ for $x$: the multiple of $Q(S)$ kills $f$, so $S^Nf=\sum_ib_iS^if$, and evaluating at index $0$ gives the method in one line,

$$f(N)=\sum_{i=0}^{37}b_i\,f(i).$$

The values $f(0),\dots,f(37)$ come from the layer recursion, so everything reduces to $x^N\bmod Q(x)$. The $2\times2$ tower shows the mechanics: there $Q$ is $x^3-3x^2-3x+1$, so $x^3\equiv3x^2+3x-1$, and one more multiplication by $x$ and substitution give $x^4\equiv12x^2+8x-3$, whence $h(4)=12\cdot9+8\cdot2-3\cdot1=121$ as listed above. Powers come from exponentiation by squaring in $\mathbb{F}_q[x]/(Q)$, which needs nothing but an associative product. ^[7] The published loop squares once for each of the $\lfloor\log_2N\rfloor+1$ bits and multiplies into the result once per one-bit; its last squaring is never used, so a sharper loop would need $\lfloor\log_2N\rfloor$ squarings. Each product multiplies two polynomials of degree below $38$ and folds $x^{38},\dots,x^{74}$ back down with the recurrence, at $O(38^2)$ modular operations.

Section 6 pays again. As $Q$ is a polynomial in $x^2$, reduction never mixes exponents of different parity, so $x^N\bmod Q$ is non-zero only at the $19$ indices congruent to $N$ modulo $2$ — confirmed for $N=40,41,97,1000$ and $1001$. The code keeps arrays of length $38$ anyway, but its zero tests cut the convolution from $38^2=1444$ products to $19^2=361$; the reduction gains less, as each non-zero high coefficient still runs over all $38$ recurrence coefficients. One instrumented product gives $361+684=1045$ modular products against $1444+1406=2850$ dense: a saving near $2.7$, not $4$.

9. Reading the bits of $10^{10000}$

The exponent has $10001$ decimal digits and fits no machine integer. In binary it has $33220$ bits, $11683$ of them $1$, so the loop performs $33220+11683=44903$ polynomial products. In C++ the bits come from halving the decimal string repeatedly: at each digit $d$ a carried remainder $r\in\{0,1\}$ makes the dividend $10r+d$, whose quotient becomes the next digit and whose remainder is carried on; the remainder left at the end of a pass is the least significant bit. Repeating until the string is $0$ yields the bits from least significant upwards, exactly the order the square-and-multiply loop wants.

How the Code Works

The C++ program ^[10] fixes MOD${}=q$, CELLS${}=9$, FULL_MASK${}=511$ and STATE_COUNT${}=512$, and works in 32-bit residues: mul_mod widens to 64 bits before the remainder, and inv_mod calls pow_mod with exponent $q-2$, the Fermat inverse Berlekamp–Massey needs.

dfs_layer_fill implements Section 3 without building the grid graph. Given the cells already filled it takes the lowest free one — __builtin_ctz of the free bits — and considers up to three legal placements: let that cell carry an upright block, which sets its bit in the outgoing mask; pair it with its right neighbour; or pair it with the neighbour below. Since that cell's role is settled before any later cell is touched, each filling is generated once, and those two flat branches suffice: a block to its left or above would reuse a decided cell. build_transitions runs the search from all $512$ incoming masks and stores each row as a list of (next, count) pairs, $3658$ in all.

generate_sequence is the walk count of Section 2: from current[0] = 1, each of its $1300$ steps scans the $512$ states, skips zero entries, pushes each value along the stored transitions into next, and records current[0] as the value of $f$ at that step. berlekamp_massey maintains the connection polynomial C, its previous version B, the order L and the gap m, and returns the $38$ coefficients in forward form; validate_recurrence checks the $1263$ recurrence equations at indices $38$ to $1300$, and the first $38$ terms become initial.

combine_polynomials is the product of Section 8: a convolution of two length-$38$ arrays into a buffer of length $76$, then a downward sweep replacing each coefficient of $x^{i}$, $i\ge38$, by its expansion through the recurrence. linear_recurrence_nth is the square-and-multiply loop over the exponent bits, and bits_from_decimal supplies them by the string halving of Section 9. nth_term short-circuits any $n\lt38$ by reading initial, and run_checkpoints, skippable with --skip-checkpoints, verifies the five statement values.

Python and Java take the same steps, with three real differences. Python works in arbitrary-precision integers, so bits_from_decimal parses the string with int() and then repeatedly reads value_int & 1 and shifts right; CPython's default limit on decimal conversions is $4300$ digits, so the main block first raises it with sys.set_int_max_str_digits(20000). ^[8] Java keeps residues in int and widens each product to long in mulMod, safe because $(q-1)^2\lt2^{63}$, and reads the exponent from a BigInteger with testBit and shiftRight. Next, berlekampMassey in Java guards its inner read with i < C.size() where C++ and Python index C[i] for every $i\le L$ unconditionally; instrumenting the run shows C always holds more than L coefficients, so the condition is true on every iteration and no out-of-range read occurs in any language, though in C++ one would be undefined behaviour rather than a wrong answer. Finally only C++ runs validate_recurrence and the checkpoints.

Complexity and Verification

Building the table runs $512$ filling searches with $3940$ leaves, scans $512^2$ candidate entries and stores $3658$ transitions, in well under a millisecond. The $1300$ layer steps inspect $665600$ states and execute $1391267$ transition updates; the product $1300\times3658=4\,755\,400$ is only the bound before zero states are skipped. Berlekamp–Massey, of cost $O(K^2)$ for $K$ terms, performs $47994$ products in its discrepancy sums, $705$ correction products and $39$ inversions. The evaluation dominates: $44903$ polynomial products of $O(L^2)$ operations each with $L=38$, measured at about $1009$ products per call in Section 8, so some $4.5\times10^{7}$ modular multiplications. With $B$ the bit length of $N$ that part is $O(L^2B)$, and converting the exponent adds $O(B^2)$ digit or word operations, each of the $B$ halvings costing time proportional to the length. The transitions, the $1301$-term prefix, the length-$38$ arrays and the $33220$ stored bits need a few hundred kilobytes, excluding each language runtime.

Measured on one machine, with GCC at -O2, OpenJDK 17 and CPython 3.11, C++ takes about $0.3$ s and Java about $0.2$ s, while Python takes about $2.6$ s, of which only $0.005$ s goes to the transitions, $0.07$ s to the layer steps and $0.002$ s to Berlekamp–Massey. The three outputs are byte-for-byte identical and the C++ run passes every checkpoint. Note what each checkpoint tests: $2$, $4$ and $10$ lie below $38$, so those three only read the stored prefix and test the transfer matrix, while $f(10^3)$ and $f(10^6)$ exercise the recurrence and the polynomial power.

The following were verified for this article, outside the published programs. A direct cell-by-cell backtracking over the actual box, using no transfer matrix, gives $f(1)=0$, $f(2)=229$, $f(3)=0$ and $f(4)=117805$; an exact big-integer layer recursion gives $f(6)=64647289$, $f(8)=35669566217$, $f(10)=19690797527709$ and $f(12)=10870506600976757$, agreeing with the published table of OEIS A028452 as far as it goes. ^[2] The entry formula of Section 3 was checked for all $512^2$ pairs against a separate perfect-matching routine, confirming the symmetry, the row sum $131$ and that $|s|+|t|$ is odd wherever the entry is non-zero. It also gives $\det T=1$: with $J$ the permutation matrix of the complement map, $(TJ)_{s,u}=\mu(u\setminus s)$ for $s\subseteq u$ and $0$ otherwise, so $TJ$ is triangular for inclusion with unit diagonal, while complementation pairs the $512$ masks into $256$ transpositions, giving $\det J=1$. The two certificates of Section 5 were computed, the integer recurrence reproduces every exact term up to $f(80)$, and all five statement checkpoints were reproduced through the published evaluation path itself.

Footnotes and References

The references below supply the background facts used above. The transfer matrix and its entry formula, the symmetry and parity identities, the shift-operator derivation and the numerical certificates are derived explicitly in this article.

Project Euler 324 — Building a Tower. The official statement defines f(n) as the number of fillings of a 3 × 3 × n tower by 2 × 1 × 1 blocks, lets the blocks rotate freely while counting rotations and reflections of the tower as distinct, and gives f(2) = 229, f(4) = 117805 and the residues of f(10), f(10³) and f(10⁶) mod 100000007.
OEIS A028452 — perfect matchings in P₃ × P₃ × P₂ₙ. The entry notes that the sequence also counts tilings of a 3 × 3 × 2n box by 1 × 1 × 2 bricks, tabulates 1, 229, 117805, 64647289, 35669566217, … — the subsequence f(2n) — and gives a recurrence of order 23 with a generating function of denominator degree 19.
Wikipedia — Adjacency matrix. In the n-th power of an adjacency matrix, "the element (i, j) gives the number of (directed or undirected) walks of length n from vertex i to vertex j", and integer entries record edge multiplicities.
Wikipedia — Cayley–Hamilton theorem. Every square matrix over a commutative ring satisfies its own characteristic polynomial det(λI − A), monic of degree n, so Aⁿ is a combination of lower powers.
Wikipedia — Constant-recursive sequence. The order is the smallest d admitting a recurrence of that length, and order at most d is characterised by a representation sₙ = u Aⁿ v with A of size d.
Wikipedia — Berlekamp–Massey algorithm. The algorithm returns the shortest linear recurrence, equivalently the minimal polynomial, of a linearly recurrent sequence; it needs every non-zero field element to be invertible, and its discrepancies vanish once twice the order has been read.
Wikipedia — Exponentiation by squaring. Powers can be formed in any semigroup, polynomials included, with about ⌊log₂ n⌋ squarings plus one multiplication per one-bit, the iterative form given there consuming the bits from least significant upwards.
Python — sys.set_int_max_str_digits documents the limit on conversions between integers and decimal strings, 4300 digits by default in CPython 3.11. Java — BigInteger documents its arbitrary-precision integers and the methods testBit(n) and shiftRight(n).
J. Propp, A Reciprocity Theorem for Domino Tilings, Electronic Journal of Combinatorics 8 (2001), #R18. For the number T(m, n) of domino tilings of an m × n rectangle the counts obey a recurrence in n for fixed m, hence extend to negative n, and satisfy T(m, −2 − n) = ±T(m, n), the sign being negative only for m ≡ 2 (mod 4) and odd n.
C++ — Euler324.cpp. The linked immutable revision contains dfs_layer_fill, build_transitions, generate_sequence, berlekamp_massey, combine_polynomials and linear_recurrence_nth. It is the implementation record for this article; the independent checks above are not packaged in it.

Problem 324 source code

C++

#include <array>
#include <cstdint>
#include <functional>
#include <iostream>
#include <string>
#include <utility>
#include <vector>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;

constexpr int MOD = 100000007;
constexpr int CELLS = 9;
constexpr int FULL_MASK = (1 << CELLS) - 1;
constexpr int STATE_COUNT = 1 << CELLS;

int add_mod(const int a, const int b) {
    int x = a + b;
    if (x >= MOD) {
        x -= MOD;
    }
    return x;
}

int sub_mod(const int a, const int b) {
    int x = a - b;
    if (x < 0) {
        x += MOD;
    }
    return x;
}

int mul_mod(const int a, const int b) {
    return static_cast<int>((static_cast<i64>(a) * static_cast<i64>(b)) % MOD);
}

int pow_mod(int base, u64 exp) {
    int result = 1;
    while (exp > 0) {
        if ((exp & 1ULL) != 0ULL) {
            result = mul_mod(result, base);
        }
        base = mul_mod(base, base);
        exp >>= 1ULL;
    }
    return result;
}

int inv_mod(const int a) {
    return pow_mod(a, static_cast<u64>(MOD - 2));
}

void dfs_layer_fill(const int filled_mask, const int next_mask, std::array<int, STATE_COUNT>& ways) {
    if (filled_mask == FULL_MASK) {
        ++ways[static_cast<std::size_t>(next_mask)];
        return;
    }

    const int free_bits = (~filled_mask) & FULL_MASK;
    const int bit = __builtin_ctz(static_cast<unsigned int>(free_bits));

    // Place a domino across layers: marks this cell in next_mask.
    dfs_layer_fill(filled_mask | (1 << bit), next_mask | (1 << bit), ways);

    const int r = bit / 3;
    const int c = bit % 3;

    // Place inside the current layer along columns.
    if (c < 2 && ((filled_mask & (1 << (bit + 1))) == 0)) {
        dfs_layer_fill(filled_mask | (1 << bit) | (1 << (bit + 1)), next_mask, ways);
    }

    // Place inside the current layer along rows.
    if (r < 2 && ((filled_mask & (1 << (bit + 3))) == 0)) {
        dfs_layer_fill(filled_mask | (1 << bit) | (1 << (bit + 3)), next_mask, ways);
    }
}

std::array<std::vector<std::pair<int, int>>, STATE_COUNT> build_transitions() {
    std::array<std::vector<std::pair<int, int>>, STATE_COUNT> transitions;

    for (int mask = 0; mask < STATE_COUNT; ++mask) {
        std::array<int, STATE_COUNT> ways{};
        dfs_layer_fill(mask, 0, ways);

        auto& edges = transitions[static_cast<std::size_t>(mask)];
        edges.reserve(64);
        for (int nxt = 0; nxt < STATE_COUNT; ++nxt) {
            const int cnt = ways[static_cast<std::size_t>(nxt)];
            if (cnt > 0) {
                edges.emplace_back(nxt, cnt % MOD);
            }
        }
    }

    return transitions;
}

std::vector<int> generate_sequence(const int terms,
                                   const std::array<std::vector<std::pair<int, int>>, STATE_COUNT>& transitions) {
    std::vector<int> seq(static_cast<std::size_t>(terms + 1), 0);

    std::array<int, STATE_COUNT> current{};
    std::array<int, STATE_COUNT> next{};
    current[0] = 1;
    seq[0] = 1;

    for (int step = 1; step <= terms; ++step) {
        next.fill(0);
        for (int mask = 0; mask < STATE_COUNT; ++mask) {
            const int value = current[static_cast<std::size_t>(mask)];
            if (value == 0) {
                continue;
            }
            for (const auto& [nxt, cnt] : transitions[static_cast<std::size_t>(mask)]) {
                const int add = mul_mod(value, cnt);
                next[static_cast<std::size_t>(nxt)] = add_mod(next[static_cast<std::size_t>(nxt)], add);
            }
        }
        current = next;
        seq[static_cast<std::size_t>(step)] = current[0];
    }

    return seq;
}

std::vector<int> berlekamp_massey(const std::vector<int>& sequence) {
    std::vector<int> C{1};
    std::vector<int> B{1};

    int L = 0;
    int m = 1;
    int b = 1;

    for (int n = 0; n < static_cast<int>(sequence.size()); ++n) {
        int d = sequence[static_cast<std::size_t>(n)];
        for (int i = 1; i <= L; ++i) {
            d = add_mod(d, mul_mod(C[static_cast<std::size_t>(i)], sequence[static_cast<std::size_t>(n - i)]));
        }

        if (d == 0) {
            ++m;
            continue;
        }

        const std::vector<int> T = C;
        const int coef = mul_mod(d, inv_mod(b));

        if (static_cast<int>(C.size()) < static_cast<int>(B.size()) + m) {
            C.resize(static_cast<std::size_t>(static_cast<int>(B.size()) + m), 0);
        }

        for (int i = 0; i < static_cast<int>(B.size()); ++i) {
            const int idx = i + m;
            C[static_cast<std::size_t>(idx)] = sub_mod(C[static_cast<std::size_t>(idx)],
                                                       mul_mod(coef, B[static_cast<std::size_t>(i)]));
        }

        if (2 * L <= n) {
            L = n + 1 - L;
            B = T;
            b = d;
            m = 1;
        } else {
            ++m;
        }
    }

    std::vector<int> recurrence(static_cast<std::size_t>(L), 0);
    for (int i = 1; i <= L; ++i) {
        recurrence[static_cast<std::size_t>(i - 1)] = (MOD - C[static_cast<std::size_t>(i)]) % MOD;
    }
    return recurrence;
}

std::vector<int> combine_polynomials(const std::vector<int>& a, const std::vector<int>& b,
                                     const std::vector<int>& recurrence) {
    const int k = static_cast<int>(recurrence.size());
    std::vector<int> prod(static_cast<std::size_t>(2 * k), 0);

    for (int i = 0; i < k; ++i) {
        if (a[static_cast<std::size_t>(i)] == 0) {
            continue;
        }
        for (int j = 0; j < k; ++j) {
            if (b[static_cast<std::size_t>(j)] == 0) {
                continue;
            }
            const int idx = i + j;
            prod[static_cast<std::size_t>(idx)] = add_mod(
                prod[static_cast<std::size_t>(idx)],
                mul_mod(a[static_cast<std::size_t>(i)], b[static_cast<std::size_t>(j)]));
        }
    }

    for (int i = 2 * k - 2; i >= k; --i) {
        const int coef = prod[static_cast<std::size_t>(i)];
        if (coef == 0) {
            continue;
        }
        for (int j = 1; j <= k; ++j) {
            const int idx = i - j;
            prod[static_cast<std::size_t>(idx)] = add_mod(
                prod[static_cast<std::size_t>(idx)], mul_mod(coef, recurrence[static_cast<std::size_t>(j - 1)]));
        }
    }

    prod.resize(static_cast<std::size_t>(k));
    return prod;
}

int linear_recurrence_nth(const std::vector<int>& initial, const std::vector<int>& recurrence,
                          const std::vector<int>& bits_lsb_first) {
    const int k = static_cast<int>(recurrence.size());

    std::vector<int> result(static_cast<std::size_t>(k), 0);
    result[0] = 1;

    std::vector<int> x_poly(static_cast<std::size_t>(k), 0);
    if (k == 1) {
        x_poly[0] = recurrence[0];
    } else {
        x_poly[1] = 1;
    }

    for (const int bit : bits_lsb_first) {
        if (bit != 0) {
            result = combine_polynomials(result, x_poly, recurrence);
        }
        x_poly = combine_polynomials(x_poly, x_poly, recurrence);
    }

    int answer = 0;
    for (int i = 0; i < k; ++i) {
        answer = add_mod(answer, mul_mod(result[static_cast<std::size_t>(i)], initial[static_cast<std::size_t>(i)]));
    }
    return answer;
}

std::vector<int> bits_from_u64(u64 n) {
    std::vector<int> bits;
    bits.reserve(64);
    while (n > 0ULL) {
        bits.push_back(static_cast<int>(n & 1ULL));
        n >>= 1ULL;
    }
    if (bits.empty()) {
        bits.push_back(0);
    }
    return bits;
}

std::vector<int> bits_from_decimal(std::string value) {
    std::vector<int> bits;
    bits.reserve(value.size() * 4);

    while (!(value.size() == 1 && value[0] == '0')) {
        int carry = 0;
        std::string next;
        next.reserve(value.size());

        for (char ch : value) {
            const int cur = carry * 10 + static_cast<int>(ch - '0');
            const int q = cur / 2;
            carry = cur % 2;
            if (!next.empty() || q != 0) {
                next.push_back(static_cast<char>('0' + q));
            }
        }

        bits.push_back(carry);
        value = next.empty() ? "0" : next;
    }

    if (bits.empty()) {
        bits.push_back(0);
    }
    return bits;
}

int nth_term(const u64 n, const std::vector<int>& initial, const std::vector<int>& recurrence) {
    if (n < initial.size()) {
        return initial[static_cast<std::size_t>(n)];
    }
    return linear_recurrence_nth(initial, recurrence, bits_from_u64(n));
}

int nth_term_decimal(const std::string& n_decimal, const std::vector<int>& initial,
                     const std::vector<int>& recurrence) {
    return linear_recurrence_nth(initial, recurrence, bits_from_decimal(n_decimal));
}

bool validate_recurrence(const std::vector<int>& sequence, const std::vector<int>& recurrence) {
    const int k = static_cast<int>(recurrence.size());
    for (int n = k; n < static_cast<int>(sequence.size()); ++n) {
        int expected = 0;
        for (int i = 1; i <= k; ++i) {
            expected = add_mod(expected, mul_mod(recurrence[static_cast<std::size_t>(i - 1)],
                                                 sequence[static_cast<std::size_t>(n - i)]));
        }
        if (expected != sequence[static_cast<std::size_t>(n)]) {
            return false;
        }
    }
    return true;
}

bool run_checkpoints(const std::vector<int>& initial, const std::vector<int>& recurrence) {
    if (nth_term(2ULL, initial, recurrence) != 229) {
        std::cerr << "Checkpoint failed: f(2)\n";
        return false;
    }
    if (nth_term(4ULL, initial, recurrence) != 117805) {
        std::cerr << "Checkpoint failed: f(4)\n";
        return false;
    }
    if (nth_term(10ULL, initial, recurrence) != 96149360) {
        std::cerr << "Checkpoint failed: f(10)\n";
        return false;
    }
    if (nth_term(1000ULL, initial, recurrence) != 24806056) {
        std::cerr << "Checkpoint failed: f(10^3)\n";
        return false;
    }
    if (nth_term(1000000ULL, initial, recurrence) != 30808124) {
        std::cerr << "Checkpoint failed: f(10^6)\n";
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    const auto transitions = build_transitions();
    const std::vector<int> sequence = generate_sequence(1300, transitions);
    const std::vector<int> recurrence = berlekamp_massey(sequence);

    if (recurrence.empty()) {
        std::cerr << "Failed to derive recurrence\n";
        return 2;
    }

    if (!validate_recurrence(sequence, recurrence)) {
        std::cerr << "Recurrence validation failed\n";
        return 3;
    }

    std::vector<int> initial(recurrence.size(), 0);
    for (std::size_t i = 0; i < recurrence.size(); ++i) {
        initial[i] = sequence[i];
    }

    if (!skip_checkpoints && !run_checkpoints(initial, recurrence)) {
        return 4;
    }

    std::string exponent = "1";
    exponent.append(10000, '0');  // 10^10000

    const int answer = nth_term_decimal(exponent, initial, recurrence);
    std::cout << answer << '\n';
    return 0;
}

Python

MOD = 100000007
CELLS = 9
FULL_MASK = (1 << CELLS) - 1
STATE_COUNT = 1 << CELLS

def add_mod(a, b):
    return (a + b) % MOD

def sub_mod(a, b):
    return (a - b) % MOD

def mul_mod(a, b):
    return (a * b) % MOD

def pow_mod(base, exp):
    return pow(base, exp, MOD)

def inv_mod(a):
    return pow_mod(a, MOD - 2)

def dfs_layer_fill(filled_mask, next_mask, ways):
    if filled_mask == FULL_MASK:
        ways[next_mask] += 1
        return
        
    free_bits = (~filled_mask) & FULL_MASK
    bit = (free_bits & -free_bits).bit_length() - 1
    
    dfs_layer_fill(filled_mask | (1 << bit), next_mask | (1 << bit), ways)
    
    r = bit // 3
    c = bit % 3
    
    if c < 2 and (filled_mask & (1 << (bit + 1))) == 0:
        dfs_layer_fill(filled_mask | (1 << bit) | (1 << (bit + 1)), next_mask, ways)
        
    if r < 2 and (filled_mask & (1 << (bit + 3))) == 0:
        dfs_layer_fill(filled_mask | (1 << bit) | (1 << (bit + 3)), next_mask, ways)

def build_transitions():
    transitions = [[] for _ in range(STATE_COUNT)]
    for mask in range(STATE_COUNT):
        ways = [0] * STATE_COUNT
        dfs_layer_fill(mask, 0, ways)
        for nxt in range(STATE_COUNT):
            if ways[nxt] > 0:
                transitions[mask].append((nxt, ways[nxt] % MOD))
    return transitions

def generate_sequence(terms, transitions):
    seq = [0] * (terms + 1)
    current = [0] * STATE_COUNT
    current[0] = 1
    seq[0] = 1
    
    for step in range(1, terms + 1):
        nxt = [0] * STATE_COUNT
        for mask in range(STATE_COUNT):
            if current[mask] == 0:
                continue
            val = current[mask]
            for nx, cnt in transitions[mask]:
                nxt[nx] = (nxt[nx] + val * cnt) % MOD
        current = nxt
        seq[step] = current[0]
        
    return seq

def berlekamp_massey(sequence):
    C = [1]
    B = [1]
    L = 0
    m = 1
    b = 1
    
    for n in range(len(sequence)):
        d = sequence[n]
        for i in range(1, L + 1):
            d = (d + C[i] * sequence[n - i]) % MOD
            
        if d == 0:
            m += 1
            continue
            
        T = list(C)
        coef = (d * inv_mod(b)) % MOD
        
        if len(C) < len(B) + m:
            C.extend([0] * (len(B) + m - len(C)))
            
        for i in range(len(B)):
            C[i + m] = (C[i + m] - coef * B[i]) % MOD
            
        if 2 * L <= n:
            L = n + 1 - L
            B = T
            b = d
            m = 1
        else:
            m += 1
            
    recurrence = [0] * L
    for i in range(1, L + 1):
        recurrence[i - 1] = (MOD - C[i]) % MOD
    return recurrence

def combine_polynomials(a, b, recurrence):
    k = len(recurrence)
    prod = [0] * (2 * k)
    
    for i in range(k):
        if a[i] == 0:
            continue
        for j in range(k):
            if b[j] == 0:
                continue
            prod[i + j] = (prod[i + j] + a[i] * b[j]) % MOD
            
    for i in range(2 * k - 2, k - 1, -1):
        coef = prod[i]
        if coef == 0:
            continue
        for j in range(1, k + 1):
            prod[i - j] = (prod[i - j] + coef * recurrence[j - 1]) % MOD
            
    return prod[:k]

def linear_recurrence_nth(initial, recurrence, bits_lsb_first):
    k = len(recurrence)
    result = [0] * k
    result[0] = 1
    
    x_poly = [0] * k
    if k == 1:
        x_poly[0] = recurrence[0]
    else:
        x_poly[1] = 1
        
    for bit in bits_lsb_first:
        if bit != 0:
            result = combine_polynomials(result, x_poly, recurrence)
        x_poly = combine_polynomials(x_poly, x_poly, recurrence)
        
    ans = 0
    for i in range(k):
        ans = (ans + result[i] * initial[i]) % MOD
    return ans

def bits_from_decimal(value_str):
    value_int = int(value_str)
    if value_int == 0:
        return [0]
    bits = []
    while value_int > 0:
        bits.append(value_int & 1)
        value_int >>= 1
    return bits

def solve():
    transitions = build_transitions()
    seq = generate_sequence(1300, transitions)
    recurrence = berlekamp_massey(seq)
    
    initial = seq[:len(recurrence)]
    
    exponent = "1" + "0" * 10000
    bits = bits_from_decimal(exponent)
    ans = linear_recurrence_nth(initial, recurrence, bits)
    return str(ans)

if __name__ == '__main__':
    import sys
    sys.set_int_max_str_digits(20000)
    print(solve())

Java

import java.util.*;
import java.math.BigInteger;

public class Euler324 {
    static final int MOD = 100000007;
    static final int CELLS = 9;
    static final int FULL_MASK = (1 << CELLS) - 1;
    static final int STATE_COUNT = 1 << CELLS;

    static int addMod(int a, int b) {
        int x = a + b;
        return x >= MOD ? x - MOD : x;
    }

    static int subMod(int a, int b) {
        int x = a - b;
        return x < 0 ? x + MOD : x;
    }

    static int mulMod(int a, int b) {
        return (int) (((long) a * b) % MOD);
    }

    static int powMod(int base, long exp) {
        int res = 1;
        while (exp > 0) {
            if ((exp & 1) != 0)
                res = mulMod(res, base);
            base = mulMod(base, base);
            exp >>= 1;
        }
        return res;
    }

    static int invMod(int a) {
        return powMod(a, MOD - 2);
    }

    static void dfsLayerFill(int filledMask, int nextMask, int[] ways) {
        if (filledMask == FULL_MASK) {
            ways[nextMask]++;
            return;
        }

        int freeBits = (~filledMask) & FULL_MASK;
        int bit = Integer.numberOfTrailingZeros(freeBits);

        dfsLayerFill(filledMask | (1 << bit), nextMask | (1 << bit), ways);

        int r = bit / 3;
        int c = bit % 3;

        if (c < 2 && (filledMask & (1 << (bit + 1))) == 0) {
            dfsLayerFill(filledMask | (1 << bit) | (1 << (bit + 1)), nextMask, ways);
        }

        if (r < 2 && (filledMask & (1 << (bit + 3))) == 0) {
            dfsLayerFill(filledMask | (1 << bit) | (1 << (bit + 3)), nextMask, ways);
        }
    }

    static class Edge {
        int nxt, cnt;

        Edge(int nxt, int cnt) {
            this.nxt = nxt;
            this.cnt = cnt;
        }
    }

    @SuppressWarnings("unchecked")
    static List<Edge>[] buildTransitions() {
        List<Edge>[] transitions = new List[STATE_COUNT];
        for (int mask = 0; mask < STATE_COUNT; ++mask) {
            int[] ways = new int[STATE_COUNT];
            dfsLayerFill(mask, 0, ways);
            transitions[mask] = new ArrayList<>();
            for (int nxt = 0; nxt < STATE_COUNT; ++nxt) {
                if (ways[nxt] > 0) {
                    transitions[mask].add(new Edge(nxt, ways[nxt] % MOD));
                }
            }
        }
        return transitions;
    }

    static int[] generateSequence(int terms, List<Edge>[] transitions) {
        int[] seq = new int[terms + 1];
        int[] current = new int[STATE_COUNT];
        int[] next = new int[STATE_COUNT];

        current[0] = 1;
        seq[0] = 1;

        for (int step = 1; step <= terms; ++step) {
            Arrays.fill(next, 0);
            for (int mask = 0; mask < STATE_COUNT; ++mask) {
                int value = current[mask];
                if (value == 0)
                    continue;
                for (Edge e : transitions[mask]) {
                    next[e.nxt] = addMod(next[e.nxt], mulMod(value, e.cnt));
                }
            }
            System.arraycopy(next, 0, current, 0, STATE_COUNT);
            seq[step] = current[0];
        }
        return seq;
    }

    static int[] berlekampMassey(int[] sequence) {
        List<Integer> C = new ArrayList<>();
        List<Integer> B = new ArrayList<>();
        C.add(1);
        B.add(1);

        int L = 0, m = 1, b = 1;

        for (int n = 0; n < sequence.length; ++n) {
            int d = sequence[n];
            for (int i = 1; i <= L; ++i) {
                if (i < C.size()) {
                    d = addMod(d, mulMod(C.get(i), sequence[n - i]));
                }
            }

            if (d == 0) {
                m++;
                continue;
            }

            List<Integer> T = new ArrayList<>(C);
            int coef = mulMod(d, invMod(b));

            while (C.size() < B.size() + m) {
                C.add(0);
            }

            for (int i = 0; i < B.size(); ++i) {
                C.set(i + m, subMod(C.get(i + m), mulMod(coef, B.get(i))));
            }

            if (2 * L <= n) {
                L = n + 1 - L;
                B = T;
                b = d;
                m = 1;
            } else {
                m++;
            }
        }

        int[] recurrence = new int[L];
        for (int i = 1; i <= L; ++i) {
            recurrence[i - 1] = (MOD - C.get(i)) % MOD;
        }
        return recurrence;
    }

    static int[] combinePolynomials(int[] a, int[] b, int[] recurrence) {
        int k = recurrence.length;
        int[] prod = new int[2 * k];

        for (int i = 0; i < k; ++i) {
            if (a[i] == 0)
                continue;
            for (int j = 0; j < k; ++j) {
                if (b[j] == 0)
                    continue;
                prod[i + j] = addMod(prod[i + j], mulMod(a[i], b[j]));
            }
        }

        for (int i = 2 * k - 2; i >= k; --i) {
            int coef = prod[i];
            if (coef == 0)
                continue;
            for (int j = 1; j <= k; ++j) {
                prod[i - j] = addMod(prod[i - j], mulMod(coef, recurrence[j - 1]));
            }
        }

        int[] res = new int[k];
        System.arraycopy(prod, 0, res, 0, k);
        return res;
    }

    static int linearRecurrenceNth(int[] initial, int[] recurrence, List<Integer> bitsLsbFirst) {
        int k = recurrence.length;
        int[] result = new int[k];
        result[0] = 1;

        int[] xPoly = new int[k];
        if (k == 1)
            xPoly[0] = recurrence[0];
        else
            xPoly[1] = 1;

        for (int bit : bitsLsbFirst) {
            if (bit != 0) {
                result = combinePolynomials(result, xPoly, recurrence);
            }
            xPoly = combinePolynomials(xPoly, xPoly, recurrence);
        }

        int answer = 0;
        for (int i = 0; i < k; ++i) {
            answer = addMod(answer, mulMod(result[i], initial[i]));
        }
        return answer;
    }

    static List<Integer> bitsFromDecimal(String valueStr) {
        BigInteger val = new BigInteger(valueStr);
        List<Integer> bits = new ArrayList<>();
        if (val.equals(BigInteger.ZERO)) {
            bits.add(0);
            return bits;
        }
        while (val.compareTo(BigInteger.ZERO) > 0) {
            bits.add(val.testBit(0) ? 1 : 0);
            val = val.shiftRight(1);
        }
        return bits;
    }

    public static String solve() {
        List<Edge>[] transitions = buildTransitions();
        int[] seq = generateSequence(1300, transitions);
        int[] recurrence = berlekampMassey(seq);

        int[] initial = new int[recurrence.length];
        System.arraycopy(seq, 0, initial, 0, recurrence.length);

        StringBuilder sb = new StringBuilder();
        sb.append("1");
        for (int i = 0; i < 10000; i++)
            sb.append("0");
        String exponentStr = sb.toString();

        List<Integer> bits = bitsFromDecimal(exponentStr);
        int ans = linearRecurrenceNth(initial, recurrence, bits);
        return String.valueOf(ans);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}