第？課——基於矩陣快速冪的遞推解法
由於中間的數論部分我自己學的很差，沒有辦法寫出清晰的部落格來，所以這裡跳過了數論部分的部落格，來到矩陣快速冪。
遞推
遞推是一個非常常用的工具。比如經典的斐波那契數列：
\[f(x)= \left\{
\begin{array}{**lr**}
1 &, 0\leq x\leq 1 \\
f(x-1)+f(x-2)&, 2 \leq x \\
\end{array}
\right.
\]
很明顯，\(f(n)\)依賴於\(f(n-1)\)和\(f(n-2)\)，所以我們需要先計算\(f(n-1)f(n-2)\)才能計算\(f(n)\)。
假設我們現在需要求f(1e9+7)，你很快發現了，這個數字非常大。所以我們要求只需要結果對\(MOD\)取模就好了，而\(MOD=1e9+7\)。問題是：我們的迭代演算法是\(O(n)\)的。那如何快速的求解這樣一個遞推問題呢？
矩陣與遞推的聯絡
讓我們站在巨人的肩膀上來看這個遞推問題的第二項以及之後：
\[\begin{pmatrix}f(n)\\? \\\end{pmatrix}=
\begin{pmatrix}1&1\\?&?\\\end{pmatrix}\dot\\
\begin{pmatrix}f(n-1)\\f(n-2)\end{pmatrix}
\]
數學家們把這個問題用矩陣的形式表現了出來。但是矩陣上還有一行是空的。打個比方，我們在求\(f(2)\)的時候，矩陣可以寫成：
\[\begin{pmatrix}f(2)\\? \\\end{pmatrix}=
\begin{pmatrix}1&1\\?&?\\\end{pmatrix}\dot\\
\begin{pmatrix}f(1)\\f(0)\end{pmatrix}
\]
那如果我們需要繼續計算\(f(3)\)呢？
\[\begin{pmatrix}f(3)\\? \\\end{pmatrix}=
\begin{pmatrix}1&1\\?&?\\\end{pmatrix}\dot\\
\begin{pmatrix}f(2)\\f(1)\end{pmatrix}
\]
我們需要知道
\[\begin{pmatrix}f(2)\\f(1)\end{pmatrix}
\]
可是\(f(2)f(1)\)哪裡來呢？不妨把我們的通項改寫一下，在計算\(f(n)\)的同時順帶計算\(f(n-1)\)？
\[\begin{pmatrix}f(n)\\f(n-1) \\\end{pmatrix}=
\begin{pmatrix}1&1\\1&0\\\end{pmatrix}\dot\\
\begin{pmatrix}f(n-1)\\f(n-2)\end{pmatrix}
\]
！！wow，我們得到了一個強大的新遞推式子。這時候有人要不樂意了，這不是複雜化了麼，我們本來只要兩個數加一加，現在還要算一個2$\times$4的矩陣乘法？可是，我們的矩陣是常數。 再稍微改寫一下這個式子：
\[\begin{pmatrix}f(n)\\f(n-1) \\\end{pmatrix}=
\begin{pmatrix}1&1\\1&0\\\end{pmatrix}^{n-1} \\ \dot\\
\begin{pmatrix}f(1)\\f(0)\end{pmatrix}
\]
相信你看到這裡已經頓悟了，因為這個矩陣的高次冪我們可以大做文章！
快速冪與矩陣快速冪
快速冪
快速冪很簡單，這裡直接給出程式碼：
int quickpow(int a,int b,int x) {
while(x>0) {
if(x&1) a=a*b;
b=b*b;
x>>=1;
}
return a;
}
快速冪原理給一個友情連結：
快速冪總結
矩陣乘法
先放一個市面上常見的通用矩乘：
for(int i=1; i<=n; i++) {
for(int j=1; j<=n; j++) {
for(int k=1; k<=n; k++) {
tmp.mat[i][j]+=(a.mat[i][k]%mod*b.mat[k][j]%mod)%mod; \\拖慢速度
tmp.mat[i][j]%=mod;
}
}
}
這個矩陣乘非常符合人的想法，但是有一個點讓它比較慢，我們矩陣快速冪只需要進行冪次乘法，可以做一些最佳化：
for(int i=1; i<=n; i++) {
for(int j=1; j<=n; j++) {
for(int k=1; k<=n; k++) {
tmp.mat[i][k]+=(a.mat[i][j]%mod*b.mat[j][k]%mod)%mod;
tmp.mat[i][k]%=mod;
}
}
}
看出不一樣的地方了嗎？唯一的區別就在下標處。把k放在第一維會大大增加定址的計算量，所以需要把k放在第二維上就能減少非常多的計算量。其實矩陣小的時候也沒什麼區別
矩陣快速冪
其實和普通快速冪沒有多大區別，只需要過載一個乘法運算子就好了。這裡給一份結構體加快速冪的程式碼：
#include <iostream>
#include <string.h>
using namespace std;
const int mod = (int)1e9+7;
struct MyMat{
MyMat(int n_=24,int m_=24):n(n_),m(m_){
memset(mat,0,sizeof(mat));
}
int mat[25][25];
int m,n;
friend MyMat operator*(const MyMat&a,const MyMat&b){
MyMat tmp;
if(a.m!=b.n)throw("Matrix size mismatch");
for(int i=1; i<=a.n; i++) {
for(int j=1; j<=b.m; j++) {
for(int k=1; k<=a.m; k++) {
tmp.mat[i][j]+=(a.mat[i][k]%mod*b.mat[k][j]%mod)%mod;
tmp.mat[i][j]%=mod;
}
}
}
tmp.n = a.n,tmp.m=b.m;
return tmp;
}
friend istream& operator>>(istream&is,MyMat& M){
for(int i=1;i<=M.n;i++){
for(int j=1;j<=M.m;j++){
is>>M.mat[i][j];
}
}
return is;
}
friend ostream& operator<<(ostream&os,const MyMat& M){
for(int i=1;i<=M.n;i++){
for(int j=1;j<=M.m;j++){
if(j!=1)os<<" ";
os<<M.mat[i][j];
}
os<<"\n";
}
return os;
}
};
MyMat quickpow(MyMat a,MyMat b,int x) {
while(x>0) {
if(x&1) a=b*a;
b=b*b;
x>>=1;
}
return a;
}
int main(){
MyMat a(2,3),b(3,2);
cin>>a>>b;
auto ans = a*b;
cout<<"a:\n"<<a<<"*\nb:\n"<<b<<"=\n"<<ans;
ans = quickpow(a,ans,5);
cout<<"a*ans^5=\n"<<ans;
}
感興趣的可以當一個參考。下面是本地執行測試輸出
g++ -o test test.cpp
./test
1 2 1 2 1 2
1 2 3 1 2 3
a:
1 2 1
2 1 2
*
b:
1 2
3 1
2 3
=
9 7
9 11
a*ans^5=
2480048 2480080 2480048
3188656 3188624 3188656
解題
斐波那契數列
讓我們回到經典的斐波那契數列，寫出這個矩陣遞推式後，加上矩陣快速冪
\[\begin{pmatrix}f(n)\\f(n-1) \\\end{pmatrix}=
\begin{pmatrix}1&1\\1&0\\\end{pmatrix}^{n-1} \\ \dot\\
\begin{pmatrix}f(1)\\f(0)\end{pmatrix}
\]
我們就能快速地寫出斐波那契數列的程式碼了，和上面不同的地方只有quickpow函式的x改為了long long：
MyMat quickpow(MyMat a,MyMat b,long long x) {
while(x>0) {
if(x&1) a=b*a;
b=b*b;
x>>=1;
}
return a;
}
int main(){
long long n;
cin>>n;
if(n<2)cout<<"1"<<endl;
else {
MyMat A(2,2);
A.mat[1][1]=1;
A.mat[1][2]=1;
A.mat[2][1]=1;
MyMat x(2,1);
x.mat[1][1]=1;
x.mat[2][1]=1;
auto ans = quickpow(x,A,n);
cout<<ans.mat[1][1]<<endl;
}
return 0;
}
A Simple Math Problem
重點：
這個矩陣我就不寫了，自己動手嘗試一下吧～
Count
重點：
\[\left\{
\begin{array}{**lr**}
f(n) = f(n-1)+2f(n-2)+n^3 \\n^3 = (n-1)^3-3n^3-3n-1\\n^2 = (n-1)^2+2n-1\\n = n-1
\end{array}
\right.
\]
\[\begin{pmatrix}f(n)\\ f(n-1) \\ n^3\\ n^2\\ n\\\ 1 \end{pmatrix}=
\begin{pmatrix}1&2&1&3&3&1\\1&0&0&0&0&0\\0&0&1&3&3&1\\0&0&0&1&2&1\\0&0&0&0&1&1\\0&0&0&0&0&1 \end{pmatrix}\\ \dot\\
\begin{pmatrix}f(n-1)\\f(n-2)\\(n-1)^3\\(n-1)^2\\n-1\\1\end{pmatrix}
\]
//AC程式碼
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=123456789;
int t;
ll n;
struct Matrix {
ll mat[25][25];
} pre,A;;
Matrix operator * (Matrix a,Matrix b) {
Matrix tmp;
memset(tmp.mat,0,sizeof(tmp.mat));
for(int i=1; i<=6; i++) {
for(int j=1; j<=6; j++) {
for(int k=1; k<=6; k++) {
tmp.mat[i][j]+=(a.mat[i][k]%mod*b.mat[k][j]%mod)%mod;
tmp.mat[i][j]%=mod;
}
}
}
return tmp;
}
Matrix quickpow(Matrix a,Matrix b,ll x) { //A*B^x
while(x>0) {
if(x&1) a=a*b;
b=b*b;
x>>=1;
}
return a;
}
void init() {
memset(pre.mat,0,sizeof(pre.mat));
memset(A.mat,0,sizeof(A.mat));
//1 2 27 9 3 1
pre.mat[1][1]=1;
pre.mat[1][2]=2;
pre.mat[1][3]=27;
pre.mat[1][4]=9;
pre.mat[1][5]=3;
pre.mat[1][6]=1;
A.mat[1][2]=2;
A.mat[2][1]=1;
A.mat[2][2]=1;
A.mat[3][2]=1;
A.mat[3][3]=1;
A.mat[4][3]=3;
A.mat[4][4]=1;
A.mat[5][3]=3;
A.mat[5][4]=2;
A.mat[5][5]=1;
A.mat[6][3]=1;
A.mat[6][4]=1;
A.mat[6][5]=1;
A.mat[6][6]=1;
}
int main() {
init();
scanf("%d",&t);
while(t--) {
scanf("%lld",&n);
Matrix ans=quickpow(pre,A,n-1);
printf("%lld\n",ans.mat[1][1]);
}
return 0;
}