Problem 497: Drunken Tower of Hanoi
View on Project EulerProject Euler Problem 497 Solution
EulerSolve provides an optimized solution for Project Euler Problem 497, Drunken Tower of Hanoi, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each \(n\), three rods are placed at positions \(3^n\), \(6^n\), and \(9^n\) on a circular track of length \(10^n\). The task is to perform the canonical three-rod Hanoi transfer from the first rod to the third rod, using the middle rod as auxiliary, and to measure the total weighted travel cost of the mover. The mover begins at the middle rod. A direct simulation would expand a Hanoi move list of exponential length, so the solution does not enumerate individual disk states. Instead, it counts how many times each directed rod-to-rod transition occurs and combines those counts with a closed-form cost for the corresponding geometric movement. The required output is $$\sum_{n=1}^{10000} E_n \pmod{10^9},$$ where \(E_n\) denotes the cost of the \(n\)-disk instance. Mathematical Approach The key observation is that the Hanoi recursion only needs a tiny fixed state: the orientation of the subproblem and the number of directed transitions between the three rods. Step 1: Encode each subproblem by its orientation Label the rods \(0,1,2\) in increasing position order. An orientation \((f,t,a)\) means: move the whole tower from source rod \(f\) to target rod \(t\), using rod \(a\) as auxiliary. Since \((f,t,a)\) is a permutation of \((0,1,2)\), there are exactly six orientations....
Detailed mathematical approach
Problem Summary
For each \(n\), three rods are placed at positions \(3^n\), \(6^n\), and \(9^n\) on a circular track of length \(10^n\). The task is to perform the canonical three-rod Hanoi transfer from the first rod to the third rod, using the middle rod as auxiliary, and to measure the total weighted travel cost of the mover.
The mover begins at the middle rod. A direct simulation would expand a Hanoi move list of exponential length, so the solution does not enumerate individual disk states. Instead, it counts how many times each directed rod-to-rod transition occurs and combines those counts with a closed-form cost for the corresponding geometric movement.
The required output is
$$\sum_{n=1}^{10000} E_n \pmod{10^9},$$
where \(E_n\) denotes the cost of the \(n\)-disk instance.
Mathematical Approach
The key observation is that the Hanoi recursion only needs a tiny fixed state: the orientation of the subproblem and the number of directed transitions between the three rods.
Step 1: Encode each subproblem by its orientation
Label the rods \(0,1,2\) in increasing position order. An orientation \((f,t,a)\) means: move the whole tower from source rod \(f\) to target rod \(t\), using rod \(a\) as auxiliary. Since \((f,t,a)\) is a permutation of \((0,1,2)\), there are exactly six orientations.
For each orientation and disk count \(n\), define a \(3\times 3\) matrix
$$T_n^{(f,t,a)}=\bigl(T_n^{(f,t,a)}(i,j)\bigr)_{0\le i,j\le 2},$$
where \(T_n^{(f,t,a)}(i,j)\) is the number of directed transitions \(i\to j\) that occur after the mover has already reached the source rod of that subproblem. This convention leaves the very first walk of the whole instance to a separate term.
Step 2: Derive the transition-count recurrence
Let \(U_{u,v}\) be the \(3\times 3\) matrix whose only nonzero entry is a \(1\) at row \(u\), column \(v\). For one disk, once the mover is already standing at the source rod, the task is just a single loaded move:
$$T_1^{(f,t,a)}=U_{f,t}.$$
For \(n\ge 2\), the standard Hanoi decomposition says:
1. move the top \(n-1\) disks from \(f\) to \(a\) using \(t\),
2. return from \(a\) to \(f\),
3. carry the largest disk from \(f\) to \(t\),
4. walk from \(t\) to \(a\),
5. move the \(n-1\) disks from \(a\) to \(t\) using \(f\).
Therefore
$$T_n^{(f,t,a)}=T_{n-1}^{(f,a,t)}+\Delta^{(f,t,a)}+T_{n-1}^{(a,t,f)},$$
with
$$\Delta^{(f,t,a)}=U_{a,f}+U_{f,t}+U_{t,a}.$$
This recurrence is exact, and it only updates six fixed-size matrices for each new \(n\).
Step 3: Derive the directed cost kernel on the circle
Now fix a circle size \(K\) and rod positions \(1\le x,y\le K\). The code uses a weighted travel rule that depends on direction.
If \(x<y\), the move goes forward through the arc from \(x\) to \(y\). The crossed segment weights are
$$2x-1,\ 2x+1,\ \dots,\ 2y-3,$$
so the cost is the arithmetic-series sum
$$h_K(x,y)=\sum_{r=x}^{y-1}(2r-1)=(y-x)(x+y-2).$$
If \(x>y\), the directed move wraps around the far side of the circle. The segment weights are then
$$2K-2y+1,\ 2K-2y-1,\ \dots,\ 2K-2x+3,$$
which gives
$$h_K(x,y)=\sum_{r=y}^{x-1}(2K-2r-1)=(x-y)(2K-x-y).$$
Combining both cases,
$$h_K(x,y)=\begin{cases} (y-x)(x+y-2), & x<y,\\ (x-y)(2K-x-y), & x>y,\\ 0, & x=y. \end{cases}$$
Step 4: Turn transition counts into the cost of one instance
Suppose the three rod positions are \(p_0<p_1<p_2\) on a circle of size \(K\), and the tower must move from rod \(0\) to rod \(2\) using rod \(1\). The mover starts at rod \(1\), so the initial walk is \(1\to 0\), with cost \(h_K(p_1,p_0)\).
After that initial placement, every directed transition counted by \(T_n^{(0,2,1)}(i,j)\) contributes the cost \(h_K(p_i,p_j)\). Hence
$$E(n,K,p_0,p_1,p_2)=h_K(p_1,p_0)+\sum_{i=0}^{2}\sum_{j=0}^{2}T_n^{(0,2,1)}(i,j)\,h_K(p_i,p_j).$$
This formula converts the abstract transition matrix into the exact weighted cost of the whole \(n\)-disk task.
Step 5: Specialize to Problem 497
For the problem sequence, the geometric parameters are
$$K_n=10^n,\qquad p_0=3^n,\qquad p_1=6^n,\qquad p_2=9^n.$$
Because
$$3^n<6^n<9^n<10^n\qquad (n\ge 1),$$
the relative rod order never changes, so the correct branch of \(h_K\) is determined once and for all by the rod indices. The total requested sum is therefore
$$S(N)=\sum_{n=1}^{N} E(n,10^n,3^n,6^n,9^n)\pmod{10^9}.$$
Since only the last nine digits are needed, every multiplication and addition may be reduced modulo \(10^9\) during the dynamic program.
Step 6: Worked Example \((n,K,p_0,p_1,p_2)=(2,5,1,3,5)\)
For \(n=2\), the canonical orientation is \((0,2,1)\). The recurrence gives
$$T_2^{(0,2,1)}=U_{0,1}+U_{1,0}+U_{0,2}+U_{2,1}+U_{1,2}.$$
So after the initial walk \(1\to 0\), the mover performs the directed transitions
$$0\to 1,\qquad 1\to 0,\qquad 0\to 2,\qquad 2\to 1,\qquad 1\to 2.$$
With \(K=5\) and positions \((1,3,5)\), the individual costs are
$$h_5(3,1)=12,\qquad h_5(1,3)=4,\qquad h_5(1,5)=16,\qquad h_5(5,3)=4,\qquad h_5(3,5)=12.$$
Therefore
$$E(2,5,1,3,5)=12+4+12+16+4+12=60,$$
which matches the small exact checkpoint used by the implementation.
How the Code Works
The C++, Python, and Java implementations first enumerate the six orientations and store one \(3\times 3\) transition-count matrix for each of them. At \(n=1\), every orientation contains exactly one source-to-target move. For each larger \(n\), every new matrix is produced from two previously stored matrices plus the three fixed bridge transitions \(a\to f\), \(f\to t\), and \(t\to a\).
In parallel, the implementations advance the four power sequences \(3^n\), \(6^n\), \(9^n\), and \(10^n\) modulo \(10^9\). After each update they evaluate the canonical orientation \((0,2,1)\), apply the cost kernel to the current positions, and add the result to a running sum modulo \(10^9\).
The C++ implementation also checks the derivation on two exact small instances before computing the long sum:
$$E(2,5,1,3,5)=60,\qquad E(3,20,4,9,17)=2358.$$
Those checkpoints confirm that both the transition recurrence and the geometric cost formula are wired correctly.
Complexity Analysis
There are always 6 orientations, and each orientation stores only 9 transition counts. Every step from \(n\) to \(n+1\) performs constant-time arithmetic on this fixed state, followed by one fixed-size cost evaluation. Hence the total running time is \(O(N)\) for summing up to \(N\), and the working memory is \(O(1)\).
Footnotes and References
- Problem page: https://projecteuler.net/problem=497
- Tower of Hanoi: Wikipedia — Tower of Hanoi
- Dynamic programming: Wikipedia — Dynamic programming
- Modular arithmetic: Wikipedia — Modular arithmetic
Problem 497 source code
C++
#include <algorithm>
#include <array>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <map>
#include <string>
#include <tuple>
#include <vector>
namespace {
using u64 = std::uint64_t;
using u128 = __uint128_t;
constexpr u64 kMod = 1'000'000'000ULL;
struct Options {
int n_max = 10'000;
bool run_checkpoints = true;
};
bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& out) {
if (arg.rfind(prefix, 0U) != 0U) {
return false;
}
const std::string tail = arg.substr(prefix.size());
if (tail.empty()) {
return false;
}
try {
out = std::stoi(tail);
} catch (...) {
return false;
}
return true;
}
bool parse_arguments(int argc, char** argv, Options& options) {
for (int i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--skip-checkpoints") {
options.run_checkpoints = false;
continue;
}
if (parse_int_after_prefix(arg, "--n-max=", options.n_max)) {
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
if (options.n_max <= 0) {
std::cerr << "--n-max must be positive.\n";
return false;
}
return true;
}
template <typename T>
using Matrix = std::array<std::array<T, 3>, 3>;
template <typename T>
Matrix<T> zero_matrix() {
Matrix<T> m{};
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
m[i][j] = static_cast<T>(0);
}
}
return m;
}
template <typename T>
Matrix<T> add_matrix(const Matrix<T>& a, const Matrix<T>& b, const T mod = 0) {
Matrix<T> out = zero_matrix<T>();
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
if (mod == 0) {
out[i][j] = a[i][j] + b[i][j];
} else {
out[i][j] = (a[i][j] + b[i][j]) % mod;
}
}
}
return out;
}
template <typename T>
Matrix<T> base_matrix_for_orientation(const int from, const int to) {
Matrix<T> m = zero_matrix<T>();
m[from][to] = static_cast<T>(1);
return m;
}
u64 mod_sub(const u64 a, const u64 b) {
return (a >= b) ? (a - b) : (a + kMod - b);
}
u64 hit_mod(const u64 k, const u64 x, const u64 y, const bool ascending) {
if (x == y) return 0ULL;
if (ascending) {
const u64 a = mod_sub(y, x);
const u64 b = (x + y + kMod - 2ULL) % kMod;
return static_cast<u64>((static_cast<u128>(a) * b) % kMod);
}
const u64 a = mod_sub(x, y);
const u64 two_k = (2ULL * k) % kMod;
const u64 b = mod_sub(two_k, (x + y) % kMod);
return static_cast<u64>((static_cast<u128>(a) * b) % kMod);
}
u128 hit_exact(const u64 k, const u64 x, const u64 y) {
if (x == y) {
return 0;
}
if (x < y) {
return static_cast<u128>(y - x) * static_cast<u128>(x + y - 2ULL);
}
return static_cast<u128>(x - y) * static_cast<u128>(2ULL * k - x - y);
}
std::string to_string_u128(u128 v) {
if (v == 0) {
return "0";
}
std::string s;
while (v > 0) {
s.push_back(static_cast<char>('0' + static_cast<int>(v % 10U)));
v /= 10U;
}
std::reverse(s.begin(), s.end());
return s;
}
struct OrientationIndex {
std::vector<std::tuple<int, int, int>> triples;
std::map<std::tuple<int, int, int>, int> index_by_tuple;
};
OrientationIndex build_orientations() {
OrientationIndex out;
const std::vector<std::tuple<int, int, int>> all = {
{0, 1, 2}, {0, 2, 1}, {1, 0, 2}, {1, 2, 0}, {2, 0, 1}, {2, 1, 0}};
out.triples = all;
for (int i = 0; i < static_cast<int>(all.size()); ++i) {
out.index_by_tuple[all[static_cast<std::size_t>(i)]] = i;
}
return out;
}
template <typename T>
std::vector<Matrix<T>> next_transition_counts(const std::vector<Matrix<T>>& prev,
const OrientationIndex& ori,
const T mod = 0) {
std::vector<Matrix<T>> next(6, zero_matrix<T>());
for (int id = 0; id < 6; ++id) {
const auto [f, t, a] = ori.triples[static_cast<std::size_t>(id)];
const int id1 = ori.index_by_tuple.at({f, a, t});
const int id2 = ori.index_by_tuple.at({a, t, f});
Matrix<T> cur = add_matrix(prev[static_cast<std::size_t>(id1)],
prev[static_cast<std::size_t>(id2)], mod);
if (mod == 0) {
cur[a][f] += 1;
cur[f][t] += 1;
cur[t][a] += 1;
} else {
cur[a][f] = (cur[a][f] + 1) % mod;
cur[f][t] = (cur[f][t] + 1) % mod;
cur[t][a] = (cur[t][a] + 1) % mod;
}
next[static_cast<std::size_t>(id)] = cur;
}
return next;
}
u64 E_mod(const Matrix<u64>& c, const u64 k, const u64 a, const u64 b, const u64 cc) {
std::array<u64, 3> pos{a, b, cc};
u64 out = hit_mod(k, b, a, false); // initial walk to first source rod
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
if (c[i][j] == 0ULL) {
continue;
}
const bool ascending = (i < j);
const u64 h = hit_mod(k, pos[static_cast<std::size_t>(i)],
pos[static_cast<std::size_t>(j)], ascending);
out = (out + static_cast<u64>((static_cast<u128>(c[i][j]) * h) % kMod)) % kMod;
}
}
return out;
}
u128 E_exact(const Matrix<u128>& c, const u64 k, const u64 a, const u64 b, const u64 cc) {
std::array<u64, 3> pos{a, b, cc};
u128 out = hit_exact(k, b, a);
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
if (c[i][j] == 0U) {
continue;
}
out += c[i][j] * hit_exact(k, pos[static_cast<std::size_t>(i)],
pos[static_cast<std::size_t>(j)]);
}
}
return out;
}
u64 solve_sum_last9(const int n_max) {
const OrientationIndex ori = build_orientations();
const int canonical = ori.index_by_tuple.at({0, 2, 1}); // move from rod 0 to rod 2 using rod 1
std::vector<Matrix<u64>> cur_mod(6, zero_matrix<u64>());
for (int id = 0; id < 6; ++id) {
const auto [f, t, a] = ori.triples[static_cast<std::size_t>(id)];
(void)a;
cur_mod[static_cast<std::size_t>(id)] = base_matrix_for_orientation<u64>(f, t);
}
u64 p3 = 1ULL;
u64 p6 = 1ULL;
u64 p9 = 1ULL;
u64 p10 = 1ULL;
u64 total = 0ULL;
for (int n = 1; n <= n_max; ++n) {
if (n > 1) {
cur_mod = next_transition_counts(cur_mod, ori, kMod);
}
p3 = static_cast<u64>((static_cast<u128>(p3) * 3ULL) % kMod);
p6 = static_cast<u64>((static_cast<u128>(p6) * 6ULL) % kMod);
p9 = static_cast<u64>((static_cast<u128>(p9) * 9ULL) % kMod);
p10 = static_cast<u64>((static_cast<u128>(p10) * 10ULL) % kMod);
const u64 e = E_mod(cur_mod[static_cast<std::size_t>(canonical)], p10, p3, p6, p9);
total += e;
total %= kMod;
}
return total;
}
bool run_checkpoints() {
const OrientationIndex ori = build_orientations();
const int canonical = ori.index_by_tuple.at({0, 2, 1});
std::vector<Matrix<u128>> cur(6, zero_matrix<u128>());
for (int id = 0; id < 6; ++id) {
const auto [f, t, a] = ori.triples[static_cast<std::size_t>(id)];
(void)a;
cur[static_cast<std::size_t>(id)] = base_matrix_for_orientation<u128>(f, t);
}
// n = 2
cur = next_transition_counts(cur, ori, static_cast<u128>(0));
if (E_exact(cur[static_cast<std::size_t>(canonical)], 5ULL, 1ULL, 3ULL, 5ULL) !=
static_cast<u128>(60ULL)) {
std::cerr << "Checkpoint failed: E(2,5,1,3,5)\n";
return false;
}
// n = 3
cur = next_transition_counts(cur, ori, static_cast<u128>(0));
if (E_exact(cur[static_cast<std::size_t>(canonical)], 20ULL, 4ULL, 9ULL, 17ULL) !=
static_cast<u128>(2358ULL)) {
std::cerr << "Checkpoint failed: E(3,20,4,9,17)\n";
return false;
}
return true;
}
} // namespace
int main(int argc, char** argv) {
Options options;
if (!parse_arguments(argc, argv, options)) {
return 1;
}
if (options.run_checkpoints && !run_checkpoints()) {
return 1;
}
const u64 ans = solve_sum_last9(options.n_max);
std::cout << std::setw(9) << std::setfill('0') << ans << '\n';
return 0;
}
Python
MOD = 1000000000
def zero_matrix():
return [[0] * 3 for _ in range(3)]
def add_matrix(a, b, mod=0):
out = zero_matrix()
for i in range(3):
for j in range(3):
if mod == 0:
out[i][j] = a[i][j] + b[i][j]
else:
out[i][j] = (a[i][j] + b[i][j]) % mod
return out
def base_matrix_for_orientation(frm, to):
m = zero_matrix()
m[frm][to] = 1
return m
def mod_sub(a, b):
return (a - b) if a >= b else (a + MOD - b)
def hit_mod(k, x, y, ascending):
if x == y: return 0
if ascending:
a = mod_sub(y, x)
b = (x + y + MOD - 2) % MOD
return (a * b) % MOD
a = mod_sub(x, y)
two_k = (2 * k) % MOD
b = mod_sub(two_k, (x + y) % MOD)
return (a * b) % MOD
class OrientationIndex:
def __init__(self):
self.triples = [
(0, 1, 2), (0, 2, 1), (1, 0, 2),
(1, 2, 0), (2, 0, 1), (2, 1, 0)
]
self.index_by_tuple = {t: i for i, t in enumerate(self.triples)}
def next_transition_counts(prev, ori, mod=0):
nxt = [zero_matrix() for _ in range(6)]
for id_ in range(6):
f, t, a = ori.triples[id_]
id1 = ori.index_by_tuple[(f, a, t)]
id2 = ori.index_by_tuple[(a, t, f)]
cur = add_matrix(prev[id1], prev[id2], mod)
if mod == 0:
cur[a][f] += 1
cur[f][t] += 1
cur[t][a] += 1
else:
cur[a][f] = (cur[a][f] + 1) % mod
cur[f][t] = (cur[f][t] + 1) % mod
cur[t][a] = (cur[t][a] + 1) % mod
nxt[id_] = cur
return nxt
def E_mod(c, k, a, b, cc):
pos = [a, b, cc]
out = hit_mod(k, b, a, False)
for i in range(3):
for j in range(3):
if c[i][j] == 0: continue
ascending = (i < j)
h = hit_mod(k, pos[i], pos[j], ascending)
out = (out + c[i][j] * h) % MOD
return out
def solve_sum_last9(n_max):
ori = OrientationIndex()
canonical = ori.index_by_tuple[(0, 2, 1)]
cur_mod = [zero_matrix() for _ in range(6)]
for id_ in range(6):
f, t, a = ori.triples[id_]
cur_mod[id_] = base_matrix_for_orientation(f, t)
p3, p6, p9, p10 = 1, 1, 1, 1
total = 0
for n in range(1, n_max + 1):
if n > 1:
cur_mod = next_transition_counts(cur_mod, ori, MOD)
p3 = (p3 * 3) % MOD
p6 = (p6 * 6) % MOD
p9 = (p9 * 9) % MOD
p10 = (p10 * 10) % MOD
e = E_mod(cur_mod[canonical], p10, p3, p6, p9)
total = (total + e) % MOD
return total
def solve():
return f"{solve_sum_last9(10000):09d}"
if __name__ == '__main__':
print(solve())
Java
import java.util.HashMap;
import java.util.Map;
import java.util.Objects;
public class Euler497 {
private static final long MOD = 1000000000L;
private static long[][] zeroMatrix() {
return new long[3][3];
}
private static long[][] addMatrix(long[][] a, long[][] b, long mod) {
long[][] out = zeroMatrix();
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
if (mod == 0) {
out[i][j] = a[i][j] + b[i][j];
} else {
out[i][j] = (a[i][j] + b[i][j]) % mod;
}
}
}
return out;
}
private static long[][] baseMatrixForOrientation(int from, int to) {
long[][] m = zeroMatrix();
m[from][to] = 1L;
return m;
}
private static long modSub(long a, long b) {
return (a >= b) ? (a - b) : (a + MOD - b);
}
private static long hitMod(long k, long x, long y, boolean ascending) {
if (x == y)
return 0L;
if (ascending) {
long a = modSub(y, x);
long b = (x + y + MOD - 2L) % MOD;
return (a * b) % MOD;
}
long a = modSub(x, y);
long twoK = (2L * k) % MOD;
long b = modSub(twoK, (x + y) % MOD);
return (a * b) % MOD;
}
static class Tuple {
int f, t, a;
Tuple(int f, int t, int a) {
this.f = f;
this.t = t;
this.a = a;
}
@Override
public int hashCode() {
return Objects.hash(f, t, a);
}
@Override
public boolean equals(Object o) {
if (!(o instanceof Tuple))
return false;
Tuple that = (Tuple) o;
return f == that.f && t == that.t && a == that.a;
}
}
static class OrientationIndex {
Tuple[] triples;
Map<Tuple, Integer> indexByTuple;
OrientationIndex() {
triples = new Tuple[] {
new Tuple(0, 1, 2), new Tuple(0, 2, 1), new Tuple(1, 0, 2),
new Tuple(1, 2, 0), new Tuple(2, 0, 1), new Tuple(2, 1, 0)
};
indexByTuple = new HashMap<>();
for (int i = 0; i < triples.length; i++) {
indexByTuple.put(triples[i], i);
}
}
}
private static long[][][] nextTransitionCounts(long[][][] prev, OrientationIndex ori, long mod) {
long[][][] next = new long[6][][];
for (int id = 0; id < 6; ++id) {
Tuple t = ori.triples[id];
int id1 = ori.indexByTuple.get(new Tuple(t.f, t.a, t.t));
int id2 = ori.indexByTuple.get(new Tuple(t.a, t.t, t.f));
long[][] cur = addMatrix(prev[id1], prev[id2], mod);
if (mod == 0) {
cur[t.a][t.f]++;
cur[t.f][t.t]++;
cur[t.t][t.a]++;
} else {
cur[t.a][t.f] = (cur[t.a][t.f] + 1) % mod;
cur[t.f][t.t] = (cur[t.f][t.t] + 1) % mod;
cur[t.t][t.a] = (cur[t.t][t.a] + 1) % mod;
}
next[id] = cur;
}
return next;
}
private static long eMod(long[][] c, long k, long a, long b, long cc) {
long[] pos = { a, b, cc };
long out = hitMod(k, b, a, false);
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
if (c[i][j] == 0L)
continue;
boolean ascending = (i < j);
long h = hitMod(k, pos[i], pos[j], ascending);
out = (out + c[i][j] * h) % MOD;
}
}
return out;
}
public static void main(String[] args) {
int nMax = 10000;
OrientationIndex ori = new OrientationIndex();
int canonical = ori.indexByTuple.get(new Tuple(0, 2, 1));
long[][][] curMod = new long[6][][];
for (int id = 0; id < 6; ++id) {
Tuple t = ori.triples[id];
curMod[id] = baseMatrixForOrientation(t.f, t.t);
}
long p3 = 1L;
long p6 = 1L;
long p9 = 1L;
long p10 = 1L;
long total = 0L;
for (int n = 1; n <= nMax; ++n) {
if (n > 1) {
curMod = nextTransitionCounts(curMod, ori, MOD);
}
p3 = (p3 * 3L) % MOD;
p6 = (p6 * 6L) % MOD;
p9 = (p9 * 9L) % MOD;
p10 = (p10 * 10L) % MOD;
long e = eMod(curMod[canonical], p10, p3, p6, p9);
total = (total + e) % MOD;
}
System.out.printf("%09d\n", total);
}
}