矩阵快速幂模板（解决斐波那契数列和HDU 5015 233 Matrix）

$a^{156}=a^{10011100_2}=a^{2^7\ast 1+2^6\ast 0+2^5\ast 0+2^4\ast 1+2^3\ast 1+2^2\ast 1+2^1\ast 0+2^0\ast 0}$

普通快速幂的模板：

int fastPow(int base,int n,int mod){
    int ans=1;//单位元
    while(n){
        if(n&1) ans=(ans*base)%mod;
        base=(base*base)%mod;//ccpc的痛
        n=n>>1;
    }
    return ans%mod;
}

base的变化: $bas{e}^{2^0}$ $bas{e}^{2^1}$ $bas{e}^{2^2}$ $bas{e}^{2^3}$ $bas{e}^{2^4}$

正好对应n的二进制各位。

矩阵快速幂模板：（矩阵大小依据题目给定）

//矩阵乘法
vector<vector<int>> matrixMult(vector<vector<int>> &a,vector<vector<int>> &b){
    int len=a.size();//这里定义的矩阵乘法其实就是方阵，大小是len
    vector<vector<int>> c(len,vector<int>(len,0));
    for(int i=0;i<len;i++){//三重循环解决普通矩阵乘法
        for(int j=0;j<len;j++){
            for(int k=0;k<len;k++){
                c[i][j]+=a[i][k]*b[k][j]%mod;
            }
        }
    }
    return c;
}
//矩阵快速幂运算
vector<vector<int>> matrixPow(vector<vector<int>> &base,int n){
    int len=base.size();
    vector<vector<int>> ans(len,vector<int>(len,0));
    for(int i=0;i<len;i++){//ans首先定义为单位阵
        for(int j=0;j< len;j++)
            if(i==j) ans[i][j]=1;
    }
    while(n){
        if(n&1) ans=matrixMult(ans,base);//和普通的快速幂的区别就在于乘法是矩阵乘法
        base=matrixMult(base,base);
        n=n>>1;
    }
    return ans;
}

在实际应用中，最关键的是寻找关系矩阵base，毕竟上述板子只是加速幂运算的，需要在实际问题中挖掘出其中的矩阵幂运算。

应用举例：

1. 解决斐波那契数列

（主要参考：https://leetcode-cn.com/problems/fibonacci-number/solution/fei-bo-na-qi-shu-by-leetcode-solution-o4ze/）

​ 斐波那契数，通常用 F(n) 表示，形成的序列称为斐波那契数列 。该数列由 0 和 1 开始，后面的每一项数字都是前面两项数字的和。也就是：

给定 n 计算 F(n) 。

构造递推关系：

因此有：

令：

那么可以化简得到：

也就是$$F(n)=bas{e}^{n-1}[0][0]$$ 所以只需要取base矩阵的n-1次幂运算完成后的第一行第一列的元素即为结果。完整代码如下：

void fib(int n){
    if(n<2) cout << n;
    vector<vector<int>> base{{1,1},{1,0}};
    vector<vector<int>> res=matrixPow(base,n-1);
    cout << res[0][0];
}
//和一般的板子的区别在于这里base是二阶方阵，定义运算的时候直接展开简洁一点
vector<vector<int>> matrixMult(vector<vector<int>> &a, vector<vector<int>> &b){
    vector<vector<int>> c{{0,0},{0,0}};
    for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
            c[i][j]=a[i][0]*b[0][j]+a[i][1]*b[1][j];
        }
    }
    return c;
}
vector<vector<int>> matrixPow(vector<vector<int>> &base , int n){
    vector<vector<int>> ans{{1,0},{0,1}};//单位阵（单位元）
    while(n){
        if(n&1) ans=matrixMult(ans,base);
        base=matrixMult(base,base);
        n=n>>1;
    }
    return ans;
}

2. 233 Matrix

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=5015

(主要参考：https://zhuanlan.zhihu.com/p/42639682)


题意就是给定一个矩阵的第一行，然后输入第一列，依据递推关系$$a_{i,j}=a_{i-1,j}+a_{i,j-1}$$

可以把整个矩阵确定，要求的是$$a_{n,m}$$

可以观察到第一行233，2333，2333的变化规律其实就是前一个数字乘10加3，即$$a_{0,j}=a_{0,j-1}\times 10+3$$, 而 $$233=23\ast 10+3$$ , 我们可以把上述矩阵改写成下面的形式：

这样第二列的元素实际可以由第一列的元素左乘一个关系矩阵得到：

而第三列的元素也可以继续由第二列的元素左乘该关系矩阵得到，于是就找到了这个递推关系，最后一列实际就是：

完整代码如下：

一个投机取巧的地方是#define int long long， 这样就不用手动把所有的int换成long long了，但是因为main()函数不能返回long long类型，因此将主函数返回类型改成了signed.

另一个需要注意的地方是在一切可能超范围的位置都要想着%一下。

res=(res + a[i]*basePow[n][i]%mod) %mod;

#include <iostream>
#include<cstdio>
#include<bits/stdc++.h>
#define int long long
#define mod 10000007
using namespace std;
void showMatrix(vector<vector<int>> &tmp){
    int len=tmp.size();
    for(int i=0;i<len;i++){
        for(int j=0;j<len;j++)
            cout << tmp[i][j] << " ";
        cout <<endl;
    }
}
void initBase(vector<vector<int>> &base,int n){//构造关系矩阵
    for(int i=0;i<=n+1;i++){
        base[i][n+1]=1;
        if(i==n+1) continue;
        for(int j=0;j<=i;j++)
            base[i][j]=1;
        base[i][0]=10;
    }
}
//定义矩阵乘法
vector<vector<int>> matrixMult(vector<vector<int>> &a,vector<vector<int>> &b){
    int len=a.size();//这里定义的矩阵乘法其实就是方阵，大小是len
    // 此题中有len=n+2,因为添加了两行两列用于辅助计算
    vector<vector<int>> c(len,vector<int>(len,0));
    for(int i=0;i<len;i++){
        for(int j=0;j<len;j++){
            for(int k=0;k<len;k++){
                c[i][j]+=a[i][k]*b[k][j]%mod;
            }
        }
    }
    //cout << "c" << endl;
    //showMatrix(c);
    return c;
}
//定义矩阵快速幂运算
vector<vector<int>> matrixPow(vector<vector<int>> &base,int n){
    int len=base.size();
    vector<vector<int>> ans(len,vector<int>(len,0));
    for(int i=0;i<len;i++){//ans首先定义为单位阵
        for(int j=0;j< len;j++)
            if(i==j) ans[i][j]=1;
    }
    //showMatrix(ans);
    while(n){
        if(n&1) ans=matrixMult(ans,base);
        base=matrixMult(base,base);
        n=n>>1;
    }
    return ans;
}

signed main() {
    int n,m;
    int a[12];
    while(scanf("%d %d",&n,&m)!=EOF){
        for(int i=0;i<n;i++){
            scanf("%d",&a[i+1]);
        }
        a[0]=23;
        a[n+1]=3;
        //构造关系矩阵
        vector<vector<int>> base(n+2,vector<int>(n+2,0));
        vector<vector<int>> basePow(n+2,vector<int>(n+2,0));
        initBase(base,n);
        //cout << "initMatrix" << endl;
        //showMatrix(base);
        basePow=matrixPow(base,m);
        //showMatrix(basePow);
        int res=0;
        for(int i=0;i<=n+1;i++){
            res=(res + a[i]*basePow[n][i]%mod) %mod;
        }
        cout << res << endl;
    }
    return 0;
}