∑MathAdvanced

Linear Recurrence

Key Points

•
A linear recurrence defines each term as a fixed linear combination of a small, fixed number of previous terms.
•
You can fast-forward a linear recurrence by turning it into a matrix multiplication and exponentiating the matrix in O( $k^{3}$ $lo g$ n) time.
•
The characteristic polynomial $x^{k}$ - $c_{1}$ $x^{k - 1}$ - $\dots$ - $c_{k}$ encodes everything about the recurrence, including closed forms and fast algorithms.
•
Kitamasa’s algorithm reduces the complexity to O( $k^{2}$ $lo g$ n) by working with polynomials modulo the characteristic polynomial.
•
Berlekamp–Massey can recover the minimal recurrence from the first few terms of any linearly recurrent sequence over a field (e.g., modulo a prime).
•
Always handle small n < k directly from initial values and be consistent with indexing and matrix orientation to avoid off-by-one errors.
•
Use modular arithmetic (and 128-bit intermediates in C++) to prevent overflow when coefficients and moduli are large.

Prerequisites

→Modular arithmetic — Linear recurrences in programming contests are almost always computed modulo a prime to keep numbers bounded and enable inverses.
→Matrices and vectors — Matrix exponentiation relies on constructing and multiplying the companion matrix with state vectors.
→Exponentiation by squaring — Both matrix and polynomial methods rely on fast powering to achieve O(log n).
→Polynomial arithmetic — Kitamasa reduces powers of x modulo the characteristic polynomial, which requires polynomial multiplication and reduction.
→Basic dynamic programming — Many recurrences come from DPs with constant-width windows; recognizing this pattern helps modeling.
→Number theory over finite fields — Berlekamp–Massey assumes field operations; working modulo a prime ensures well-defined inverses.

Detailed Explanation

Tap terms for definitions

01Overview

Hook: Imagine pressing a “skip ahead 10^12 songs” button on a playlist that moves to the next song by a simple rule: each song’s score depends on the last k scores in a fixed linear way. How could you possibly jump that far instantly? Concept: That is exactly what linear recurrences let us do. A linear recurrence of order k defines a sequence (a_n) where each new term is a fixed linear combination of the previous k terms. Famous examples include the Fibonacci sequence and many counting problems in combinatorics and DP. Example: For Fibonacci, a_n = a_{n-1} + a_{n-2} with a_0 = 0, a_1 = 1. We can represent the update rule as multiplying a fixed k×k matrix (the companion matrix) by a k-dimensional state vector of the last k terms. Then a_n can be computed by fast matrix exponentiation in O(k^3 \log n). For even faster solutions, we can use polynomial-based methods (Kitamasa) in O(k^2 \log n), and when the recurrence is unknown, Berlekamp–Massey can infer it from observed terms. These tools are essential in competitive programming when n is huge (e.g., 10^18) but k is modest (≤ 500).

02Intuition & Analogies

Hook: Think of a vending machine with k memory slots. To produce the next snack, the machine blends fixed portions of the last k snacks. If you know the blending recipe (coefficients) and the initial k snacks, you can simulate output snack by snack. But to jump to the billionth snack, you need a fast-forward button. Concept: The state of the machine at step n is just the last k outputs stacked in a vector. Pushing the machine one step forward applies the same linear recipe each time. Linear operations compose nicely: applying the recipe twice is equivalent to multiplying by the same matrix twice. Therefore, jumping ahead many steps equals raising the transition (companion) matrix to a high power and applying it once. Exponentiation by squaring lets us do this in logarithmic time with respect to n. Example: For Fibonacci, store [a_n, a_{n-1}] as the state. The update is [a_{n+1}, a_n] = [[1,1],[1,0]] · [a_n, a_{n-1}]. To get a_{10^{12}}, compute [[1,1],[1,0]]^{10^{12}} using fast exponentiation instead of iterating 10^{12} times. For higher-order recurrences, the matrix just gets larger with ones on the subdiagonal and the coefficients in the top row. Alternatively, think polynomially: the recurrence says x^n equals a combination of the previous k powers modulo the characteristic polynomial. Exponentiating x by n modulo that polynomial gives a compact linear combination of the first k terms—this is the heart of Kitamasa’s algorithm.

