矩阵快速幂

看高中信息竞赛 数学部分

https://www.cnblogs.com/dx123/p/16669615.html

题单

做前三题：
http://218.5.5.242:9018/JudgeOnline/status.php?user_id=24B13&cid=1973

题目:

有一段长为n的阶梯，最多可以往上跨k步，请求出走到第n级阶梯的总方案数mod 7777777的值。输入
两个正整数n,k(n<=2^31-1,k<=10)

输出
总方案数 mod 7777777的值

样例输入
4
2
样例输出
5

算法思路

我们令 \(f[i]\) 表示走到第 \(i\) 级台阶的方案数。规则是每一步可以跨 \(1\) 到 \(k\) 级，所以满足递推：

\[f[i]=f[i-1]+f[i-2]+\dots+f[i-k] \]

并且 \(f[0]=1\) （不动也算一种方法），对于 \(i<0\) 认为 \(f[i]=0\) 。

直接递推到 \(n\) 会太慢（ \(n\) 可达约 \(2^{31}-1\) ），但这是线性常系数递推，所以可以用矩阵快速幂求解。构造 \(k\times k\) 的转移矩阵 \(M\) ，使得向量

\[V_t=[f[t],f[t-1],\dots,f[t-k+1]]^T \]

满足 \(V_{t+1}=M\cdot V_t\) 。

1. 当 \(n < k\) ，直接递推计算。
2. 当 \(n \ge k\) ，用矩阵快速幂：
   \[V_t = [f[t], f[t-1], \dots, f[t-k+1]]^T \]
   \[V_{t+1} = M \cdot V_t \]
   其中转移矩阵：
   \[M = \begin{bmatrix} 1 & 1 & 1 & \dots & 1\\ 1 & 0 & 0 & \dots & 0\\ 0 & 1 & 0 & \dots & 0\\ \vdots & & \ddots & & \vdots\\ 0 & 0 & \dots & 1 & 0 \end{bmatrix} \]

矩阵 \(M\) 的第一行全为 1，其余为向下移位（第 \(i\) 行第 \(i-1\) 列为1，其余为0）。于是

\[V_n = M^{\,n-(k-1)} \cdot V_{k-1} \]

（当 \(n\ge k-1\) 时）。如果 \(n<k\) ，我们直接用前面的小范围 dp 结果返回 \(f[n]\) 。

时间复杂度：矩阵乘法 \(O(k^3)\) ，二分幂要做 \(O(\log n)\) 次，合计 \(O(k^3\log n)\) 。空间 \(O(k^2)\) 。模数是 7777777，注意用 64 位中间乘法。

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

const int MAXK = 12;     // k ≤ 10
const ll MOD = 7777777;

// ---------------------------
// 矩阵结构体（封装在 struct 中）
// ---------------------------
struct Matrix {
    int n;                // 矩阵维度
    ll a[MAXK][MAXK];    

    Matrix(int _n = 0, bool ident = false) {
        n = _n;
        memset(a, 0, sizeof(a));
        if (ident) {
            for (int i = 0; i < n; ++i) a[i][i] = 1;
        }
    }

    // 矩阵乘法（重载 * 运算符）
    Matrix operator*(const Matrix &B) const {
        Matrix C(n);
        for (int i = 0; i < n; ++i) 
            for (int j = 0; j < n; ++j)
                for (int k = 0; k < n; ++k)
                    C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;    
        return C;
    }
};

// ---------------------------
// 矩阵快速幂
// ---------------------------
Matrix mat_pow(Matrix base, long long e) {
    Matrix res(base.n, true);  // 单位矩阵
    while (e > 0) {
        if (e & 1) res = res * base;
        base = base * base;
        e >>= 1;
    }
    return res;
}


int main() {
    long long n;
    int k;
    cin >> n >> k;

    ll f[MAXK]; // 存储前 k 项
    f[0] = 1;
    for (int i = 1; i < k; ++i) {
        f[i] = 0;
        for (int j = 1; j <= k && i - j >= 0; ++j)
            f[i] = (f[i] + f[i - j]) % MOD;
    }

    if (n < k) {
        cout << f[n] % MOD << '\n';
        return 0;
    }

    // 构造转移矩阵 M
    Matrix M(k);
    for (int j = 0; j < k; ++j) M.a[0][j] = 1;   // 第一行全为 1
    for (int i = 1; i < k; ++i) M.a[i][i - 1] = 1; // 下移部分

    // 快速幂求 M^(n - (k-1))
    Matrix P = mat_pow(M, n - (k - 1));

    // 初始列向量 V = [f[k-1], f[k-2], ..., f[0]]
    reverse(f, f+k);
    ll ans = 0;
    for (int j = 0; j < k; ++j)
        ans = (ans + P.a[0][j] * f[j]) % MOD;

    cout << ans % MOD << '\n';
    return 0;
}

