Problem 312: Cyclic Paths on Sierpiński Graphs

Project Euler Problem 312 Solution

EulerSolve provides an optimized solution for Project Euler Problem 312, Cyclic Paths on Sierpiński Graphs, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $C(n)$ be the number of Hamiltonian cycles on the Sierpinski graph $S_n$ described in the problem. We are asked to compute $$C(C(C(10000))) \pmod{13^8}.$$ The direct values explode immediately, so the only viable route is: $$\text{graph recurrence} \Longrightarrow \text{closed form} \Longrightarrow \text{modular exponent reduction}.$$ Mathematical Approach 1) Decompose $S_{n+1}$ into three copies of $S_n$. The graph $S_{n+1}$ consists of three corner-sharing copies of $S_n$. A Hamiltonian cycle on $S_{n+1}$ cannot stay closed inside one copy; it must enter each copy through one shared corner and leave through the other shared corner. Therefore each copy contributes a Hamiltonian path between two corner vertices . Let $P(n)$ denote the number of such corner-to-corner Hamiltonian paths in $S_n$ for a fixed ordered copy position. Then the three copies are independent once their corner roles are fixed, so $$C(n+1)=P(n)^3.$$ 2) Why $P(n)=3C(n)$. Inside one copy of $S_n$, a Hamiltonian path is obtained from a Hamiltonian cycle by choosing which one of the three outer corners plays the role of the “unused outer corner” when the copy is embedded into the next level. There are exactly three symmetric choices, and each choice converts bijectively between the cycle state and the required corner-to-corner path state....

Detailed mathematical approach

Problem Summary

A Sierpiński graph of order $1$, written $S_1$, is a triangle: three vertices, three edges. The graph $S_{n+1}$ is obtained by positioning three copies of $S_n$ so that every pair of copies has one common corner. Let $C(n)$ be the number of cycles that pass exactly once through all the vertices of $S_n$: its Hamiltonian cycles. A path through all the vertices is a Hamiltonian path. A cycle counts as a set of edges: reversing it or starting elsewhere gives the same cycle, while cycles related by a symmetry of the graph count separately. ^[2] The statement gives $C(1)=C(2)=1$, $C(3)=8$ (the eight cycles are drawn there), $C(5)=71328803586048$, $C(10\,000)\bmod 10^8=37652224$ and $C(10\,000)\bmod 13^8=617720485$, and asks for

$$C(C(C(10\,000)))\bmod 13^8.$$

The whole solution in one sentence: for $n\ge3$ the Hamiltonian cycles of $S_n$ are counted exactly by $8\cdot 12^{(3^{n-2}-3)/2}$, and modulo $13^8$ such a power depends on $n$ only through the residue of $n$ modulo $12\cdot 13^6$, so the triple nesting collapses into a handful of modular exponentiations. ^[1]

Mathematical Approach

1. The graphs, their size and their degrees

Each copy of $S_n$ inside $S_{n+1}$ has three corners: two are glued to the neighbouring copies, the third stays a corner of $S_{n+1}$. Gluing identifies three pairs of vertices, so $V_{n+1}=3V_n-3$ and $E_{n+1}=3E_n$; with $V_1=E_1=3$,

$$V_n=\frac{3^n+3}{2},\qquad E_n=3^n,$$

that is $3,6,15,42,123$ vertices and $3,9,27,81,243$ edges for $n=1,\dots,5$. ^[3] Every vertex of $S_1$ has degree $2$, and a glued vertex inherits degree $2$ from each of its copies, so in $S_n$ exactly the three corners have degree $2$ and the other $V_n-3$ vertices have degree $4$. These graphs are the Sierpiński gasket graphs of the literature; the name Sierpiński graph there denotes a different family $S(n,k)$ with $k^n$ vertices. ^[3]

2. Three quantities, and the smallest cases by hand

Counting cycles alone does not recurse: two kinds of paths must be counted alongside them. For the corners of $S_n$ define

$C(n)$: the Hamiltonian cycles of $S_n$;
$P(n)$: the Hamiltonian paths whose two endpoints are two given corners;
$Q(n)$: the paths between two given corners that cover every vertex except the third corner.

The six symmetries of the triangle act on $S_n$, so $P(n)$ and $Q(n)$ do not depend on which corners are chosen. In $S_1$ with corners $a,b,c$ the only Hamiltonian path from $a$ to $b$ is $a\to c\to b$ and the only path avoiding $c$ is the edge itself, so $P(1)=Q(1)=C(1)=1$.