03Formal Definition

A homogeneous linear recurrence of order k over a field (or ring) is a sequence (

a_{n}

) satisfying, for n

\geq

k, \[

a_{n} = c_{1}

a_{n - 1}

c_{2}

a_{n - 2}

\dots

c_{k}

a_{n - k}

, \] with fixed coefficients

c_{1}

\dots

c_{k}

and initial terms

a_{0}

\dots

a_{k - 1}

. Its characteristic polynomial is \[ P(x) =

x^{k}

c_{1}

x^{k - 1}

c_{2}

x^{k - 2}

\dots

c_{k}

. \] Let A denote the k

\times

k companion matrix with first row [

c_{1}

\dots

c_{k}

], ones on the subdiagonal, and zeros elsewhere. Define the state vector

s_{n}

= [

a_{n}

a_{n - 1}

\dots

a_{n - k + 1}

]^T. Then

s_{n} = A\

s_{n - 1}

for n

\geq

k-1, hence

s_{n} = A^{n - k + 1}

s_{k - 1}

. By the Cayley–Hamilton theorem, P(A)=0, ensuring that powers of A can be reduced to a linear combination of I, A,

\dots

A^{k - 1}

, consistent with the recurrence. If P has distinct roots

r_{1}

\dots

r_{k}

in an extension field, then the closed form is \[

a_{n}

\sum_{i = 1}^{k}

α_{i}

i^{n}

, \] for suitable constants

α_{i}

determined by the initial conditions. With repeated roots, polynomial factors in n appear, e.g.,

α_{i, 0}

i^{n}

α_{i, 1}

n r_

i^{n}

\dots

04When to Use

Computing a_n for very large n (up to 10^{18} or more) when k is small/moderate and the coefficients are constant.
Speeding up DPs with constant-width windows where dp[n] is a fixed linear combination of the last k dp values (e.g., tilings with fixed-size tiles, constrained walks with finite memory).
Counting paths in automata or Markov-like systems with a fixed state transition that’s linear.
Sequences with periodic or modular behavior where you want a_n modulo a prime (e.g., 10^9+7).
Recovering hidden recurrences from the first 2k terms using Berlekamp–Massey, then predicting future terms quickly.
Improving from O(nk) straightforward DP to O(k^3 \log n) via matrix exponentiation or O(k^2 \log n) via Kitamasa when n is huge. Concrete examples: Fibonacci-like sequences, linear homogeneous recurrences arising in linear DP transitions, counting domino tilings of 2\times n boards, or evaluating linear feedback shift register (LFSR) sequences modulo a prime.

⚠️Common Mistakes

Off-by-one indexing: Mixing 0-based and 1-based a_n definitions or misaligning the state vector with the companion matrix. Always define clearly whether your top row maps [a_{n},...,a_{n-k+1}] to [a_{n+1},...,a_{n-k+2}].
Forgetting base cases: For n < k, you must return the provided initial terms. Don’t exponentiate needlessly.
Wrong companion matrix orientation: The standard form uses coefficients in the first row and ones on the subdiagonal. Verify with a small example (e.g., Fibonacci) before generalizing.
Modular overflow: In C++, intermediate products can exceed 64-bit when summing many terms. Use 128-bit temporaries (__int128) and reduce modulo after each multiply-add.
Using O(nk) simulation for huge n: Switch to matrix exponentiation or Kitamasa; otherwise you time out.
Misusing Berlekamp–Massey: It requires field arithmetic; use a prime modulus (or work over rationals with care). Ensure you feed enough terms (typically 2k) and handle zero discrepancies properly.
Numerical instability: Avoid floating-point root finding for closed forms in programming contests; stick to exact modular methods.

Key Formulas

Linear Recurrence

a_{n} = c_{1} a_{n - 1} + c_{2} a_{n - 2} + \dots + c_{k} a_{n - k}

