Problem 670: Colouring a Strip

Project Euler Problem 670 Solution

EulerSolve provides an optimized solution for Project Euler Problem 670, Colouring a Strip, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $A(n)$ denote the number of valid colourings of a $2\times n$ strip using four colours. The implementations work modulo $$M=10^6+4321=1{,}000{,}004{,}321.$$ The key difficulty is that the legality of a new column is not determined only by the colours visible in that column. The algorithm therefore converts the strip into a finite-state process that remembers just enough local information to decide every future extension. Mathematical Approach The counting problem is turned into a transfer-matrix problem. Each state describes the current colours on the top and bottom rows, how many more columns the current horizontal segments must continue, and whether the previous move was a special monochromatic vertical column. Step 1: Encode the Necessary Local Memory After processing a column, the algorithm stores a state $$s=(r_{\text{top}},c_{\text{top}},r_{\text{bot}},c_{\text{bot}},\nu),$$ where $c_{\text{top}},c_{\text{bot}}\in\{1,2,3,4\}$ are the current colours, $r_{\text{top}},r_{\text{bot}}\in\{0,1,2\}$ are the remaining numbers of columns that the current top and bottom horizontal segments must still occupy, and $\nu\in\{0,1\}$ is a one-bit memory flag. If a segment has total length $1$, its stored remainder is $0$; if its total length is $2$, the stored remainder is $1$; if its total length is $3$, the stored remainder is $2$....

Detailed mathematical approach

Problem Summary

Let $A(n)$ denote the number of valid colourings of a $2\times n$ strip using four colours. The implementations work modulo

$$M=10^6+4321=1{,}000{,}004{,}321.$$

The key difficulty is that the legality of a new column is not determined only by the colours visible in that column. The algorithm therefore converts the strip into a finite-state process that remembers just enough local information to decide every future extension.

Mathematical Approach

The counting problem is turned into a transfer-matrix problem. Each state describes the current colours on the top and bottom rows, how many more columns the current horizontal segments must continue, and whether the previous move was a special monochromatic vertical column.

Step 1: Encode the Necessary Local Memory

After processing a column, the algorithm stores a state

$$s=(r_{\text{top}},c_{\text{top}},r_{\text{bot}},c_{\text{bot}},\nu),$$

where $c_{\text{top}},c_{\text{bot}}\in\{1,2,3,4\}$ are the current colours, $r_{\text{top}},r_{\text{bot}}\in\{0,1,2\}$ are the remaining numbers of columns that the current top and bottom horizontal segments must still occupy, and $\nu\in\{0,1\}$ is a one-bit memory flag.

If a segment has total length $1$, its stored remainder is $0$; if its total length is $2$, the stored remainder is $1$; if its total length is $3$, the stored remainder is $2$. Thus only segment lengths $1,2,3$ can occur.

The raw state space is bounded by

$$3\cdot 4\cdot 3\cdot 4\cdot 2=288,$$

and a reachability search later removes the impossible ones.

Step 2: Describe the First Column

The first column can begin in two qualitatively different ways.

First, it may be a monochromatic vertical column. In that case top and bottom have the same colour, both stored remainders are $0$, and the flag is set to $\nu=1$. There are exactly $4$ such initial states.

Second, the first column may use different top and bottom colours. Then one chooses ordered colours $c_{\text{top}}\neq c_{\text{bot}}$ and independent segment lengths $\ell_{\text{top}},\ell_{\text{bot}}\in\{1,2,3\}$. The stored remainders are $\ell_{\text{top}}-1$ and $\ell_{\text{bot}}-1$, and the flag is $\nu=0$.

This gives the initial row vector $v_1$, whose entries count all admissible descriptions of a strip of length $1$.

Step 3: Generate All Legal One-Column Extensions

Suppose the current state is $s=(r_{\text{top}},c_{\text{top}},r_{\text{bot}},c_{\text{bot}},\nu)$.

If $r_{\text{top}}>0$, then the next top cell is forced to keep colour $c_{\text{top}}$, and the new remainder becomes $r_{\text{top}}-1$. If $r_{\text{top}}=0$, then a new top segment begins: choose a new colour different from $c_{\text{top}}$, choose a new length $\ell\in\{1,2,3\}$, and store the remainder $\ell-1$. The bottom row behaves symmetrically.

