#include <bits/stdc++.h>
using namespace std;

using ll = long long;
using Matrix = vector<vector<ll>>;

// ===================================================================
// This file contains a collection of Matrix Exponentiation algorithms.
// Each function is ready to be used as a "black box".
// Read the comments above each one to understand:
//   - What it solves
//   - What input it expects
//   - What it returns
//   - Time complexity
//   - Important constraints / assumptions
// ===================================================================

// ===================================================================
// SECTION 1: BASIC MATRIX OPERATIONS
// ===================================================================

// -------------------------------------------------------------------
// matMul(A, B, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Multiplies two matrices A and B.
//
// INPUT:
//   A : a matrix of size (n x p), B : a matrix of size (p x m)
//   mod : modulo value (e.g., 1e9+7)
//
// OUTPUT:
//   Returns a new matrix C = A * B of size (n x m), with all entries
//   reduced modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(n * p * m)  (naive triple loop)
//
// CONSTRAINTS / PRECONDITIONS:
//   - The number of columns of A must equal the number of rows of B.
//   - mod should be positive.
//   - Entries should be non‑negative or already reduced; multiplication
//     may overflow 64‑bit if mod is large, so 'mod' is typically < 1e9.
//   - For larger matrices, consider using __int128 if needed.
Matrix matMul(const Matrix& A, const Matrix& B, ll mod) {
    int n = (int)A.size();
    int p = (int)A[0].size();
    int m = (int)B[0].size();
    Matrix C(n, vector<ll>(m, 0));
    for (int i = 0; i < n; ++i) {
        for (int k = 0; k < p; ++k) {
            if (A[i][k] == 0) continue;        // skip zeros for speed
            ll aik = A[i][k];
            for (int j = 0; j < m; ++j) {
                C[i][j] = (C[i][j] + aik * B[k][j]) % mod;
            }
        }
    }
    return C;
}

// -------------------------------------------------------------------
// matAdd(A, B, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Adds two matrices of the same size.
//
// INPUT:
//   A, B : matrices of size (n x m)
//   mod  : modulo value
//
// OUTPUT:
//   Returns a new matrix C = (A + B) modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(n * m)
Matrix matAdd(const Matrix& A, const Matrix& B, ll mod) {
    int n = (int)A.size();
    int m = (int)A[0].size();
    Matrix C(n, vector<ll>(m, 0));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < m; ++j)
            C[i][j] = (A[i][j] + B[i][j]) % mod;
    return C;
}

// -------------------------------------------------------------------
// matPow(base, exp, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Raises a square matrix 'base' to the power 'exp' using binary
//   exponentiation. This is the core function for matrix exponentiation.
//
// INPUT:
//   base : a square matrix of size (n x n)
//   exp  : exponent (non‑negative integer, can be up to 1e18)
//   mod  : modulo value
//
// OUTPUT:
//   Returns the matrix base^exp, with all entries reduced modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(n^3 * log(exp))  (each multiplication is O(n^3), repeated log exp times)
//
// CONSTRAINTS / PRECONDITIONS:
//   - base must be square.
//   - exp >= 0.
//   - mod > 0.
Matrix matPow(Matrix base, ll exp, ll mod) {
    int n = (int)base.size();
    // Initialize result as identity matrix
    Matrix res(n, vector<ll>(n, 0));
    for (int i = 0; i < n; ++i) res[i][i] = 1 % mod;

    while (exp > 0) {
        if (exp & 1) res = matMul(res, base, mod);
        base = matMul(base, base, mod);
        exp >>= 1;
    }
    return res;
}

// -------------------------------------------------------------------
// matVecMul(M, v, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Multiplies a matrix M (size n x m) by a column vector v (size m x 1).
//
// INPUT:
//   M : matrix of size (n x m)
//   v : vector of length m
//   mod : modulo value
//
// OUTPUT:
//   Returns a vector of length n = M * v (mod 'mod').
//
// TIME COMPLEXITY:
//   O(n * m)
vector<ll> matVecMul(const Matrix& M, const vector<ll>& v, ll mod) {
    int n = (int)M.size();
    int m = (int)M[0].size();
    vector<ll> res(n, 0);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            res[i] = (res[i] + M[i][j] * v[j]) % mod;
        }
    }
    return res;
}