Explanation: Each term is a fixed linear combination of the previous k terms. The coefficients $c_{i}$ are constants that define the sequence's rule.

Characteristic Polynomial

P (x) = x^{k} - c_{1} x^{k - 1} - c_{2} x^{k - 2} - \dots - c_{k}

Explanation: This polynomial encodes the recurrence. Its roots determine the closed form and enable polynomial reduction methods like Kitamasa.

Matrix Power Form

s_{n} = A^{n - k + 1} s_{k - 1}

Explanation: If $s_{n}$ is the state vector [ $a_{n}$ , $a_{n - 1}$ , ..., $a_{n - k + 1}$ ]^T and A is the companion matrix, then advancing by n−k+1 steps is multiplying by A that many times.

Companion Matrix

a_{n} a_{n - 1} ⋮ a_{n - k + 1} = c_{1} 10 ⋮ 0 c_{2} 01 ⋮ 0 \dots \dots \dots ⋱ \dots c_{k} 00 ⋮ 1 a_{n - 1} a_{n - 2} ⋮ a_{n - k}

Explanation: This shows how the next state is a fixed linear transformation of the previous state. The top row holds the coefficients and the subdiagonal shifts the window.

Cayley–Hamilton

P (A) = 0

Explanation: Substituting the companion matrix A into its characteristic polynomial yields the zero matrix. This allows reducing $A^{t}$ to a combination of I, A, ..., $A^{k - 1}$ .

Closed Form with Multiplicities

a_{n} = i = 1 \sum m j = 0 \sum m_{i} - 1 α_{i, j} n^{j} r_{i}^{n}

Explanation: If P(x) has roots $r_{i}$ with multiplicities $m_{i}$ , the general solution is a linear combination of terms $r_{i}$ ^n times polynomials in n of degree $m_{i}$ −1.

Matrix Exponentiation Complexity

T_{mat} (n, k) = O (k^{3} lo g n)

Explanation: Raising a k×k matrix to the nth power by repeated squaring takes O(log n) matrix multiplications, each O( $k^{3}$ ) using the naive algorithm.

Kitamasa Complexity

T_{kit} (n, k) = O (k^{2} lo g n)

Explanation: Exponentiating x modulo P(x) uses O(log n) polynomial multiplications and reductions, each O( $k^{2}$ ) with naive convolution.

Summation Reference

i = 0 \sum n i = \frac{n ( n + 1 )}{2}

Explanation: Commonly appears when analyzing cost accumulations. Here it reminds that adding up linear costs yields quadratic growth in k if nested.

Complexity Analysis

We consider three main approaches. 1) Naive DP simulation takes O(nk) time and O(k) space, since each term uses k previous terms. This is only feasible for n up to around 10^7 when k is small. 2) Matrix exponentiation converts the k-th order recurrence into a k×k companion matrix. Exponentiation by squaring requires O(log n) matrix multiplications; with naive O(

k^{3}

) multiplication, the total time is O(

k^{3}

log n). Space is O(

k^{2}

) to store matrices. For modest k (≤ 60–100), this is typically fine, and even larger k can pass with careful optimization. 3) Kitamasa’s algorithm works directly with polynomials modulo the characteristic polynomial, maintaining coefficient vectors of length k. Polynomial multiplication and reduction each cost O(

k^{2}

), and we need O(log n) such steps, yielding O(

k^{2}

log n) time and O(k) space. It is usually faster than matrix exponentiation for larger k. If FFT/NTT is applicable, polynomial multiplications can be improved to O(k log k), giving O(k log k log n), but constant factors are higher. When the recurrence is unknown but the sequence is linearly recurrent, Berlekamp–Massey recovers a minimal recurrence in O(

L^{2}

), where L is the length of the input prefix (often ≈ 2k). After recovery, compute

a_{n}

via Kitamasa in O(

k^{2}

log n). All methods use O(k) or O(

k^{2}

) memory. Over modular arithmetic, use 128-bit intermediates in C++ to avoid overflow of partial sums. Always handle

n < k

in O(1) by returning initial values.

