神奇的矩阵乘法

2021-01-02

神奇的矩阵乘法

矩阵快速幂

啊其实非常简单只需要在实数的快速幂基础上，重新定义针对矩阵的乘法即可。看代码。

#include <iostream>

using namespace std;
long long n, k;
const int mod = 1000000007;
struct Matrix{
    long long m[105][105];
}A, I; // I表示单位阵
//针对 struct Matrix 重载运算符*
Matrix operator*(const Matrix &a, const Matrix &b){
    Matrix tmp;
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            tmp.m[i][j] = 0;
        }
    }
    //很normal的矩阵相乘
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            for(int k = 1; k <= n; k++){
                tmp.m[i][j] += a.m[i][k] * b.m[k][j] % mod;
                tmp.m[i][j] %= mod;
            }
        }
    }
    return tmp;
}
int main()
{
    cin >> n >> k;
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++)
            cin >> A.m[i][j];
    }
    for(int i = 1; i <= n; i++){
        I.m[i][i] = 1;
    }

    while(k > 0){
        if(k % 2){
            I = I * A;
        }
        A = A * A;
        k = k >> 1;
    }
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++){
            cout << I.m[i][j] << " ";
        }
        cout << endl;
    }
    return 0;
}

矩阵乘法的大用处

跳台阶很熟悉了，一般无非给定一个台阶数 $n$，一次最多能跳 $k$ 次，然后递推方程dp[n] = dp[n - 1] + dp[n - 2] + ... + dp[n - k]，当 $n$ 不是很大的时候我们可以直接开数组算出来，但是当 $n$ 很大（$1 to 2^{31}-1$），直接开数组直接内存爆炸，而且 $k$ 也不是固定的时候，该怎么求？这就是这个题：https://vijos.org/p/1067

参考链接：http://www.matrix67.com/blog/archives/276

举个例子当 $k = 3$ 的时候，$n \le3$是初始状态，当$n > 3$时候初始状态

$\pmb x = (1,2,3)^T$

继续递推：

$\pmb x = (2,3,5)^T,(3,5,8)^T,...$

向量的最后一个元素代表到达$n$的方案数，我们需要构造一个矩阵 $M$ ，使得每 $\pmb x$ 与 $M$ 相乘一次，就相当于一次递推。我们观察上面发现，与矩阵 $M$ 相乘后，第$k$个元素是向量前$k-1$个元素的和，然后$1-(k-1)$个元素是原向量的 $2 - k$ 个元素。由此得出当$k = 3$时，

$M = \begin{pmatrix}0&1&0\\0&0&1\\1&1&1\end{pmatrix}$

当$k = 4$时，

$M = \begin{pmatrix}0&1&0&0\\0&0&1&0\\0&0&0&1\\1&1&1&1\end{pmatrix}$

可知最终答案为：

$\pmb y_k,\pmb y = M^{n-k}\pmb x$

#include <iostream>

using namespace std;
const int mod = 7777777;
long long k,n,dp[15], ans[15];
struct Matrix{
    long long mat[12][12];
}A,I;
// 重载运算符
Matrix operator*(const Matrix &a, const Matrix &b){
    Matrix tmp;
    for(int i = 1; i <= k; i++){
        for(int j = 1; j <= k; j++)
            tmp.mat[i][j] = 0;
    }
    // 很normal的矩阵相乘
    for(int i = 1; i <= k; i++){
        for(int j = 1; j <= k; j++){
            for(int q = 1; q <= k; q++){
                tmp.mat[i][j] += a.mat[i][q]*b.mat[q][j] % mod;
                tmp.mat[i][j] %= mod;
            }
        }
    }
    return tmp;
}

int main()
{
    cin >> k >> n;
    // 构造初始矩阵,大小:k*k
    // 最后一行最为1 右上角(k-1)*(k-1)的单位矩阵
    for(int i = 1; i <= k; i++){
        A.mat[k][i] = 1;
        A.mat[i][i + 1] = 1;
        I.mat[i][i] = 1;
    }
    // 初始化dp数组
    dp[0] = 1;
    int sum = 1;
    for(int i = 1; i <= k; i++){
        dp[i] = sum;
        sum += dp[i];
    }
    //矩阵快速幂 n-k 次幂,
    int p = n - k;
    while(p > 0){
        if(p % 2){
            I = I * A;
        }
        A = A * A;
        p = p >> 1;
    }

    //然后再与dp向量相乘
    int row = 1;
    for(int i = 1; i <= k; i++){
        for(int j = 1; j <= k; j++){
           ans[row] += I.mat[i][j] * dp[j] % mod;
           ans[row] %= mod;
        }
        row++;
    }
    cout << ans[k] << endl;
    return 0;
}