In $S_2$ call the glued vertices $p,q,r$ and the corners $u,v,w$, so that the three triangles are $\{u,p,r\}$, $\{v,p,q\}$ and $\{w,q,r\}$. The Hamiltonian paths from $u$ to $v$ are $u\to r\to w\to q\to p\to v$ and $u\to p\to r\to w\to q\to v$, so $P(2)=2$; those avoiding $w$ are $u\to r\to q\to p\to v$, $u\to r\to p\to q\to v$ and $u\to p\to r\to q\to v$, so $Q(2)=3$; and $C(2)=1$.

3. What a Hamiltonian cycle looks like inside one copy

Write the three copies of $S_n$ inside $S_{n+1}$ as $X,Y,Z$, with $p$ shared by $X$ and $Y$, $q$ by $Y$ and $Z$, $r$ by $Z$ and $X$, and let $u\in X$, $v\in Y$, $w\in Z$ be the corners of $S_{n+1}$: so $X$ has corners $u,p,r$, $Y$ has $v,p,q$ and $Z$ has $w,q,r$. For a subgraph $H$ write $d_K(z)$ for the number of edges of $H$ at $z$ lying in the copy $K$.

Let $H$ be a Hamiltonian cycle of $S_{n+1}$. Every vertex $z$ of $X$ other than $p$ and $r$ has all its neighbours inside $X$, hence $d_X(z)=2$. Being a proper subgraph of a cycle, $H\cap X$ is a disjoint union of paths, isolated vertices allowed. The degrees of any graph sum to twice its number of edges, so the number of odd degrees is even, and here only $p$ and $r$ can be odd. They are not both even: a path with an edge has two ends of degree $1$, so with every degree even $H\cap X$ would have no edge at all — impossible, since both edges of $H$ at $u$ lie in $X$. So both are odd, and neither exceeds $2$:

$$d_X(p)=d_X(r)=1 .$$

Then $p$ and $r$ are the only ends and every other vertex has degree $2$: $H\cap X$ is a Hamiltonian path of $X$ from $p$ to $r$, and likewise in $Y$ (from $p$ to $q$) and in $Z$ (from $q$ to $r$). Conversely three such paths meet pairwise in exactly one endpoint, so they join end to end into one cycle through every vertex, and the three choices are independent. Hence

$$C(n+1)=P(n)^3 .$$

In particular $C(3)=P(2)^3=2^3=8$, the eight pictures of the statement.

4. The two recursions for the paths

Let $H$ be a Hamiltonian path of $S_{n+1}$ from $u$ to $v$; $P(n+1)$ counts these. In $Z$ the corner $w$ is an interior vertex with both edges inside $Z$, so Section 3 applies word for word: $H\cap Z$ is a Hamiltonian path of $Z$ from $q$ to $r$, in $P(n)$ ways. Each of $q$ and $r$ thus spends one of its two $H$-edges in $Z$, leaving $d_X(r)=1$ and $d_Y(q)=1$.

In $X$ the vertex $u$ is an end of $H$, so $d_X(u)=1$. Since the number of odd degrees is even and only $u,p,r$ can be odd, $d_X(p)$ is even, so it is $0$ or $2$.

$d_X(p)=2$: the ends in $X$ are $u$ and $r$ and everything else has degree $2$, so $H\cap X$ is a Hamiltonian path from $u$ to $r$ — $P(n)$ ways.
$d_X(p)=0$: both edges of $H$ at $p$ lie in $Y$, so $H\cap X$ covers $V(X)$ without $p$ and runs from $u$ to $r$ — $Q(n)$ ways.

The copy $Y$ behaves the same way with $v$ for $u$, and the two edges at $p$ are split between $X$ and $Y$, so exactly one copy takes $d(p)=2$ and the other $0$. Conversely each allowed combination joins into one path $u\to\dots\to r\to\dots\to q\to\dots\to v$ through every vertex once, and restricting it recovers the same choices, so nothing is counted twice. That gives two cases worth $P(n)Q(n)$ each, times the $P(n)$ choices in $Z$:

$$P(n+1)=2\,P(n)^2Q(n)\qquad (n\ge1).$$