Code Examples

Compute a_n with Matrix Exponentiation (General k, modulo prime)

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 using int64 = long long;
5 const int64 MOD = 1000000007LL; // prime modulus
6 
7 // Multiply two k x k matrices modulo MOD
8 vector<vector<int64>> matMul(const vector<vector<int64>>& A, const vector<vector<int64>>& B) {
9     int k = (int)A.size();
10     vector<vector<int64>> C(k, vector<int64>(k, 0));
11     for (int i = 0; i < k; ++i) {
12         for (int t = 0; t < k; ++t) {
13             __int128 a = A[i][t];
14             if (a == 0) continue;
15             for (int j = 0; j < k; ++j) {
16                 C[i][j] = (C[i][j] + (int64)((a * B[t][j]) % MOD)) % MOD;
17             }
18         }
19     }
20     return C;
21 }
22 
23 // Multiply matrix A (k x k) with vector v (k), result is k-vector
24 vector<int64> matVecMul(const vector<vector<int64>>& A, const vector<int64>& v) {
25     int k = (int)A.size();
26     vector<int64> r(k, 0);
27     for (int i = 0; i < k; ++i) {
28         __int128 sum = 0;
29         for (int j = 0; j < k; ++j) {
30             sum += (__int128)A[i][j] * v[j];
31         }
32         r[i] = (int64)(sum % MOD);
33     }
34     return r;
35 }
36 
37 // Fast exponentiation of matrix A to power e
38 vector<vector<int64>> matPow(vector<vector<int64>> A, long long e) {
39     int k = (int)A.size();
40     vector<vector<int64>> R(k, vector<int64>(k, 0));
41     for (int i = 0; i < k; ++i) R[i][i] = 1; // identity
42     while (e > 0) {
43         if (e & 1) R = matMul(R, A);
44         A = matMul(A, A);
45         e >>= 1;
46     }
47     return R;
48 }
49 
50 // Computes a_n for a k-th order recurrence: a_n = sum_{i=1..k} c[i-1] * a_{n-i}
51 // Inputs:
52 //   c: coefficients of size k (0-based), c[0] multiplies a_{n-1}
53 //   init: initial values a_0..a_{k-1}
54 //   n: index to compute
55 int64 linearRecurrenceMatrix(const vector<int64>& c, const vector<int64>& init, long long n) {
56     int k = (int)c.size();
57     if (n < (long long)k) return (init[(size_t)n] % MOD + MOD) % MOD;
58 
59     // Build companion matrix A (k x k)
60     // State vector is [a_{t}, a_{t-1}, ..., a_{t-k+1}]^T
61     vector<vector<int64>> A(k, vector<int64>(k, 0));
62     // top row: coefficients
63     for (int j = 0; j < k; ++j) A[0][j] = (c[j] % MOD + MOD) % MOD;
64     // subdiagonal ones
65     for (int i = 1; i < k; ++i) A[i][i-1] = 1;
66 
67     // Raise A to the power (n - (k - 1))
68     long long e = n - (k - 1);
69     vector<vector<int64>> Ae = matPow(A, e);
70 
71     // Build initial state s_{k-1} = [a_{k-1}, a_{k-2}, ..., a_0]^T
72     vector<int64> s(k);
73     for (int i = 0; i < k; ++i) s[i] = (init[(size_t)(k - 1 - i)] % MOD + MOD) % MOD;
74 
75     // s_n = Ae * s_{k-1}
76     vector<int64> sn = matVecMul(Ae, s);
77     // a_n is the first entry of s_n
78     return sn[0];
79 }
80 
81 int main() {
82     ios::sync_with_stdio(false);
83     cin.tie(nullptr);
84 
85     // Example: Fibonacci a_n = a_{n-1} + a_{n-2}, a_0=0, a_1=1
86     vector<int64> c = {1, 1};
87     vector<int64> init = {0, 1};
88     long long n = 1000000000000LL; // 1e12
89     cout << linearRecurrenceMatrix(c, init, n) << "\n";
90 
91     // Another example: a_n = 2 a_{n-1} - 3 a_{n-2} + 5 a_{n-3}
92     vector<int64> c2 = {2, -3, 5};
93     vector<int64> init2 = {1, 4, 9}; // a_0, a_1, a_2
94     long long n2 = 1000000; // 1e6
95     cout << linearRecurrenceMatrix(c2, init2, n2) << "\n";
96 
97     return 0;
98 }
99

