矩阵加速递推学习笔记

· · 个人记录

之前写的太naive了。

因此今日重新写一篇。

前置知识：矩阵乘法

矩阵的一系列操作

直接贴代码吧..

只是一个大概的模板。具体根据题目而定。

class Matrix{
  private:
    int row,col,s[N][N];    //行，列，元素
  public:
    Matrix() {      //无参数构造函数
        row=col=0;
        memset(s,0,sizeof(s));
    }
    Matrix(int n,int m):row(n),col(m){memset(s,0,sizeof(s));}       //重载构造函数
    void readp() {      //输入一个矩阵
        for(int i=1;i<=row;i++) {
            for(int j=1;j<=col;j++) cin>>s[i][j];
        }
    }
    void printp() {     //输出一个矩阵
        for(int i=1;i<=row;i++) {
            for(int j=1;j<=col;j++) cout<<s[i][j]<<" ";
            cout<<endl;
        }
    }
    void transpose() {      //矩阵的转置
        Matrix y(col,row);
        for(int i=1;i<=row;i++) {
            for(int j=1;j<=col;j++) y.s[j][i]=s[i][j];
        }
        *this=y;
    }
    void identity() {       //将该矩阵改成单位矩阵
        for(int i=1;i<=row;i++) s[i][i]=1;
    }
    Matrix operator + (Matrix x) {      //重载+
        Matrix y(row,col);
        for(int i=1;i<=row;i++) {
            for(int j=1;j<=col;j++) y.s[i][j]+=s[i][j]+x.s[i][j];
        }
        return y;
    }
    Matrix operator * (Matrix x) {      //重载*
        Matrix y(row,x.col);
        for(int i=1;i<=row;i++) {
            for(int j=1;j<=x.col;j++) {
                for(int k=1;k<=col;k++) {
                    y.s[i][j]+=s[i][k]*x.s[k][j];
                }
            }
        }
        return y;
    }
    Matrix operator *= (Matrix x) {     //重载*=
        return *this=*this*x;
    }
};

矩阵乘法

重新温故一下矩阵乘法。

举个例子

A= \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ \end{bmatrix}, B= \begin{bmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \\ b_{31} & b_{32} \end{bmatrix}

于是

C=AB=\begin{bmatrix} a_{11}b_{11}+a_{12}b_{21}+a_{13}b_{31} & a_{11}b_{12}+a_{12}b_{22}+a_{13}b_{32} \\ a_{21}b_{11}+a_{22}b_{21}+a_{23}b_{31} & a_{21}b_{12}+a_{22}b_{22}+a_{23}b_{32} \end{bmatrix}

注意，只有满足 A 为 n\times m 的矩阵，B 为 m\times p 的矩阵，矩阵 A 与矩阵 B 才能相乘（AB）。所得的矩阵 C 为 n\times p 的矩阵。

$$ C_{ij}=\sum_{k=1}^ma_{ik}b_{kj} $$ 矩阵乘法满足结合律，即 $(AB)C=A(BC)$. 矩阵乘法满足分配率，即 $(A+B)C=AC+BC$，$C(A+B)=CA+CB$. 矩阵乘法对数乘满足结合律，即 $\lambda(AB)=(\lambda A)B=A(\lambda B)$. 对于转置矩阵也有性质：$(AB)^T=B^TA^T$. 但是要注意的是矩阵乘法**不满足交换律**，即有时可以 $AB$ 却无法 $BA$，因为要求**第一个矩阵的列数与第二个矩阵的行数相等**。 **矩阵满足结合律**，这是矩阵加速递推的一个关键。 ## [洛谷P3390]【模板】矩阵快速幂 给定 $n\times n$ 的矩阵 $A$，求 $A^k$。 ## 数据范围 $1\le n\le 100$，$0\le k\le 10^12$，$|A_{ij}|\le 1000$。 ## 链接 [P3390 【模板】矩阵快速幂](https://www.luogu.com.cn/problem/P3390) ## 思路 因为是 $n\times n$ 的，且矩阵满足结合律，因此可以进行快速幂。 ## 代码 ```cpp #include<bits/stdc++.h> #define N 110 #define mod 1000000007 #define ll long long using namespace std; int n; ll k; class Matrix{ private: int len; ll s[N][N]; public: Matrix() { len=0; memset(s,0,sizeof(0)); } Matrix(int _len):len(_len){memset(s,0,sizeof(s));} Matrix operator * (Matrix x) { Matrix y(len); for(int i=1;i<=len;i++) { for(int j=1;j<=len;j++) { for(int k=1;k<=len;k++) y.s[i][j]=(y.s[i][j]+s[i][k]*x.s[k][j])%mod; } } return y; } Matrix operator *= (Matrix x) { return *this=*this*x; } void idenitity() { for(int i=1;i<=len;i++) s[i][i]=1; } void readp() { for(int i=1;i<=len;i++) { for(int j=1;j<=len;j++) cin>>s[i][j]; } } void printp() { for(int i=1;i<=len;i++) { for(int j=1;j<=len;j++) cout<<s[i][j]<<" "; cout<<endl; } } }; Matrix qpow(Matrix x,ll y) { Matrix a(n); Matrix b(n); a=x; b.idenitity(); while(y) { if(y&1) b*=a; a*=a; y>>=1; } return b; } int main() { cin>>n>>k; Matrix A(n); A.readp(); qpow(A,k).printp(); return 0; } ``` ------ # 矩阵加速递推 斐波那契数列，即 $f_1=1,f_2=1,f_n=f_{n-1}+f_{n-2},n\ge 2$。 $O(n)$ 的递推是很显然的。 那如果求的项数很大呢？ ## [洛谷P1962] 斐波那契数列 大家都知道，斐波那契数列是满足如下性质的一个数列： $$ F_n = \left\{\begin{aligned} 1 \space (n \le 2) \\ F_{n-1}+F_{n-2} \space (n\ge 3) \end{aligned}\right. $$ 请你求出 $F_n\bmod 10^9+7$ 的值。 ## 数据范围 $1\le n\le 2^63

链接

P1962 斐波那契数列

思路

考虑这样一个矩阵

\begin{bmatrix} F_{n-1} \\ F_{n-2} \end{bmatrix}

如何变成

\begin{bmatrix} F_n \\ F_{n-1} \end{bmatrix}

呢？

显然

\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} F_{n-1} \\ F_{n-2} \end{bmatrix} = \begin{bmatrix} F_{n} \\ F_{n-1} \end{bmatrix}