After those local choices are made, the next column must satisfy two filters.

First, an ordinary next column must use different colours on the two rows:

$$c'_{\text{top}}\neq c'_{\text{bot}}.$$

Second, if both current remainders are $0$ and the flag is $\nu=0$, then both rows are not allowed to restart simultaneously into an ordinary split column. In other words, a synchronized restart of both rows is permitted only immediately after the special vertical situation.

There is also one extra branch. When both remainders are $0$, the next column may be a monochromatic vertical column of colour $v$, provided

$$v\neq c_{\text{top}},\qquad v\neq c_{\text{bot}}.$$

This moves the automaton to the state $(0,v,0,v,1)$.

Step 4: Build the Reachable Transfer System

Starting from all initial states, a breadth-first search enumerates the reachable state set $\mathcal S$. This is an exact pruning step: unreachable raw states never contribute and are discarded before any matrix algebra is done.

Now define the transfer matrix $T$ by

$$T_{ij}=\#\{\text{legal one-column extensions from state }s_i\text{ to state }s_j\} \pmod M.$$

For this problem the entries are effectively $0$ or $1$, but writing them as counts keeps the formulation clean. If $v_n$ is the row vector of counts after $n$ columns, then

$$v_n=v_1T^{\,n-1} \pmod M.$$

Step 5: Identify the Terminal States

Not every state after the $n$-th column corresponds to a complete strip of length $n$. If a remainder is still positive, then some horizontal segment would need extra columns beyond the end of the strip.

Therefore the valid terminal set is

$$\mathcal T=\{\,s\in\mathcal S: r_{\text{top}}=0 \text{ and } r_{\text{bot}}=0\,\}.$$

The final answer is

$$A(n)=\sum_{s\in\mathcal T} v_n(s)\pmod M.$$

Because the target value uses $n=10^{16}$, the power $T^{n-1}$ is computed by binary exponentiation.

Worked Example: Why $A(2)=120$

The implementations include the checkpoint $A(2)=120$. This can be derived directly from the state model.

If the first column is monochromatic, there are $4$ choices. The second column can then be another monochromatic vertical column in one of the $3$ other colours, or an ordinary split column with two different colours chosen from the remaining $3$ colours, which gives $3\cdot 2=6$ choices. This contributes

$$4\cdot(3+6)=36.$$

If the first column is split with different top and bottom colours, there are $12$ ordered colour pairs. Now inspect the possible initial segment lengths.

If the lengths are $(1,1)$, both rows end immediately. Since the flag is then $0$, the second column cannot restart both rows into another split column, so it must be a monochromatic vertical column using one of the two unused colours. Contribution:

$$12\cdot 2=24.$$

If the lengths are $(2,1)$, the top row continues while the bottom row restarts, and the restarted bottom colour must differ from both the continuing top colour and the old bottom colour. That gives $2$ choices, hence another $24$. By symmetry, the pattern $(1,2)$ contributes another $24$.

If the lengths are $(2,2)$, both rows simply continue and the strip ends exactly after the second column, so the contribution is $12$.

Any pattern involving a segment of length $3$ cannot finish within two columns, so it contributes $0$. Altogether,

$$36+24+24+24+12=120.$$

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. They first enumerate all admissible first-column states, then run a reachability search to collect and sort the finite state set. After that, they build the transition matrix by explicitly generating every legal successor of every reachable state.

Once the matrix is ready, the implementation keeps a count vector for the current strip length and raises the matrix to the needed power with exponentiation by squaring. Whenever a binary digit of $n-1$ is $1$, the vector is multiplied by the current matrix power. At the end, only states with both remainders equal to $0$ are summed. One implementation also includes small checkpoints such as $A(2)=120$, $A(5)=45876$, and $A(100)=53275818$.

Complexity Analysis

Let $S=|\mathcal S|$. The reachability phase is $O(S)$ up to a small constant branching factor, and the raw state space never exceeds $288$ states. The dominant cost is matrix exponentiation: dense matrix multiplication is $O(S^3)$, so the full power computation costs $O(S^3\log n)$. The vector-matrix multiplications add $O(S^2\log n)$, and the memory usage is $O(S^2)$ for the matrix plus $O(S)$ for the state lists and vectors. The implementations skip zero entries in the inner loops, which improves constants without changing the asymptotic bound.