// ===================================================================
// SECTION 2: FIBONACCI AND LINEAR RECURRENCES (using Matrix)
// ===================================================================

// -------------------------------------------------------------------
// fib(n, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Computes the n‑th Fibonacci number (F(0)=0, F(1)=1) modulo 'mod'.
//
// INPUT:
//   n   : index (n >= 0, can be up to 1e18)
//   mod : modulo value
//
// OUTPUT:
//   Returns F(n) modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(log n)  (2x2 matrix exponentiation)
//
// CONSTRAINTS / PRECONDITIONS:
//   - mod > 0.
//   - Works for any non‑negative integer n.
ll fib(ll n, ll mod) {
    if (n == 0) return 0;
    if (n == 1) return 1 % mod;
    // Fibonacci transition matrix: [[1,1],[1,0]]
    Matrix base = {{1 % mod, 1 % mod}, {1 % mod, 0}};
    Matrix res = matPow(base, n - 1, mod);
    // F(n) = res[0][0] * F(1) + res[0][1] * F(0) = res[0][0]
    return res[0][0];
}

// -------------------------------------------------------------------
// linearRecurrence(init, coeff, n, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Computes the n‑th term of a linear recurrence.
//
//   The recurrence is defined as:
//     f[n] = coeff[0] * f[n-1] + coeff[1] * f[n-2] + ... + coeff[k-1] * f[n-k]
//   for n >= k, with initial terms f[0], f[1], ..., f[k-1] given in 'init'.
//
// INPUT:
//   init  : vector of length k, containing f[0] ... f[k-1]
//   coeff : vector of length k, coefficients in the order shown above
//   n     : index of the term to compute (0‑based, n >= 0)
//   mod   : modulo value
//
// OUTPUT:
//   Returns f[n] modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(k^3 * log n) using matrix exponentiation.
//   For k up to ~50 it is acceptable; for larger k use Kitamasa (see below).
//
// CONSTRAINTS / PRECONDITIONS:
//   - k >= 1.
//   - n >= 0.
//   - If n < k, the function directly returns init[n].
//   - The recurrence must hold for n >= k.
//   - All values should be reduced modulo 'mod'.
ll linearRecurrence(const vector<ll>& init, const vector<ll>& coeff, ll n, ll mod) {
    int k = (int)init.size();
    if (n < k) return init[n] % mod;

    // Build companion matrix of size k x k.
    // State vector: [f[t], f[t-1], ..., f[t-k+1]]^T
    Matrix T(k, vector<ll>(k, 0));
    for (int i = 0; i < k; ++i) T[0][i] = coeff[i] % mod;       // first row
    for (int i = 1; i < k; ++i) T[i][i-1] = 1;                  // sub-diagonal

    // We need T^(n - k + 1) because we start from state at t = k-1.
    Matrix Tpow = matPow(T, n - k + 1, mod);

    // Initial state at t = k-1: [f[k-1], f[k-2], ..., f[0]]^T
    vector<ll> state(k);
    for (int i = 0; i < k; ++i) state[i] = init[k-1-i] % mod;

    vector<ll> res = matVecMul(Tpow, state, mod);
    return res[0] % mod;
}

// ===================================================================
// SECTION 3: ADVANCED LINEAR RECURRENCE (Kitamasa / Polynomial)
// ===================================================================
// These functions compute the n‑th term of a linear recurrence in
// O(k^2 log n) instead of O(k^3 log n), which is useful when k is large.
//
// TERMINOLOGY:
//   - Characteristic polynomial: derived from the recurrence.
//     For recurrence f[n] = c0*f[n-1] + c1*f[n-2] + ... + c[k-1]*f[n-k],
//     the characteristic polynomial is:
//       P(x) = x^k - c0*x^(k-1) - c1*x^(k-2) - ... - c[k-1].
//   - We compute x^n mod P(x) using binary exponentiation of polynomials.
//   - Then f[n] = sum_{i=0}^{k-1} r[i] * f[i], where r is the remainder
//     polynomial.
// -------------------------------------------------------------------