We model the k-th order recurrence with a k×k companion matrix that maps the state [a_t, a_{t-1}, ..., a_{t-k+1}] to [a_{t+1}, ..., a_{t-k+2}]. Then a_n is obtained by raising this matrix to the appropriate power and multiplying it by the initial state. Modular arithmetic ensures values stay bounded. The code handles n < k directly and uses __int128 to prevent overflow during accumulation.

Time: O(k^3 log n)Space: O(k^2)

Kitamasa’s Algorithm: O(k^2 log n) computation of a_n

1 #include <bits/stdc++.h>
2 using namespace std;
3 using int64 = long long;
4 const int64 MOD = 1000000007LL;
5 
6 // We want to compute a_n where a_n = sum_{i=1..k} c[i-1]*a_{n-i}
7 // Kitamasa computes coefficients beta[0..k-1] such that
8 // a_n = sum_{i=0..k-1} beta[i] * a_i.
9 // It does so by computing x^n modulo P(x) = x^k - c_1 x^{k-1} - ... - c_k
10 // represented as a vector of length k (reduced basis).
11 
12 // Multiply two polynomials A and B (degree < k) modulo P using naive O(k^2)
13 vector<int64> polyMulMod(const vector<int64>& A, const vector<int64>& B, const vector<int64>& c) {
14     int k = (int)c.size();
15     vector<__int128> tmp(2*k - 1, 0);
16     for (int i = 0; i < k; ++i) if (A[i]) {
17         for (int j = 0; j < k; ++j) if (B[j]) {
18             tmp[i + j] += (__int128)A[i] * B[j];
19         }
20     }
21     // Reduce modulo P(x) = x^k - c1 x^{k-1} - ... - ck
22     for (int deg = 2*k - 2; deg >= k; --deg) {
23         if (tmp[deg] == 0) continue;
24         __int128 coef = tmp[deg];
25         int shift = deg - k;
26         // x^{deg} = (c1 x^{deg-1} + ... + ck x^{deg-k}) under modulo P
27         for (int j = 1; j <= k; ++j) {
28             tmp[shift + (k - j)] += coef * ( (__int128)((c[j-1] % MOD + MOD) % MOD) );
29         }
30         tmp[deg] = 0;
31     }
32     vector<int64> res(k, 0);
33     for (int i = 0; i < k; ++i) res[i] = (int64)(tmp[i] % MOD);
34     return res;
35 }
36 
37 // Exponentiate x to the n modulo P, where x is represented by vecX = [0,1,0,0,...] (i.e., x)
38 // We maintain polynomials of degree < k as vectors of length k.
39 vector<int64> kitamasaCoeff(long long n, const vector<int64>& c) {
40     int k = (int)c.size();
41     // base: x^0 mod P = 1 => [1,0,0,...]
42     vector<int64> res(k, 0); res[0] = 1; // represents 1
43     // base x^1 mod P = x => [0,1,0,...]
44     vector<int64> x(k, 0); if (k >= 2) x[1] = 1; else x[0] = c[0]; // handle k=1 carefully
45 
46     long long e = n;
47     vector<int64> base = x;
48     while (e > 0) {
49         if (e & 1) res = polyMulMod(res, base, c);
50         base = polyMulMod(base, base, c);
51         e >>= 1;
52     }
53     return res; // coefficients beta s.t. a_n = sum beta[i] * a_i
54 }
55 
56 // Compute a_n using Kitamasa
57 int64 linearRecurrenceKitamasa(const vector<int64>& c, const vector<int64>& init, long long n) {
58     int k = (int)c.size();
59     if (n < (long long)k) return (init[(size_t)n] % MOD + MOD) % MOD;
60     vector<int64> beta = kitamasaCoeff(n, c);
61     __int128 ans = 0;
62     for (int i = 0; i < k; ++i) {
63         ans += (__int128)beta[i] * ((init[(size_t)i] % MOD + MOD) % MOD);
64     }
65     return (int64)(ans % MOD);
66 }
67 
68 int main(){
69     ios::sync_with_stdio(false);
70     cin.tie(nullptr);
71 
72     // Example: Fibonacci
73     vector<int64> c = {1, 1};
74     vector<int64> init = {0, 1};
75     long long n = 1000000000000LL; // 1e12
76     cout << linearRecurrenceKitamasa(c, init, n) << "\n";
77 
78     // Example: 3rd-order recurrence
79     vector<int64> c2 = {2, -3, 5};
80     vector<int64> init2 = {1, 4, 9};
81     long long n2 = 1000000LL;
82     cout << linearRecurrenceKitamasa(c2, init2, n2) << "\n";
83     return 0;
84 }
85