题目描述

a[1]=a[2]=a[3]=1

a[x]=a[x-3]+a[x-1]
求a数列的第n项对1000000007（10^9+7）取余的值。

输入
第一行一个整数T，表示询问个数。

以下T行，每行一个正整数n。

输出
每行输出一个非负整数表示答案。

样例输入
3
6
8
10
样例输出
4
9
19

思路

定义数列满足：
\(a_1=a_2=a_3=1\) ，且对 \(n\ge4\) 有

\[a_n = a_{n-1} + a_{n-3}. \]

把它写成向量形式：令

\[v_n = \begin{bmatrix} a_n \\ a_{n-1} \\ a_{n-2} \end{bmatrix}. \]

可以找到一个 \(3\times3\) 的转移矩阵 \(A\) ，使得

\[v_n = A \cdot v_{n-1}, \]

其中

\[A = \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}. \]

因此对任意 \(n\ge3\) ：

\[v_n = A^{\,n-3} \cdot v_3, \qquad v_3=\begin{bmatrix}1\\1\\1\end{bmatrix}. \]

我们只需用矩阵快速幂计算 \(A^{\,n-3}\) ，然后乘以 \(v_3\) ，取第一分量就是 \(a_n\) （对 \(10^9+7\) 取模）。

时间复杂度：每次查询 \(O(3^3\log n)=O(\log n)\) 。

---

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const ll MOD = 1000000007LL;

// OI 风格矩阵快速幂：处理 3x3 矩阵
struct Matrix {
    ll a[3][3];
    // 构造为零矩阵
    Matrix() { memset(a, 0, sizeof(a)); }
    // 生成单位矩阵
    static Matrix identity() {
        Matrix I;
        for (int i = 0; i < 3; ++i) I.a[i][i] = 1;
        return I;
    }
    // 矩阵乘法（模 MOD）
    Matrix operator*(const Matrix &other) const {
        Matrix res;
        // 3x3 矩阵乘法，常数较小，写成三重循环
        for (int i = 0; i < 3; ++i) {
            for (int k = 0; k < 3; ++k) if (a[i][k]) {
                ll aik = a[i][k];
                for (int j = 0; j < 3; ++j) {
                    res.a[i][j] = (res.a[i][j] + aik * other.a[k][j]) % MOD;
                }
            }
        }
        return res;
    }
};

// 快速幂计算 base^exp
Matrix matrix_pow(Matrix base, long long exp) {
    Matrix res = Matrix::identity();
    while (exp > 0) {
        if (exp & 1LL) res = res * base;
        base = base * base;
        exp >>= 1LL;
    }
    return res;
}

int main() {
    int T;
    cin >> T;
    while (T--) {
        long long n;
        cin >> n;
        // 边界：n = 1,2,3
        if (n <= 3) {
            cout << 1 << '\n';
            continue;
        }

        // 定义转移矩阵 A
        Matrix A;
        // A = [[1,0,1],
        //      [1,0,0],
        //      [0,1,0]]
        A.a[0][0] = 1; A.a[0][1] = 0; A.a[0][2] = 1;
        A.a[1][0] = 1; A.a[1][1] = 0; A.a[1][2] = 0;
        A.a[2][0] = 0; A.a[2][1] = 1; A.a[2][2] = 0;

        // 计算 A^(n-3)
        Matrix P = matrix_pow(A, n - 3);

        // v3 = [a3, a2, a1] = [1,1,1]
        // v_n = P * v3，我们只需第一个分量：
        // a_n = P.row0 * v3 = P.a[0][0]*1 + P.a[0][1]*1 + P.a[0][2]*1
        ll an = (P.a[0][0] + P.a[0][1] + P.a[0][2]) % MOD;
        cout << an << '\n';
    }
    return 0;
}