【笔记】矩阵快速幂

前置芝士

快速幂。

什么是矩阵？

矩阵，是由 \(\begin{bmatrix}\end{bmatrix}\) 组成的一个方阵（就这么理解好啦）。

比如：\(\begin{bmatrix}1& 2\\3&4\end{bmatrix}\) 是一个 \(2\times 2\) 的矩阵。

矩阵乘法

矩阵乘法的条件：仅当第 \(1\) 个矩阵的列数 \(=\) 第 \(2\) 个矩阵的行数才有意义。

U401562 矩阵乘法（站外题）

矩阵乘法的计算方法：

\[c_{i,j} = \sum_{k=1}^{n} a_{i,k}\times b_{k,j} \]

了解了这些以后，代码就迎刃而解：

第一层循环枚举 \(k\)，第二、三层循环分别枚举 \(i\)，\(j\)。

for(k=1;k<=n;k++)
    for(i=1;i<=m;i++)
        for(j=1;j<=p;j++)
            c[i][j]+=a[i][k]*b[k][j];

就是 \(A\) 的第 \(i\) 行分别乘以 \(B\) 的第 \(j\) 列的和就是 \(C_{i,j}\)。

这样就可以得到一个 \(m\times p\) 的矩阵（第一个矩阵的行 \(\times\) 第二个矩阵的列）。

---

矩阵快速幂

不懂快速幂的可以见顶部。

矩阵快速幂，用于解决 \(n\) 很大的情况下求解某些东西。

比如求斐波那契数列的第 \(n\) 项。我们知道，斐波那契数列有一个通项公式：

\[\begin{cases}f_i = 1&i\le 3\\f_i = f_{i-1} + f_{i-2}&i\ge 3\\\end{cases} \]

• 首先特判 \(n=1,2\) 的情况。

将上面的递推公式变为矩阵：

\[\begin{bmatrix}?&?\\?&?\end{bmatrix}\times \begin{bmatrix}f_{i-1}\\f_{i-2}\\\end{bmatrix} = \begin{bmatrix}f_i\\f_{i-1}\end{bmatrix} \]

矩阵 \(2\) 是一个 \(2\times 1\) 的矩阵，矩阵 \(3\) 也是一个 \(2\times 1\) 的矩阵，因此矩阵 \(1\) 是一个 \(2\times 2\) 的矩阵。

然后通过递推公式 \(f_i = f_{i-1}\times 1+f_{i-2}\times 1\) 和 \(f_{i-1} = f_{i-1}\times 1\)，可构造矩阵：

\[\begin{bmatrix}1&1\\1&0\end{bmatrix}\times \begin{bmatrix}f_{i-1}\\f_{i-2}\\\end{bmatrix} = \begin{bmatrix}f_i\\f_{i-1}\end{bmatrix} \]

通过如上公式，不停迭代，如果求 \(f_n\)，就是乘以 \(n-2\) 个 \(\begin{bmatrix}1&1\\1&0\end{bmatrix}\)。

因此可以得到如下公式：

\[\begin{bmatrix}1&1\\1&0\end{bmatrix}^{n-2}\times \begin{bmatrix}f_{i-1}\\f_{i-2}\\\end{bmatrix} = \begin{bmatrix}f_n\\f_{n-1}\end{bmatrix} \]

咦，这不就是快速幂吗？只不过把乘法改成了矩阵乘法而已嘛。

