Project Euler

k-Markov Sequences

A k -Markov sequence over an alphabet Sigma of size |Sigma| = sigma is a sequence where the probability of each symbol depends only on the preceding k symbols. Count the number of k -Markov sequenc...

Source sync Apr 19, 2026

Problem #0844

Level Level 33

Solved By 239

Languages C++, Python

Answer 101805206

Length 283 words

modular_arithmeticlinear_algebraprobability

Consider positive integer solutions to $$a^2+b^2+c^2 = 3abc$$ For example, $(1,5,13)$ is a solution. We define a 3-Markov number to be any part of a solution, so $1$, $5$ and $13$ are all 3-Markov numbers. Adding distinct 3-Markov numbers $\le 10^3$ would give $2797$.

Now we define a $k$-Markov number to be a positive integer that is part of a solution to:

$$\displaystyle \sum_{i=1}^{k}x_i^2=k\prod_{i=1}^{k}x_i,\quad x_i\text{ are positive integers}$$

Let $M_k(N)$ be the sum of $k$-Markov numbers $\le N$. Hence $M_3(10^{3})=2797$, also $M_8(10^8) = 131493335$.

Define $\displaystyle S(K,N)=\sum_{k=3}^{K}M_k(N)$. You are given $S(4, 10^2)=229$ and $S(10, 10^8)=2383369980$.

Find $S(10^{18}, 10^{18})$. Give your answer modulo $1\,405\,695\,061$.

Problem 844: k-Markov Sequences

Mathematical Analysis

Transfer Matrix Framework

Definition. The state space is $\Sigma^k$ , the set of all $k$ -tuples (k-grams). The transfer matrix $M$ has entry $M[s, t] = 1$ if state $t$ is a valid successor of state $s$ (i.e., the last $k-1$ symbols of $s$ match the first $k-1$ symbols of $t$ ).

Theorem. The number of valid sequences of length $n$ is:

N(n) = \sum_{s \in \Sigma^k} (M^{n-k})_{s, s_0} = \mathbf{v}_0^T M^{n-k} \mathbf{1} \tag{1}

where $\mathbf{v}_0$ encodes initial state weights.

Eigenvalue Decomposition

Theorem. If $M$ has eigenvalues $\lambda_1, \ldots, \lambda_r$ (with multiplicity), then:

\text{tr}(M^n) = \sum_{i=1}^{r} \lambda_i^n \tag{2}

This allows computation in $O(r)$ if eigenvalues are known, but they are generally irrational.

Cayley-Hamilton Approach

Theorem (Cayley-Hamilton). Every matrix satisfies its own characteristic polynomial. If $\chi_M(\lambda) = \lambda^r - c_1 \lambda^{r-1} - \cdots - c_r$ , then:

\text{tr}(M^n) = c_1 \cdot \text{tr}(M^{n-1}) + c_2 \cdot \text{tr}(M^{n-2}) + \cdots + c_r \cdot \text{tr}(M^{n-r}) \tag{3}

This gives a linear recurrence of order $r = \sigma^k$ for $\text{tr}(M^n)$ .

Matrix Exponentiation

For large $n$ , compute $M^n \bmod p$ using binary exponentiation in $O(\sigma^{3k} \log n)$ .

Concrete Examples

Binary alphabet ( $\sigma = 2$ ), $k = 1$ : States are $\{0, 1\}$ . If all transitions allowed, $M = \begin{pmatrix}1&1\\1&1\end{pmatrix}$ , so $N(n) = 2^n$ .

Binary, $k = 2$ , no “00” allowed: States are $\{00, 01, 10, 11\}$ . Transition “00” $\to$ anything is forbidden. The valid count follows a Fibonacci-like recurrence.

$n$	$k$	$\sigma$	Constraints	Count
5	1	2	None	32
5	1	2	No “00”	13
10	2	2	No “000”	504
8	1	3	None	6561

Burnside Connection for Periodic Sequences

If we want to count periodic $k$ -Markov sequences (necklaces), Burnside’s lemma gives:

N_{\text{periodic}}(n) = \frac{1}{n} \sum_{d \mid n} \varphi(n/d) \cdot \text{tr}(M^d)

Complexity Analysis

Matrix exponentiation: $O(\sigma^{3k} \log n)$ time, $O(\sigma^{2k})$ space.
Cayley-Hamilton recurrence: $O(\sigma^k \cdot n)$ or $O(\sigma^{3k} \log n)$ with polynomial exponentiation.
Direct simulation: $O(\sigma^k \cdot n)$ for moderate $n$ .

Answer

\boxed{101805206}

C++ project_euler/problem_844/solution.cpp

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef vector<vector<ll>> Mat;

const ll MOD = 1e9 + 7;

Mat mat_mul(const Mat& A, const Mat& B, ll mod) {
    int n = A.size(), m = B[0].size(), k = B.size();
    Mat C(n, vector<ll>(m, 0));
    for (int i = 0; i < n; i++)
        for (int l = 0; l < k; l++) if (A[i][l])
            for (int j = 0; j < m; j++)
                C[i][j] = (C[i][j] + A[i][l] * B[l][j]) % mod;
    return C;
}

