矩阵快速幂与矩阵递推加速
1 矩阵
1.1 介绍
由m×n个数a_{ij}排成的m行n列的长方形数表称为m行n列的矩阵，简称m×n矩阵，其用中括号括起。下面这三个都是矩阵：
可以看出，以上矩阵分别是2x3矩阵、3x3矩阵和4x4矩阵。行数与列数都等于n的矩阵称为n阶矩阵或n阶方阵，因此后两个矩阵也可称为3阶矩阵和4阶方阵。最后一个矩阵还可称为单位矩阵，因为其正对角线上的数字为1，其余数字为0。
1.2 矩阵的基本运算
矩阵的基本运算包括加法，减法，数乘，乘法等，下面依次介绍。
矩阵加法就是矩阵对应位置上的数字相加，这要求相加的两个矩阵必须有相同的大小（同样的m和n），下面是加法的一个例子：
矩阵加法满足以下运算律：
矩阵减法与矩阵加法类似。
矩阵数乘就是用一个常数乘以矩阵所有位置上的数，下面是一个例子：
矩阵数乘满足以下运算律：
1.3 矩阵乘法
矩阵乘法只有在两个矩阵A=(a_{ij})中左矩阵的列数和右矩阵B=(b_{ij})的行数相等时才能进行。如果左矩阵A是m×n矩阵，右矩阵B是n×p矩阵，它们的乘积矩阵是m×p矩阵C=(c_{ij})，矩阵乘法的式子记为C=AB，矩阵内数之间的关系如下:
下面是一个例子：
矩阵的乘法满足以下运算律：
注意，矩阵乘法不满足交换律。
2 矩阵快速幂
2.1 快速幂原理
一般来说，求一个数的n次方幂，我们会通过连乘n-1次来求得。但只需做个简单的改进就能减少连乘的次数。
幂运算有个法则：a^b\times a^c=a^{(b+c)}，也就是说幂运算是可以拆开来算的。假设是计算2^8，对半拆也就是两个2^4相乘，而计算2^4也可以用相同的方法，这样最终只需要计算最简单的2^2=2\times 2，这样的思路就是快速幂，通过利用现有的结果进行运算，这样对于n次幂，只需要计算\log n次。下面是快速幂的代码。
long long fast_power(int a,int b,int mod) //mod用于取余
{
long long ans=1,base=a;
while(b)
{
if(b&1) //判定是否是奇数
ans=ans*base%mod;
base=base*base%mod;
b>>=1;
}
return ans%mod; //非多余步，如果没有b=0这种数据可以去掉取余
}在快速幂的代码的基础上将int类型换成矩阵的类型，就是矩阵快速幂了。要注意的是，只有n阶方阵才会有次幂。
2.2 矩阵与矩阵乘法的实现
现在来实现矩阵与矩阵乘法。矩阵可以直接用二维数组表示，矩阵乘法通过重载X运算符实现，请看代码。
int mod=int(1e9+7);
struct matrix
{
long long data[101][101];
int n,m;
matrix(int _n,int _m)
{
memset(data,0,sizeof(data));
n=_n;
m=_m;
}
matrix operator*(const matrix &b) const //矩阵乘法
{
matrix t(n,b.m);
if (m==b.n)
{
for (int i=1;i<=n;i++) //基本枚举，下同
for (int j=1;j<=b.m;j++)
for (int k=1;k<=m;k++)
t.data[i][j]=(t.data[i][j]+data[i][k]*b.data[k][j])%mod; //取模运算不能丢
}
return t;
}
matrix& operator*=(const matrix b) //乘赋值
{
return *this=*this*b;
}
void show() //显示
{
for (int i=1;i<=n;i++)
{
for (int j=1;j<=m;j++)
cout<<data[i][j]<<' ';
cout<<'\n';
}
}
};2.3 最终整体
写完矩阵和矩阵乘法之后，只需要小小修改快速幂的函数就能实现矩阵快速幂了，只不过这里还得讨论一下ans=1怎么改。1在计算中可以使非0的数与之相乘仍是该数，同理矩阵中也有一个矩阵有类似的作用，它就是上面提到的单位矩阵。所以ans=1的部分用单位矩阵替代1即可，请看代码。
matrix& fast_power(int b)
{
if (m!=n) return *this;
matrix ans(n,m);
ans.id();
while(b)
{
if (b&1)
ans*=*this;
(*this)*=*this;
b>>=1;
}
*this=ans;
return *this;
}
void id() //单位化
{
if (m!=n) return;
for (int i=1;i<=n;i++)
data[i][i]=1;
}矩阵快速幂的完整代码如下。
int mod=int(1e9+7);
struct matrix
{
long long data[101][101];
int n,m;
matrix(int _n,int _m)
{
memset(data,0,sizeof(data));
n=_n;
m=_m;
}
matrix operator*(const matrix &b) const //矩阵乘法
{
matrix t(n,b.m);
if (m==b.n)
{
for (int i=1;i<=n;i++) //基本枚举，下同
for (int j=1;j<=b.m;j++)
for (int k=1;k<=m;k++)
t.data[i][j]=(t.data[i][j]+data[i][k]*b.data[k][j])%mod; //取模运算不能丢
}
return t;
}
matrix& operator*=(const matrix b) //乘赋值
{
return *this=*this*b;
}
matrix& fast_power(int b)
{
if (m!=n) return *this;
matrix ans(n,m);
ans.id();
while(b)
{
if (b&1)
ans*=*this;
(*this)*=*this;
b>>=1;
}
*this=ans;
return *this;
}
void id() //单位化
{
if (m!=n) return;
for (int i=1;i<=n;i++)
data[i][i]=1;
}
void show() //显示
{
for (int i=1;i<=n;i++)
{
for (int j=1;j<=m;j++)
cout<<data[i][j]<<' ';
cout<<'\n';
}
}
};3 矩阵加速
3.1 线性方程组
线性方程组是指各个方程中的未知量均为一次的方程组。它长下面的样子：
这是一个三元线性方程组，其中x_i表未知量，a_{ij}称为系数，b_j称为常数项。将系数用矩阵来表示就是下面的系数矩阵：
由矩阵乘法可以得到关系：
这就是矩阵加速的基本原理。
3.2 矩阵加速
求斐波那契数列的基本方法是递归或递推，写起来很简单，但是这些写法是不能快速计算出斐波那契数列中下标比较大的项的。而现在，我们可以通过构造系数矩阵并使用矩阵快速幂来快速计算下标比较大的项。这一方法称为被矩阵加速。
对于斐波那契数列，将f[n-1]和f[n-2]作为未知量，递推式可以转化为为线性方程组：
那么其矩阵乘法的关系式就是：
你可能会问：为什么是两个方程而不是一个，一个方程看起来也是够的？用两个方程这是因为递推式f[n]=f[n-1]+f[n-2]决定的矩阵乘法关系式是：
系数矩阵不是方阵，没法进行幂运算，并且如果将系数矩阵看作为一个特殊的对应法则g。那么是这样：
你可以从左推导到右，但是你不能从右推导到左，或者说你不能求出对应法则g^{-1}，也就是逆矩阵，因为有逆矩阵的一个充分条件是原矩阵是方阵。这就是在说，有多少个未知量就需要构造多少个方程。
回到正题，得到矩阵乘法关系式后，只要有初始条件，就可以通过不断的将系数矩阵乘以由未知量构成的矩阵来得到结果，又因为矩阵乘法的结合律，只需要先求出系数矩阵的次幂然后再乘以由初始条件构成的矩阵就可以得到结果。
3.3 矩阵构造法
矩阵加速的核心就是系数矩阵的构造，所以下面将举一些系数矩阵构造的例子来讲解。
递推式f[n]=f[n-1]+f[n-2]，初始条件f[1]=f[2]=1。求第n项
斐波那契数列，这个已经在上面提到过。不再过多赘述。
递推式f[n]=f[n-1]+f[n-2]+1，初始条件f[1]=f[2]=1。求第n项
注意递推式中的1也要算作为未知量，但因为它是常数，所以需要保证其值不变。矩阵乘法关系式为：
递推式f[n]=f[n-1]+f[n-2]+n+1，初始条件f[1]=f[2]=1。求第n项
跟第二个同理，不过要注意对于n，递推一遍后应该是变成n+1。矩阵乘法关系式为：
递推式f[n]=f[n-1]+f[n-2]，初始条件f[1]=f[2]=1。求前n项和s[n]
由于s[n]是我们要求的答案，所以s[n]也要当成未知量。可以获得一个额外的递推式s[n-1]=s[n-2]+f[n-1]。那么矩阵乘法关系式为：
下面再来实操一道题目。
首先根据递推式子得出矩阵关系式，这里初始条件有三个，那应该使用三个未知量，即：
代码如下：
#include<iostream>
#include<cstring>
using namespace std;
int mod=int(1e9+7);
struct matrix
{
long long data[101][101];
int n,m;
matrix(int _n,int _m)
{
memset(data,0,sizeof(data));
n=_n;
m=_m;
}
matrix operator*(const matrix &b) const //矩阵乘法
{
matrix t(n,b.m);
if (m==b.n)
{
for (int i=1;i<=n;i++) //基本枚举，下同
for (int j=1;j<=b.m;j++)
for (int k=1;k<=m;k++)
t.data[i][j]=(t.data[i][j]+data[i][k]*b.data[k][j])%mod; //取模运算不能丢
}
return t;
}
matrix& operator*=(const matrix b) //乘赋值
{
return *this=*this*b;
}
matrix& fast_power(int b)
{
if (m!=n) return *this;
matrix ans(n,m);
ans.id();
while(b)
{
if (b&1)
ans*=*this;
(*this)*=*this;
b>>=1;
}
*this=ans;
return *this;
}
void id() //单位化
{
if (m!=n) return;
for (int i=1;i<=n;i++)
data[i][i]=1;
}
void show() //显示
{
for (int i=1;i<=n;i++)
{
for (int j=1;j<=m;j++)
cout<<data[i][j]<<' ';
cout<<'\n';
}
}
};
int main()
{
int t,n;
matrix X(3,1);
X.data[1][1]=X.data[2][1]=X.data[3][1]=1; //设定初始条件
matrix A(3,3);
A.data[1][1]=A.data[1][3]=A.data[2][1]=A.data[3][2]=1; //设定系数矩阵
cin>>t;
while (t--)
{
matrix temp=A;
cin>>n;
temp.fast_power(n-1); //a1到an需要n-1次幂
temp*=X;
cout<<temp.data[3][1]<<'\n';
}
return 0;
}更新日志：
2024.8.6 - 文章全面更新