矩陣快速冪

用來求矩陣的n次方，方法與快速冪相似，所以我們先寫個快速冪

求A^b次方:









int powf(int a,int b){
    int ans=1;
    while(b){
        if(b & 1) ans*=a;
        a*=a;
        b>>=1;
    }
    return ans;
}

矩陣快速冪

矩陣快速冪只是把整數換成矩陣而已，其概念是一樣的，差別在於要寫一個矩陣相乘的函式。

求(矩陣a)^b:

先副上Code後面再來解釋




























#define matrix vector<vector<int>> 
#define N 矩陣大小(N*N)

matrix matrix_multiply(matrix &a,matrix &b){  //矩陣相乘函式
    matrix ans(N,vector<int>(N));
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            for(int k=0;k<N;k++){
                ans[i][j]+=a[i][k]*b[k][j];
            }
        }
    }
    return ans;
}

matrix matrix_powf(matrix a,int b){
    matrix ans(N,vector<int>(N));
    for(int i=0;i<N;i++){
        ans[i][i]=1;
    }
    
    while(b){
        if(b & 1) ans=matrix_multiply(ans,a);
        a=matrix_multiply(a,a);
        b>>=1;
    }
    return ans;
}

在matrix_powf()中，有3點被更改，分別是

初始化ans

在快速冪中，ans初始化為1，因為要把b在二進位表示法中最後一位為1時，把a乘到答案裡，但在矩陣快速冪裡，初始化的ans就要換乘矩陣裡的1，也就是[[1,0];[0,1]] (以下都稱作A)，驗算一下就可得知



1  2   1  0    1*1+2*0  1*0+2*1   1  2
     x       =                  =
3  4   0  1    3*1+4*0  3*0+4*1   3  4

所以任何2*2的矩陣乘上A都會是自己，藉由觀察A的1是由左上往右下斜的，所以在初始化可先把A設為0，再用for跑一遍把斜線上給的點都設成1。




matrix ans(N,vector<int>(N,0));
for(int i=0;i<N;i++){
    ans[i][i]=1;
}

把答案乘到ans裡

在原本快速冪的if(b & 1) ans*=a中，要把a乘到答案裡很簡單，但在矩陣中，要寫一個函式把矩陣相乘的結果傳回來，再把它賦予給原本的矩陣。


if(b & 1) ans=matrix_multiply(ans,a);

矩陣平方

同上的概念,也是把矩陣丟到matrix_multiply，再把回傳的結果賦予給b。


b=matrix_multiply(b,b);

費式數列

\(f(0)=0,f(1)=1,f(2)=1,f(3)=2,f(4)=3,f(5)=5\)

如果要求費式數列有很多方法：

top down




int f(int n){
    if(n==0 || n==1) return 1;
    return f(n-1)+f(n-2);
}

botton up










int f(int n){
    vector<int> fi(n+1);
    fi[0]=1;
    fi[1]=1;
    for(int i=2;i<=n;i++){
        fi[i]=fi[i-1]+fi[i-2];
    }
    return fi[n]
}

但如果要求費式數列第1e9個mod m呢? 如果用上述兩種方法時間複雜度肯定會爆掉，這時候矩陣快速冪就派上用場了。

我們可以先把f1和f2寫成矩陣的形式(以下稱為初始矩陣)


[f(1) f(2)]

這時當我們乘上矩陣[[0,1];[1,1]] (以下都稱為A)的話



              0  1   
[f(1) f(2)] x      = [f(1)*0+f(2)*1 f(1)*1+f(2)*1] = [f(2) f(3)]
              1  1

如果再乘以A



              0  1   
[f(2) f(3)] x      = [f(2)*0+f(3)*1 f(2)*1+f(3)*1] = [f(3) f(4)]
              1  1

如果要得到[f(3) f(4)]就得把[f(1) f[2]]乘上A的2次方



              0  1   0  1  
[f(1) f(2)] x      x       = [f(3) f(4)]
              1  1   1  1

每多乘A一次，數列都會往後移一格
結論 : 如果想得到[f(n) f(n+1)]就必須把[f(1) f(2)]乘上A^(n-1)

最後就回到了上面的矩陣快速冪，但為了方便撰寫，初始矩陣可以從費式數列第0項開始，A就可以很直觀的乘上n次而非n-1，這樣最後初始矩陣乘上A的時候，A[1][0]就是答案了。










//R=A的n次方
                              R[0][0] R[0][1]
[f(n) f(n+1)] = [f(0) f(1)] x                  
                              R[1][0] R[1][1]
                              
= [f(0)*R[0][0]+f(1)*R[1][0] f(0)*R[0][1]+f(1)*R[1][1]]

= [0*R[0][0]+1*R[1][0] 0*R[0][1]+1*R[1][1]]

= [R[1][0] R[1][1]]

完整Code












































#include <bits/stdc++.h>
#define matrix vector<vector<int>> 
using namespace std;

matrix matrix_multiply(matrix &a,matrix &b){  
    matrix ans(2,vector<int>(2));
    for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
            for(int k=0;k<2;k++){
                ans[i][j]+=a[i][k]*b[k][j];
            }
        }
    }
    return ans;
}

matrix matrix_powf(matrix a,int b){
    matrix ans(2,vector<int>(2));
    for(int i=0;i<2;i++){
        ans[i][i]=1;
    }
    while(b){
        if(b & 1) ans=matrix_multiply(ans,a);
        a=matrix_multiply(a,a);
        b>>=1;
    }
    return ans;
}