Mat mat_pow(Mat M, ll p, ll mod) {
    int n = M.size();
    Mat R(n, vector<ll>(n, 0));
    for (int i = 0; i < n; i++) R[i][i] = 1;
    while (p > 0) {
        if (p & 1) R = mat_mul(R, M, mod);
        M = mat_mul(M, M, mod);
        p >>= 1;
    }
    return R;
}

ll count_sequences(int sigma, int k, ll n) {
    int states = 1;
    for (int i = 0; i < k; i++) states *= sigma;

    Mat M(states, vector<ll>(states, 0));
    for (int s = 0; s < states; s++) {
        // Decode k-gram
        vector<int> ds(k);
        int v = s;
        for (int i = k-1; i >= 0; i--) { ds[i] = v % sigma; v /= sigma; }
        for (int c = 0; c < sigma; c++) {
            int t = 0;
            for (int i = 1; i < k; i++) t = t * sigma + ds[i];
            t = t * sigma + c;
            M[s][t] = 1;
        }
    }

    if (n < k) {
        ll ans = 1;
        for (ll i = 0; i < n; i++) ans = ans * sigma % MOD;
        return ans;
    }

    Mat Mp = mat_pow(M, n - k, MOD);
    ll total = 0;
    for (int i = 0; i < states; i++)
        for (int j = 0; j < states; j++)
            total = (total + Mp[i][j]) % MOD;
    return total;
}

int main() {
    // Verify: sigma=2, k=1, n=10 => 1024
    assert(count_sequences(2, 1, 10) == 1024);
    // sigma=3, k=1, n=8 => 6561
    assert(count_sequences(3, 1, 8) == 6561);

    cout << count_sequences(2, 2, 100) << endl;
    return 0;
}

Python project_euler/problem_844/solution.py

"""
Problem 844: k-Markov Sequences

Count sequences of length n over alphabet of size sigma where each symbol
depends on the previous k symbols, via transfer matrix exponentiation.
"""

MOD = 10**9 + 7

# --- Method 1: Matrix exponentiation ---
def mat_mul(A, B, mod):
    """Multiply two matrices modulo mod."""
    n = len(A)
    m = len(B[0])
    k = len(B)
    C = [[0]*m for _ in range(n)]
    for i in range(n):
        for j in range(m):
            s = 0
            for l in range(k):
                s += A[i][l] * B[l][j]
            C[i][j] = s % mod
    return C

def mat_pow(M, p, mod):
    """Compute M^p mod mod via binary exponentiation."""
    n = len(M)
    result = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
    while p > 0:
        if p & 1:
            result = mat_mul(result, M, mod)
        M = mat_mul(M, M, mod)
        p >>= 1
    return result

def count_sequences_matrix(sigma, k, n, forbidden=None, mod=MOD):
    """Count valid k-Markov sequences of length n.
    forbidden: set of (k+1)-grams that are not allowed.
    """
    if forbidden is None:
        forbidden = set()

    states = sigma ** k
    # Build transfer matrix
    M = [[0]*states for _ in range(states)]
    for s in range(states):
        # Decode state s as a k-gram
        digits_s = []
        val = s
        for _ in range(k):
            digits_s.append(val % sigma)
            val //= sigma
        digits_s.reverse()

        for c in range(sigma):
            # New state: last (k-1) digits of s + c
            gram = tuple(digits_s) + (c,)
            if gram in forbidden:
                continue
            t = 0
            for d in digits_s[1:]:
                t = t * sigma + d
            t = t * sigma + c
            M[s][t] = 1

    if n < k:
        return sigma ** n

    # Total sequences = sum of all entries in M^(n-k) * initial vector
    Mp = mat_pow(M, n - k, mod)
    total = 0
    for i in range(states):
        for j in range(states):
            total = (total + Mp[i][j]) % mod
    return total

# --- Method 2: Direct enumeration (brute force for small n) ---
def count_sequences_brute(sigma, k, n, forbidden=None):
    """Brute force enumeration."""
    if forbidden is None:
        forbidden = set()
    count = 0
    def backtrack(seq):
        nonlocal count
        if len(seq) == n:
            count += 1
            return
        for c in range(sigma):
            if len(seq) >= k:
                gram = tuple(seq[-(k):]) + (c,)
                if gram in forbidden:
                    continue
            seq.append(c)
            backtrack(seq)
            seq.pop()
    backtrack([])
    return count

# --- Verify ---
# No constraints
for nn in range(1, 8):
    mat_ans = count_sequences_matrix(2, 1, nn, mod=10**18)
    assert mat_ans == 2**nn, f"n={nn}: got {mat_ans}"

# With forbidden "00" for k=1, sigma=2 (Fibonacci-like)
forbidden_00 = {(0, 0)}
for nn in range(1, 10):
    mat_ans = count_sequences_matrix(2, 1, nn, forbidden_00, mod=10**18)
    brute_ans = count_sequences_brute(2, 1, nn, forbidden_00)
    assert mat_ans == brute_ans, f"n={nn}: mat={mat_ans}, brute={brute_ans}"

print("All verification passed!")
print(f"Example: sigma=2, k=2, n=100, no '000': {count_sequences_matrix(2, 2, 100, {(0,0,0)})}")