Problem 490: Jumping Frog
View on Project EulerProject Euler Problem 490 Solution
EulerSolve provides an optimized solution for Project Euler Problem 490, Jumping Frog, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each \(n\), let \(f(n)\) be the number of routes in which a frog starts on stone \(1\), ends on stone \(n\), visits every stone \(1,2,\dots,n\) exactly once, and every jump has length at most \(3\). Equivalently, \(f(n)\) counts Hamiltonian paths from \(1\) to \(n\) in the graph $$G_n=(\{1,\dots,n\},E_n),\qquad \{i,j\}\in E_n \iff 1\le |i-j|\le 3.$$ The quantity required by the problem is $$S(L)=\sum_{n=1}^{L} f(n)^3 \pmod{10^9}.$$ Since the main input has \(L=10^{14}\), it is useless to compute the values one by one. The solution instead converts the route count into a finite-state transfer process and then sums cubes by binary lifting in the third tensor power. Mathematical Approach The crucial observation is local: when the stones are processed from left to right, a future stone can only connect to the previous three stones. That means the entire history collapses to a small frontier state. Step 1: Rewrite the Problem as a Path in a Banded Graph A legal frog route is a permutation $$a_1,a_2,\dots,a_n$$ of \(\{1,\dots,n\}\) such that $$a_1=1,\qquad a_n=n,\qquad |a_{t+1}-a_t|\le 3 \text{ for every } t.$$ Viewed as an undirected graph, consecutive jumps form a path that uses every vertex exactly once. So we are counting Hamiltonian paths in \(G_n\) whose endpoints are \(1\) and \(n\)....
Detailed mathematical approach
Problem Summary
For each \(n\), let \(f(n)\) be the number of routes in which a frog starts on stone \(1\), ends on stone \(n\), visits every stone \(1,2,\dots,n\) exactly once, and every jump has length at most \(3\).
Equivalently, \(f(n)\) counts Hamiltonian paths from \(1\) to \(n\) in the graph
$$G_n=(\{1,\dots,n\},E_n),\qquad \{i,j\}\in E_n \iff 1\le |i-j|\le 3.$$
The quantity required by the problem is
$$S(L)=\sum_{n=1}^{L} f(n)^3 \pmod{10^9}.$$
Since the main input has \(L=10^{14}\), it is useless to compute the values one by one. The solution instead converts the route count into a finite-state transfer process and then sums cubes by binary lifting in the third tensor power.
Mathematical Approach
The crucial observation is local: when the stones are processed from left to right, a future stone can only connect to the previous three stones. That means the entire history collapses to a small frontier state.
Step 1: Rewrite the Problem as a Path in a Banded Graph
A legal frog route is a permutation
$$a_1,a_2,\dots,a_n$$
of \(\{1,\dots,n\}\) such that
$$a_1=1,\qquad a_n=n,\qquad |a_{t+1}-a_t|\le 3 \text{ for every } t.$$
Viewed as an undirected graph, consecutive jumps form a path that uses every vertex exactly once. So we are counting Hamiltonian paths in \(G_n\) whose endpoints are \(1\) and \(n\).
Because edges only join vertices whose indices differ by at most \(3\), once we have passed far enough to the right, an old vertex can no longer receive any new edge. This is what makes a transfer-matrix approach possible.
Step 2: Describe a Partial Construction by a Frontier State
Suppose we have already processed stones \(1,\dots,t\). Only the stones
$$\max(1,t-2),\max(1,t-1),t$$
can still connect to a future stone, so only those active stones need to be remembered.
For each active stone we record two pieces of information:
$$\text{its current degree } \in \{0,1,2\},\qquad \text{its connected-component label.}$$
Degree \(2\) is the maximum possible degree in a path. Component labels tell us whether two active stones are already linked through the partial route built so far.
Different label names do not matter; only the partition matters. Therefore component labels are canonically renumbered in first-occurrence order. After this normalization, the set of reachable states is finite.
Step 3: Legal Transitions When a New Stone Is Added
When stone \(t+1\) is introduced, it may connect to zero, one, or two currently active stones.
The local restrictions are exactly the restrictions for extending a path:
$$\deg(v)\le 2 \text{ for every active stone } v,$$
and if the new stone connects to two active stones, those two stones must lie in different components; otherwise the new edges would close a cycle before the route is complete.
Once \(t+1\ge 4\), the oldest active stone leaves the frontier forever, because later stones are too far away to reach it. At that moment its degree must already be final:
$$\deg(1)=1 \text{ when stone } 1 \text{ leaves},$$
$$\deg(v)=2 \text{ for every later stone } v \text{ when it leaves}.$$
The first condition fixes stone \(1\) as the starting endpoint. The second condition forces every interior stone to be an interior vertex of the Hamiltonian path.
There is one more connectivity condition. If removing the oldest active stone would make its entire connected component disappear from the frontier, then that component is now sealed off from every future stone and can never merge with the rest of the route. Such a transition is impossible and must be discarded.
Step 4: The Transfer Matrix After the Warm-Up Phase
The first few lengths are exceptional because the frontier has not yet stabilized. Direct inspection gives
$$f(1)=f(2)=f(3)=1.$$
After the first four stones have been processed, the frontier always has size \(3\). Let the reachable canonical frontier states in this stable regime be
$$q_1,q_2,\dots,q_d.$$
Let \(v_n\in (\mathbb Z/10^9\mathbb Z)^d\) be the column vector whose \(i\)-th entry counts partial constructions of size \(n\) currently in state \(q_i\). Then for every \(n\ge 4\),
$$v_{n+1}=A v_n,$$
where \(A_{ij}\) is the number of legal local extensions from state \(q_j\) to state \(q_i\).
Now define a selector vector \(c\) by
$$c_i=\begin{cases}1,&\text{if } q_i \text{ represents one connected path with endpoint } n,\\ 0,&\text{otherwise.}\end{cases}$$
Then
$$f(n)=c^{\mathsf T} v_n \qquad (n\ge 4).$$
Step 5: Why the Cube Sum Lives in the Third Tensor Power
For any state vector \(x\),
$$\bigl(c^{\mathsf T}x\bigr)^3=\langle c^{\otimes 3},x^{\otimes 3}\rangle,$$
where \(c^{\otimes 3}=c\otimes c\otimes c\) and \(\langle \cdot,\cdot\rangle\) is the natural dot product on the third tensor power.
Therefore block sums of cubes can be encoded by third-order tensors. For \(r\ge 0\), define
$$\Phi_r(x)=\sum_{m=0}^{2^r-1} \bigl(c^{\mathsf T}A^m x\bigr)^3.$$
There exists a tensor \(T_r\) such that
$$\Phi_r(x)=\langle T_r,x^{\otimes 3}\rangle.$$
The base case is immediate:
$$T_0=c^{\otimes 3}.$$
Step 6: Doubling Formula for Block Sums
Let
$$A_r=A^{2^r}.$$
A block of length \(2^{r+1}\) splits into two consecutive blocks of length \(2^r\), so
$$\Phi_{r+1}(x)=\Phi_r(x)+\Phi_r(A_r x).$$
Hence the tensors satisfy
$$T_{r+1}=T_r+\mathcal T_{A_r}(T_r),$$
where \(\mathcal T_M\) means: apply the matrix \(M\) in each of the three tensor slots. By definition,
$$\langle \mathcal T_M(T),x^{\otimes 3}\rangle=\langle T,(Mx)^{\otimes 3}\rangle.$$
This is the tensor analogue of the usual doubling identity for prefix sums.
Worked Example
For \(n=4\), every pair of stones is at distance at most \(3\), so \(G_4\) is the complete graph \(K_4\).
To count routes from \(1\) to \(4\), we only need to place stones \(2\) and \(3\) in the middle while preserving a Hamiltonian path. The two valid routes are
$$1,2,3,4 \qquad \text{and} \qquad 1,3,2,4,$$
so
$$f(4)=2.$$
The implementations also verify the checkpoints
$$f(6)=14,\qquad f(10)=254,$$
$$S(10)=18230635,$$
and
$$S(1000)\equiv 225031475 \pmod{10^9}.$$
These values are strong consistency checks for both the transfer matrix and the tensor-doubling sum.
How the Code Works
The C++, Python, and Java implementations all follow the same structure. First they expand the first few stones explicitly, because the frontier size is still growing and the left endpoint has special degree requirements. After this warm-up they enumerate every reachable canonical frontier state and build the stationary transition matrix.
Next they assemble the state vector for length \(4\) and the selector that recognizes a completed route. They then precompute matrix powers \(A^{2^r}\) together with the block tensors \(T_r\) from the doubling formula above.
Finally, they read the binary expansion of \(L-3\). Whenever bit \(r\) is set, they add the contribution of the whole block represented by \(T_r\) for the current state vector and then advance that vector by \(A^{2^r}\). The three trivial initial terms \(f(1)^3,f(2)^3,f(3)^3\) are added separately. All arithmetic is reduced modulo \(10^9\).
Complexity Analysis
Let \(d\) be the number of reachable canonical frontier states. Once the local rule is fixed, \(d\) is a small constant independent of \(L\), so the main cost is logarithmic in \(L\).
Building matrix powers costs \(O(d^3\log L)\). Each tensor transform applies the matrix in three slots, which costs \(O(d^4)\) per level, so all tensor block precomputations cost \(O(d^4\log L)\). Evaluating the selected block tensors during the binary walk costs \(O(d^3\log L)\).
The memory usage is \(O(d^3\log L)\) for the stored third-order tensors and \(O(d^2\log L)\) for the matrix powers.
Footnotes and References
- Problem page: https://projecteuler.net/problem=490
- Hamiltonian path: Wikipedia — Hamiltonian path
- Transfer-matrix method: Wikipedia — Transfer-matrix method
- Kronecker product and tensor powers: Wikipedia — Kronecker product
- Exponentiation by squaring: Wikipedia — Exponentiation by squaring
Problem 490 source code
C++
#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <queue>
#include <string>
#include <unordered_map>
#include <vector>
namespace {
constexpr uint64_t MOD = 1000000000ULL;
struct State {
int k = 0;
std::array<int, 4> deg{0, 0, 0, 0};
std::array<int, 4> comp{0, 0, 0, 0};
};
struct StateHash {
size_t operator()(const State& s) const noexcept {
size_t h = static_cast<size_t>(s.k);
for (int i = 0; i < s.k; ++i) {
h = h * 3 + static_cast<size_t>(s.deg[i]);
}
for (int i = 0; i < s.k; ++i) {
h = h * 3 + static_cast<size_t>(s.comp[i]);
}
return h;
}
};
struct StateEq {
bool operator()(const State& a, const State& b) const noexcept {
if (a.k != b.k) return false;
for (int i = 0; i < a.k; ++i) {
if (a.deg[i] != b.deg[i]) return false;
}
for (int i = 0; i < a.k; ++i) {
if (a.comp[i] != b.comp[i]) return false;
}
return true;
}
};
State canonicalize(State s) {
int map[5];
for (int& v : map) v = -1;
int next = 0;
for (int i = 0; i < s.k; ++i) {
int c = s.comp[i];
if (map[c] == -1) map[c] = next++;
s.comp[i] = map[c];
}
return s;
}
bool comps_all_same(const State& s) {
for (int i = 1; i < s.k; ++i) {
if (s.comp[i] != s.comp[0]) return false;
}
return true;
}
std::vector<std::pair<State, uint64_t>> transitions(const State& s, bool remove, bool remove_is_one) {
std::vector<std::pair<State, uint64_t>> out;
int k = s.k;
std::vector<std::vector<int>> subsets;
subsets.push_back({});
for (int i = 0; i < k; ++i) subsets.push_back({i});
for (int i = 0; i < k; ++i) {
for (int j = i + 1; j < k; ++j) {
subsets.push_back({i, j});
}
}
auto add_state = [&](const State& ns) {
for (auto& entry : out) {
if (StateEq{}(entry.first, ns)) {
entry.second += 1;
return;
}
}
out.push_back({ns, 1});
};
for (const auto& subset : subsets) {
if (subset.size() == 2 && s.comp[subset[0]] == s.comp[subset[1]]) continue;
State ns;
ns.k = k + 1;
ns.deg = s.deg;
ns.comp = s.comp;
bool ok = true;
for (int idx : subset) {
if (ns.deg[idx] >= 2) {
ok = false;
break;
}
}
if (!ok) continue;
for (int idx : subset) ns.deg[idx] += 1;
int new_deg = static_cast<int>(subset.size());
int new_comp = 0;
if (subset.empty()) {
int maxc = -1;
for (int i = 0; i < k; ++i) maxc = std::max(maxc, ns.comp[i]);
new_comp = maxc + 1;
} else if (subset.size() == 1) {
new_comp = ns.comp[subset[0]];
} else {
int comp_a = ns.comp[subset[0]];
int comp_b = ns.comp[subset[1]];
new_comp = std::min(comp_a, comp_b);
for (int i = 0; i < k; ++i) {
if (ns.comp[i] == comp_a || ns.comp[i] == comp_b) ns.comp[i] = new_comp;
}
}
ns.deg[k] = new_deg;
ns.comp[k] = new_comp;
if (remove) {
int leaving_deg = ns.deg[0];
if (remove_is_one) {
if (leaving_deg != 1) continue;
} else {
if (leaving_deg != 2) continue;
}
int leaving_comp = ns.comp[0];
for (int i = 1; i <= k; ++i) {
ns.deg[i - 1] = ns.deg[i];
ns.comp[i - 1] = ns.comp[i];
}
ns.k = k;
bool present = false;
for (int i = 0; i < ns.k; ++i) {
if (ns.comp[i] == leaving_comp) {
present = true;
break;
}
}
if (!present) continue;
}
ns = canonicalize(ns);
add_state(ns);
}
return out;
}
bool is_final_small(const State& s, int n) {
if (!comps_all_same(s)) return false;
if (n == 1) {
return s.k == 1 && s.deg[0] == 0;
}
if (n == 2) {
return s.k == 2 && s.deg[0] == 1 && s.deg[1] == 1;
}
if (n == 3) {
return s.k == 3 && s.deg[0] == 1 && s.deg[1] == 2 && s.deg[2] == 1;
}
return s.k == 3 && s.deg[0] == 2 && s.deg[1] == 2 && s.deg[2] == 1;
}
bool is_final_generic(const State& s) {
return s.k == 3 && comps_all_same(s) && s.deg[0] == 2 && s.deg[1] == 2 && s.deg[2] == 1;
}
std::vector<uint64_t> compute_small_f(int n_max) {
std::unordered_map<State, uint64_t, StateHash, StateEq> dp;
State start;
start.k = 1;
start.deg = {0, 0, 0, 0};
start.comp = {0, 0, 0, 0};
dp[start] = 1;
std::vector<uint64_t> f(n_max + 1, 0);
for (int t = 1; t <= n_max; ++t) {
uint64_t total = 0;
for (const auto& kv : dp) {
if (is_final_small(kv.first, t)) total += kv.second;
}
f[t] = total;
if (t == n_max) break;
bool remove = (t + 1 >= 4);
bool remove_is_one = remove && (t - 2 == 1);
std::unordered_map<State, uint64_t, StateHash, StateEq> next;
for (const auto& kv : dp) {
for (const auto& nxt : transitions(kv.first, remove, remove_is_one)) {
next[nxt.first] += kv.second * nxt.second;
}
}
dp.swap(next);
}
return f;
}
struct Model {
int d = 0;
std::vector<State> states;
std::unordered_map<State, int, StateHash, StateEq> index;
std::vector<std::vector<uint64_t>> A;
std::vector<uint64_t> v4;
std::vector<uint64_t> c;
};
Model build_model() {
std::unordered_map<State, uint64_t, StateHash, StateEq> dp;
State start;
start.k = 1;
start.deg = {0, 0, 0, 0};
start.comp = {0, 0, 0, 0};
dp[start] = 1;
for (int t = 1; t <= 3; ++t) {
bool remove = (t + 1 >= 4);
bool remove_is_one = remove && (t - 2 == 1);
std::unordered_map<State, uint64_t, StateHash, StateEq> next;
for (const auto& kv : dp) {
for (const auto& nxt : transitions(kv.first, remove, remove_is_one)) {
uint64_t add = (kv.second * nxt.second) % MOD;
uint64_t& slot = next[nxt.first];
slot += add;
if (slot >= MOD) slot %= MOD;
}
}
dp.swap(next);
}
Model model;
std::queue<State> q;
for (const auto& kv : dp) {
if (model.index.emplace(kv.first, model.d).second) {
model.states.push_back(kv.first);
q.push(kv.first);
++model.d;
}
}
while (!q.empty()) {
State cur = q.front();
q.pop();
for (const auto& nxt : transitions(cur, true, false)) {
if (model.index.emplace(nxt.first, model.d).second) {
model.states.push_back(nxt.first);
q.push(nxt.first);
++model.d;
}
}
}
model.A.assign(model.d, std::vector<uint64_t>(model.d, 0));
for (int j = 0; j < model.d; ++j) {
const State& cur = model.states[j];
for (const auto& nxt : transitions(cur, true, false)) {
int i = model.index[nxt.first];
model.A[i][j] = (model.A[i][j] + nxt.second) % MOD;
}
}
model.v4.assign(model.d, 0);
for (const auto& kv : dp) {
int id = model.index[kv.first];
model.v4[id] = (model.v4[id] + kv.second) % MOD;
}
model.c.assign(model.d, 0);
for (int i = 0; i < model.d; ++i) {
if (is_final_generic(model.states[i])) model.c[i] = 1;
}
return model;
}
using Matrix = std::vector<std::vector<uint64_t>>;
Matrix mul_mat(const Matrix& A, const Matrix& B) {
int d = static_cast<int>(A.size());
Matrix C(d, std::vector<uint64_t>(d, 0));
for (int i = 0; i < d; ++i) {
for (int k = 0; k < d; ++k) {
if (A[i][k] == 0) continue;
for (int j = 0; j < d; ++j) {
__int128 term = static_cast<__int128>(A[i][k]) * B[k][j];
C[i][j] = static_cast<uint64_t>((C[i][j] + term) % MOD);
}
}
}
return C;
}
std::vector<uint64_t> mul_mat_vec(const Matrix& A, const std::vector<uint64_t>& v) {
int d = static_cast<int>(A.size());
std::vector<uint64_t> out(d, 0);
for (int i = 0; i < d; ++i) {
__int128 sum = 0;
for (int j = 0; j < d; ++j) {
sum += static_cast<__int128>(A[i][j]) * v[j];
}
out[i] = static_cast<uint64_t>(sum % MOD);
}
return out;
}
int idx3(int i, int j, int k, int d) {
return (i * d + j) * d + k;
}
std::vector<uint64_t> transform_tensor(const std::vector<uint64_t>& T, const Matrix& A) {
int d = static_cast<int>(A.size());
int size = d * d * d;
std::vector<uint64_t> U(size, 0), V(size, 0), W(size, 0);
for (int a = 0; a < d; ++a) {
for (int j = 0; j < d; ++j) {
for (int k = 0; k < d; ++k) {
__int128 sum = 0;
for (int i = 0; i < d; ++i) {
sum += static_cast<__int128>(A[i][a]) * T[idx3(i, j, k, d)];
}
U[idx3(a, j, k, d)] = static_cast<uint64_t>(sum % MOD);
}
}
}
for (int a = 0; a < d; ++a) {
for (int b = 0; b < d; ++b) {
for (int k = 0; k < d; ++k) {
__int128 sum = 0;
for (int j = 0; j < d; ++j) {
sum += static_cast<__int128>(A[j][b]) * U[idx3(a, j, k, d)];
}
V[idx3(a, b, k, d)] = static_cast<uint64_t>(sum % MOD);
}
}
}
for (int a = 0; a < d; ++a) {
for (int b = 0; b < d; ++b) {
for (int c = 0; c < d; ++c) {
__int128 sum = 0;
for (int k = 0; k < d; ++k) {
sum += static_cast<__int128>(A[k][c]) * V[idx3(a, b, k, d)];
}
W[idx3(a, b, c, d)] = static_cast<uint64_t>(sum % MOD);
}
}
}
return W;
}
uint64_t eval_tensor(const std::vector<uint64_t>& T, const std::vector<uint64_t>& v) {
int d = static_cast<int>(v.size());
uint64_t total = 0;
for (int i = 0; i < d; ++i) {
if (v[i] == 0) continue;
for (int j = 0; j < d; ++j) {
if (v[j] == 0) continue;
uint64_t vij = (v[i] * v[j]) % MOD;
for (int k = 0; k < d; ++k) {
if (v[k] == 0) continue;
uint64_t t = T[idx3(i, j, k, d)];
if (t == 0) continue;
__int128 term = static_cast<__int128>(t) * vij % MOD;
term = term * v[k] % MOD;
total += static_cast<uint64_t>(term);
if (total >= MOD) total %= MOD;
}
}
}
return total % MOD;
}
struct Powers {
std::vector<Matrix> A_pow;
std::vector<std::vector<uint64_t>> T_pow;
};
Powers build_powers(const Matrix& A, const std::vector<uint64_t>& c, int max_bits) {
int d = static_cast<int>(A.size());
Powers powers;
powers.A_pow.resize(max_bits);
powers.T_pow.resize(max_bits);
powers.A_pow[0] = A;
std::vector<uint64_t> T1(d * d * d, 0);
for (int i = 0; i < d; ++i) {
if (c[i] == 0) continue;
for (int j = 0; j < d; ++j) {
if (c[j] == 0) continue;
uint64_t cij = (c[i] * c[j]) % MOD;
for (int k = 0; k < d; ++k) {
if (c[k] == 0) continue;
T1[idx3(i, j, k, d)] = (cij * c[k]) % MOD;
}
}
}
powers.T_pow[0] = std::move(T1);
for (int bit = 1; bit < max_bits; ++bit) {
powers.A_pow[bit] = mul_mat(powers.A_pow[bit - 1], powers.A_pow[bit - 1]);
std::vector<uint64_t> transformed = transform_tensor(powers.T_pow[bit - 1], powers.A_pow[bit - 1]);
std::vector<uint64_t> merged(transformed.size(), 0);
for (size_t i = 0; i < transformed.size(); ++i) {
uint64_t val = powers.T_pow[bit - 1][i] + transformed[i];
if (val >= MOD) val %= MOD;
merged[i] = val;
}
powers.T_pow[bit] = std::move(merged);
}
return powers;
}
uint64_t sum_from_4(uint64_t L, const Model& model, const Powers& powers) {
if (L < 4) return 0;
uint64_t len = L - 3;
std::vector<uint64_t> v = model.v4;
uint64_t sum = 0;
int bit = 0;
while (len > 0) {
if (len & 1ULL) {
sum = (sum + eval_tensor(powers.T_pow[bit], v)) % MOD;
v = mul_mat_vec(powers.A_pow[bit], v);
}
len >>= 1ULL;
++bit;
}
return sum;
}
uint64_t compute_S_mod(uint64_t L, const Model& model, const Powers& powers) {
uint64_t total = 0;
if (L >= 1) total = (total + 1) % MOD;
if (L >= 2) total = (total + 1) % MOD;
if (L >= 3) total = (total + 1) % MOD;
if (L >= 4) total = (total + sum_from_4(L, model, powers)) % MOD;
return total;
}
bool validate(const Model& model, const Powers& powers) {
auto f = compute_small_f(40);
if (f[6] != 14) {
std::cerr << "Validation failed: f(6) = " << f[6] << "\n";
return false;
}
if (f[10] != 254) {
std::cerr << "Validation failed: f(10) = " << f[10] << "\n";
return false;
}
if (f[40] != 1439682432976ULL) {
std::cerr << "Validation failed: f(40) = " << f[40] << "\n";
return false;
}
__int128 sum10 = 0;
for (int i = 1; i <= 10; ++i) {
__int128 v = f[i];
sum10 += v * v * v;
}
if (static_cast<uint64_t>(sum10) != 18230635ULL) {
std::cerr << "Validation failed: S(10) = " << static_cast<uint64_t>(sum10) << "\n";
return false;
}
__int128 sum20 = 0;
for (int i = 1; i <= 20; ++i) {
__int128 v = f[i];
sum20 += v * v * v;
}
if (static_cast<uint64_t>(sum20) != 104207881192114219ULL) {
std::cerr << "Validation failed: S(20) = " << static_cast<uint64_t>(sum20) << "\n";
return false;
}
uint64_t s1000 = compute_S_mod(1000, model, powers);
if (s1000 != 225031475ULL) {
std::cerr << "Validation failed: S(1000) mod 1e9 = " << s1000 << "\n";
return false;
}
uint64_t s1e6 = compute_S_mod(1000000ULL, model, powers);
if (s1e6 != 363486179ULL) {
std::cerr << "Validation failed: S(1e6) mod 1e9 = " << s1e6 << "\n";
return false;
}
return true;
}
} // namespace
int main(int argc, char** argv) {
uint64_t L = 100000000000000ULL;
if (argc > 1) {
L = std::stoull(argv[1]);
}
Model model = build_model();
uint64_t max_L = L < 1000000ULL ? 1000000ULL : L;
uint64_t max_len = (max_L >= 4) ? (max_L - 3) : 1;
int max_bits = 0;
while (max_bits < 63 && (1ULL << max_bits) <= max_len) ++max_bits;
if (max_bits == 0) max_bits = 1;
Powers powers = build_powers(model.A, model.c, max_bits);
if (!validate(model, powers)) return 1;
uint64_t answer = compute_S_mod(L, model, powers);
std::cout << answer << '\n';
return 0;
}
Python
MOD = 1000000000
def canonicalize(k, deg, comp):
m = {}
nxt = 0
new_comp = list(comp)
for i in range(k):
c = comp[i]
if c not in m:
m[c] = nxt
nxt += 1
new_comp[i] = m[c]
return tuple(new_comp)
def comps_all_same(k, comp):
for i in range(1, k):
if comp[i] != comp[0]: return False
return True
def transitions(k, deg, comp, remove, remove_is_one):
out = {}
subsets = [[]]
for i in range(k): subsets.append([i])
for i in range(k):
for j in range(i + 1, k):
subsets.append([i, j])
for subset in subsets:
if len(subset) == 2 and comp[subset[0]] == comp[subset[1]]:
continue
ok = True
for idx in subset:
if deg[idx] >= 2:
ok = False
break
if not ok: continue
ndeg = list(deg) + [len(subset)]
ncomp = list(comp) + [0]
nk = k + 1
for idx in subset:
ndeg[idx] += 1
if len(subset) == 0:
maxc = -1
for i in range(k):
maxc = max(maxc, ncomp[i])
new_comp = maxc + 1
elif len(subset) == 1:
new_comp = ncomp[subset[0]]
else:
comp_a, comp_b = ncomp[subset[0]], ncomp[subset[1]]
new_comp = min(comp_a, comp_b)
for i in range(k):
if ncomp[i] in (comp_a, comp_b):
ncomp[i] = new_comp
ncomp[k] = new_comp
if remove:
leaving_deg = ndeg[0]
if remove_is_one:
if leaving_deg != 1: continue
else:
if leaving_deg != 2: continue
leaving_comp = ncomp[0]
ndeg = ndeg[1:]
ncomp = ncomp[1:]
nk = k
present = False
for i in range(nk):
if ncomp[i] == leaving_comp:
present = True
break
if not present: continue
ncomp = canonicalize(nk, ndeg, ncomp)
st = (nk, tuple(ndeg), ncomp)
out[st] = out.get(st, 0) + 1
return out
def is_final_generic(k, deg, comp):
return k == 3 and comps_all_same(k, comp) and deg[0] == 2 and deg[1] == 2 and deg[2] == 1
def build_model():
dp = {(1, (0,), (0,)): 1}
for t in range(1, 4):
remove = (t + 1 >= 4)
remove_is_one = remove and (t - 2 == 1)
nxt = {}
for st, count in dp.items():
for nst, ncount in transitions(st[0], st[1], st[2], remove, remove_is_one).items():
nxt[nst] = (nxt.get(nst, 0) + count * ncount) % MOD
dp = nxt
states = []
idx_map = {}
d = 0
q = []
for st in dp:
if st not in idx_map:
idx_map[st] = d
states.append(st)
q.append(st)
d += 1
while q:
cur = q.pop(0)
for nst in transitions(cur[0], cur[1], cur[2], True, False):
if nst not in idx_map:
idx_map[nst] = d
states.append(nst)
q.append(nst)
d += 1
A = [[0]*d for _ in range(d)]
for j, cur in enumerate(states):
for nst, count in transitions(cur[0], cur[1], cur[2], True, False).items():
i = idx_map[nst]
A[i][j] = (A[i][j] + count) % MOD
v4 = [0]*d
for st, count in dp.items():
v4[idx_map[st]] = (v4[idx_map[st]] + count) % MOD
c = [0]*d
for i, st in enumerate(states):
if is_final_generic(st[0], st[1], st[2]):
c[i] = 1
return d, states, A, v4, c
def mul_mat(A, B, d):
C = [[0]*d for _ in range(d)]
for i in range(d):
for k in range(d):
if A[i][k] == 0: continue
for j in range(d):
C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD
return C
def mul_mat_vec(A, v, d):
out = [0]*d
for i in range(d):
out[i] = sum(A[i][j] * v[j] for j in range(d)) % MOD
return out
def idx3(i, j, k, d):
return (i * d + j) * d + k
def transform_tensor(T, A, d):
U = [0] * (d*d*d)
for a in range(d):
for j in range(d):
for k in range(d):
s = 0
for i in range(d):
s += A[i][a] * T[idx3(i, j, k, d)]
U[idx3(a, j, k, d)] = s % MOD
V = [0] * (d*d*d)
for a in range(d):
for b in range(d):
for k in range(d):
s = 0
for j in range(d):
s += A[j][b] * U[idx3(a, j, k, d)]
V[idx3(a, b, k, d)] = s % MOD
W = [0] * (d*d*d)
for a in range(d):
for b in range(d):
for c in range(d):
s = 0
for k in range(d):
s += A[k][c] * V[idx3(a, b, k, d)]
W[idx3(a, b, c, d)] = s % MOD
return W
def eval_tensor(T, v, d):
total = 0
for i in range(d):
if v[i] == 0: continue
for j in range(d):
if v[j] == 0: continue
vij = (v[i] * v[j]) % MOD
for k in range(d):
if v[k] == 0: continue
t = T[idx3(i, j, k, d)]
if t == 0: continue
term = (t * vij) % MOD
term = (term * v[k]) % MOD
total = (total + term) % MOD
return total
def build_powers(A, c, max_bits, d):
A_pow = [None] * max_bits
T_pow = [None] * max_bits
A_pow[0] = A
T1 = [0] * (d*d*d)
for i in range(d):
if c[i] == 0: continue
for j in range(d):
if c[j] == 0: continue
cij = (c[i] * c[j]) % MOD
for k in range(d):
if c[k] == 0: continue
T1[idx3(i, j, k, d)] = (cij * c[k]) % MOD
T_pow[0] = T1
for bit in range(1, max_bits):
A_pow[bit] = mul_mat(A_pow[bit - 1], A_pow[bit - 1], d)
transformed = transform_tensor(T_pow[bit - 1], A_pow[bit - 1], d)
merged = [0] * (d*d*d)
for i in range(d*d*d):
merged[i] = (T_pow[bit-1][i] + transformed[i]) % MOD
T_pow[bit] = merged
return A_pow, T_pow
def sum_from_4(L, d, v4, A_pow, T_pow):
if L < 4: return 0
length = L - 3
v = list(v4)
sum_val = 0
bit = 0
while length > 0:
if length & 1:
sum_val = (sum_val + eval_tensor(T_pow[bit], v, d)) % MOD
v = mul_mat_vec(A_pow[bit], v, d)
length >>= 1
bit += 1
return sum_val
def compute_S_mod(L, d, v4, A_pow, T_pow):
total = 0
if L >= 1: total = (total + 1) % MOD
if L >= 2: total = (total + 1) % MOD
if L >= 3: total = (total + 1) % MOD
if L >= 4: total = (total + sum_from_4(L, d, v4, A_pow, T_pow)) % MOD
return total
def solve():
L = 100000000000000
d, states, A, v4, c = build_model()
max_bits = 0
max_len = L - 3 if L >= 4 else 1
while max_bits < 63 and (1 << max_bits) <= max_len:
max_bits += 1
if max_bits == 0: max_bits = 1
A_pow, T_pow = build_powers(A, c, max_bits, d)
answer = compute_S_mod(L, d, v4, A_pow, T_pow)
return str(answer)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.LinkedList;
import java.util.List;
import java.util.Map;
import java.util.Queue;
public class Euler490 {
private static final long MOD = 1000000000L;
static class State {
int k;
int[] deg;
int[] comp;
State(int k, int[] deg, int[] comp) {
this.k = k;
this.deg = deg;
this.comp = comp;
}
@Override
public int hashCode() {
int h = k;
for (int i = 0; i < k; ++i)
h = h * 3 + deg[i];
for (int i = 0; i < k; ++i)
h = h * 3 + comp[i];
return h;
}
@Override
public boolean equals(Object obj) {
if (!(obj instanceof State))
return false;
State s = (State) obj;
if (k != s.k)
return false;
for (int i = 0; i < k; ++i)
if (deg[i] != s.deg[i])
return false;
for (int i = 0; i < k; ++i)
if (comp[i] != s.comp[i])
return false;
return true;
}
}
private static int[] canonicalize(int k, int[] comp) {
int[] map = new int[comp.length + 5];
Arrays.fill(map, -1);
int nxt = 0;
int[] newComp = Arrays.copyOf(comp, comp.length);
for (int i = 0; i < k; ++i) {
int c = comp[i];
if (map[c] == -1)
map[c] = nxt++;
newComp[i] = map[c];
}
return newComp;
}
private static boolean compsAllSame(int k, int[] comp) {
for (int i = 1; i < k; ++i) {
if (comp[i] != comp[0])
return false;
}
return true;
}
private static Map<State, Integer> transitions(State s, boolean remove, boolean removeIsOne) {
Map<State, Integer> out = new HashMap<>();
int k = s.k;
List<int[]> subsets = new ArrayList<>();
subsets.add(new int[0]);
for (int i = 0; i < k; ++i)
subsets.add(new int[] { i });
for (int i = 0; i < k; ++i) {
for (int j = i + 1; j < k; ++j) {
subsets.add(new int[] { i, j });
}
}
for (int[] subset : subsets) {
if (subset.length == 2 && s.comp[subset[0]] == s.comp[subset[1]])
continue;
boolean ok = true;
for (int idx : subset) {
if (s.deg[idx] >= 2) {
ok = false;
break;
}
}
if (!ok)
continue;
int[] ndeg = Arrays.copyOf(s.deg, k + 1);
int[] ncomp = Arrays.copyOf(s.comp, k + 1);
int nk = k + 1;
for (int idx : subset)
ndeg[idx]++;
int newComp = 0;
if (subset.length == 0) {
int maxc = -1;
for (int i = 0; i < k; ++i)
maxc = Math.max(maxc, ncomp[i]);
newComp = maxc + 1;
} else if (subset.length == 1) {
newComp = ncomp[subset[0]];
} else {
int compA = ncomp[subset[0]];
int compB = ncomp[subset[1]];
newComp = Math.min(compA, compB);
for (int i = 0; i < k; ++i) {
if (ncomp[i] == compA || ncomp[i] == compB)
ncomp[i] = newComp;
}
}
ndeg[k] = subset.length;
ncomp[k] = newComp;
if (remove) {
int leavingDeg = ndeg[0];
if (removeIsOne) {
if (leavingDeg != 1)
continue;
} else {
if (leavingDeg != 2)
continue;
}
int leavingComp = ncomp[0];
for (int i = 1; i <= k; ++i) {
ndeg[i - 1] = ndeg[i];
ncomp[i - 1] = ncomp[i];
}
ndeg[k] = 0;
ncomp[k] = 0;
nk = k;
boolean present = false;
for (int i = 0; i < nk; ++i) {
if (ncomp[i] == leavingComp) {
present = true;
break;
}
}
if (!present)
continue;
}
ncomp = canonicalize(nk, ncomp);
State nst = new State(nk, ndeg, ncomp);
out.put(nst, out.getOrDefault(nst, 0) + 1);
}
return out;
}
private static boolean isFinalGeneric(State s) {
return s.k == 3 && compsAllSame(s.k, s.comp) && s.deg[0] == 2 && s.deg[1] == 2 && s.deg[2] == 1;
}
static class Model {
int d;
List<State> states;
long[][] A;
long[] v4;
long[] c;
}
private static Model buildModel() {
Map<State, Long> dp = new HashMap<>();
State start = new State(1, new int[] { 0 }, new int[] { 0 });
dp.put(start, 1L);
for (int t = 1; t <= 3; ++t) {
boolean remove = (t + 1 >= 4);
boolean removeIsOne = remove && (t - 2 == 1);
Map<State, Long> next = new HashMap<>();
for (Map.Entry<State, Long> kv : dp.entrySet()) {
State st = kv.getKey();
long count = kv.getValue();
for (Map.Entry<State, Integer> tr : transitions(st, remove, removeIsOne).entrySet()) {
long add = (count * tr.getValue()) % MOD;
next.put(tr.getKey(), (next.getOrDefault(tr.getKey(), 0L) + add) % MOD);
}
}
dp = next;
}
Model model = new Model();
model.states = new ArrayList<>();
Map<State, Integer> idxMap = new HashMap<>();
model.d = 0;
Queue<State> q = new LinkedList<>();
for (State st : dp.keySet()) {
if (!idxMap.containsKey(st)) {
idxMap.put(st, model.d++);
model.states.add(st);
q.add(st);
}
}
while (!q.isEmpty()) {
State cur = q.poll();
for (Map.Entry<State, Integer> tr : transitions(cur, true, false).entrySet()) {
State nst = tr.getKey();
if (!idxMap.containsKey(nst)) {
idxMap.put(nst, model.d++);
model.states.add(nst);
q.add(nst);
}
}
}
int d = model.d;
model.A = new long[d][d];
for (int j = 0; j < d; ++j) {
State cur = model.states.get(j);
for (Map.Entry<State, Integer> tr : transitions(cur, true, false).entrySet()) {
int i = idxMap.get(tr.getKey());
model.A[i][j] = (model.A[i][j] + tr.getValue()) % MOD;
}
}
model.v4 = new long[d];
for (Map.Entry<State, Long> kv : dp.entrySet()) {
int id = idxMap.get(kv.getKey());
model.v4[id] = (model.v4[id] + kv.getValue()) % MOD;
}
model.c = new long[d];
for (int i = 0; i < d; ++i) {
if (isFinalGeneric(model.states.get(i)))
model.c[i] = 1;
}
return model;
}
private static long[][] mulMat(long[][] A, long[][] B) {
int d = A.length;
long[][] C = new long[d][d];
for (int i = 0; i < d; ++i) {
for (int k = 0; k < d; ++k) {
if (A[i][k] == 0)
continue;
for (int j = 0; j < d; ++j) {
long term = (A[i][k] * B[k][j]) % MOD;
C[i][j] = (C[i][j] + term) % MOD;
}
}
}
return C;
}
private static long[] mulMatVec(long[][] A, long[] v) {
int d = A.length;
long[] out = new long[d];
for (int i = 0; i < d; ++i) {
long sum = 0;
for (int j = 0; j < d; ++j) {
sum = (sum + A[i][j] * v[j]) % MOD;
}
out[i] = sum;
}
return out;
}
private static int idx3(int i, int j, int k, int d) {
return (i * d + j) * d + k;
}
private static long[] transformTensor(long[] T, long[][] A) {
int d = A.length;
long[] U = new long[d * d * d];
for (int a = 0; a < d; ++a) {
for (int j = 0; j < d; ++j) {
for (int k = 0; k < d; ++k) {
long sum = 0;
for (int i = 0; i < d; ++i) {
sum = (sum + A[i][a] * T[idx3(i, j, k, d)]) % MOD;
}
U[idx3(a, j, k, d)] = sum;
}
}
}
long[] V = new long[d * d * d];
for (int a = 0; a < d; ++a) {
for (int b = 0; b < d; ++b) {
for (int k = 0; k < d; ++k) {
long sum = 0;
for (int j = 0; j < d; ++j) {
sum = (sum + A[j][b] * U[idx3(a, j, k, d)]) % MOD;
}
V[idx3(a, b, k, d)] = sum;
}
}
}
long[] W = new long[d * d * d];
for (int a = 0; a < d; ++a) {
for (int b = 0; b < d; ++b) {
for (int c = 0; c < d; ++c) {
long sum = 0;
for (int k = 0; k < d; ++k) {
sum = (sum + A[k][c] * V[idx3(a, b, k, d)]) % MOD;
}
W[idx3(a, b, c, d)] = sum;
}
}
}
return W;
}
private static long evalTensor(long[] T, long[] v) {
int d = v.length;
long total = 0;
for (int i = 0; i < d; ++i) {
if (v[i] == 0)
continue;
for (int j = 0; j < d; ++j) {
if (v[j] == 0)
continue;
long vij = (v[i] * v[j]) % MOD;
for (int k = 0; k < d; ++k) {
if (v[k] == 0)
continue;
long t = T[idx3(i, j, k, d)];
if (t == 0)
continue;
long term = (t * vij) % MOD;
term = (term * v[k]) % MOD;
total = (total + term) % MOD;
}
}
}
return total;
}
static class Powers {
long[][][] APow;
long[][] TPow;
}
private static Powers buildPowers(long[][] A, long[] c, int maxBits) {
int d = A.length;
Powers p = new Powers();
p.APow = new long[maxBits][][];
p.TPow = new long[maxBits][];
p.APow[0] = A;
long[] T1 = new long[d * d * d];
for (int i = 0; i < d; ++i) {
if (c[i] == 0)
continue;
for (int j = 0; j < d; ++j) {
if (c[j] == 0)
continue;
long cij = (c[i] * c[j]) % MOD;
for (int k = 0; k < d; ++k) {
if (c[k] == 0)
continue;
T1[idx3(i, j, k, d)] = (cij * c[k]) % MOD;
}
}
}
p.TPow[0] = T1;
for (int bit = 1; bit < maxBits; ++bit) {
p.APow[bit] = mulMat(p.APow[bit - 1], p.APow[bit - 1]);
long[] transformed = transformTensor(p.TPow[bit - 1], p.APow[bit - 1]);
long[] merged = new long[d * d * d];
for (int i = 0; i < d * d * d; ++i) {
merged[i] = (p.TPow[bit - 1][i] + transformed[i]) % MOD;
}
p.TPow[bit] = merged;
}
return p;
}
private static long sumFrom4(long L, int d, long[] v4, Powers powers) {
if (L < 4)
return 0;
long len = L - 3;
long[] v = Arrays.copyOf(v4, v4.length);
long sum = 0;
int bit = 0;
while (len > 0) {
if ((len & 1) != 0) {
sum = (sum + evalTensor(powers.TPow[bit], v)) % MOD;
v = mulMatVec(powers.APow[bit], v);
}
len >>= 1;
bit++;
}
return sum;
}
private static long computeSMod(long L, Model model, Powers powers) {
long total = 0;
if (L >= 1)
total = (total + 1) % MOD;
if (L >= 2)
total = (total + 1) % MOD;
if (L >= 3)
total = (total + 1) % MOD;
if (L >= 4)
total = (total + sumFrom4(L, model.d, model.v4, powers)) % MOD;
return total;
}
public static void main(String[] args) {
long L = 100000000000000L;
Model model = buildModel();
long maxL = Math.max(L, 1000000L);
long maxLen = (maxL >= 4) ? (maxL - 3) : 1;
int maxBits = 0;
while (maxBits < 63 && (1L << maxBits) <= maxLen)
maxBits++;
if (maxBits == 0)
maxBits = 1;
Powers powers = buildPowers(model.A, model.c, maxBits);
long answer = computeSMod(L, model, powers);
System.out.println(answer);
}
}