For $Q(n+1)$ the corner $w$ is deleted. Only $Z$ changes: $H\cap Z$ covers $V(Z)$ without $w$ and still runs from $q$ to $r$, which is $Q(n)$ ways, while $X$ and $Y$ contribute $P(n)Q(n)$ in two ways as before. Hence

$$Q(n+1)=2\,P(n)Q(n)^2\qquad (n\ge2).$$

From $P(2)=2$ and $Q(2)=3$ these give $P(3)=24$, $Q(3)=36$, $P(4)=41472$ and $Q(4)=62208$. An exhaustive search on $S_3$ and $S_4$ returns exactly these four numbers, and $C(4)=P(3)^3=24^3=13824$.

5. A pitfall at the bottom of the recursion

The argument for $P$ survives at $n=1$, giving $P(2)=2$, the right value; the one for $Q$ does not. For $n=1$ the copy $Z$ is only its three corners, so deleting $w$ leaves just $q$ and $r$; then $H\cap Z$ may be empty, with $q$ and $r$ covered by the other copies and $X$, $Y$ each carrying a Hamiltonian path ending at $p$. That extra case contributes $P(1)^2=1$, so

$$Q(2)=2\,P(1)Q(1)^2+P(1)^2=3 ,$$

not $2$. Using $Q(2)=2$ would give $P(3)=16$ and $C(4)=4096$ instead of $13824$, and every later value wrong. For $n\ge2$ the copy $Z$ has vertices besides its three corners; if $d_Z(q)$ and $d_Z(r)$ were both even, every degree in $H\cap Z$ would be even and $H\cap Z$ would have no edge, leaving those vertices off $H$ altogether. So the degenerate case cannot return.

6. One recursion, then a closed form

Dividing the two recursions gives $Q(n+1)/P(n+1)=Q(n)/P(n)$ for $n\ge2$, and the ratio starts at $Q(2)/P(2)=3/2$. So $Q(n)=\tfrac32P(n)$ for every $n\ge2$, and substituting removes $Q$ entirely:

$$P(n+1)=3P(n)^3\qquad (n\ge2),\quad P(2)=2 .$$

Comparing with $C(n+1)=P(n)^3$ shows $P(n+1)=3\,C(n+1)$, so $P(m)=3C(m)$ for $m\ge3$ and therefore

$$C(m+1)=\bigl(3C(m)\bigr)^3\qquad (m\ge3),$$

starting from $C(3)=8$. Writing $C(n)=8\cdot 12^{e_n}$ turns this into $e_{n+1}=3e_n+3$ with $e_3=0$, since $(3\cdot 8\cdot 12^{e})^3=8\cdot 12^{3e+3}$. As $e_n+\tfrac32$ triples at each step, $e_n=(3^{n-2}-3)/2$ and

$$\begin{aligned}C(n)&=8\cdot 12^{(3^{n-2}-3)/2}\\&=2^{3^{n-2}}\,3^{(3^{n-2}-3)/2}\end{aligned}$$

for $n\ge3$. The exponents $e_4=3$, $e_5=12$, $e_6=39$ give $C(4)=13824$, $C(5)=8\cdot 12^{12}=71328803586048$ — the value in the statement — and a $43$-digit $C(6)$, all three matching the published sequence of these counts. ^[4]

7. Why residues are unavoidable

The number of decimal digits of $C(n)$ is about $3^{n-2}\log_{10}(2\sqrt3)$, roughly $0.54\cdot 3^{n-2}$: already $C(10)$ has $3540$ digits, and $C(10\,000)$ has roughly $0.54\cdot 3^{9998}$ of them, a digit count that itself needs $4770$ digits. The inner values can never be formed, so the task is to compute $C(N)\bmod 13^8$ knowing $N$ only through a residue.

8. The order of 12 modulo a power of 13

Since $13$ does not divide $12$, the powers $12^e\bmod 13^k$ repeat with period $\operatorname{ord}_{13^k}(12)$, the least $t\ge1$ with $12^t\equiv1\pmod{13^k}$, so the exponent matters only modulo that order. ^[6] Because $12\equiv-1\pmod{13}$, an odd power of $12$ is $\equiv-1\pmod{13}$ and cannot be $1$ modulo $13^k$, so the order is even. For an even exponent $2m$ write $12^{2m}-1=144^m-1^m$. Here $13$ divides $144-1=143=11\cdot13$ exactly once, so the lifting-the-exponent lemma ^[5] gives

