矩阵快速幂
矩阵乘法是一种矩阵运算,满足交换律,可以写成幂的形式,所以我们可以使用矩阵快速幂来解决一些问题
先来看普通的快速幂:
快速幂
ll qpow(int a,int p){
if(p == 0)return 1;
ll temp = (qpow(a,p >> 1)) % M;
if(p % 2 == 0){
return (temp * temp) % M;
}
else{
return (temp * temp * a) % M;
}
}
把幂拆成 n * 2 或者n * 2 + 1的形式,避免重复运算,提高效率。因为结果通常很大,记得不要弄错MOD
矩阵乘法
设A为 的矩阵,B为 的矩阵,那么称 的矩阵C为矩阵A与B的乘积,记作 ,其中矩阵C中的第 行第 列元素可以表示为:
如下所示:
(摘自百度百科)
所以我们说:矩阵能相乘是有条件的:A矩阵的列数必须等于B的行数,得到的矩阵行数等于A的行数,列数等于B的列数
以下是矩阵乘法的代码
struct Mtrix{
ll a[19][19];
ll n,m;
};
Mtrix X(Mtrix a,Mtrix b){
Mtrix ans;
for(ll i = 1;i <= a.n;i++)
for(ll j = 1;j <= b.m;j++)
ans.a[i][j] = 0;
ans.n = a.n;
ans.m = b.m;
for(ll i = 1;i <= ans.n;i++){
for(ll j = 1;j <= ans.m;j++){
for(ll k = 1;k <= a.m;k++){
ans.a[i][j] += (a.a[i][k] * b.a[k][j]) % M;
ans.a[i][j] %= M;
}
}
}
return ans;
}
矩阵快速幂
模仿一下一般快速幂,我们就可以得到矩阵快速幂的代码
Mtrix M_qpow(Mtrix a,ll p){
if(p == 1)return a;
Mtrix temp = M_qpow(a,p >> 1);
if(p % 2 == 0)return X(temp,temp);
else return X(a,X(temp,temp));
}
我们不知道当p == 0时的矩阵值是多少,就干脆返回p == 1时为原矩阵就好
模板题:P3390 【模板】矩阵快速幂
利用矩阵快速幂优化递推问题
先看这题:
P1962 斐波那契数列
题目背景
大家都知道,斐波那契数列是满足如下性质的一个数列:
• f(1) = 1
• f(2) = 1
• f(n) = f(n-1) + f(n-2) (n ≥ 2 且 n 为整数)
题目描述
请你求出 f(n) mod 1000000007 的值。
输入输出格式
输入格式:
·第 1 行:一个整数 n
输出格式:
第 1 行: f(n) mod 1000000007 的值
说明
对于 60% 的数据: n ≤ 92
对于 100% 的数据: n在long long(INT64)范围内。
第一眼:水题!!
然后
瞟了一眼数据范围n在long long(INT64)范围内,吓个半死。
普通记忆化所需要的大量空间是肯定支持不了了,这时候就要用的矩阵快速幂了
首先依据斐波那契的定义式:
那么对于一个斐波那契矩阵:
第几项 | 值 |
---|---|
n | f(n) |
n + 1 | f(n + 1) |
我们可以这样写:
第几项 | 值 |
---|---|
n | f(n) = f(n) |
n + 1 | f(n + 1) = f(n) + f(n - 1) |
然后是不是可以由这个得到:
第几项 | 值 |
---|---|
n - 1 | f(n - 1) |
n | f(n) |
具体:
f(n)=上面那个矩阵第一行×0+第二行×1
f(n+1)=上面矩阵第一行×1+第二行×1
怎么样,是不是特别熟悉?没错,矩阵运算!
它长这样:
0 | 1 |
---|---|
1 | 1 |
这就是斐波那契的递推矩阵
再观察我们可以发现,f(n)其实就是f(1);f(2)【纵向排列】乘以递推矩阵的(n - 1)次方,那么问题就解决了
代码:
#include<iostream>
#include<algorithm>
#include<cstdio>
#include<queue>
#define ll long long
using namespace std;
ll RD(){
ll flag = 1,out = 0;char c = getchar();
while(c < '0' || c > '9'){if(c == '-')flag = -1;c = getchar();}
while(c >= '0' && c <= '9'){out = out * 10 + c - '0';c = getchar();}
return flag * out;
}
const ll M = 1000000007;
ll n;
struct Mtrix{
ll a[19][19];
ll n,m;
};
Mtrix X(Mtrix a,Mtrix b){
Mtrix ans;
for(ll i = 1;i <= a.n;i++)
for(ll j = 1;j <= b.m;j++)
ans.a[i][j] = 0;
ans.n = a.n;
ans.m = b.m;
for(ll i = 1;i <= ans.n;i++){
for(ll j = 1;j <= ans.m;j++){
for(ll k = 1;k <= a.m;k++){
ans.a[i][j] += (a.a[i][k] * b.a[k][j]) % M;
ans.a[i][j] %= M;
}
}
}
return ans;
}
Mtrix M_qpow(Mtrix a,ll p){
if(p == 1)return a;
Mtrix temp = M_qpow(a,p >> 1);
if(p % 2 == 0)return X(temp,temp);
else return X(a,X(temp,temp));
}
Mtrix a,ori;
int main(){
n = RD();
if(n == 1){
cout<<1<<endl;
return 0;
}
a.n = 2,a.m = 1;
a.a[1][1] = 1,a.a[2][1] = 1;
ori.n = 2,ori.m = 2;
ori.a[1][1] = 0,ori.a[1][2] = 1,ori.a[2][1] = 1,ori.a[2][2] = 1;
Mtrix last = M_qpow(ori,n - 1);
Mtrix ans = X(last,a);
cout<<ans.a[1][1] % M<<endl;
return 0;
}
总结
利用矩阵快速幂来求解递推是递推的一大法宝,要熟练运用