從快速冪運算到矩陣快速冪
快速冪運算
HDU2035 求
http://acm.hdu.edu.cn/showproblem.php?pid=2035
題目是很簡單的，因為b也不大所以時間複雜度為n的演算法也能ac
#include <iostream> using namespace std; int a, b; int pow(int a,int b,int mod){ int ans=1; for(int i=0;i<b;i++) ans=ans*a%mod; return ans; } int main(){ while(cin>>a>>b){ if(!a&&!b)return 0; cout<<pow(a,b,1000)<<endl; } return 0; }
但是如果b大一點的話，可能就會TLE了。
時間複雜度log(n)的演算法：快速冪運算
#include <iostream> using namespace std; int a, b; int qkpow(int a, int b, int mod){ int ans = 1; while (b){ if(b&1) ans =ans*a%mod; //如果b是奇數，手動乘一次 a=a*a%mod; //a=a的平方 b>>=1; //此處等效於b/=2 }//也就是說 迴圈的出口就是當b=0 上一次b=1的時候 return ans; }//這樣做的時間複雜度是log(n) //本質是把指數拆分為二進位制數之和 int main(){ while(cin>>a>>b){ if(!a&&!b)return 0; cout<<qkpow(a,b,1000)<<endl; } return 0; }
大概這也是求逆元的基礎吧。
矩陣快速冪
此處建議自行了解一點線性代數的相關知識：矩陣乘法。
簡單來說，矩陣快速冪就是在快速冪的基礎上，其乘法用矩陣乘法來代替。
因為是一個矩陣乘自己，根據矩陣乘法的定義，該矩陣只能是方陣。
矩陣提供了一個手段實現用自身乘自身來代替遞推，而且是快速遞推，時間複雜度是log(n)。
理論上說，遞推均可用矩陣快速冪。
#include<bits/stdc++.h> using namespace std; typedef long long ll; const int N=1e2+5; int a[N][N],b[N][N],c[N][N]; void input(int p[][N],int x){ for(int i=0;i<x;i++) for(int j=0;j<x;j++) scanf("%d",&p[i][j]); } void output(int p[][N],int x){ for(int i=0;i<x;i++){ for(int j=0;j<x;j++) printf("%d ",p[i][j]); printf("\n"); } } void multi(int p1[][N],int p2[][N],int x){ memset(p,0,sizeof(p)); for(int i=0;i<x;i++) for(int j=0;j<x;j++) for(int k=0;k<x;k++) p[i][j]+=p1[i][k]*p2[k][j]; for(int i=0; i<n; i++) for(int j=0; j<n; j++) p1[i][j]=p[i][j]; } int res[N][N]; void Pow(int p[][N],int n,int x) { memset(res,0,sizeof(res)); for(int i=0; i<x; i++) res[i][i]=1; while(n) { if(n&1) multi(res,p,N);//res=res*a;複製直接在multi裡面實現了； multi(p,p,N);//a=a*a n>>=1; } } int main(){ }
斐波那契數列
http://poj.org/problem?id=3070
題目：斐波那契數列f(n),給一個n，求f(n)%10000,n<=1e9;
其實我感覺我這種寫法比那種結構體模板的乾淨啊。
斐波那契數列的遞推式：
斐波那契數列的矩陣遞推式：
注意遞推式的預設，此處感謝青竹大佬
3070也是一遍過了，程式碼如下：
#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N=2;
const int mod=10000;
void multi(ll p1[][N],ll p2[][N],int x) {
ll p[N][N];
memset(p,0,sizeof(p));
for(int i=0; i<x; i++)
for(int j=0; j<x; j++)
for(int k=0; k<x; k++)
p[i][j]+=p1[i][k]*p2[k][j]%mod;
for(int i=0; i<x; i++)
for(int j=0; j<x; j++)
p1[i][j]=p[i][j];
}
ll res[N][N];
void Pow(ll p[][N],ll n) {
memset(res,0,sizeof(res));
for(int i=0; i<N; i++) res[i][i]=1;
while(n) {
if(n&1) multi(res,p,N);
multi(p,p,N);
n>>=1;
}
}
int main() {
ll n;
while(cin>>n) {
if(n==-1) break;
ll a[N][N]= {1,1,1,0};
Pow(a,n-1);
cout<<res[0][0]%mod<<endl;
}
return 0;
}
總而言之，利用矩陣快速冪實現快速遞推的關鍵就是找到遞推式，構建初始矩陣
常用矩陣遞推式
搬運自https://blog.csdn.net/zhangxiaoduoduo/article/details/81807338
其他例題
http://poj.org/problem?id=3233
http://acm.hdu.edu.cn/showproblem.php?pid=2276
數三角形
題目
Problem Description
小朋友們都知道將一個三角形三條邊的中點連線會劃分成四個小的三角形，如果一直這樣劃分下去將會出現無數個三角形，如下所示：
小朋友發現裡面有兩種三角形，小朋友將 △ 命名為“帥氣”的三角形 ，將 ▽ 命名為“漂亮”的三角形，小朋友想知道經過 N 次劃分之後有多少個“帥氣”的三角形，初始的時候只有一個“帥氣”的三角形，比如上面示例中劃分一次之後就有3個“帥氣”的三角形了，劃分兩次就有10個“帥氣”的三角形了。
Input
多組輸入，一行包含一個整數 N（0≤N≤10e18） ,表示劃分的次數。
Output
輸出一個整數表示劃分 N次後 “帥氣”的三角形的數量，結果對 1000000007（10e9 +7）求餘數。
Sample Input
0
9
999
Sample Output
1
131328
265758117
solution
通過觀察我們可以發現：
每個上三角完成一次劃分後，生成3個上三角和1個下三角；
每個下三角完成一次劃分後，生成1個下三角和3個上三角。
即，設第i次劃分後上三角個數為$x_i$，下三角個數為$y_i$，有
$
x_{i}=3x_{i-1}+y_{i-1}
$
$
y_{i}=x_{i-1}+3y_{i-1}
$
由此我們可以完成遞推。
#include<iostream>
#include<cstdio>
using namespace std;
const int N=1e3+7;
const int mod=1e9+7;
int main() {
long long x[N],y[N];
x[0]=1,y[0]=0;
for(int i=1;i<=1000;++i){
x[i]=(3*x[i-1]+y[i-1])%mod;
y[i]=(x[i-1]+3*y[i-1])%mod;
}
int n;
while(cin>>n) cout<<x[n]<<endl;
return 0;
}
但是，這樣的做法的時間複雜度是O(n)，這也就意味著在（0≤N≤10e18）的前提下一定會超時。
那麼O(logn)的做法是什麼呢？
矩陣快速冪。
要使用矩陣快速冪，我們就要推匯出來轉移矩陣。
解得 轉移矩陣為
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod = 1e9 + 7;
struct node { //矩陣類
ll matrix[2][2];
node() { memset(matrix, 0, sizeof(matrix)); } //初始化
void one() { //單位矩陣E
matrix[0][0] = 1;
matrix[1][1] = 1;
}
node operator*(node other) { //對*過載 定義矩陣乘法
node ans; //記錄乘積
for (int i = 0; i < 2; i++)
for (int j = 0; j < 2; j++)
for (int k = 0; k < 2; k++) {
ans.matrix[i][j] += matrix[i][k] * other.matrix[k][j];
ans.matrix[i][j] %= mod;
}
return ans;
}
};
node power(node a, ll b) { //快速冪 矩陣a的b次方
node res;
res.one(); //單位矩陣
while (b) {
if (b & 1) res = res * a;
a = a * a;
b >>= 1;
}
return res;
}
int main() {
node temp;
temp.matrix[0][0] = 3;
temp.matrix[0][1] = 1;
temp.matrix[1][0] = 1;
temp.matrix[1][1] = 3;
ll n, flag = 0;
while (cin >> n) {
if (n == 0)
cout << 1 << endl;
else
cout << power(temp, n).matrix[0][0] << endl;
}
return 0;
}
這個程式碼是隊友的，我稍微改了一點。
然後格式我也改了一下，這個格式是正常的格式（PE警告
我在比賽的時候卡住了，沒有想出來轉移矩陣怎麼解……
我沒想出來轉移矩陣怎麼解，又要求O(logn)的話怎麼辦呢？
快速冪。
因為
$
x_{i}=3x_{i-1}+y_{i-1}
$
$
y_{i}=x_{i-1}+3y_{i-1}
$
所以
$
x_{i}+y_{i}=4x_{i-1}+4y_{i-1}=4(x_{i-1}+y_{i-1})
$
$
x_{i}-y_{i}=2x_{i-1}-2y_{i-1}=2(x_{i-1}-y_{i-1})
$
所以
最後注意模1e9+7不能用除法，/2等效於*它的逆元。這裡涉及到費馬小定理，不展開，有興趣自己去了解。
#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
ll mod=1e9+7;
ll inv=500000004;
ll qkpow(ll a,ll b) {
ll ans=1;
while(b) {
if(b&1) ans=ans*a%mod;
a=a*a%mod;
b>>=1;
}
return ans;
}
int main() {
ll n;
while(cin>>n) cout<<(qkpow(4,n)+qkpow(2,n))*inv%mod<<endl;
return 0;
}
模板
qkpow & inv
ll qkpow(ll a,ll b) {
ll ans=1;
while(b) {
if(b&1) ans=ans*a%mod;
a=a*a%mod;
b>>=1;
}
return ans;
}
ll getInv(ll a){return qkpow(a,mod-2);} //求一個數的逆元
矩陣快速冪
楊文冠的2*2
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod = 1e9 + 7;
struct node { //矩陣類
ll matrix[2][2];
node() { memset(matrix, 0, sizeof(matrix)); } //初始化
void one() { //單位矩陣E
matrix[0][0] = 1;
matrix[1][1] = 1;
}
node operator*(node other) { //對*過載 定義矩陣乘法
node ans; //記錄乘積
for (int i = 0; i < 2; i++)
for (int j = 0; j < 2; j++)
for (int k = 0; k < 2; k++) {
ans.matrix[i][j] += matrix[i][k] * other.matrix[k][j];
ans.matrix[i][j] %= mod;
}
return ans;
}
};
node power(node a, ll b) { //快速冪 矩陣a的b次方
node res;
res.one(); //單位矩陣
while (b) {
if (b & 1) res = res * a;
a = a * a;
b >>= 1;
}
return res;
}
int main() {
node temp;
temp.matrix[0][0] = 3;
temp.matrix[0][1] = 1;
temp.matrix[1][0] = 1;
temp.matrix[1][1] = 3;
ll n, flag = 0;
while (cin >> n) {
if (n == 0)
cout << 1 << endl;
else
cout << power(temp, n).matrix[0][0] << endl;
}
return 0;
}
我自己寫的
#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N=2;
const int mod=10000;
void multi(ll p1[][N],ll p2[][N],int x) {
ll p[N][N];
memset(p,0,sizeof(p));
for(int i=0; i<x; i++)
for(int j=0; j<x; j++)
for(int k=0; k<x; k++)
p[i][j]+=p1[i][k]*p2[k][j]%mod;
for(int i=0; i<x; i++)
for(int j=0; j<x; j++)
p1[i][j]=p[i][j];
}
ll res[N][N];
void Pow(ll p[][N],ll n) {
memset(res,0,sizeof(res));
for(int i=0; i<N; i++) res[i][i]=1;
while(n) {
if(n&1) multi(res,p,N);
multi(p,p,N);
n>>=1;
}
}
int main() {
ll n;
while(cin>>n) {
if(n==-1) break;
ll a[N][N]= {1,1,1,0};
Pow(a,n-1);
cout<<res[0][0]%mod<<endl;
}
return 0;
}