// Helper: Multiply two polynomials modulo the characteristic polynomial.
// Both polynomials have degree < k. The product is reduced using
// the relation: x^k = c0*x^(k-1) + c1*x^(k-2) + ... + c[k-1].
vector<ll> polyMulMod(const vector<ll>& a, const vector<ll>& b,
                      const vector<ll>& coeff, ll mod) {
    int k = (int)coeff.size();
    vector<ll> res(2 * k - 1, 0);
    // Multiply
    for (int i = 0; i < k; ++i) {
        if (a[i] == 0) continue;
        for (int j = 0; j < k; ++j) {
            res[i+j] = (res[i+j] + a[i] * b[j]) % mod;
        }
    }
    // Reduce terms with degree >= k using the recurrence.
    // We iterate from high degree down to k.
    for (int deg = 2*k - 2; deg >= k; --deg) {
        if (res[deg] == 0) continue;
        ll coef = res[deg];
        // For each i from 0 to k-1, x^deg = coeff[i] * x^(deg-1-i)
        // Because x^k = coeff[0]*x^(k-1) + coeff[1]*x^(k-2) + ... + coeff[k-1]
        // So x^deg = x^(deg-k) * x^k = sum_{i=0}^{k-1} coeff[i] * x^(deg-1-i)
        for (int i = 0; i < k; ++i) {
            res[deg - 1 - i] = (res[deg - 1 - i] + coef * coeff[i]) % mod;
        }
        // The term res[deg] is now eliminated.
    }
    // Return only the first k coefficients.
    vector<ll> reduced(k);
    for (int i = 0; i < k; ++i) reduced[i] = res[i] % mod;
    return reduced;
}

// -------------------------------------------------------------------
// linearRecurrenceKitamasa(init, coeff, n, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Computes the n‑th term of a linear recurrence using Kitamasa's
//   algorithm (polynomial exponentiation) – faster than matrix method
//   for large k.
//
// INPUT:
//   init  : vector of length k, f[0] ... f[k-1]
//   coeff : vector of length k, recurrence coefficients (as defined above)
//   n     : index to compute (0‑based)
//   mod   : modulo value
//
// OUTPUT:
//   Returns f[n] modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(k^2 * log n)  (due to polynomial multiplication each step)
//
// CONSTRAINTS / PRECONDITIONS:
//   - k >= 1.
//   - n >= 0.
//   - Works for large k (e.g., k up to 500) within time limits.
//   - All arithmetic is modulo 'mod'.
ll linearRecurrenceKitamasa(const vector<ll>& init, const vector<ll>& coeff,
                            ll n, ll mod) {
    int k = (int)init.size();
    if (n < k) return init[n] % mod;

    // We want to compute x^n modulo characteristic polynomial.
    // Start with polynomial representing x^1.
    vector<ll> pol(k, 0);
    if (k == 1) {
        // For k=1, characteristic polynomial is x - c0, so x ≡ c0 (mod P).
        pol[0] = coeff[0] % mod;
    } else {
        pol[1] = 1; // x
    }
    // result polynomial starts as 1 (x^0)
    vector<ll> res(k, 0);
    res[0] = 1;

    ll exp = n;
    while (exp > 0) {
        if (exp & 1) {
            res = polyMulMod(res, pol, coeff, mod);
        }
        pol = polyMulMod(pol, pol, coeff, mod);
        exp >>= 1;
    }

    // Now res represents x^n mod P(x), i.e., res[i] = coefficient of x^i.
    // f[n] = sum_{i=0}^{k-1} res[i] * f[i]
    ll ans = 0;
    for (int i = 0; i < k; ++i) {
        ans = (ans + res[i] * (init[i] % mod)) % mod;
    }
    return ans;
}

// ===================================================================
// SECTION 4: GRAPH APPLICATIONS
// ===================================================================

// -------------------------------------------------------------------
// countWalks(adj, k, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Given an adjacency matrix of a directed/undirected graph, computes
//   the matrix (adj^k) where entry (i,j) is the number of walks of
//   exactly length k from node i to node j.
//
// INPUT:
//   adj : square matrix (n x n), where adj[i][j] = number of edges i->j
//   k   : length of walks (non‑negative integer)
//   mod : modulo value
//
// OUTPUT:
//   Returns the matrix adj^k modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(n^3 * log k)
//
// CONSTRAINTS / PRECONDITIONS:
//   - adj must be square.
//   - k >= 0. For k=0, the result is the identity matrix (walk of length 0).
Matrix countWalks(const Matrix& adj, ll k, ll mod) {
    return matPow(adj, k, mod);
}

