Problem 324: Building a Tower
View on Project EulerProject 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
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. While doing that, a domino can be placed in exactly three ways:
1) entirely inside the current layer, horizontally along columns,
2) entirely inside the current layer, vertically along rows,
3) across the layer boundary, meaning one cube is in the current layer and the other cube will occupy the same cell in the next layer.
The third type creates a 1 in the outgoing mask \(t\).
4) DFS generates the transfer counts.
The code uses a depth-first search on the current layer. At each step it picks the first free cell and branches over the legal placements listed above. When the layer becomes fully covered, we have produced exactly one transition
$$s\to t.$$
Let \(T_{s,t}\) be the number of ways this can happen.
Because the layer is only \(3\times 3\), this DFS is finite and can be done once for all 512 incoming masks.
5) Transfer-matrix dynamic programming.
Let \(v_k[s]\) be the number of ways to build the first \(k\) layers and finish with outgoing mask \(s\). Then
$$v_{k+1}[t]=\sum_{s=0}^{511} v_k[s]\,T_{s,t}\pmod M.$$
The initial condition is
$$v_0[0]=1,\qquad v_0[s]=0\ \text{for }s\ne 0,$$
because before building anything there is no pending occupancy from a previous layer.
The quantity we actually want is the fully closed tower count
$$f(k)=v_k[0],$$
since at the top of a valid tower no domino may stick out into an imaginary next layer.
6) Why a linear recurrence must exist.
The vector \(v_k\) lives in a 512-dimensional vector space over \(\mathbb Z_M\), and it evolves by multiplication with the same fixed \(512\times 512\) matrix \(T\). Therefore the scalar sequence \(f(k)=v_k[0]\) is automatically a linear recurrence sequence. In principle its recurrence length is at most 512, although the minimal recurrence can be smaller.
7) Recover the minimal recurrence with Berlekamp-Massey.
The program generates the first 1301 terms
$$f(0),f(1),\dots,f(1300)$$
using the 512-state DP. Then Berlekamp-Massey reconstructs the shortest recurrence
$$f(n)=\sum_{j=1}^{L} c_j f(n-j)\pmod M.$$
This is the crucial compression step: after that, we never need the 512-state DP again.
8) Small checkpoints.
The C++ code validates the recovered recurrence against known terms:
$$f(2)=229,\qquad f(4)=117805,\qquad f(10)=96149360,$$
and also much later values
$$f(10^3)=24806056,\qquad f(10^6)=30808124.$$
These checkpoints are important: they verify both the transfer-DP part and the recurrence-evaluation part.
9) Why giant matrix exponentiation is replaced by Kitamasa-style polynomial arithmetic.
Once we know the recurrence, computing \(f(N)\) is equivalent to reducing the monomial \(x^N\) modulo the characteristic polynomial
$$P(x)=x^L-\sum_{j=1}^{L} c_j x^{L-j}.$$
If
$$x^N \equiv a_0+a_1x+\cdots+a_{L-1}x^{L-1}\pmod{P(x)},$$
then
$$f(N)=a_0f(0)+a_1f(1)+\cdots+a_{L-1}f(L-1)\pmod M.$$
The code computes these coefficients by repeated squaring and polynomial combination, which is the standard Kitamasa viewpoint.
10) Handling \(N=10^{10000}\).
The exponent does not fit into any machine integer type. So the code stores it as the decimal string
$$1000\ldots 0$$
with 10000 zeros, repeatedly divides that string by 2, and records the remainders. This produces the bits of \(N\) in least-significant-bit-first order, which is exactly what binary exponentiation needs.
Algorithm
1) Precompute all layer transitions \(T_{s,t}\) by DFS over a \(3\times 3\) slice.
2) Run the 512-state DP long enough to obtain a prefix of \(f(n)\).
3) Apply Berlekamp-Massey to extract the minimal linear recurrence.
4) Validate the recurrence on additional known terms.
5) Convert \(N=10^{10000}\) from decimal to binary bits.
6) Evaluate \(f(N)\) modulo \(M\) using Kitamasa-style polynomial squaring.
Complexity Analysis
The transition precomputation is finite and bounded by the 512 masks. The prefix DP is also small. If the minimal recurrence length is \(L\), the final evaluation costs roughly
$$O(L^2\log N)$$
modular operations and
$$O(L)$$
memory. That is what makes the astronomical index \(10^{10000}\) manageable.
Checks And Final Result
The code checks
$$f(2)=229,\quad f(4)=117805,\quad f(10)=96149360,\quad f(10^3)=24806056,\quad f(10^6)=30808124.$$
For the required exponent
$$N=10^{10000},$$
the final answer is
$$f(N)=96972774.$$
Further Reading
- Problem page: https://projecteuler.net/problem=324
- Berlekamp-Massey algorithm: https://en.wikipedia.org/wiki/Berlekamp-Massey_algorithm
- Linear recurrences / Kitamasa: https://cp-algorithms.com/algebra/linear-recurrence.html
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());
}
}