前言: 矩阵快速幂就是缝合怪
快速幂
这篇文章写的真心好啊,%%%.其实在理解其它的算法时,用纸笔推导每一步在干什么,也能达到差不多的理解程度,所以一定要耐心.
首先给一个幂运算:
a
p
a^p
ap
在代码实现中,实现p次方的复杂度为
O
(
p
)
O(p)
O(p).我们考虑缩减指数,若p是偶数,则
a
p
=
(
a
∗
a
)
p
2
a^p=(a*a)^{frac p 2}
ap=(a∗a)2p.只进行了一次乘法运算,却将指数缩减到了原来的二分之一(这也是为什么快速幂的复杂度是
l
o
g
n
log n
log n).
有了想法,我们再看能不能实现.由唯一分解定理得:整数 m = ∏ a i p i = ∏ a i 2 k + n = ∏ ( a i 2 k ) 1 ( a i n ) 1 m=prod a_i^{p_i}=prod a_i^{2k+n}=prod (a_i^{2k})^1(a_i^{n})^1 m=∏aipi=∏ai2k+n=∏(ai2k)1(ain)1,也就是说,m可分解为若干个幂次为1的数的乘积(好像是废话),那么,在缩减指数的时候,若 p m o d 2 = = 1 p mod 2==1 p mod 2==1,ans*=a即可.
初步代码：
ll fpm(ll a,ll p,ll mod)
{
ll ans=1;
while(p>1)
{
if(p%2==1)
{
ans=(ans*a)%mod;
}
p/=2;
a=(a*a)%mod;
}
ans=(ans*a)%mod;
return ans;
}
优化代码:
ll fpm(ll a,ll p,ll mod)
{
ll ans=1;
while(p>=1)
{
if(p&1) (ans*=a)%=mod;
p>>=1;
(a*=a)%=mod;
}
return ans;
}
分析while循环,每次做的事为p>>=1和(a*=a)%=mod,结束条件为p>=1,可将其改写为for循环节省行数
ll fpm(ll a,ll p,ll mod)
{
ll ans=1;
for(;p;p>>=1,(a*=a)%=mod)
if(p&1) (ans*=a)%=mod;
return ans;
}
注意(a*=a)%=mod)是在if(p&1) (ans*=a)%=mod;后执行的.
矩阵快速幂
矩阵乘法已经在线性代数中提到了.顾名思义,矩阵快速幂用于解决诸如 A p A^p Ap的问题,为什么会用到一个矩阵的幂呢?我们先来回顾斐波那契数列 f n = f n − 1 + f n − 2 f_n=f_{n-1}+f_{n-2} fn=fn−1+fn−2,这个简单的递推式在数据规模为1e7时能轻松通过,但到了1e8便很危险了.这时候矩阵便派上用场了.
斐波那契数列的第n项有通项公式 F n = ( 5 + 1 2 ) n − ( 5 − 1 2 ) n 5 F_n= frac {(frac {sqrt{5}+1} {2})^n-(frac{sqrt{5}-1}{2})^n } {sqrt 5} Fn=5(25+1)n−(25−1)n,但对精度要求极高,不适用.但我们能不能构造一个类似的通项公式,只不过求出的是一个矩阵呢?
可能你会疑惑,为什么要求一个矩阵呢?请看
[
f
n
f
n
−
1
]
egin{bmatrix}f_n\f_{n-1}end{bmatrix}
[fnfn−1],这是不是一个矩阵?求出的是不是斐波那契数列?这样,我们的矩阵、快速幂与递推便有机统一了.现在到了在做题时最有思维难度的部分:构造矩阵.
思考:
[
f
[
i
]
f
[
i
−
1
]
]
=
[
a
b
c
d
]
×
[
f
[
i
−
1
]
f
[
i
−
2
]
]
egin{bmatrix}f[i]\f[i-1]end{bmatrix}=egin{bmatrix}a&b\c&dend{bmatrix} imesegin{bmatrix}f[i-1]\f[i-2]end{bmatrix}
[f[i]f[i−1]]=[acbd]×[f[i−1]f[i−2]]
第一个矩阵该填什么?首先它肯定是长这个样子的,也就是要找出a,b,c,d满足
a
×
f
[
i
−
1
]
+
b
×
f
[
i
−
2
]
=
f
[
i
]
a imes f[i-1]+b imes f[i-2]=f[i]
a×f[i−1]+b×f[i−2]=f[i]
c × f [ i − 1 ] + d × f [ i − 2 ] = f [ i − 1 ] c imes f[i-1]+d imes f[i-2]=f[i-1] c×f[i−1]+d×f[i−2]=f[i−1]
c=1,d=2是显然的,a和b等于多少呢?但观察形式,它不就是最基本的斐波那契递推式吗?
于是就有了:
f
[
i
]
=
1
×
f
[
i
−
1
]
+
1
×
f
[
i
−
2
]
(1)
ag 1f[i]=1 imes f[i-1]+1 imes f[i-2]
f[i]=1×f[i−1]+1×f[i−2](1)
f
[
i
−
1
]
=
1
×
f
[
i
−
1
]
+
0
×
f
[
i
−
2
]
(2)
ag 2f[i-1]=1 imes f[i-1]+0 imes f[i-2]
f[i−1]=1×f[i−1]+0×f[i−2](2)
由(1)(2)发现式(3)
[
f
[
i
]
f
[
i
−
1
]
]
=
[
1
1
1
0
]
×
[
f
[
i
−
1
]
f
[
i
−
2
]
]
(3)
ag 3 egin{bmatrix}f[i]\f[i-1]end{bmatrix}=egin{bmatrix}1&1\1&0end{bmatrix} imes egin{bmatrix}f[i-1]\f[i-2]end{bmatrix}
[f[i]f[i−1]]=[1110]×[f[i−1]f[i−2]](3)
ohhhhhhh!构造出来了!剩下的不就简单多了吗,根据式(3)归纳出式(4)
[
f
[
n
+
1
]
f
[
n
]
]
=
[
1
1
1
0
]
n
−
1
×
[
f
[
2
]
f
[
1
]
]
(4)
ag 4 egin{bmatrix}f[n+1]\f[n]end{bmatrix}=egin{bmatrix}1&1\1&0end{bmatrix}^{n-1} imes egin{bmatrix}f[2]\f[1]end{bmatrix}
[f[n+1]f[n]]=[1110]n−1×[f[2]f[1]](4)
至此,我们已经完成了"构造一个矩阵的通项公式"的任务了,可以编写代码了.
接下来该考虑的事情是如何将矩阵乘法运算套入快速幂,在快速幂的求解中,即将 a p a^p ap拆分成若干个幂次为1的数的乘积,只用到了结合律,显然其对矩阵乘法的运算也是成立的.(其实在这里也满足交换律,毕竟全都是同一个矩阵A).那么我们直接将快速幂中涉及到乘法的部分换为矩阵乘法即可.
void fpm_mat(matrix &ans,matrix a,ll p)
{
for(;p;p>>=1,a.mul(a))
if(p&1) ans.mul(a);
}
封装了一个矩阵板子,支持乘法操作(我居然自己造了一个板子,太不可思议了)
struct matrix
{
ll g[N][N];
void mul(matrix b)
{
matrix a;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
a.g[i][j]=g[i][j],g[i][j]=0;
for(int i=1;i<=n;i++)
for(int k=1;k<=n;k++)
{
int s=a.g[i][k];
for(int j=1;j<=n;j++)
(g[i][j]+=s*b.g[k][j])%=mod;
}
};
索然无味的AC代码:
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=105,mod=1e9+7;
ll n,k;
struct matrix
{
ll g[N][N];
void mul(matrix b)
{
matrix a;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
a.g[i][j]=g[i][j],g[i][j]=0;
for(int i=1;i<=n;i++)
for(int k=1;k<=n;k++)
{
int s=a.g[i][k];
for(int j=1;j<=n;j++)
(g[i][j]+=s*b.g[k][j])%=mod;
}
};
void fpm_mat(matrix &ans,matrix a,ll p)
{
for(;p;p>>=1,a.mul(a))
if(p&1) ans.mul(a);
}
int main()
{
cin>>n>>k;
matrix a;
for(int i=1;i<=n;i++)
{
for(int j=1;j<=n;j++)
{
cin>>a.g[i][j];
}
}
matrix ans;
//构造单位矩阵
for(int i=1;i<=n;i++)
{
for(int j=1;j<=n;j++)
{
if(i==j) ans.g[i][j]=1;
else ans.g[i][j]=0;
}
}
fpm_mat(ans,a,k);
for(int i=1;i<=n;i++)
{
for(int j=1;j<=n;j++)
{
cout<<ans.g[i][j]<<' ';
}
cout<<endl;
}
return 0;
}
这里涉及到一个叫做单位矩阵(identity matrix)的东西,就是主对角线上元素全为1的矩阵. I × A = A I imes A=A I×A=A.
好了,回到我们的斐波那契数列
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=105,mod=1e9+7;
ll n,k;
struct matrix
{
ll g[N][N];
void mul(matrix b,int n,int m)
{
matrix a;
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)
a.g[i][j]=g[i][j],g[i][j]=0;
for(int i=1;i<=n;i++)
for(int k=1;k<=n;k++)
{
int s=a.g[i][k];
for(int j=1;j<=m;j++)
(g[i][j]+=s*b.g[k][j])%=mod;
}
}
void idt()
{
for(int i=1;i<=2;i++)
for(int j=1;j<=2;j++)
if(i==j) g[i][j]=1;
else g[i][j]=0;
}
};
void fpm_mat(matrix &ans,matrix a,ll p)
{
for(;p;p>>=1,a.mul(a,2,2))
if(p&1) ans.mul(a,2,2);
}
int main()
{
cin>>n;
if(n==1)
{
cout<<1;
return 0;
}
matrix a,idt;
idt.idt();
a.g[1][1]=1,a.g[1][2]=1,a.g[2][1]=1,a.g[2][2]=0;
fpm_mat(idt,a,n);//结果存在idt里
matrix f;
f.g[1][1]=1,f.g[2][1]=1;
idt.mul(f,2,1);
cout<<idt.g[2][1];
return 0;
}
代码虽然AC了,但是因为状态不好,写的极度糟糕.
改进了一下板子,增加了idt构造单位矩阵功能:
struct matrix
{
ll g[N][N];
void mul(matrix b,int n,int m)
{
matrix a;
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)
a.g[i][j]=g[i][j],g[i][j]=0;
for(int i=1;i<=n;i++)
for(int k=1;k<=n;k++)
{
int s=a.g[i][k];
for(int j=1;j<=m;j++)
(g[i][j]+=s*b.g[k][j])%=mod;
}
}
void idt(int n)
{
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
if(i==j) g[i][j]=1;
else g[i][j]=0;
}
};