// -------------------------------------------------------------------
// walksBetween(adj, u, v, k, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Returns the number of walks of exactly length k from node u to node v
//   in the graph described by adjacency matrix 'adj'.
//
// INPUT:
//   adj : square matrix (n x n)
//   u, v: 0‑based indices of nodes
//   k   : length of walk
//   mod : modulo
//
// OUTPUT:
//   Returns adj^k[u][v] modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(n^3 * log k)  (dominated by matPow)
//
// CONSTRAINTS / PRECONDITIONS:
//   - 0 <= u,v < n.
ll walksBetween(const Matrix& adj, int u, int v, ll k, ll mod) {
    Matrix p = matPow(adj, k, mod);
    return p[u][v] % mod;
}

// ===================================================================
// SECTION 5: DP WITH MATRIX EXPONENTIATION (Generic Transition)
// ===================================================================

// -------------------------------------------------------------------
// applyTransition(T, init, steps, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Applies a linear transition 'steps' times to an initial state vector.
//   State evolves as: state_{t+1} = T * state_t   (mod 'mod').
//
// INPUT:
//   T     : transition matrix of size (m x m)
//   init  : initial state vector of length m
//   steps : number of transitions to apply (non‑negative integer)
//   mod   : modulo value
//
// OUTPUT:
//   Returns the state vector after 'steps' applications: state_steps = T^steps * init.
//
// TIME COMPLEXITY:
//   O(m^3 * log steps)
//
// CONSTRAINTS / PRECONDITIONS:
//   - T must be square.
//   - init length must equal m.
//   - steps >= 0.
vector<ll> applyTransition(const Matrix& T, const vector<ll>& init,
                            ll steps, ll mod) {
    int m = (int)T.size();
    Matrix Tpow = matPow(T, steps, mod);
    return matVecMul(Tpow, init, mod);
}

// ===================================================================
// SECTION 6: SUM OF FIRST n TERMS OF LINEAR RECURRENCE
// ===================================================================

// -------------------------------------------------------------------
// sumLinearRecurrence(init, coeff, n, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Computes the sum of the first n terms of a linear recurrence:
//     S(n) = sum_{i=0}^{n-1} f[i]
//   where f follows the same recurrence as defined in linearRecurrence().
//
// INPUT:
//   init  : initial terms f[0] ... f[k-1] (length k)
//   coeff : recurrence coefficients (length k)
//   n     : number of terms to sum (n >= 0)
//   mod   : modulo value
//
// OUTPUT:
//   Returns S(n) modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(k^3 * log n)  using an augmented matrix of size (k+1)
//
// CONSTRAINTS / PRECONDITIONS:
//   - n >= 0.
//   - If n <= k, it computes the sum directly.
//   - All arithmetic is modulo 'mod'.
ll sumLinearRecurrence(const vector<ll>& init, const vector<ll>& coeff,
                       ll n, ll mod) {
    int k = (int)init.size();
    if (n == 0) return 0;
    if (n <= k) {
        ll s = 0;
        for (int i = 0; i < n; ++i) s = (s + init[i]) % mod;
        return s;
    }

    // Build augmented transition matrix of size (k+1) x (k+1)
    // State: [f[t], f[t-1], ..., f[t-k+1], S(t)]^T
    // where S(t) = sum_{i=0}^{t-1} f[i]
    Matrix T(k+1, vector<ll>(k+1, 0));
    // first row for f[t+1]
    for (int i = 0; i < k; ++i) T[0][i] = coeff[i] % mod;
    // shift rows
    for (int i = 1; i < k; ++i) T[i][i-1] = 1;
    // last row: S(t+1) = S(t) + f[t]  =>  T[k][0] = 1, T[k][k] = 1
    T[k][0] = 1;
    T[k][k] = 1;

    // Initial state at t = k-1:
    // state[0..k-1] = f[k-1], f[k-2], ..., f[0]
    vector<ll> state(k+1);
    for (int i = 0; i < k; ++i) state[i] = init[k-1-i] % mod;
    // S(k-1) = sum_{i=0}^{k-2} f[i]
    ll sum_init = 0;
    for (int i = 0; i < k-1; ++i) sum_init = (sum_init + init[i]) % mod;
    state[k] = sum_init;

    // exponent = n - (k-1)
    Matrix Tpow = matPow(T, n - k + 1, mod);
    vector<ll> res = matVecMul(Tpow, state, mod);
    return res[k] % mod;
}