$$v_{13}\bigl(144^m-1\bigr)=1+v_{13}(m),$$

where $v_{13}$ is the exponent of $13$ in a factorisation. So $13^k$ divides $12^{2m}-1$ exactly when $13^{k-1}$ divides $m$, and

$$\operatorname{ord}_{13^k}(12)=2\cdot 13^{k-1},$$

which is $4394$ for $k=4$, $742586$ for $k=6$ and $125497034$ for $k=8$. Direct computation confirms these orders for all $k\le8$.

9. Halving the exponent safely

By Section 8 we need $e_N=(3^{N-2}-3)/2$ modulo $M=2\cdot 13^{k-1}$. Halving is not well defined modulo $M$, so the doubled exponent is reduced instead: compute $2e_N=3^{N-2}-3$ modulo $2M=4\cdot 13^{k-1}$. If $s$ is that residue then $s$ is even, because $3^{N-2}-3$ is even and the modulus $2M$ is even, and $2e_N\equiv s \pmod{2M}$ gives $e_N\equiv s/2\pmod{M}$. In short: take $x=3^{N-2}\bmod(4\cdot 13^{k-1})$, which is odd, set $s=(x-3)\bmod(4\cdot 13^{k-1})$, and halve.

10. How much of $N$ is needed

What remains is $3^{N-2}$ modulo $4\cdot 13^{k-1}$ for an unknown huge $N$. As $3$ is coprime to that modulus, the exponent $N-2$ may be reduced modulo any multiple of the order of $3$. The Carmichael function $\lambda(m)$ is the smallest exponent with $x^{\lambda(m)}\equiv1\pmod m$ for every $x$ coprime to $m$; it is a multiple of every such order, and for coprime factors it is the lcm of their values: ^[6]

$$\begin{aligned}\lambda(4\cdot 13^{k-1})&=\operatorname{lcm}(2,12\cdot 13^{k-2})\\&=12\cdot 13^{k-2}\end{aligned}$$

for $k\ge2$. Hence, for $N\ge3$, the value $C(N)\bmod 13^k$ is determined by $N\bmod(12\cdot 13^{k-2})$, which is $57921708$ for $k=8$. The true order of $3$ is only $6\cdot 13^{k-2}$, half of $\lambda$, for every $2\le k\le8$, but any multiple of the order serves.

An example with $k=3$, checkable by hand: $13^3=2197$, the order of $12$ is $338$, the doubled modulus $676$ and $\lambda=156$. For $N=10$ the reduced exponent is $(10-2)\bmod156=8$, so $x=3^8\bmod676=477$, $s=474$ and $e=237$ — and the exact $e_{10}=(3^8-3)/2=3279$ is indeed $\equiv237\pmod{338}$. Then $8\cdot 12^{237}\equiv811\pmod{2197}$, which is what the recurrence $C(m+1)=(3C(m))^3$ from $C(3)=8$ gives modulo $2197$, for $N=10$ and $N=166$ alike, as the period $156$ predicts.

11. The inner values are multiples of 12

Knowing $N$ modulo $12\cdot 13^{j}$ means knowing it modulo $12$ and modulo $13^{j}$, and the first half is free: for $N\ge4$ the exponent $e_N\ge3$, so $C(N)=8\cdot 12^{e_N}$ is a multiple of $12$, indeed of $96$. Every inner value of the tower is therefore $\equiv0\pmod{12}$, only its residue modulo a power of $13$ has to be computed, and since $\gcd(12,13^j)=1$ the Chinese remainder theorem ^[7] combines the two congruences into a unique residue modulo $12\cdot 13^j$.

12. The chain of three levels

Let $a=C(10\,000)$ and $b=C(a)$; the answer is $C(b)\bmod 13^8$. Reading Section 10 backwards fixes the precision of each level:

$C(b)\bmod 13^8$ needs $b$ modulo $12\cdot 13^6=57921708$;
$b\bmod 13^6$ needs $a$ modulo $12\cdot 13^4=342732$;
$a\bmod 13^4$ needs only the exact $N=10\,000$.