于是若要求 F_n

\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n-2} \begin{bmatrix} F_{2} \\ F_{1} \end{bmatrix} = \begin{bmatrix} F_{n} \\ F_{n-1} \end{bmatrix}

于是便可以矩阵快速幂+矩阵乘法解决。

时间复杂度 O(2^3 \log n)，可以通过。

代码

#include<bits/stdc++.h>
#define N 5
#define ll long long 
#define mod 1000000007

using namespace std;

ll n;
class Matrix{
  private:
    int row,col;
    ll s[N][N];
  public:
    Matrix() {
        row=col=0;
        memset(s,0,sizeof(s));
    }
    Matrix(int _row,int _col):row(_row),col(_col){memset(s,0,sizeof(s));}
    Matrix(int n):row(n),col(n){memset(s,0,sizeof(s));}
    Matrix operator * (Matrix x) {
        Matrix y(row,x.col);
        for(int i=1;i<=row;i++) {
            for(int j=1;j<=x.col;j++) {
                for(int k=1;k<=col;k++) y.s[i][j]=(y.s[i][j]+s[i][k]*x.s[k][j])%mod;
            }
        }
        return y;
    }
    Matrix operator *= (Matrix x) {
        return *this=*this*x;
    }
    void idenitity() {
        for(int i=1;i<=row;i++) s[i][i]=1;
    }
    void printp() {
        cout<<s[1][1]<<endl;
    }
    void init(int x[][N]) {
        for(int i=1;i<=row;i++) {
            for(int j=1;j<=col;j++) s[i][j]=x[i][j];
        }
    }
}Init(2),seq(2,1);

void init() {
    int I[N][N],S[N][N];
    memset(I,0,sizeof(I));
    memset(S,0,sizeof(S));
    I[1][1]=I[1][2]=I[2][1]=S[1][1]=S[2][1]=1;
    Init.init(I);
    seq.init(S);
}

Matrix qpow(Matrix x,ll y) {
    Matrix a;
    Matrix b(2);
    a=x;
    b.idenitity();
    while(y) {
        if(y&1) b*=a;
        a*=a;
        y>>=1;
    }
    return b;
}