// ===================================================================
// SECTION 7: ADDITIONAL TRICKS / PATTERNS
// ===================================================================

// -------------------------------------------------------------------
// (Note) Binary exponentiation for scalars:
//   Already available via std::pow? But not needed, we can use
//   our matPow with 1x1 matrix, or implement a fastPow function.
//   We'll provide a simple fastPow for completeness.
// -------------------------------------------------------------------

// fastPow(base, exp, mod): returns (base^exp) % mod.
// Use this when you need scalar exponentiation.
ll fastPow(ll base, ll exp, ll mod) {
    base %= mod;
    ll res = 1 % mod;
    while (exp > 0) {
        if (exp & 1) res = (res * base) % mod;
        base = (base * base) % mod;
        exp >>= 1;
    }
    return res;
}

// -------------------------------------------------------------------
// linearRecurrenceWithConstant(init, a, b, n, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Computes the n‑th term of a recurrence of the form:
//     f(0) = init
//     f(n) = a * f(n-1) + b   for n >= 1
//
// INPUT:
//   init : the initial value f(0)
//   a    : coefficient of f(n-1)
//   b    : constant term added each step
//   n    : index to compute (n >= 0)
//   mod  : modulo value
//
// OUTPUT:
//   Returns f(n) modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(log n)  (using 2x2 matrix exponentiation)
//
// CONSTRAINTS / PRECONDITIONS:
//   - n >= 0.
//   - All values should be non‑negative or reduced modulo 'mod'.
ll linearRecurrenceWithConstant(ll init, ll a, ll b, ll n, ll mod) {
    if (n == 0) return init % mod;

    // State vector: [f(t), 1]^T
    // Transition: [f(t+1)]   = [a  b] * [f(t)]
    //             [   1   ]     [0  1]   [  1 ]
    Matrix T = {{a % mod, b % mod}, {0, 1}};
    vector<ll> state = {init % mod, 1};

    vector<ll> res = applyTransition(T, state, n, mod);
    return res[0];
}

// -------------------------------------------------------------------
// twoSequences(a, b, n, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Computes the n‑th terms of two sequences defined as:
//     x(0) = 1, y(0) = 0
//     x(n) = p * x(n-1) + q * y(n-1)
//     y(n) = r * x(n-1) + s * y(n-1)
//   This is a generic example; you can change the coefficients and
//   initial values as needed.
//
// INPUT:
//   p, q, r, s : coefficients of the recurrence
//   n           : index to compute (n >= 0)
//   mod         : modulo value
//
// OUTPUT:
//   Returns a pair {x(n), y(n)} modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(log n)  (using 2x2 matrix exponentiation)
//
// CONSTRAINTS / PRECONDITIONS:
//   - n >= 0.
//   - Works for any integer coefficients.
pair<ll, ll> twoSequences(ll p, ll q, ll r, ll s, ll n, ll mod) {
    if (n == 0) return {1 % mod, 0};

    // State: [x(t), y(t)]^T
    // Transition: [x(t+1)]   = [p  q] * [x(t)]
    //              [y(t+1)]     [r  s]   [y(t)]
    Matrix T = {{p % mod, q % mod}, {r % mod, s % mod}};
    vector<ll> state = {1 % mod, 0};

    vector<ll> res = applyTransition(T, state, n, mod);
    return {res[0], res[1]};
}