The first level gives $a\equiv3177\pmod{28561}$, so $a=3177+28561t$ for some $t$. As $3177\equiv9$ and $28561\equiv1$ modulo $12$, the condition $a\equiv0\pmod{12}$ forces $t\equiv3\pmod{12}$, and $t=3$ gives $a\equiv88860\pmod{342732}$. The published checkpoint $C(10\,000)\bmod 13^8=617720485$ is a free cross-check: it reduces to the same $3177$ modulo $13^4$. The second level feeds $88860$ into the closed form with $k=6$ and lifts with $b\equiv0$ modulo $12$; the third evaluates the closed form with $k=8$. Both $a$ and $b$ are far larger than $4$, so the special values $C(1)=C(2)=1$ never interfere.

How the Code Works

The C++ program ^[8] is a direct transcription of Sections 8 to 12. mul_mod multiplies two residues, pow_mod is square-and-multiply exponentiation on top of it and pow_u64 builds the small powers $13^k$. inverse_mod inverts a unit modulo $13^k$ as $a^{\varphi-1}$ with $\varphi=13^k-13^{k-1}$, by Euler's theorem rather than the extended Euclidean algorithm — correct only because the modulus is a power of $13$. crt_coprime is Section 11: from $x\equiv a_1\pmod{m_1}$ and $x\equiv a_2\pmod{m_2}$ it returns $a_1+m_1t$ with $t=(a_2-a_1)m_1^{-1}\bmod m_2$.

The heart is c_mod_13_power_from_n_residue, with mod $=13^k$, ord12 $=2\cdot 13^{k-1}$, two_ord twice that and lambda $=12\cdot 13^{k-2}$ (or $2$ for $k=1$, a branch solve never needs). The exponent (n_residue + lambda - 2) % lambda stands for $n-2$, the added lambda keeping the subtraction non-negative in unsigned arithmetic; the helper raises $3$ to it modulo two_ord, forms numer = (x + two_ord - 3) % two_ord, the non-negative even $s$ of Section 9, halves it, reduces modulo ord12 and returns $8\cdot 12^{e}\bmod 13^k$ — Sections 10, 9 and 8 in that order. The flag n_is_large says whether n_residue is the true $n$, when $n\le2$ must return $1$, or only a residue of a huge value, where the closed form always applies. solve walks the three levels of Section 12, calling crt_coprime(0, 12, …) twice because the inner values are multiples of $12$.

c_mod_iterative is the independent route: it iterates $c\mapsto(3c)^3$ from $C(3)=8$, the recurrence of Section 6, and serves only the checks. run_checkpoints verifies $C(1)=C(2)=1$ and $C(5)=71328803586048$ modulo $10^{15}$, above $C(5)$ itself, so these compare exact integers; it then checks both published residues of $C(10\,000)$ and that the closed form agrees with the iterative route at $n=10\,000$ modulo $13^8$. main runs the checks unless --skip-checkpoints is given.

C++ keeps residues in std::uint64_t and widens each product to unsigned __int128 in mul_mod, so it is correct for any modulus below $2^{64}$; Java keeps them in long and lets BigInteger widen, with multiply then mod in mul_mod and modPow in pow_mod; Python needs neither, its integers being unbounded and its three-argument pow already modular. ^[9] The widening is headroom, not necessity: the largest modulus anywhere is $13^8=815730721$, so no product exceeds $7\cdot10^{17}$, within signed $64$-bit range, and the checkpoint modulus $10^{15}$ merely has to exceed $C(5)$. Only the C++ file carries the checkpoints; the two ports evaluate the chain alone.

Complexity and Verification

The chain evaluates the closed form three times — one exponentiation modulo $4\cdot 13^{k-1}$ and one modulo $13^k$ each time — plus an inversion in each of the two Chinese remainder steps. The exponent of $3$ stays below $\lambda=12\cdot 13^{k-2}$, that of $12$ below $2\cdot 13^{k-1}$ and the inverses use $\varphi-1$; the largest of them is $2\cdot 13^7=125497034$, of $27$ bits, so one exponentiation costs at most $27$ squarings and $27$ multiplications: $O(\log 13^k)$ per level, fewer than a thousand modular multiplications in all, in $O(1)$ memory. The checkpoints dominate: the two long recurrence checks run $9997$ iterations each with two mul_mod calls per iteration, $39988$ in total. Everything finishes instantly.