int f(int n){
    //建立A模板
    matrix A(2,vector<int>(2,1));
    A[0][0]=0;
    
    matrix result=matrix_powf(A,n);
    return result[1][0];
}
int main(){
    int n;
    while(cin >> n){
        cout << f(n) << '\n';
    }
    return 0;
}

題目:Zerojudge : b525 先別管這個了，你聽過turtlebee嗎？

費氏數列大數實作

當第n項沒有到很大時，矩陣快速冪反而比較慢

Code













































































































#include <bits/stdc++.h>
#define matrix vector<vector<string>>
using namespace std;

string add(string a,string b){
    reverse(a.begin(),a.end());
    reverse(b.begin(),b.end());
    
    string ans="";
    int rest=0;
    for(int i=0;i<min(a.size(),b.size());i++){
        int cnt=a[i]-'0'+b[i]-'0'+rest;
        rest=0;
        if(cnt>9){
            rest=cnt/10;
            cnt=cnt%10;
        }
        ans+=(cnt+'0');
    }
    for(int i=min(a.size(),b.size());i<max(a.size(),b.size());i++){
        int cnt=rest;
        if(a.size() > b.size()) cnt+=a[i]-'0';
        else cnt+=b[i]-'0'; 
        rest=0;
        if(cnt>9){
            rest=cnt/10;
            cnt=cnt%10;
        }
        ans+=(cnt+'0');
    }
    if(rest!=0) ans+=(rest+'0');
    reverse(ans.begin(),ans.end());
    return ans;
}

string mul(string a,string b){
    reverse(a.begin(),a.end());
    reverse(b.begin(),b.end());

    string ans;
    vector<int> ad(a.size()+b.size()+1,0);
    
    for(int i=0;i<a.size();i++){
        for(int j=0;j<b.size();j++){
            int cnt=(a[i]-'0')*(b[j]-'0');
            ad[i+j]+=cnt;
        }
    }

    for(int i=0,rest=0;i<ad.size();i++){
        int cnt=ad[i]+rest;
        rest=0;
        if(cnt>9){
            rest=cnt/10;
        }
        ans+=(cnt%10)+'0';
    }
    int i=ans.size()-1;
    for(;i>=0;i--){
        if(ans[i]!='0') break;
    }
    ans=ans.substr(0,i+1);

    reverse(ans.begin(),ans.end());
    return ans;
}

matrix maxtrix_mul(matrix& a,matrix& b){
    matrix ans(2,vector<string>(2,"0"));
    for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
            for(int k=0;k<2;k++){
                ans[i][j]=add(ans[i][j],mul(a[i][k],b[k][j]));
            }
        }
    }
    return ans;
}


matrix matrix_pow(matrix& A,int n){
    matrix ans(2,vector<string>(2,"0"));
    ans[0][0]="1";ans[1][1]="1";
    while(n){
        if(n & 1){
            ans=maxtrix_mul(ans,A);
        }
        n >>= 1;
        A=maxtrix_mul(A,A);
    }
    return ans;
}

string fi(int n){
    if(n==0) return "0";
    
    matrix A(2,vector<string>(2,"1"));
    A[0][0]="0";
    matrix result=matrix_pow(A,n);
    return result[1][0];
}

int main(){
    int n;
    while(cin >> n){
        cout << "The Fibonacci number for " << n <<  " is " << fi(n) << '\n';
    }    
    return 0;
}

費氏數列小知識

當題目要取F(i) mod 1e9+9時是有周期的，週期為333333336，f(333333336)=0
費氏數列有公式解，如下

進階運用

題目:e811: 3. 密碼產生器 (Password)

三階矩陣快速冪應用
根據公式:A(N)＝P×A(N-1)＋Q×A(N-2)＋R

  A2 0 0       P Q R       A1 0 0
[ A1 0 0 ] = [ 1 0 0 ] * [ A0 0 0 ]
  1  0 0       0 0 1       1  0 0
  
                      ^(n-1)
  An   0 0       P Q R       A1 0 0
[ An-1 0 0 ] = [ 1 0 0 ] * [ A0 0 0 ]
  1    0 0       0 0 1       1  0 0

Code:





































#include <bits/stdc++.h>
#define int long long
#define matrix vector<vector<int>>
using namespace std;

matrix mul(matrix &a,matrix &b){
    matrix ans(3,vector<int>(3));
    for(int i=0;i<3;i++){
        for(int j=0;j<3;j++){
            for(int k=0;k<3;k++){
                ans[i][j]=(ans[i][j]+a[i][k]*b[k][j])%(long long)1e9;
            }
        }
    }
    return ans;
}

signed main(){
    int P,Q,R,A0,A1,N;
    cin >> P >> Q >> R >> A0 >> A1 >> N;
    matrix ans(3,vector<int>(3));
    ans[0][0]=A1,ans[1][0]=A0,ans[2][0]=1;    
    
    matrix p(3,vector<int>(3));
    p[0][0]=P,p[0][1]=Q,p[0][2]=R,p[1][0]=1,p[2][2]=1;

    N--;
    while(N){
        if(N & 1){
            ans=mul(p,ans);
        }
        p=mul(p,p);
        N >>= 1;
    }
    cout << setw(4) << setfill('0') << ans[0][0] << '\n';
    return 0;
}

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`	在筆記中貼入程式碼
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.