// -------------------------------------------------------------------
// sumFirstNFibonacci(n, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Computes the sum of the first n Fibonacci numbers:
//     S(n) = F(0) + F(1) + ... + F(n-1)
//   where F(0)=0, F(1)=1.
//
// INPUT:
//   n   : number of terms to sum (n >= 0)
//   mod : modulo value
//
// OUTPUT:
//   Returns S(n) modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(log n)  (using 3x3 matrix exponentiation)
//
// CONSTRAINTS / PRECONDITIONS:
//   - n >= 0.
//   - If n == 0, returns 0.
ll sumFirstNFibonacci(ll n, ll mod) {
    if (n == 0) return 0;

    // State: [F(t), F(t-1), S(t)]^T  where S(t) = sum_{i=0}^{t-1} F(i)
    // Transition for t >= 1:
    // [F(t+1)]   = [1  1  0] * [F(t)]
    // [F(t)  ]     [1  0  0]   [F(t-1)]
    // [S(t+1)]     [1  1  1]   [S(t)  ]
    Matrix T = {{1, 1, 0}, {1, 0, 0}, {1, 1, 1}};
    // Initial state at t = 1: F(1)=1, F(0)=0, S(1)=sum F(0)=0
    vector<ll> state = {1 % mod, 0, 0};

    vector<ll> res = applyTransition(T, state, n - 1, mod);
    return res[2];
}

// -------------------------------------------------------------------
// countArraysNoThreeConsecutiveEqual(M, n, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Counts the number of arrays of length n, where each element is
//   between 1 and M (inclusive), and no three consecutive elements
//   are equal.
//
// INPUT:
//   M   : maximum value of each element (M >= 1)
//   n   : length of the array (n >= 1)
//   mod : modulo value
//
// OUTPUT:
//   Returns the number of valid arrays modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(M^3 * log n)  (matrix size is M x M)
//
// CONSTRAINTS / PRECONDITIONS:
//   - Works for small M (e.g., M <= 50) because matrix is M x M.
//   - If n is small (1 or 2), the answer is computed directly.
//   - This is a classic example; for larger M, you need Kitamasa.
ll countArraysNoThreeConsecutiveEqual(int M, ll n, ll mod) {
    if (n == 1) return M % mod;
    if (n == 2) return (M * M) % mod;

    // State: dp[len][last][prev] but we can compress.
    // Better approach: Use a 2x2 matrix for this specific problem
    // because the state can be: (number of ways ending with two equal,
    // number of ways ending with two different).
    // But to demonstrate large matrix, we'll keep it general.

    // For M up to 50, we can use a 2x2 matrix:
    // State: [ways where last two are equal, ways where last two are different]
    // Transition:
    // new_equal = old_different * 1  (choose the same as last)
    // new_different = old_equal * (M-1) + old_different * (M-2)
    // Because from equal state, you must choose a different element (M-1 choices)
    // from different state, you choose an element different from the last (M-2 choices)

    Matrix T = {{0, 1}, {M-1, M-2}};
    // Initial state for length 2:
    // equal ways = M (pairs like (1,1), (2,2), ...)
    // different ways = M * (M-1)
    vector<ll> state = {M % mod, (M * (M-1)) % mod};

    vector<ll> res = applyTransition(T, state, n - 2, mod);
    return (res[0] + res[1]) % mod;
}

// -------------------------------------------------------------------
// matMulMinPlus(A, B)
// -------------------------------------------------------------------
// PURPOSE:
//   Multiplies two matrices using the min-plus (or tropical) semiring:
//     C[i][j] = min_k ( A[i][k] + B[k][j] )
//   This is used to compute shortest paths after k steps, where
//   the graph has edge weights and you want the minimum total weight
//   of a walk of exactly k edges.
//
// INPUT:
//   A, B : square matrices of the same size, containing edge weights.
//          Use INF (a large number) for no edge.
//
// OUTPUT:
//   Returns the min-plus product matrix C.
//
// TIME COMPLEXITY:
//   O(n^3)
//
// CONSTRAINTS / PRECONDITIONS:
//   - Matrices must be square and of the same size.
//   - INF should be large enough (e.g., 4e18) to avoid overflow.
const ll INF = 4e18;

Matrix matMulMinPlus(const Matrix& A, const Matrix& B) {
    int n = (int)A.size();
    Matrix C(n, vector<ll>(n, INF));
    for (int i = 0; i < n; ++i) {
        for (int k = 0; k < n; ++k) {
            if (A[i][k] == INF) continue;
            for (int j = 0; j < n; ++j) {
                if (B[k][j] == INF) continue;
                C[i][j] = min(C[i][j], A[i][k] + B[k][j]);
            }
        }
    }
    return C;
}