The following were verified for this article. An exhaustive search, pruned by degrees and connectivity so that the $42$-vertex case stays fast, counted the Hamiltonian structures of $S_1,\dots,S_4$: $C=1,1,8,13824$, $P=1,2,24,41472$ and $Q=1,3,36,62208$, equal for all three corner pairs as Section 2 assumes and matching $C(3)=8$ and the recursions of Sections 3 to 5. The closed form agrees with the recurrence for all $3\le n\le60$ modulo $10^8$ and every $13^k$, $k\le8$; $C(5)$ and $C(6)$ match the exact integers and the published sequence. ^[4] The order of $12$ modulo $13^k$ comes out as $2\cdot 13^{k-1}$ for $k\le8$, and the period of Section 10 was checked over a full period of $N$ for $k=2,3,4$.

The residue chain was tested three ways: carrying more precision than Section 12 requires, the inner values modulo $12\cdot 13^{10}$ and $12\cdot 13^{12}$, leaves the final value unchanged, as does replacing a residue by another lift; towers starting from $n=5,6,7$, where the first inner value is an exact integer usable as a true exponent, give the same result reduced or not; and $a\bmod(12\cdot 13^4)=88860$ agrees between the closed form, the recurrence modulo $342732$ and the published checkpoint. The three programs print the same value, the C++ checkpoints pass, and an independent reimplementation agrees with them.

Footnotes and References

The references below supply the background facts used above. The recursions for $C$, $P$ and $Q$, the closed form, the order computation and the precision analysis of the three levels are derived explicitly in this article.

Project Euler 312 — Cyclic Paths on Sierpiński Graphs. The official statement builds S(n+1) from three copies of S(n) sharing one corner per pair, defines C(n) as the number of cycles through every vertex of S(n) exactly once, gives C(1) = C(2) = 1, C(3) = 8, C(5) = 71328803586048 and the two residues of C(10000) quoted above, and asks for C(C(C(10000))) mod 13⁸.
Wikipedia — Hamiltonian path. A Hamiltonian path visits each vertex exactly once and a Hamiltonian cycle is a cycle that does so; deciding whether either exists is NP-complete in general.
MathWorld — Sierpiński Gasket Graph gives the vertex count 3(3ⁿ⁻¹ + 1)/2 and the edge count 3ⁿ of Section 1, and reports that these graphs are Hamiltonian and pancyclic (Teguia and Godbole, 2006). Wikipedia — Sierpiński graph describes the other family S(n, k) of n-tuples over k symbols, with kⁿ vertices.
OEIS A246959 — "Numbers of (undirected) Hamiltonian cycles in the n-Sierpiński gasket graph" — lists 1, 1, 8, 13824, 71328803586048 and the 43-digit sixth term, and states both a(n) = 8 · 12^((3^(n−2) − 3)/2) for n ≥ 3 and a(n) = (3a(n−1))³ for n ≥ 4, the formulas proved in Sections 3 to 6.
Wikipedia — Lifting-the-exponent lemma. For an odd prime p and integers x, y not divisible by p with p dividing x − y, the lemma states v_p(xⁿ − yⁿ) = v_p(x − y) + v_p(n); Section 8 applies it with x = 144, y = 1 and p = 13.
Wikipedia — Multiplicative order defines the order of a unit a modulo n as the smallest positive k with a^k ≡ 1 (mod n) and notes that it divides λ(n). Wikipedia — Carmichael function gives λ(n) = φ(n) for 1, 2, 4 and odd prime powers, hence λ(4) = 2 and λ(13^r) = 12 · 13^(r−1), together with the lcm rule used in Section 10.
Wikipedia — Chinese remainder theorem. For coprime moduli the simultaneous congruences have a solution, unique modulo the product, constructed from Bézout coefficients; Sections 11 and 12 use the moduli 12 and a power of 13.
C++ — Euler312.cpp. The linked immutable revision contains the functions walked through above and accepts --skip-checkpoints. It is the implementation record for this article.
GCC manual — 128-bit integers documents the unsigned __int128 extension used by mul_mod, Java — BigInteger the multiply, mod and modPow used by the Java port, and Python — pow the three-argument form that computes a modular power directly.

Problem 312 source code

C++

#include <cstdint>
#include <iostream>
#include <string>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

u64 mul_mod(const u64 a, const u64 b, const u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * static_cast<u128>(b)) % static_cast<u128>(mod));
}

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

u64 pow_u64(u64 base, int exp) {
    u64 result = 1ULL;
    while (exp-- > 0) {
        result *= base;
    }
    return result;
}