Footnotes and References

Problem page: Project Euler 670
Transfer-matrix method: Wikipedia - Transfer-matrix method
Finite-state machine: Wikipedia - Finite-state machine
Breadth-first search: Wikipedia - Breadth-first search
Exponentiation by squaring: Wikipedia - Exponentiation by squaring

Problem 670 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <queue>
#include <unordered_map>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u64 MOD = 1'000'004'321ULL;
constexpr int COLORS = 4;

struct State {
    int rt;
    int ct;
    int rb;
    int cb;
    int prev_vert;

    bool operator==(const State& other) const {
        return rt == other.rt && ct == other.ct && rb == other.rb && cb == other.cb &&
               prev_vert == other.prev_vert;
    }

    bool operator<(const State& other) const {
        if (rt != other.rt) return rt < other.rt;
        if (ct != other.ct) return ct < other.ct;
        if (rb != other.rb) return rb < other.rb;
        if (cb != other.cb) return cb < other.cb;
        return prev_vert < other.prev_vert;
    }
};

struct StateHash {
    std::size_t operator()(const State& s) const {
        std::size_t h = static_cast<std::size_t>(s.rt);
        h = h * 1315423911u + static_cast<std::size_t>(s.ct + 7);
        h = h * 1315423911u + static_cast<std::size_t>(s.rb + 11);
        h = h * 1315423911u + static_cast<std::size_t>(s.cb + 13);
        h = h * 1315423911u + static_cast<std::size_t>(s.prev_vert + 17);
        return h;
    }
};

std::vector<State> initial_states() {
    std::vector<State> init;
    init.reserve(4 + 3 * 3 * 4 * 3);

    for (int v = 0; v < COLORS; ++v) {
        init.push_back({0, v, 0, v, 1});
    }

    for (int lt = 1; lt <= 3; ++lt) {
        for (int lb = 1; lb <= 3; ++lb) {
            for (int a = 0; a < COLORS; ++a) {
                for (int b = 0; b < COLORS; ++b) {
                    if (a == b) {
                        continue;
                    }
                    init.push_back({lt - 1, a, lb - 1, b, 0});
                }
            }
        }
    }

    return init;
}

std::vector<State> next_states(const State& s) {
    std::vector<State> out;

    const bool top_change = (s.rt == 0);
    const bool bottom_change = (s.rb == 0);

    if (s.rt == 0 && s.rb == 0) {
        for (int v = 0; v < COLORS; ++v) {
            if (v == s.ct || v == s.cb) {
                continue;
            }
            out.push_back({0, v, 0, v, 1});
        }
    }

    std::vector<std::pair<int, int>> top_options;
    std::vector<std::pair<int, int>> bottom_options;

    if (s.rt > 0) {
        top_options.push_back({s.rt - 1, s.ct});
    } else {
        for (int len = 1; len <= 3; ++len) {
            for (int color = 0; color < COLORS; ++color) {
                if (color == s.ct) {
                    continue;
                }
                top_options.push_back({len - 1, color});
            }
        }
    }

    if (s.rb > 0) {
        bottom_options.push_back({s.rb - 1, s.cb});
    } else {
        for (int len = 1; len <= 3; ++len) {
            for (int color = 0; color < COLORS; ++color) {
                if (color == s.cb) {
                    continue;
                }
                bottom_options.push_back({len - 1, color});
            }
        }
    }

    for (const auto& top : top_options) {
        for (const auto& bottom : bottom_options) {
            const int nt_rt = top.first;
            const int nt_ct = top.second;
            const int nt_rb = bottom.first;
            const int nt_cb = bottom.second;

            if (nt_ct == nt_cb) {
                continue;
            }

            if (top_change && bottom_change && s.prev_vert == 0) {
                continue;
            }

            out.push_back({nt_rt, nt_ct, nt_rb, nt_cb, 0});
        }
    }

    return out;
}

struct Matrix {
    int n;
    std::vector<u64> a;