我们可以定义一个结构体，用重载操作符代替矩阵乘法，于是代码如下：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll MOD=1e9+7;
const ll N=3;
ll n,k;
struct node{ 
    ll x[N][N];
    node(){
        memset(x,0,sizeof(x));
    }
    inline void bases(){
        memset(x,0,sizeof(x));
        x[1][1]=x[1][2]=x[2][1]=1;
    }
    inline void anses(){
        memset(x,0,sizeof(x));
        x[1][1]=x[2][1]=1;      
    }
    node operator *(const node &b)const{//重载操作符。
        node anss;
        for(ll k=1;k<=2;k++)//矩阵乘法模板。
            for(ll i=1;i<=2;i++)
                for(ll j=1;j<=2;j++)
                    anss.x[i][j]=(anss.x[i][j]+x[i][k]*b.x[k][j])%MOD;
        return anss;
    }
}base,ans;
inline void init(){
    base.bases();
    ans.anses();
}
inline void quick_pow(ll n){//矩阵快速幂。
    while(n){
        if(n&1) ans=base*ans;
        base=base*base;
        n>>=1;
    }
}
int main(){
	cin>>n;
	if(n<=2){
		cout<<1<<endl;
		exit(0);
	}
    init();
    quick_pow(n-2);
    cout<<ans.x[1][1]<<endl;
	return 0;
}

经典例题

• 【模板】矩阵快速幂

就是一道矩阵快速幂模板啦。

\(base\) 矩阵就是读入的 \(a\) 矩阵，\(ans\) 矩阵在图中已经给出，为 \(\begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}\)，也就是第 \(i\) 行第 \(i\) 列为 \(1\)，其余都是
\(0\)。

\(code\)

#define ll long long
const ll MOD=1e9+7;
const ll N=1e2+10;
ll n,k;
struct node{ 
  ll x[N][N];
  node(){memset(x,0,sizeof(x));} 
  inline void init(){//定义 init() 函数初始化，只有该结构体才允许该初始化。
      for(ll i=1;i<=n;i++) x[i][i]=1;
  }
//说句闲话：还可以把矩阵乘法也搬进来。
}a;
node operator * (const node &a,const node &b){//矩阵乘法。
    node ans;
    for(ll k=1;k<=n;k++)
        for(ll i=1;i<=n;i++)
            for(ll j=1;j<=n;j++){
                ans.x[i][j]=(ans.x[i][j]+a.x[i][k]*b.x[k][j]%MOD)%MOD;
            }
    return ans;
}
int main(){
	cin>>n>>k;
    for(ll i=1;i<=n;i++)
        for(ll j=1;j<=n;j++) cin>>a.x[i][j];
    node ans;
    ans.init();//初始化。
    while(k){//矩阵快速幂。
        if(k&1) ans=ans*a;
        a=a*a;
        k>>=1;
    }
    for(ll i=1;i<=n;i++){
        for(ll j=1;j<=n;j++)
            cout<<ans.x[i][j]<<" ";
        cout<<endl;
    }
	return 0;
}

• 广义斐波那契数列

题目给出 \(a_1\)，\(a_2\)，以及递推公式 \(a_i = p\times a_{i-1} + q\times a_{i-2},i\ge 3\)。

让我们求 \(a_n\bmod m\) 的值。

我们平时的斐波那契数列里，\(p=q=1\)。

这里，\(p\)、\(q\) 是个不定值。

我们试着规划矩阵：

\[\begin{bmatrix}?&?\\?&?\end{bmatrix}\times \begin{bmatrix}f_{i-1}\\f_{i-2}\\\end{bmatrix} = \begin{bmatrix}f_i\\f_{i-1}\\\end{bmatrix} \]

明显，矩阵为 \(\begin{bmatrix}p&q\\1&0\\\end{bmatrix}\)。

验证一下，\(p\times f_{i-1} + q\times f_{i-2} = f_i\)，\(1\times f_{i-2} + 0\times f_{i-2} = f_{i-1}\)，矩阵规划正确。

\(\texttt{PS}\)：有些题解是 \(1\times 2\) 的矩阵乘以 \(2\times 2\) 的矩阵，而我是 \(2\times 2\) 的矩阵乘以 \(2\times 1\) 的矩阵，因此规划结果可能不一样，不过没关系。

最后注意取余，初始化不要弄错位置。

代码就不展示了，反正和上面的几乎一样，改一下就行啦。