Kitamasa computes coefficients of x^n modulo the characteristic polynomial P(x). This produces a vector β such that a_n = Σ β_i a_i. Exponentiation by squaring makes it O(k^2 log n). The code uses polynomial multiplication with reduction according to x^k ≡ c_1 x^{k-1} + ... + c_k, consistent with the chosen coefficient orientation.

Time: O(k^2 log n)Space: O(k)

Berlekamp–Massey + Kitamasa: Infer recurrence and predict a_n

1 #include <bits/stdc++.h>
2 using namespace std;
3 using int64 = long long;
4 const int64 MOD = 1000000007LL; // must be prime for field inverses
5 
6 int64 modPow(int64 a, int64 e){
7     int64 r = 1%MOD; a%=MOD; if(a<0) a+=MOD;
8     while(e){ if(e&1) r = (__int128)r*a%MOD; a = (__int128)a*a%MOD; e>>=1; }
9     return r;
10 }
11 int64 modInv(int64 a){ return modPow(a, MOD-2); }
12 
13 // Berlekamp–Massey: given sequence s[0..m-1], find minimal recurrence
14 // Returns coefficients c[0..L-1] such that s_n = sum_{i=1..L} c[i-1]*s_{n-i}
15 vector<int64> berlekampMassey(const vector<int64>& s){
16     vector<int64> C = {1}, B = {1}; // connection polynomials
17     int64 b = 1; // last discrepancy
18     int L = 0; // current length of recurrence
19     int m = 0; // steps since last update of B
20     for(size_t n = 0; n < s.size(); ++n){
21         // compute discrepancy d = s[n] + sum_{i=1..L} C[i]*s[n-i]
22         int64 d = s[n] % MOD;
23         for(int i = 1; i <= L; ++i){ d = (d + (__int128)C[i]*s[n - i]) % MOD; }
24         if(d < 0) d += MOD;
25         if(d == 0){ ++m; continue; }
26         vector<int64> T = C; // backup C
27         int64 coef = (MOD - d) * modInv(b) % MOD; // -d / b
28         // C(x) <- C(x) + coef * x^m * B(x)
29         if((int)C.size() < (int)B.size() + m) C.resize(B.size() + m, 0);
30         for(size_t i = 0; i < B.size(); ++i){
31             C[i + m] = (C[i + m] + (__int128)coef * B[i]) % MOD;
32         }
33         if(2*L <= (int)n){
34             L = (int)n + 1 - L;
35             B = T;
36             b = d;
37             m = 1;
38         } else {
39             ++m;
40         }
41     }
42     // C is 1 + c1 x + c2 x^2 + ...; we need c[0..L-1] with sign flipped
43     vector<int64> c(L, 0);
44     for(int i = 0; i < L; ++i){ c[i] = (MOD - C[i+1]) % MOD; }
45     return c;
46 }
47 
48 // Kitamasa helper: multiply polynomials modulo P defined by c
49 vector<int64> polyMulMod(const vector<int64>& A, const vector<int64>& B, const vector<int64>& c) {
50     int k = (int)c.size();
51     vector<__int128> tmp(2*k - 1, 0);
52     for (int i = 0; i < k; ++i) if (A[i]) {
53         for (int j = 0; j < k; ++j) if (B[j]) tmp[i+j] += (__int128)A[i]*B[j];
54     }
55     for (int deg = 2*k - 2; deg >= k; --deg) {
56         if (tmp[deg] == 0) continue;
57         __int128 coef = tmp[deg];
58         int shift = deg - k;
59         for (int j = 1; j <= k; ++j) tmp[shift + (k - j)] += coef * ( (__int128)((c[j-1]%MOD+MOD)%MOD) );
60         tmp[deg] = 0;
61     }
62     vector<int64> res(k, 0);
63     for (int i = 0; i < k; ++i) res[i] = (int64)(tmp[i] % MOD);
64     return res;
65 }
66 
67 vector<int64> kitamasaCoeff(long long n, const vector<int64>& c) {
68     int k = (int)c.size();
69     vector<int64> res(k, 0); res[0] = 1;
70     vector<int64> x(k, 0); if (k >= 2) x[1] = 1; else x[0] = c[0];
71     vector<int64> base = x;
72     while(n){
73         if(n&1) res = polyMulMod(res, base, c);
74         base = polyMulMod(base, base, c);
75         n >>= 1;
76     }
77     return res;
78 }
79 
80 int64 nthFromRec(const vector<int64>& c, const vector<int64>& init, long long n){
81     int k = (int)c.size();
82     if (n < (long long)k) return (init[(size_t)n] % MOD + MOD) % MOD;
83     vector<int64> beta = kitamasaCoeff(n, c);
84     __int128 ans = 0;
85     for (int i = 0; i < k; ++i) ans += (__int128)beta[i] * ((init[(size_t)i]%MOD+MOD)%MOD);
86     return (int64)(ans % MOD);
87 }
88 
89 int main(){
90     ios::sync_with_stdio(false);
91     cin.tie(nullptr);
92 
93     // Suppose we only know the first terms of a sequence (mod MOD)
94     // We'll fabricate a sequence from a known recurrence and see if BM recovers it.
95     vector<int64> true_c = {3, MOD-1, 2}; // a_n = 3 a_{n-1} - a_{n-2} + 2 a_{n-3}
96     vector<int64> init = {5, 7, 11};
97 
98     // Generate first M terms by naive O(nk) (for demonstration)
99     int M = 100; int k = (int)true_c.size();
100     vector<int64> s(max(M, k), 0);
101     for(int i=0;i<k;++i) s[i] = (init[i]%MOD+MOD)%MOD;
102     for(int n=k;n<M;++n){
103         __int128 v = 0;
104         for(int j=1;j<=k;++j){ v += (__int128)true_c[j-1]*s[n-j]; }
105         s[n] = (int64)(v % MOD);
106     }
107 
108     // Recover recurrence using BM from first 2k..M terms
109     vector<int64> c = berlekampMassey(s);
110 
111     // Predict a very large index using recovered recurrence
112     long long nQuery = 1000000000000LL; // 1e12
113     int64 ans = nthFromRec(c, vector<int64>(s.begin(), s.begin()+ (int)c.size()), nQuery);
114     cout << ans << "\n";
115 
116     return 0;
117 }
118