    explicit Matrix(int n_) : n(n_), a(static_cast<std::size_t>(n_) * static_cast<std::size_t>(n_), 0ULL) {}

    u64& at(int i, int j) { return a[static_cast<std::size_t>(i) * static_cast<std::size_t>(n) + static_cast<std::size_t>(j)]; }
    u64 at(int i, int j) const {
        return a[static_cast<std::size_t>(i) * static_cast<std::size_t>(n) + static_cast<std::size_t>(j)];
    }
};

Matrix multiply(const Matrix& x, const Matrix& y) {
    const int n = x.n;
    Matrix z(n);

    for (int i = 0; i < n; ++i) {
        for (int k = 0; k < n; ++k) {
            const u64 xik = x.at(i, k);
            if (xik == 0ULL) {
                continue;
            }
            for (int j = 0; j < n; ++j) {
                const u64 ykj = y.at(k, j);
                if (ykj == 0ULL) {
                    continue;
                }
                u64& cell = z.at(i, j);
                cell = static_cast<u64>((static_cast<u128>(cell) + static_cast<u128>(xik) * static_cast<u128>(ykj)) % MOD);
            }
        }
    }

    return z;
}

std::vector<u64> multiply(const std::vector<u64>& v, const Matrix& m) {
    const int n = m.n;
    std::vector<u64> out(static_cast<std::size_t>(n), 0ULL);

    for (int i = 0; i < n; ++i) {
        const u64 vi = v[static_cast<std::size_t>(i)];
        if (vi == 0ULL) {
            continue;
        }
        for (int j = 0; j < n; ++j) {
            const u64 mij = m.at(i, j);
            if (mij == 0ULL) {
                continue;
            }
            u64& cell = out[static_cast<std::size_t>(j)];
            cell = static_cast<u64>((static_cast<u128>(cell) + static_cast<u128>(vi) * static_cast<u128>(mij)) % MOD);
        }
    }

    return out;
}

class Solver {
public:
    Solver() {
        init_states_ = initial_states();
        build_reachable();
        build_transition();
        build_initial_vector();
    }

    u64 count(u64 n) const {
        assert(n >= 1ULL);

        std::vector<u64> vec = init_vector_;
        Matrix p = transition_;
        u64 exp = n - 1ULL;

        while (exp > 0ULL) {
            if ((exp & 1ULL) != 0ULL) {
                vec = multiply(vec, p);
            }
            exp >>= 1ULL;
            if (exp > 0ULL) {
                p = multiply(p, p);
            }
        }

        u64 ans = 0ULL;
        for (int i = 0; i < static_cast<int>(states_.size()); ++i) {
            const State& s = states_[static_cast<std::size_t>(i)];
            if (s.rt == 0 && s.rb == 0) {
                ans += vec[static_cast<std::size_t>(i)];
                if (ans >= MOD) {
                    ans -= MOD;
                }
            }
        }
        return ans;
    }

private:
    std::vector<State> init_states_;
    std::vector<State> states_;
    std::unordered_map<State, int, StateHash> index_;
    Matrix transition_{1};
    std::vector<u64> init_vector_;

    void build_reachable() {
        std::unordered_map<State, int, StateHash> seen;
        std::queue<State> q;

        for (const State& s : init_states_) {
            if (seen.emplace(s, 1).second) {
                q.push(s);
            }
        }

        while (!q.empty()) {
            const State s = q.front();
            q.pop();

            const std::vector<State> ns = next_states(s);
            for (const State& t : ns) {
                if (seen.emplace(t, 1).second) {
                    q.push(t);
                }
            }
        }

        states_.reserve(seen.size());
        for (const auto& kv : seen) {
            states_.push_back(kv.first);
        }
        std::sort(states_.begin(), states_.end());

        index_.reserve(states_.size() * 2);
        for (int i = 0; i < static_cast<int>(states_.size()); ++i) {
            index_[states_[static_cast<std::size_t>(i)]] = i;
        }
    }

    void build_transition() {
        const int n = static_cast<int>(states_.size());
        transition_ = Matrix(n);

        for (int i = 0; i < n; ++i) {
            const std::vector<State> ns = next_states(states_[static_cast<std::size_t>(i)]);
            for (const State& t : ns) {
                const int j = index_.at(t);
                u64& cell = transition_.at(i, j);
                ++cell;
                if (cell >= MOD) {
                    cell -= MOD;
                }
            }
        }
    }