// -------------------------------------------------------------------
// matPowMinPlus(base, exp)
// -------------------------------------------------------------------
// PURPOSE:
//   Raises a square matrix to the power 'exp' using min-plus
//   multiplication. This computes the minimum weight of a walk of
//   exactly 'exp' edges between any two nodes.
//
// INPUT:
//   base : square matrix of edge weights (INF for no edge)
//   exp  : number of edges (exp >= 0)
//
// OUTPUT:
//   Returns base^exp under min-plus multiplication.
//
// TIME COMPLEXITY:
//   O(n^3 * log exp)
//
// CONSTRAINTS / PRECONDITIONS:
//   - exp >= 0.
//   - For exp = 0, the result is the identity matrix for min-plus:
//     C[i][i] = 0, C[i][j] = INF for i != j.
Matrix matPowMinPlus(Matrix base, ll exp) {
    int n = (int)base.size();
    Matrix res(n, vector<ll>(n, INF));
    for (int i = 0; i < n; ++i) res[i][i] = 0; // identity for min-plus

    while (exp > 0) {
        if (exp & 1) res = matMulMinPlus(res, base);
        base = matMulMinPlus(base, base);
        exp >>= 1;
    }
    return res;
}

// -------------------------------------------------------------------
// PrecomputedPowers
// -------------------------------------------------------------------
// PURPOSE:
//   Precomputes powers of a matrix: P[i] = base^(2^i).
//   Then, for any exponent k, you can compute base^k by multiplying
//   the relevant P[i] matrices.
//
// INPUT:
//   base : square matrix
//   maxExp : maximum exponent you will query (so we precompute up to log2(maxExp))
//   mod : modulo value
//
// OUTPUT:
//   The class provides a method query(exp) that returns base^exp.
//
// TIME COMPLEXITY:
//   Precomputation: O(n^3 * log maxExp)
//   Each query: O(n^3 * popcount(exp)) which is O(n^3 * log maxExp) in worst case.
//
// CONSTRAINTS / PRECONDITIONS:
//   - base must be square.
//   - maxExp >= 0.
class PrecomputedPowers {
private:
    vector<Matrix> powers;
    ll mod;

public:
    PrecomputedPowers(const Matrix& base, ll maxExp, ll mod) : mod(mod) {
        powers.push_back(base);
        for (ll e = 2; e <= maxExp; e <<= 1) {
            powers.push_back(matMul(powers.back(), powers.back(), mod));
        }
    }

    Matrix query(ll exp) {
        int n = (int)powers[0].size();
        Matrix res(n, vector<ll>(n, 0));
        for (int i = 0; i < n; ++i) res[i][i] = 1 % mod;

        int bit = 0;
        while (exp > 0) {
            if (exp & 1) {
                res = matMul(res, powers[bit], mod);
            }
            exp >>= 1;
            bit++;
        }
        return res;
    }
};

// -------------------------------------------------------------------
// matMulSparse(A, B, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Multiplies two sparse matrices efficiently by skipping zeros.
//   A and B are represented as vector of vectors of pairs (col, value).
//
// INPUT:
//   A, B : sparse matrices, each a vector of size n, where A[i] is a
//          vector of pairs {j, value} for non-zero entries in row i.
//   mod  : modulo value
//
// OUTPUT:
//   Returns the product matrix C as a dense matrix (n x n).
//
// TIME COMPLEXITY:
//   O( (number of non-zero in A) * (average non-zero per row in B) )
//   Much faster than O(n^3) if matrices are sparse.
//
// CONSTRAINTS / PRECONDITIONS:
//   - Matrices must be square.
//   - The input format is sparse; the output is dense.
Matrix matMulSparse(const vector<vector<pair<int, ll>>>& A,
                    const vector<vector<pair<int, ll>>>& B, ll mod) {
    int n = (int)A.size();
    Matrix C(n, vector<ll>(n, 0));

    for (int i = 0; i < n; ++i) {
        for (auto &p : A[i]) {
            int k = p.first;
            ll aik = p.second;
            if (aik == 0) continue;
            for (auto &q : B[k]) {
                int j = q.first;
                C[i][j] = (C[i][j] + aik * q.second) % mod;
            }
        }
    }
    return C;
}