u64 inverse_mod(u64 a, const u64 mod) {
    a %= mod;
    // mod is prime-power of 13 in this problem, so a^(phi-1) gives inverse for gcd(a,mod)=1.
    const u64 phi = mod - mod / 13ULL;
    return pow_mod(a, phi - 1ULL, mod);
}

u64 crt_coprime(const u64 a1, const u64 m1, const u64 a2, const u64 m2) {
    const u64 inv_m1 = inverse_mod(m1 % m2, m2);
    const u64 delta = (a2 >= a1 % m2) ? (a2 - (a1 % m2)) : (a2 + m2 - (a1 % m2));
    const u64 t = mul_mod(delta, inv_m1, m2);
    return a1 + m1 * t;
}

u64 c_mod_iterative(const u64 n, const u64 mod) {
    if (mod == 1ULL) {
        return 0ULL;
    }
    if (n <= 2ULL) {
        return 1ULL % mod;
    }

    u64 c = 8ULL % mod;  // C(3)
    for (u64 i = 4ULL; i <= n; ++i) {
        const u64 v = (3ULL * c) % mod;
        c = mul_mod(mul_mod(v, v, mod), v, mod);
    }
    return c;
}

u64 c_mod_13_power_from_n_residue(const u64 n_residue, const bool n_is_large, const int k) {
    const u64 mod = pow_u64(13ULL, k);

    if (!n_is_large && n_residue <= 2ULL) {
        return 1ULL % mod;
    }

    const u64 ord12 = 2ULL * pow_u64(13ULL, k - 1);        // ord_{13^k}(12)
    const u64 two_ord = 2ULL * ord12;                      // 4 * 13^{k-1}
    const u64 lambda = (k >= 2) ? (12ULL * pow_u64(13ULL, k - 2)) : 2ULL;

    const u64 exp = (n_residue + lambda - 2ULL) % lambda;  // n - 2 mod lambda(two_ord)
    const u64 x = pow_mod(3ULL, exp, two_ord);              // x = 3^{n-2} mod (2*ord)

    const u64 numer = (x + two_ord - 3ULL) % two_ord;       // x - 3 is even in this modulus
    const u64 e = (numer / 2ULL) % ord12;                   // ((3^{n-2}-3)/2) mod ord

    return mul_mod(8ULL % mod, pow_mod(12ULL, e, mod), mod);
}

u64 solve() {
    const u64 n0 = 10000ULL;

    // a = C(10000). For n >= 4, C(n) is divisible by 12.
    const u64 a_mod_13_4 = c_mod_13_power_from_n_residue(n0, false, 4);
    const u64 a_mod_12x13_4 = crt_coprime(0ULL, 12ULL, a_mod_13_4, pow_u64(13ULL, 4));

    // b = C(a). Need b mod (12 * 13^6) to evaluate C(b) mod 13^8.
    const u64 b_mod_13_6 = c_mod_13_power_from_n_residue(a_mod_12x13_4, true, 6);
    const u64 b_mod_12x13_6 = crt_coprime(0ULL, 12ULL, b_mod_13_6, pow_u64(13ULL, 6));

    // c = C(b) mod 13^8.
    return c_mod_13_power_from_n_residue(b_mod_12x13_6, true, 8);
}

