Problem 671: Colouring a Loop
View on Project EulerProject Euler Problem 671 Solution
EulerSolve provides an optimized solution for Project Euler Problem 671, Colouring a Loop, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary The task is to count rotationally distinct valid colourings of a two-row loop of length \(n\) using \(k\) colours, modulo \(M=1000004321\). A direct search over all full loops is hopeless for the target input, so the implementation scans the loop one column at a time and remembers only the local boundary information needed to continue a legal colouring. The admissible local patterns are exactly those encoded by the transition system: a colour region may continue horizontally on the top row or the bottom row for at most three columns, or a column may be a single vertical monochromatic domino. Ordinary columns have different top and bottom colours, while a vertical domino has equal top and bottom colours. A simultaneous restart of both rows is allowed only immediately after a vertical domino, which is how the automaton excludes forbidden four-corner junctions. Mathematical Approach Let \(F_k(n)\) denote the number of rotational classes of valid colourings of length \(n\). The solution has three layers: encode local legality with a finite automaton, turn rooted colourings into matrix traces, and then quotient by rotation with Burnside's lemma....
Detailed mathematical approach
Problem Summary
The task is to count rotationally distinct valid colourings of a two-row loop of length \(n\) using \(k\) colours, modulo \(M=1000004321\). A direct search over all full loops is hopeless for the target input, so the implementation scans the loop one column at a time and remembers only the local boundary information needed to continue a legal colouring.
The admissible local patterns are exactly those encoded by the transition system: a colour region may continue horizontally on the top row or the bottom row for at most three columns, or a column may be a single vertical monochromatic domino. Ordinary columns have different top and bottom colours, while a vertical domino has equal top and bottom colours. A simultaneous restart of both rows is allowed only immediately after a vertical domino, which is how the automaton excludes forbidden four-corner junctions.
Mathematical Approach
Let \(F_k(n)\) denote the number of rotational classes of valid colourings of length \(n\). The solution has three layers: encode local legality with a finite automaton, turn rooted colourings into matrix traces, and then quotient by rotation with Burnside's lemma.
Step 1: Encode the Boundary State
Walk around the loop column by column and record a state
$$\bigl(r_t,r_b,p,c_t,c_b\bigr).$$
Here \(r_t,r_b\in\{0,1,2\}\) are the numbers of future columns still occupied by the current top and bottom horizontal runs after the present column, \(c_t,c_b\) are the current colours, and \(p\in\{0,1\}\) marks whether the current column is a vertical monochromatic domino.
The validity conditions are
$$p=0 \implies c_t\ne c_b,$$
$$p=1 \implies r_t=r_b=0 \land c_t=c_b.$$
So the state space is finite, with
$$S=9k(k-1)+k=k(9k-8)$$
valid states.
Step 2: Build the Transition Matrix
From each state, the next column is determined by simple local continuation rules.
If \(r_t>0\) and \(r_b>0\), both horizontal runs continue, so the next state is obtained by decrementing both remainders. If exactly one row has ended, the unfinished row keeps its colour while the finished row starts a new run of length \(1\), \(2\), or \(3\), with a colour different from both its left neighbour and its vertical neighbour in the new column.
If both rows have ended, a fresh colour may create a vertical domino. Two independent new horizontal runs may also start, but only when the previous column was vertical; this is the local rule that removes four-corner meetings. These legal moves form a directed graph. Let \(A\) be its adjacency matrix.
Step 3: Count Rooted Colourings by Matrix Traces
A rooted colouring of length \(m\) that closes consistently after exactly \(m\) columns corresponds to a closed walk of length \(m\) in the state graph. Therefore the rooted count is
$$T_m=\operatorname{tr}(A^m).$$
This trace counts closed walks because the diagonal entries of \(A^m\) count returns to the same state after \(m\) steps. Since \(A\) is a finite matrix, the sequence \((T_m)\) satisfies a linear recurrence modulo \(M\).
Step 4: Recover the Linear Recurrence
The implementation first computes a prefix
$$T_0,T_1,T_2,\dots$$
by repeatedly multiplying the current matrix power by the sparse transition graph and taking the trace after each step. Once enough terms are available, Berlekamp-Massey reconstructs the shortest recurrence
$$T_{m+L}=c_1T_{m+L-1}+c_2T_{m+L-2}+\cdots+c_LT_m \pmod M.$$
After that, very large indices can be evaluated quickly using binary exponentiation in the recurrence algebra, so the huge target value of \(n\) is reduced to a logarithmic-time query in \(n\).
Step 5: Quotient by Rotation with Burnside's Lemma
The loop is cyclic, so rooted colourings that differ only by a rotation represent the same final object. Burnside's lemma gives
$$F_k(n)=\frac{1}{n}\sum_{s=0}^{n-1}\operatorname{Fix}(s),$$
where \(\operatorname{Fix}(s)\) is the number of colourings fixed by rotation through \(s\) columns. A colouring fixed by that rotation has period \(\gcd(n,s)\), hence
$$\operatorname{Fix}(s)=T_{\gcd(n,s)}.$$
Grouping rotations by their order yields the divisor form used in the program:
$$F_k(n)=\frac{1}{n}\sum_{d\mid n}\varphi(d)\,T_{n/d}\pmod M.$$
Worked Example: The Case \(n=3\)
For a loop of length \(3\), there are exactly three rotations. The identity fixes all rooted colourings of period \(3\), contributing \(T_3\). The other two rotations both force period \(1\), so each contributes \(T_1\). Therefore
$$F_k(3)=\frac{T_3+2T_1}{3}.$$
This is a clean small example of the divisor formula, because \(3\) is prime. For \(k=4\), the implementation checks that
$$F_4(3)=104,$$
which is exactly the expected Burnside average for that case.
How the Code Works
The C++, Python, and Java implementations all follow the same pipeline. First they enumerate every valid boundary state and build the directed adjacency list of legal next states. This adjacency list is the sparse representation of the transfer matrix.
Next they start from the identity matrix, repeatedly apply one more transition step, and record the traces \(T_0,T_1,\dots\). Once the trace sequence is long enough, they recover a minimal linear recurrence modulo \(M\) and keep the initial values needed to seed that recurrence.
Finally they factor \(n\), enumerate all divisors together with their Euler totient weights, evaluate \(T_{n/d}\) from the recurrence for each divisor \(d\), sum
$$\varphi(d)\,T_{n/d},$$
and multiply by the modular inverse of \(n\). The published answer is that Burnside average modulo \(M\).
Complexity Analysis
Let \(S\) be the number of states, \(E\) the number of directed transitions, and \(L\) the recovered recurrence order. Building the automaton takes \(O(S+E)\) time and memory.
If a prefix of \(m\) trace values is generated, the repeated dense-by-sparse matrix updates cost roughly \(O(mSE)\) time and \(O(S^2+E)\) memory, because the current matrix power is stored explicitly while the transition graph is sparse.
Recovering a recurrence from a prefix of length \(O(L)\) is quadratic in \(L\) in the standard Berlekamp-Massey formulation. After that, each large-index query \(T_x\) costs \(O(L^2\log x)\), and the final Burnside sum needs one such query for each divisor of \(n\). Thus the cyclic count after recurrence discovery is
$$O\bigl(\tau(n)L^2\log n\bigr),$$
plus the comparatively small cost of factoring \(n\) and enumerating its divisors.
Footnotes and References
- Problem page: Project Euler 671
- Transfer-matrix method: Wikipedia - Transfer-matrix method
- Burnside's lemma: Wikipedia - Burnside's lemma
- Euler's totient function: Wikipedia - Euler's totient function
- Berlekamp-Massey algorithm: Wikipedia - Berlekamp-Massey algorithm
- Cayley-Hamilton theorem: Wikipedia - Cayley-Hamilton theorem
Problem 671 source code
C++
#include <algorithm>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>
using namespace std;
using ll = long long;
static const ll MOD = 1000004321LL;
struct State {
int rt;
int rb;
int pv;
int ct;
int cb;
};
struct MatrixData {
int k;
int S;
vector<State> states;
vector<int> id;
vector<vector<int>> adj;
};
struct RecurrenceData {
vector<ll> coef;
vector<ll> init;
};
int encode_state(int rt, int rb, int pv, int ct, int cb, int k) {
return ((((rt * 3 + rb) * 2 + pv) * k + ct) * k + cb);
}
MatrixData build_matrix(int k) {
MatrixData md;
md.k = k;
int total = 3 * 3 * 2 * k * k;
md.id.assign(total, -1);
md.states.reserve(total);
for (int rt = 0; rt <= 2; ++rt) {
for (int rb = 0; rb <= 2; ++rb) {
for (int pv = 0; pv <= 1; ++pv) {
for (int ct = 0; ct < k; ++ct) {
for (int cb = 0; cb < k; ++cb) {
bool ok = false;
if (pv == 1) {
ok = (rt == 0 && rb == 0 && ct == cb);
} else {
ok = (ct != cb);
}
if (!ok) continue;
int idx = static_cast<int>(md.states.size());
md.states.push_back({rt, rb, pv, ct, cb});
md.id[encode_state(rt, rb, pv, ct, cb, k)] = idx;
}
}
}
}
}
md.S = static_cast<int>(md.states.size());
md.adj.assign(md.S, {});
auto idx_of = [&](int rt, int rb, int pv, int ct, int cb) {
return md.id[encode_state(rt, rb, pv, ct, cb, k)];
};
for (int i = 0; i < md.S; ++i) {
const State& s = md.states[i];
if (s.rt > 0 && s.rb > 0) {
md.adj[i].push_back(idx_of(s.rt - 1, s.rb - 1, 0, s.ct, s.cb));
continue;
}
if (s.rt > 0 && s.rb == 0) {
for (int lb = 1; lb <= 3; ++lb) {
for (int cb = 0; cb < k; ++cb) {
if (cb == s.ct || cb == s.cb) continue;
md.adj[i].push_back(idx_of(s.rt - 1, lb - 1, 0, s.ct, cb));
}
}
continue;
}
if (s.rt == 0 && s.rb > 0) {
for (int lt = 1; lt <= 3; ++lt) {
for (int ct = 0; ct < k; ++ct) {
if (ct == s.ct || ct == s.cb) continue;
md.adj[i].push_back(idx_of(lt - 1, s.rb - 1, 0, ct, s.cb));
}
}
continue;
}
for (int c = 0; c < k; ++c) {
if (c == s.ct || c == s.cb) continue;
md.adj[i].push_back(idx_of(0, 0, 1, c, c));
}
if (s.pv == 1) {
// Avoid four-corners: separate starts only after a vertical domino.
for (int lt = 1; lt <= 3; ++lt) {
for (int lb = 1; lb <= 3; ++lb) {
for (int ct = 0; ct < k; ++ct) {
if (ct == s.ct) continue;
for (int cb = 0; cb < k; ++cb) {
if (cb == s.cb || cb == ct) continue;
md.adj[i].push_back(idx_of(lt - 1, lb - 1, 0, ct, cb));
}
}
}
}
}
}
return md;
}
using Matrix = vector<vector<uint64_t>>;
int choose_threads(int S) {
unsigned hc = thread::hardware_concurrency();
if (hc <= 1 || S < 200) return 1;
int threads = static_cast<int>(hc);
threads = min(threads, 8);
threads = min(threads, S);
return max(threads, 1);
}
Matrix multiply_dense_sparse(const Matrix& P, const vector<vector<int>>& adj, int threads) {
int S = static_cast<int>(P.size());
Matrix next(S, vector<uint64_t>(S, 0));
auto worker = [&](int start, int end) {
for (int i = start; i < end; ++i) {
const auto& row = P[i];
auto& out = next[i];
for (int k = 0; k < S; ++k) {
uint64_t val = row[k];
if (val == 0) continue;
const auto& edges = adj[k];
for (int to : edges) {
out[to] += val;
}
}
for (int j = 0; j < S; ++j) {
out[j] %= MOD;
}
}
};
if (threads <= 1) {
worker(0, S);
return next;
}
vector<thread> pool;
int block = (S + threads - 1) / threads;
for (int t = 0; t < threads; ++t) {
int start = t * block;
int end = min(S, start + block);
if (start >= end) break;
pool.emplace_back(worker, start, end);
}
for (auto& th : pool) th.join();
return next;
}
ll mod_pow(ll base, ll exp) {
ll res = 1 % MOD;
ll cur = base % MOD;
while (exp > 0) {
if (exp & 1) res = static_cast<ll>((__int128)res * cur % MOD);
cur = static_cast<ll>((__int128)cur * cur % MOD);
exp >>= 1;
}
return res;
}
ll mod_inv(ll x) {
return mod_pow((x % MOD + MOD) % MOD, MOD - 2);
}
vector<ll> berlekamp_massey(const vector<ll>& s) {
vector<ll> C(1, 1), B(1, 1);
ll L = 0;
ll m = 1;
ll b = 1;
for (int n = 0; n < static_cast<int>(s.size()); ++n) {
ll d = 0;
for (int i = 0; i <= L; ++i) {
d = (d + C[i] * s[n - i]) % MOD;
}
if (d == 0) {
++m;
continue;
}
vector<ll> T = C;
ll coef = d * mod_inv(b) % MOD;
if (C.size() < B.size() + m) C.resize(B.size() + m, 0);
for (int i = 0; i < static_cast<int>(B.size()); ++i) {
ll val = (C[i + m] + MOD - coef * B[i] % MOD) % MOD;
C[i + m] = val;
}
if (2 * L <= n) {
L = n + 1 - L;
B = T;
b = d;
m = 1;
} else {
++m;
}
}
C.erase(C.begin());
for (ll& x : C) x = (MOD - x) % MOD;
return C;
}
vector<ll> combine_poly(const vector<ll>& a, const vector<ll>& b, const vector<ll>& coef) {
int L = static_cast<int>(coef.size());
vector<ll> res(2 * L, 0);
for (int i = 0; i < L; ++i) {
if (a[i] == 0) continue;
for (int j = 0; j < L; ++j) {
if (b[j] == 0) continue;
res[i + j] = (res[i + j] + a[i] * b[j]) % MOD;
}
}
for (int i = 2 * L - 2; i >= L; --i) {
if (res[i] == 0) continue;
ll val = res[i];
for (int j = 1; j <= L; ++j) {
res[i - j] = (res[i - j] + val * coef[j - 1]) % MOD;
}
}
res.resize(L);
return res;
}
ll linear_rec(const vector<ll>& init, const vector<ll>& coef, long long n) {
int L = static_cast<int>(coef.size());
if (L == 0) return 0;
if (n < static_cast<long long>(init.size())) return init[n];
vector<ll> pol(L, 0), e(L, 0);
pol[0] = 1;
if (L == 1) {
e[0] = coef[0];
} else {
e[1] = 1;
}
long long m = n;
while (m > 0) {
if (m & 1) pol = combine_poly(pol, e, coef);
e = combine_poly(e, e, coef);
m >>= 1;
}
ll res = 0;
for (int i = 0; i < L; ++i) {
res = (res + pol[i] * init[i]) % MOD;
}
return res;
}
RecurrenceData build_recurrence(int k) {
MatrixData md = build_matrix(k);
int S = md.S;
int threads = choose_threads(S);
Matrix P(S, vector<uint64_t>(S, 0));
for (int i = 0; i < S; ++i) P[i][i] = 1;
vector<ll> seq;
seq.reserve(512);
seq.push_back(S % MOD);
vector<ll> coef;
int last_L = 0;
int stable = 0;
int max_steps = 2000;
for (int step = 1; step <= max_steps; ++step) {
P = multiply_dense_sparse(P, md.adj, threads);
ll tr = 0;
for (int i = 0; i < S; ++i) {
tr += static_cast<ll>(P[i][i]);
if (tr >= MOD) tr %= MOD;
}
seq.push_back(tr % MOD);
coef = berlekamp_massey(seq);
int L = static_cast<int>(coef.size());
if (L == last_L) {
++stable;
} else {
stable = 0;
last_L = L;
}
if (L > 0 && static_cast<int>(seq.size()) >= 2 * L + 5 && stable >= 5) {
break;
}
}
if (coef.empty()) {
coef.push_back(0);
}
vector<ll> init(coef.size(), 0);
for (size_t i = 0; i < init.size() && i < seq.size(); ++i) {
init[i] = seq[i];
}
return {coef, init};
}
vector<pair<ll, int>> factorize(ll n) {
vector<pair<ll, int>> factors;
for (ll p = 2; p * p <= n; p += (p == 2 ? 1 : 2)) {
if (n % p != 0) continue;
int cnt = 0;
while (n % p == 0) {
n /= p;
++cnt;
}
factors.push_back({p, cnt});
}
if (n > 1) factors.push_back({n, 1});
return factors;
}
void gen_divisors(int idx, ll cur_d, ll cur_phi,
const vector<pair<ll, int>>& factors,
vector<pair<ll, ll>>& out) {
if (idx == static_cast<int>(factors.size())) {
out.push_back({cur_d, cur_phi % MOD});
return;
}
ll p = factors[idx].first;
int e = factors[idx].second;
gen_divisors(idx + 1, cur_d, cur_phi, factors, out);
ll p_pow = 1;
ll phi_pow = 1;
for (int i = 1; i <= e; ++i) {
p_pow *= p;
if (i == 1) {
phi_pow = p - 1;
} else {
phi_pow *= p;
}
gen_divisors(idx + 1, cur_d * p_pow, cur_phi * phi_pow, factors, out);
}
}
ll solve_problem(int k, ll n) {
RecurrenceData rec = build_recurrence(k);
auto factors = factorize(n);
vector<pair<ll, ll>> divisors;
gen_divisors(0, 1, 1, factors, divisors);
ll total = 0;
for (const auto& dv : divisors) {
ll d = dv.first;
ll phi_mod = dv.second % MOD;
ll term = linear_rec(rec.init, rec.coef, n / d);
total = (total + phi_mod * term) % MOD;
}
ll inv_n = mod_inv(n % MOD);
return total * inv_n % MOD;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
if (solve_problem(4, 3) != 104) {
cerr << "Validation failed: F4(3) != 104\n";
return 1;
}
if (solve_problem(5, 7) != 3327300) {
cerr << "Validation failed: F5(7) != 3327300\n";
return 1;
}
if (solve_problem(6, 101) != 75309980) {
cerr << "Validation failed: F6(101) != 75309980\n";
return 1;
}
const ll n = 10004003002001LL;
cout << solve_problem(10, n) << "\n";
return 0;
}
Python
def solve():
MOD = 1000004321
def pm(b,e):
r=1; b%=MOD
while e>0:
if e&1: r=r*b%MOD
b=b*b%MOD; e>>=1
return r
def mi(x): return pm((x%MOD+MOD)%MOD,MOD-2)
def build_matrix(k):
states=[]; idx={}
for rt in range(3):
for rb in range(3):
for pv in range(2):
for ct in range(k):
for cb in range(k):
ok=False
if pv==1: ok=(rt==0 and rb==0 and ct==cb)
else: ok=(ct!=cb)
if not ok: continue
i=len(states); states.append((rt,rb,pv,ct,cb))
idx[(rt,rb,pv,ct,cb)]=i
S=len(states); adj=[[] for _ in range(S)]
for i,(rt,rb,pv,ct,cb) in enumerate(states):
if rt>0 and rb>0: adj[i].append(idx[(rt-1,rb-1,0,ct,cb)]); continue
if rt>0 and rb==0:
for lb in range(1,4):
for cb2 in range(k):
if cb2==ct or cb2==cb: continue
adj[i].append(idx[(rt-1,lb-1,0,ct,cb2)])
continue
if rt==0 and rb>0:
for lt in range(1,4):
for ct2 in range(k):
if ct2==ct or ct2==cb: continue
adj[i].append(idx[(lt-1,rb-1,0,ct2,cb)])
continue
for c in range(k):
if c==ct or c==cb: continue
adj[i].append(idx[(0,0,1,c,c)])
if pv==1:
for lt in range(1,4):
for lb in range(1,4):
for ct2 in range(k):
if ct2==ct: continue
for cb2 in range(k):
if cb2==cb or cb2==ct2: continue
adj[i].append(idx[(lt-1,lb-1,0,ct2,cb2)])
return S,adj
def bm(s):
C=[1]; B=[1]; L=0; m=1; b=1
for n in range(len(s)):
d=0
for i in range(L+1): d=(d+C[i]*s[n-i])%MOD
if d==0: m+=1; continue
T=C[:]; coef=d*mi(b)%MOD
while len(C)<len(B)+m: C.append(0)
for i in range(len(B)): C[i+m]=(C[i+m]+MOD-coef*B[i]%MOD)%MOD
if 2*L<=n: L=n+1-L; B=T; b=d; m=1
else: m+=1
C.pop(0); return [(MOD-x)%MOD for x in C]
def lr(init,coef,n):
L=len(coef)
if L==0: return 0
if n<len(init): return init[n]
def cp(a,b):
r=[0]*(2*L)
for i in range(L):
if a[i]==0: continue
for j in range(L):
if b[j]==0: continue
r[i+j]=(r[i+j]+a[i]*b[j])%MOD
for i in range(2*L-2,L-1,-1):
if r[i]==0: continue
for j in range(L): r[i-1-j]=(r[i-1-j]+r[i]*coef[j])%MOD
return r[:L]
pol=[0]*L; pol[0]=1; e=[0]*L
if L==1: e[0]=coef[0]
else: e[1]=1
m=n
while m>0:
if m&1: pol=cp(pol,e)
e=cp(e,e); m>>=1
return sum(pol[i]*init[i] for i in range(L))%MOD
def build_rec(k):
S,adj=build_matrix(k)
# Matrix power via sparse mult, trace sequence
P=[[0]*S for _ in range(S)]
for i in range(S): P[i][i]=1
seq=[S%MOD]; last_L=0; stable=0
for step in range(1,2001):
nP=[[0]*S for _ in range(S)]
for i in range(S):
for kk in range(S):
if P[i][kk]==0: continue
for to in adj[kk]: nP[i][to]=(nP[i][to]+P[i][kk])%MOD
P=nP; tr=sum(P[i][i] for i in range(S))%MOD; seq.append(tr)
coef=bm(seq); L=len(coef)
if L==last_L: stable+=1
else: stable=0; last_L=L
if L>0 and len(seq)>=2*L+5 and stable>=5: break
if not coef: coef=[0]
init=seq[:len(coef)]
return init,coef
def factorize(n):
f=[]; p=2
while p*p<=n:
if n%p==0:
c=0
while n%p==0: n//=p; c+=1
f.append((p,c))
p+=1 if p==2 else 2
if n>1: f.append((n,1))
return f
def gen_divs(facs):
divs=[(1,1)]
for p,e in facs:
nd=[]
for d,phi in divs:
pp=1; ph=1
nd.append((d,phi))
for i in range(1,e+1):
pp*=p; ph=p-1 if i==1 else ph*p
nd.append((d*pp,phi*ph))
divs=nd
return divs
k=10; n=10004003002001
init,coef=build_rec(k)
facs=factorize(n)
divs=gen_divs(facs)
total=0
for d,phi in divs:
total=(total+phi%MOD*lr(init,coef,n//d))%MOD
return str(total*mi(n%MOD)%MOD)
if __name__=='__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
public class Euler671 {
static final long MOD = 1000004321L;
static class State {
int rt, rb, pv, ct, cb;
State(int rt, int rb, int pv, int ct, int cb) {
this.rt = rt;
this.rb = rb;
this.pv = pv;
this.ct = ct;
this.cb = cb;
}
}
static int encodeState(int rt, int rb, int pv, int ct, int cb, int k) {
return ((((rt * 3 + rb) * 2 + pv) * k + ct) * k + cb);
}
static class MatrixData {
int S;
List<List<Integer>> adj;
}
static MatrixData buildMatrix(int k) {
MatrixData md = new MatrixData();
List<State> states = new ArrayList<>();
HashMap<Integer, Integer> idMap = new HashMap<>();
for (int rt = 0; rt <= 2; ++rt) {
for (int rb = 0; rb <= 2; ++rb) {
for (int pv = 0; pv <= 1; ++pv) {
for (int ct = 0; ct < k; ++ct) {
for (int cb = 0; cb < k; ++cb) {
boolean ok;
if (pv == 1) {
ok = (rt == 0 && rb == 0 && ct == cb);
} else {
ok = (ct != cb);
}
if (!ok)
continue;
int idx = states.size();
states.add(new State(rt, rb, pv, ct, cb));
idMap.put(encodeState(rt, rb, pv, ct, cb, k), idx);
}
}
}
}
}
md.S = states.size();
md.adj = new ArrayList<>();
for (int i = 0; i < md.S; ++i)
md.adj.add(new ArrayList<>());
for (int i = 0; i < md.S; ++i) {
State s = states.get(i);
if (s.rt > 0 && s.rb > 0) {
md.adj.get(i).add(idMap.get(encodeState(s.rt - 1, s.rb - 1, 0, s.ct, s.cb, k)));
continue;
}
if (s.rt > 0 && s.rb == 0) {
for (int lb = 1; lb <= 3; ++lb) {
for (int cb = 0; cb < k; ++cb) {
if (cb == s.ct || cb == s.cb)
continue;
md.adj.get(i).add(idMap.get(encodeState(s.rt - 1, lb - 1, 0, s.ct, cb, k)));
}
}
continue;
}
if (s.rt == 0 && s.rb > 0) {
for (int lt = 1; lt <= 3; ++lt) {
for (int ct = 0; ct < k; ++ct) {
if (ct == s.ct || ct == s.cb)
continue;
md.adj.get(i).add(idMap.get(encodeState(lt - 1, s.rb - 1, 0, ct, s.cb, k)));
}
}
continue;
}
for (int c = 0; c < k; ++c) {
if (c == s.ct || c == s.cb)
continue;
md.adj.get(i).add(idMap.get(encodeState(0, 0, 1, c, c, k)));
}
if (s.pv == 1) {
for (int lt = 1; lt <= 3; ++lt) {
for (int lb = 1; lb <= 3; ++lb) {
for (int ct = 0; ct < k; ++ct) {
if (ct == s.ct)
continue;
for (int cb = 0; cb < k; ++cb) {
if (cb == s.cb || cb == ct)
continue;
md.adj.get(i).add(idMap.get(encodeState(lt - 1, lb - 1, 0, ct, cb, k)));
}
}
}
}
}
}
return md;
}
static long modPow(long base, long exp) {
long res = 1;
long cur = base % MOD;
while (exp > 0) {
if ((exp & 1) != 0)
res = (res * cur) % MOD;
cur = (cur * cur) % MOD;
exp >>= 1;
}
return res;
}
static long modInv(long x) {
return modPow((x % MOD + MOD) % MOD, MOD - 2);
}
static List<Long> berlekampMassey(List<Long> s) {
List<Long> C = new ArrayList<>();
C.add(1L);
List<Long> B = new ArrayList<>();
B.add(1L);
int L = 0, m = 1;
long b = 1;
for (int n = 0; n < s.size(); ++n) {
long d = 0;
for (int i = 0; i <= L; ++i) {
d = (d + C.get(i) * s.get(n - i)) % MOD;
}
if (d == 0) {
m++;
continue;
}
List<Long> T = new ArrayList<>(C);
long coef = (d * modInv(b)) % MOD;
while (C.size() < B.size() + m)
C.add(0L);
for (int i = 0; i < B.size(); ++i) {
long val = (C.get(i + m) + MOD - (coef * B.get(i)) % MOD) % MOD;
C.set(i + m, val);
}
if (2 * L <= n) {
L = n + 1 - L;
B = T;
b = d;
m = 1;
} else {
m++;
}
}
C.remove(0);
for (int i = 0; i < C.size(); ++i)
C.set(i, (MOD - C.get(i)) % MOD);
return C;
}
static long[] combinePoly(long[] a, long[] b, long[] coef) {
int L = coef.length;
long[] res = new long[2 * L];
for (int i = 0; i < L; ++i) {
if (a[i] == 0)
continue;
for (int j = 0; j < L; ++j) {
if (b[j] == 0)
continue;
res[i + j] = (res[i + j] + a[i] * b[j]) % MOD;
}
}
for (int i = 2 * L - 2; i >= L; --i) {
if (res[i] == 0)
continue;
long val = res[i];
for (int j = 1; j <= L; ++j) {
res[i - j] = (res[i - j] + val * coef[j - 1]) % MOD;
}
}
long[] out = new long[L];
System.arraycopy(res, 0, out, 0, L);
return out;
}
static long linearRec(List<Long> init, List<Long> coefList, long n) {
int L = coefList.size();
if (L == 0)
return 0;
if (n < init.size())
return init.get((int) n);
long[] coef = new long[L];
for (int i = 0; i < L; ++i)
coef[i] = coefList.get(i);
long[] pol = new long[L];
pol[0] = 1;
long[] e = new long[L];
if (L == 1) {
e[0] = coef[0];
} else {
e[1] = 1;
}
long m = n;
while (m > 0) {
if ((m & 1) != 0)
pol = combinePoly(pol, e, coef);
e = combinePoly(e, e, coef);
m >>= 1;
}
long res = 0;
for (int i = 0; i < L; ++i) {
res = (res + pol[i] * init.get(i)) % MOD;
}
return res;
}
static class Factor {
long p;
int e;
Factor(long p, int e) {
this.p = p;
this.e = e;
}
}
static List<Factor> factorize(long n) {
List<Factor> res = new ArrayList<>();
for (long p = 2; p * p <= n; p += (p == 2 ? 1 : 2)) {
if (n % p == 0) {
int cnt = 0;
while (n % p == 0) {
n /= p;
cnt++;
}
res.add(new Factor(p, cnt));
}
}
if (n > 1)
res.add(new Factor(n, 1));
return res;
}
static class Divisor {
long d, phi;
Divisor(long d, long phi) {
this.d = d;
this.phi = phi;
}
}
static List<Divisor> genDivisors(List<Factor> factors) {
List<Divisor> divs = new ArrayList<>();
divs.add(new Divisor(1, 1));
for (Factor f : factors) {
List<Divisor> nextDivs = new ArrayList<>();
for (Divisor dv : divs) {
nextDivs.add(dv);
long pPow = 1;
long phiPow = 1;
for (int i = 1; i <= f.e; ++i) {
pPow *= f.p;
phiPow = (i == 1) ? (f.p - 1) : (phiPow * f.p);
nextDivs.add(new Divisor(dv.d * pPow, dv.phi * phiPow));
}
}
divs = nextDivs;
}
return divs;
}
static long solveProblem(int k, long n) {
MatrixData md = buildMatrix(k);
int S = md.S;
long[][] P = new long[S][S];
for (int i = 0; i < S; ++i)
P[i][i] = 1;
long[][] nextP = new long[S][S];
List<Long> seq = new ArrayList<>();
seq.add((long) S % MOD);
int lastL = 0;
int stable = 0;
List<Long> coef = new ArrayList<>();
for (int step = 1; step <= 2000; ++step) {
for (int i = 0; i < S; ++i) {
for (int j = 0; j < S; ++j)
nextP[i][j] = 0;
for (int kIdx = 0; kIdx < S; ++kIdx) {
long val = P[i][kIdx];
if (val > 0) {
for (int to : md.adj.get(kIdx)) {
nextP[i][to] += val;
}
}
}
for (int j = 0; j < S; ++j)
nextP[i][j] %= MOD;
}
long[][] tmp = P;
P = nextP;
nextP = tmp;
long tr = 0;
for (int i = 0; i < S; ++i)
tr = (tr + P[i][i]) % MOD;
seq.add(tr);
coef = berlekampMassey(seq);
int L = coef.size();
if (L == lastL)
stable++;
else {
stable = 0;
lastL = L;
}
if (L > 0 && seq.size() >= 2 * L + 5 && stable >= 5)
break;
}
if (coef.isEmpty())
coef.add(0L);
List<Long> init = new ArrayList<>();
for (int i = 0; i < coef.size(); ++i)
init.add(seq.get(i));
List<Divisor> divs = genDivisors(factorize(n));
long total = 0;
for (Divisor dv : divs) {
long term = linearRec(init, coef, n / dv.d);
total = (total + (dv.phi % MOD) * term) % MOD;
}
long invN = modInv(n % MOD);
return (total * invN) % MOD;
}
public static String solve() {
long n = 10004003002001L;
return Long.toString(solveProblem(10, n));
}
public static void main(String[] args) {
System.out.println(solve());
}
}