// -------------------------------------------------------------------
// matrixGeometricSum(A, k, mod)
// -------------------------------------------------------------------
// PURPOSE:
//   Computes the sum: S = A + A^2 + A^3 + ... + A^k
//
// INPUT:
//   A   : square matrix
//   k   : number of terms (k >= 1)
//   mod : modulo value
//
// OUTPUT:
//   Returns the matrix S modulo 'mod'.
//
// TIME COMPLEXITY:
//   O(n^3 * log k)
//
// CONSTRAINTS / PRECONDITIONS:
//   - A must be square.
//   - k >= 1.
Matrix matrixGeometricSum(const Matrix& A, ll k, ll mod) {
    int n = (int)A.size();
    // Build augmented matrix of size 2n x 2n:
    // [A  I]
    // [0  I]
    Matrix T(2*n, vector<ll>(2*n, 0));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            T[i][j] = A[i][j] % mod;          // top-left: A
            T[i][j+n] = (i == j) ? 1 : 0;     // top-right: I
        }
        T[i+n][i+n] = 1;                      // bottom-right: I
    }

    // We need sum_{i=1}^{k} A^i = (sum_{i=0}^{k} A^i) - I.
    // T^(k+1) top-right block = sum_{i=0}^{k} A^i.
    Matrix Tk = matPow(T, k + 1, mod);

    Matrix S(n, vector<ll>(n, 0));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            S[i][j] = Tk[i][j+n] % mod;
            if (i == j) {
                S[i][j] = (S[i][j] - 1 + mod) % mod; // subtract I
            }
        }
    }
    return S;
}


// -------------------------------------------------------------------
// (Note) Exponentiating a matrix with non‑integer exponents? Not used.
// -------------------------------------------------------------------

// ===================================================================
// SECTION 8: COMMON TERMS EXPLAINED
// ===================================================================
//
// 1. Matrix Exponentiation:
//    Raising a square matrix to a power using binary exponentiation.
//    Used to accelerate linear recurrences and DP transitions that can
//    be expressed as repeated multiplication by a fixed matrix.
//
// 2. Linear Recurrence:
//    A sequence where each term is a linear combination of previous terms.
//    Example: Fibonacci, Tribonacci, etc.
//
// 3. Transition Matrix / Companion Matrix:
//    A matrix that transforms the state vector from time t to t+1.
//    For a recurrence of order k, the companion matrix is k x k.
//
// 4. State Vector:
//    A column vector containing the current values needed to compute the
//    next values (e.g., last k terms).
//
// 5. Kitamasa / Polynomial Exponentiation:
//    A technique to compute the n‑th term of a linear recurrence without
//    building the full matrix, by computing x^n modulo the characteristic
//    polynomial. Runs in O(k^2 log n) and is better for large k.
//
// 6. Characteristic Polynomial:
//    For recurrence f[n] = c0 f[n-1] + c1 f[n-2] + ... + c[k-1] f[n-k],
//    the characteristic polynomial is P(x) = x^k - c0 x^(k-1) - ... - c[k-1].
//
// 7. Modulo:
//    All operations are performed modulo a given number to prevent overflow
//    and keep numbers small. In competitive programming, MOD is often 1e9+7
//    or 998244353.
//
// ===================================================================

// ===================================================================
// EXAMPLE USAGE (you can ignore or uncomment to test)
// ===================================================================

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    const ll MOD = 1000000007;

    // Fibonacci
    cout << fib(10, MOD) << "\n";            // 55

    // Linear recurrence: Fibonacci (init [0,1], coeff [1,1])
    vector<ll> init = {0, 1};
    vector<ll> coeff = {1, 1};
    cout << linearRecurrence(init, coeff, 10, MOD) << "\n";  // 55
    cout << linearRecurrenceKitamasa(init, coeff, 10, MOD) << "\n"; // 55

    // Sum of first 5 Fibonacci numbers: 0+1+1+2+3 = 7
    cout << sumLinearRecurrence(init, coeff, 5, MOD) << "\n"; // 7

    // Graph walks: simple directed graph with 2 nodes, edge 0->1 and 1->0
    Matrix adj = {{0,1},{1,0}};
    Matrix p = countWalks(adj, 3, MOD);
    // p[0][1] = number of walks of length 3 from 0 to 1.
    cout << p[0][1] << "\n"; // 1? Actually for two nodes, length 3 from 0 to 1: 0->1->0->1 => 1, 0->1->0->? Only one.

    return 0;
}