bool run_checkpoints() {
    // Given checks from statement.
    if (c_mod_iterative(1ULL, 1000000000000000ULL) != 1ULL) {
        std::cerr << "Checkpoint failed: C(1)\n";
        return false;
    }
    if (c_mod_iterative(2ULL, 1000000000000000ULL) != 1ULL) {
        std::cerr << "Checkpoint failed: C(2)\n";
        return false;
    }
    if (c_mod_iterative(5ULL, 1000000000000000ULL) != 71328803586048ULL) {
        std::cerr << "Checkpoint failed: C(5)\n";
        return false;
    }
    if (c_mod_iterative(10000ULL, 100000000ULL) != 37652224ULL) {
        std::cerr << "Checkpoint failed: C(10000) mod 1e8\n";
        return false;
    }
    if (c_mod_iterative(10000ULL, pow_u64(13ULL, 8)) != 617720485ULL) {
        std::cerr << "Checkpoint failed: C(10000) mod 13^8\n";
        return false;
    }

    // Internal consistency: closed-form modulo 13^8 matches iterative value for n=10000.
    if (c_mod_13_power_from_n_residue(10000ULL, false, 8) != 617720485ULL) {
        std::cerr << "Checkpoint failed: closed-form/iterative mismatch\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;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve() << '\n';
    return 0;
}

Python

def inverse_mod(a, mod):
    a %= mod
    phi = mod - mod // 13
    return pow(a, phi - 1, mod)

def crt_coprime(a1, m1, a2, m2):
    inv_m1 = inverse_mod(m1 % m2, m2)
    delta = a2 - (a1 % m2)
    t = (delta * inv_m1) % m2
    return a1 + m1 * t

def c_mod_13_power_from_n_residue(n_residue, n_is_large, k):
    mod = 13 ** k
    
    if not n_is_large and n_residue <= 2:
        return 1 % mod
        
    ord12 = 2 * (13 ** (k - 1))
    two_ord = 2 * ord12
    lambda_val = 12 * (13 ** (k - 2)) if k >= 2 else 2
    
    exp = (n_residue + lambda_val - 2) % lambda_val
    x = pow(3, exp, two_ord)
    
    numer = (x + two_ord - 3) % two_ord
    e = (numer // 2) % ord12
    
    return ((8 % mod) * pow(12, e, mod)) % mod

def solve():
    n0 = 10000
    
    a_mod_13_4 = c_mod_13_power_from_n_residue(n0, False, 4)
    a_mod_12x13_4 = crt_coprime(0, 12, a_mod_13_4, 13 ** 4)
    
    b_mod_13_6 = c_mod_13_power_from_n_residue(a_mod_12x13_4, True, 6)
    b_mod_12x13_6 = crt_coprime(0, 12, b_mod_13_6, 13 ** 6)
    
    result = c_mod_13_power_from_n_residue(b_mod_12x13_6, True, 8)
    return str(result)

if __name__ == '__main__':
    print(solve())

Java

import java.math.BigInteger;

public class Euler312 {
    static long pow_u64(long base, int exp) {
        long result = 1;
        while (exp-- > 0)
            result *= base;
        return result;
    }

    static long mul_mod(long a, long b, long mod) {
        return BigInteger.valueOf(a).multiply(BigInteger.valueOf(b)).mod(BigInteger.valueOf(mod)).longValue();
    }

    static long pow_mod(long base, long exp, long mod) {
        if (mod == 1)
            return 0;
        return BigInteger.valueOf(base).modPow(BigInteger.valueOf(exp), BigInteger.valueOf(mod)).longValue();
    }

    static long inverse_mod(long a, long mod) {
        a %= mod;
        long phi = mod - mod / 13;
        return pow_mod(a, phi - 1, mod);
    }

    static long crt_coprime(long a1, long m1, long a2, long m2) {
        long inv_m1 = inverse_mod(m1 % m2, m2);
        long delta = (a2 >= a1 % m2) ? (a2 - (a1 % m2)) : (a2 + m2 - (a1 % m2));
        long t = mul_mod(delta, inv_m1, m2);
        return a1 + m1 * t;
    }

    static long c_mod_13_power_from_n_residue(long n_residue, boolean n_is_large, int k) {
        long mod = pow_u64(13, k);

        if (!n_is_large && n_residue <= 2) {
            return 1 % mod;
        }

        long ord12 = 2 * pow_u64(13, k - 1);
        long two_ord = 2 * ord12;
        long lambda_val = (k >= 2) ? (12 * pow_u64(13, k - 2)) : 2;

        long exp = (n_residue + lambda_val - 2) % lambda_val;
        long x = pow_mod(3, exp, two_ord);

        long numer = (x + two_ord - 3) % two_ord;
        long e = (numer / 2) % ord12;

        return mul_mod(8 % mod, pow_mod(12, e, mod), mod);
    }

    public static String solve() {
        long n0 = 10000;

        long a_mod_13_4 = c_mod_13_power_from_n_residue(n0, false, 4);
        long a_mod_12x13_4 = crt_coprime(0, 12, a_mod_13_4, pow_u64(13, 4));

        long b_mod_13_6 = c_mod_13_power_from_n_residue(a_mod_12x13_4, true, 6);
        long b_mod_12x13_6 = crt_coprime(0, 12, b_mod_13_6, pow_u64(13, 6));

        long result = c_mod_13_power_from_n_residue(b_mod_12x13_6, true, 8);
        return String.valueOf(result);
    }

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