We first generate a test sequence from a known recurrence. Berlekamp–Massey recovers the minimal recurrence over a prime modulus. Then we compute a_n for huge n using Kitamasa with the recovered coefficients. In real problems, you are given the first 2k to a few hundred terms; BM infers the recurrence, and Kitamasa fast-forwards it.

Time: O(L^2 + k^2 log n), where L is the number of known terms (typically ≈ 2k)Space: O(L + k)

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	using int64 = long long;
5	const int64 MOD = 1000000007LL; // prime modulus
6
7	// Multiply two k x k matrices modulo MOD
8	vector<vector<int64>> matMul(const vector<vector<int64>>& A, const vector<vector<int64>>& B) {
9	int k = (int)A.size();
10	vector<vector<int64>> C(k, vector<int64>(k, 0));
11	for (int i = 0; i < k; ++i) {
12	for (int t = 0; t < k; ++t) {
13	__int128 a = A[i][t];
14	if (a == 0) continue;
15	for (int j = 0; j < k; ++j) {
16	C[i][j] = (C[i][j] + (int64)((a * B[t][j]) % MOD)) % MOD;
17	}
18	}
19	}
20	return C;
21	}
22
23	// Multiply matrix A (k x k) with vector v (k), result is k-vector
24	vector<int64> matVecMul(const vector<vector<int64>>& A, const vector<int64>& v) {
25	int k = (int)A.size();
26	vector<int64> r(k, 0);
27	for (int i = 0; i < k; ++i) {
28	__int128 sum = 0;
29	for (int j = 0; j < k; ++j) {
30	sum += (__int128)A[i][j] * v[j];
31	}
32	r[i] = (int64)(sum % MOD);
33	}
34	return r;
35	}
36
37	// Fast exponentiation of matrix A to power e
38	vector<vector<int64>> matPow(vector<vector<int64>> A, long long e) {
39	int k = (int)A.size();
40	vector<vector<int64>> R(k, vector<int64>(k, 0));
41	for (int i = 0; i < k; ++i) R[i][i] = 1; // identity
42	while (e > 0) {
43	if (e & 1) R = matMul(R, A);
44	A = matMul(A, A);
45	e >>= 1;
46	}
47	return R;
48	}
49
50	// Computes a_n for a k-th order recurrence: a_n = sum_{i=1..k} c[i-1] * a_{n-i}
51	// Inputs:
52	// c: coefficients of size k (0-based), c[0] multiplies a_{n-1}
53	// init: initial values a_0..a_{k-1}
54	// n: index to compute
55	int64 linearRecurrenceMatrix(const vector<int64>& c, const vector<int64>& init, long long n) {
56	int k = (int)c.size();
57	if (n < (long long)k) return (init[(size_t)n] % MOD + MOD) % MOD;
58
59	// Build companion matrix A (k x k)
60	// State vector is [a_{t}, a_{t-1}, ..., a_{t-k+1}]^T
61	vector<vector<int64>> A(k, vector<int64>(k, 0));
62	// top row: coefficients
63	for (int j = 0; j < k; ++j) A[0][j] = (c[j] % MOD + MOD) % MOD;
64	// subdiagonal ones
65	for (int i = 1; i < k; ++i) A[i][i-1] = 1;
66
67	// Raise A to the power (n - (k - 1))
68	long long e = n - (k - 1);
69	vector<vector<int64>> Ae = matPow(A, e);
70
71	// Build initial state s_{k-1} = [a_{k-1}, a_{k-2}, ..., a_0]^T
72	vector<int64> s(k);
73	for (int i = 0; i < k; ++i) s[i] = (init[(size_t)(k - 1 - i)] % MOD + MOD) % MOD;
74
75	// s_n = Ae * s_{k-1}
76	vector<int64> sn = matVecMul(Ae, s);
77	// a_n is the first entry of s_n
78	return sn[0];
79	}
80
81	int main() {
82	ios::sync_with_stdio(false);
83	cin.tie(nullptr);
84
85	// Example: Fibonacci a_n = a_{n-1} + a_{n-2}, a_0=0, a_1=1
86	vector<int64> c = {1, 1};
87	vector<int64> init = {0, 1};
88	long long n = 1000000000000LL; // 1e12
89	cout << linearRecurrenceMatrix(c, init, n) << "\n";
90
91	// Another example: a_n = 2 a_{n-1} - 3 a_{n-2} + 5 a_{n-3}
92	vector<int64> c2 = {2, -3, 5};
93	vector<int64> init2 = {1, 4, 9}; // a_0, a_1, a_2
94	long long n2 = 1000000; // 1e6
95	cout << linearRecurrenceMatrix(c2, init2, n2) << "\n";
96
97	return 0;
98	}
99