int main() {
    cin>>n;

    if(n<=2) cout<<1<<endl;
    else {
        init();
        Init=qpow(Init,n-2);
        seq=Init*seq;
        seq.printp();
    }

    return 0;
}

如果不是斐波那契数列呢？

[洛谷P1939]【模板】矩阵加速（数列）

已知一个数列 a，它满足：

a_x= \begin{cases} 1 & x \in\{1,2,3\}\\ a_{x-1}+a_{x-3} & x \geq 4 \end{cases}

求 a 数列的第 n 项对 10^9+7 取余的值。

数据范围

数据组数 1\le T\le 100

1\le n\le 2\times 10^9

链接

P1939 【模板】矩阵加速（数列）

思路

考虑这样一个矩阵

\begin{bmatrix} a_{n-1} \\ a_{n-2} \\ a_{n-3} \end{bmatrix}

显然

\begin{bmatrix} 1 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} a_{n-1} \\ a_{n-2} \\ a_{n-3} \end{bmatrix} = \begin{bmatrix} a_{n} \\ a_{n-1} \\ a_{n-2} \end{bmatrix}

于是可以矩阵快速幂加速。即

\begin{bmatrix} 1 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{bmatrix}^{n-3} \begin{bmatrix} a_{3} \\ a_{2} \\ a_{1} \end{bmatrix} = \begin{bmatrix} a_{n} \\ a_{n-1} \\ a_{n-2} \end{bmatrix}

代码

#include<bits/stdc++.h>
#define N 5
#define ll long long 
#define mod 1000000007

using namespace std;

int T;
class Matrix{
  private:
    int row,col;
    ll s[N][N];
  public:
    Matrix() {
        row=col=0;
        memset(s,0,sizeof(s));
    }
    Matrix(int _row,int _col):row(_row),col(_col){memset(s,0,sizeof(s));}
    Matrix(int n):row(n),col(n){memset(s,0,sizeof(s));}
    Matrix operator * (Matrix x) {
        Matrix y(row,x.col);
        for(int i=1;i<=row;i++) {
            for(int j=1;j<=x.col;j++) {
                for(int k=1;k<=col;k++) y.s[i][j]=(y.s[i][j]+s[i][k]*x.s[k][j])%mod;
            }
        }
        return y;
    }
    Matrix operator *= (Matrix x) {
        return *this=*this*x;
    }
    void idenitity() {
        for(int i=1;i<=row;i++) s[i][i]=1;
    }
    void printp() {
        cout<<s[1][1]<<endl;
    }
    void init(int x[][N]) {
        for(int i=1;i<=row;i++) {
            for(int j=1;j<=col;j++) s[i][j]=x[i][j];
        }
    }
}Init(3),seq(3,1);

void init() {
    int I[N][N],S[N][N];
    memset(I,0,sizeof(I));
    memset(S,0,sizeof(S));
    I[1][1]=I[1][3]=I[2][1]=I[3][2]=S[1][1]=S[2][1]=S[3][1]=1;
    Init.init(I);
    seq.init(S);
}

Matrix qpow(Matrix x,int y) {
    Matrix a;
    Matrix b(3);
    a=x;
    b.idenitity();
    while(y) {
        if(y&1) b*=a;
        a*=a;
        y>>=1;
    }
    return b;
}

int main() {
    cin>>T;

    init();
    while(T--) {
        int n;
        cin>>n;
        if(n<=3) cout<<1<<endl;
        else {
            Matrix a,b;
            a=qpow(Init,n-3);
            b=a*seq;
            b.printp();
        }
    }

    return 0;
}

推广

如果是更一般的数列呢？

考虑这样一个数列

f_n=af_{n-1}+bf_{n-2}+cn^3+dn^2+en+f

可以构造出一个矩阵

\begin{bmatrix} f_{n-1} \\ f_{n-2} \\ n^3 \\ n^2 \\ n \\ 1 \end{bmatrix}

显然

\begin{bmatrix} a & b & c & d & e & f \\ 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{bmatrix} \begin{bmatrix} f_{n-1} \\ f_{n-2} \\ n^3 \\ n^2 \\ n \\ 1 \end{bmatrix} = \begin{bmatrix} f_{n} \\ f_{n-1} \\ (n+1)^3 \\ (n+1)^2 \\ n+1 \\ 1 \end{bmatrix}

其实很好构造出来的。

右下角甚至是个杨辉三角

完。