    void build_initial_vector() {
        init_vector_.assign(states_.size(), 0ULL);
        for (const State& s : init_states_) {
            const int i = index_.at(s);
            u64& cell = init_vector_[static_cast<std::size_t>(i)];
            ++cell;
            if (cell >= MOD) {
                cell -= MOD;
            }
        }
    }
};

}  // namespace

int main() {
    Solver solver;

    assert(solver.count(2ULL) == 120ULL);
    assert(solver.count(5ULL) == 45'876ULL);
    assert(solver.count(100ULL) == 53'275'818ULL);

    std::cout << solver.count(10'000'000'000'000'000ULL) << "\n";
    return 0;
}

Python

MOD = 1000004321
COLORS = 4

class State:
    def __init__(self, rt, ct, rb, cb, prev_vert):
        self.rt = rt
        self.ct = ct
        self.rb = rb
        self.cb = cb
        self.prev_vert = prev_vert

    def __hash__(self):
        return hash((self.rt, self.ct, self.rb, self.cb, self.prev_vert))

    def __eq__(self, other):
        return (self.rt == other.rt and self.ct == other.ct and 
                self.rb == other.rb and self.cb == other.cb and 
                self.prev_vert == other.prev_vert)

    def __lt__(self, other):
        return (self.rt, self.ct, self.rb, self.cb, self.prev_vert) < (other.rt, other.ct, other.rb, other.cb, other.prev_vert)

def initial_states():
    init = []
    for v in range(COLORS):
        init.append(State(0, v, 0, v, 1))

    for lt in range(1, 4):
        for lb in range(1, 4):
            for a in range(COLORS):
                for b in range(COLORS):
                    if a == b:
                        continue
                    init.append(State(lt - 1, a, lb - 1, b, 0))
    return init

def next_states(s):
    out = []
    top_change = (s.rt == 0)
    bottom_change = (s.rb == 0)

    if s.rt == 0 and s.rb == 0:
        for v in range(COLORS):
            if v == s.ct or v == s.cb:
                continue
            out.append(State(0, v, 0, v, 1))

    top_options = []
    bottom_options = []

    if s.rt > 0:
        top_options.append((s.rt - 1, s.ct))
    else:
        for length in range(1, 4):
            for color in range(COLORS):
                if color == s.ct:
                    continue
                top_options.append((length - 1, color))

    if s.rb > 0:
        bottom_options.append((s.rb - 1, s.cb))
    else:
        for length in range(1, 4):
            for color in range(COLORS):
                if color == s.cb:
                    continue
                bottom_options.append((length - 1, color))

    for top in top_options:
        for bottom in bottom_options:
            nt_rt, nt_ct = top
            nt_rb, nt_cb = bottom

            if nt_ct == nt_cb:
                continue

            if top_change and bottom_change and s.prev_vert == 0:
                continue

            out.append(State(nt_rt, nt_ct, nt_rb, nt_cb, 0))

    return out

def multiply(x, y):
    n = len(x)
    z = [[0] * n for _ in range(n)]
    for i in range(n):
        x_row = x[i]
        z_row = z[i]
        non_zero_x = [(k, v) for k, v in enumerate(x_row) if v]
        for k, xik in non_zero_x:
            y_row = y[k]
            for j in range(n):
                if y_row[j]:
                    z_row[j] += xik * y_row[j]
        for j in range(n):
            z_row[j] %= MOD
    return z

def multiply_vec(v, m):
    n = len(m)
    out = [0] * n
    non_zero_v = [(i, val) for i, val in enumerate(v) if val]
    for i, vi in non_zero_v:
        for j in range(n):
            if m[i][j]:
                out[j] += vi * m[i][j]
    for j in range(n):
        out[j] %= MOD
    return out

class Solver:
    def __init__(self):
        self.init_states = initial_states()
        self.build_reachable()
        self.build_transition()
        self.build_initial_vector()

    def build_reachable(self):
        seen = set()
        q = []

        for s in self.init_states:
            if s not in seen:
                seen.add(s)
                q.append(s)

        while q:
            s = q.pop(0)
            for t in next_states(s):
                if t not in seen:
                    seen.add(t)
                    q.append(t)

        self.states = sorted(list(seen))
        self.index = {s: i for i, s in enumerate(self.states)}

    def build_transition(self):
        n = len(self.states)
        self.transition = [[0] * n for _ in range(n)]

        for i in range(n):
            ns = next_states(self.states[i])
            for t in ns:
                j = self.index[t]
                self.transition[i][j] = (self.transition[i][j] + 1) % MOD

    def build_initial_vector(self):
        self.init_vector = [0] * len(self.states)
        for s in self.init_states:
            i = self.index[s]
            self.init_vector[i] = (self.init_vector[i] + 1) % MOD

    def count(self, n):
        vec = list(self.init_vector)
        p = [list(row) for row in self.transition]
        exp = n - 1

        while exp > 0:
            if exp % 2 != 0:
                vec = multiply_vec(vec, p)
            exp //= 2
            if exp > 0:
                p = multiply(p, p)

        ans = 0
        for i in range(len(self.states)):
            s = self.states[i]
            if s.rt == 0 and s.rb == 0:
                ans = (ans + vec[i]) % MOD

        return ans

def solve():
    solver = Solver()
    ans = solver.count(10000000000000000)
    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.HashSet;
import java.util.LinkedList;
import java.util.List;
import java.util.Objects;
import java.util.Queue;

public class Euler670 {

    static final long MOD = 1000004321L;
    static final int COLORS = 4;

    static class State implements Comparable<State> {
        int rt, ct, rb, cb, prevVert;

        State(int rt, int ct, int rb, int cb, int prevVert) {
            this.rt = rt;
            this.ct = ct;
            this.rb = rb;
            this.cb = cb;
            this.prevVert = prevVert;
        }

        @Override
        public boolean equals(Object o) {
            if (this == o)
                return true;
            if (!(o instanceof State))
                return false;
            State s = (State) o;
            return rt == s.rt && ct == s.ct && rb == s.rb && cb == s.cb && prevVert == s.prevVert;
        }

        @Override
        public int hashCode() {
            int h = rt;
            h = h * 1315423911 + (ct + 7);
            h = h * 1315423911 + (rb + 11);
            h = h * 1315423911 + (cb + 13);
            h = h * 1315423911 + (prevVert + 17);
            return h;
        }

        @Override
        public int compareTo(State o) {
            if (rt != o.rt)
                return Integer.compare(rt, o.rt);
            if (ct != o.ct)
                return Integer.compare(ct, o.ct);
            if (rb != o.rb)
                return Integer.compare(rb, o.rb);
            if (cb != o.cb)
                return Integer.compare(cb, o.cb);
            return Integer.compare(prevVert, o.prevVert);
        }
    }

    static List<State> initialStates() {
        List<State> init = new ArrayList<>();

        for (int v = 0; v < COLORS; ++v) {
            init.add(new State(0, v, 0, v, 1));
        }

        for (int lt = 1; lt <= 3; ++lt) {
            for (int lb = 1; lb <= 3; ++lb) {
                for (int a = 0; a < COLORS; ++a) {
                    for (int b = 0; b < COLORS; ++b) {
                        if (a == b)
                            continue;
                        init.add(new State(lt - 1, a, lb - 1, b, 0));
                    }
                }
            }
        }
        return init;
    }

    static class Pair {
        int first, second;

        Pair(int f, int s) {
            first = f;
            second = s;
        }
    }

    static List<State> nextStates(State s) {
        List<State> out = new ArrayList<>();

        boolean topChange = (s.rt == 0);
        boolean bottomChange = (s.rb == 0);

        if (s.rt == 0 && s.rb == 0) {
            for (int v = 0; v < COLORS; ++v) {
                if (v == s.ct || v == s.cb)
                    continue;
                out.add(new State(0, v, 0, v, 1));
            }
        }

        List<Pair> topOptions = new ArrayList<>();
        List<Pair> bottomOptions = new ArrayList<>();

        if (s.rt > 0) {
            topOptions.add(new Pair(s.rt - 1, s.ct));
        } else {
            for (int len = 1; len <= 3; ++len) {
                for (int color = 0; color < COLORS; ++color) {
                    if (color == s.ct)
                        continue;
                    topOptions.add(new Pair(len - 1, color));
                }
            }
        }

        if (s.rb > 0) {
            bottomOptions.add(new Pair(s.rb - 1, s.cb));
        } else {
            for (int len = 1; len <= 3; ++len) {
                for (int color = 0; color < COLORS; ++color) {
                    if (color == s.cb)
                        continue;
                    bottomOptions.add(new Pair(len - 1, color));
                }
            }
        }

        for (Pair top : topOptions) {
            for (Pair bottom : bottomOptions) {
                int ntRt = top.first;
                int ntCt = top.second;
                int ntRb = bottom.first;
                int ntCb = bottom.second;

                if (ntCt == ntCb)
                    continue;

                if (topChange && bottomChange && s.prevVert == 0)
                    continue;

                out.add(new State(ntRt, ntCt, ntRb, ntCb, 0));
            }
        }

        return out;
    }

    static long[][] multiply(long[][] x, long[][] y) {
        int n = x.length;
        long[][] z = new long[n][n];

        for (int i = 0; i < n; ++i) {
            for (int k = 0; k < n; ++k) {
                long xik = x[i][k];
                if (xik == 0)
                    continue;
                for (int j = 0; j < n; ++j) {
                    long ykj = y[k][j];
                    if (ykj == 0)
                        continue;
                    z[i][j] = (z[i][j] + xik * ykj) % MOD;
                }
            }
        }
        return z;
    }

    static long[] multiplyVec(long[] v, long[][] m) {
        int n = m.length;
        long[] out = new long[n];

        for (int i = 0; i < n; ++i) {
            long vi = v[i];
            if (vi == 0)
                continue;
            for (int j = 0; j < n; ++j) {
                long mij = m[i][j];
                if (mij == 0)
                    continue;
                out[j] = (out[j] + vi * mij) % MOD;
            }
        }
        return out;
    }

    static class Solver {
        List<State> initStates;
        List<State> states;
        HashMap<State, Integer> index;
        long[][] transition;
        long[] initVector;

        Solver() {
            initStates = initialStates();
            buildReachable();
            buildTransition();
            buildInitialVector();
        }

        void buildReachable() {
            HashSet<State> seen = new HashSet<>();
            Queue<State> q = new LinkedList<>();

            for (State s : initStates) {
                if (seen.add(s)) {
                    q.offer(s);
                }
            }

            while (!q.isEmpty()) {
                State s = q.poll();
                for (State t : nextStates(s)) {
                    if (seen.add(t)) {
                        q.offer(t);
                    }
                }
            }

            states = new ArrayList<>(seen);
            Collections.sort(states);

            index = new HashMap<>();
            for (int i = 0; i < states.size(); ++i) {
                index.put(states.get(i), i);
            }
        }

        void buildTransition() {
            int n = states.size();
            transition = new long[n][n];

            for (int i = 0; i < n; ++i) {
                for (State t : nextStates(states.get(i))) {
                    int j = index.get(t);
                    transition[i][j] = (transition[i][j] + 1) % MOD;
                }
            }
        }

        void buildInitialVector() {
            int n = states.size();
            initVector = new long[n];
            for (State s : initStates) {
                int i = index.get(s);
                initVector[i] = (initVector[i] + 1) % MOD;
            }
        }

        long count(long n) {
            long[] vec = initVector.clone();
            long[][] p = new long[transition.length][];
            for (int i = 0; i < transition.length; i++) {
                p[i] = transition[i].clone();
            }
            long exp = n - 1;

            while (exp > 0) {
                if ((exp & 1) != 0) {
                    vec = multiplyVec(vec, p);
                }
                exp >>= 1;
                if (exp > 0) {
                    p = multiply(p, p);
                }
            }

            long ans = 0;
            for (int i = 0; i < states.size(); ++i) {
                State s = states.get(i);
                if (s.rt == 0 && s.rb == 0) {
                    ans = (ans + vec[i]) % MOD;
                }
            }
            return ans;
        }
    }

    public static String solve() {
        Solver solver = new Solver();
        long ans = solver.count(10000000000000000L);
        return Long.toString(ans);
    }

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