矩阵快速幂（二分）

矩阵的定义

一个$n*m$矩阵是由$n*m$个实数排列成n行m列构成的，一般用一对[]括起来，通常用大写字母表示，比如下面$3 * 3$的矩阵简记为$A=(a_{ij})_{3*3}$

---

矩阵的加减

矩阵之间也可以加减，但是矩阵的加减要求是两个大小相同的矩阵，只有当两个矩阵的行列数相同时，才能进行加减，比如一个n行m列的A矩阵和一个n行m列的B矩阵的加减，就是对应位置元素的加减。

$A+B=(a_{ij}+b_{ij})_{n*m}$

$A-B=(a_{ij}-b_{ij})_{n*m}$

具体如下：

矩阵的加减满足交换律和结合律

1. 交换律：$A+B=B+A$
2. 结合律：$(A+B)+C=A+(B+C)$

---

矩阵的乘法

矩阵和数相乘

我们定义用实数k右乘矩阵A或者左乘矩阵A，其积等于矩阵A中每个元素乘上k。具体如下

$kA=Ak=\begin{bmatrix} ka_{11} &ka_{12} &...... &ka_{1m} \\ ka_{21} &ka_{22} & ......& ka_{2m}\\ ......& ...... & ...... &...... \\ ka_{n1}& ka_{n2}& ...... & ka_{nm} \end{bmatrix}$

矩阵和矩阵相乘

矩阵A乘上B当且仅当A矩阵的列数和矩阵B行数相同的时候才有定义，即只有$A_{nt}·B_{tm}$才有定义。他们相乘得到一个n行m列的矩阵$C_{n*m}$。C矩阵中每个元素的计算方法为$C_{ij}=a_{i1}b_{1j}+a_{i2}b_{2j}+...+a_{it}b_{tj}=\sum_{r=1}^{t}a_{ir}b_{rj}$，即$C_{ij}$等于A的第i行和B的第j列的每个元素对应相乘的和。比如$A=\begin{bmatrix}2 &3 \\ 1&5 \\ 4&6 \end{bmatrix} 　 B= \begin{bmatrix}1 &2 \\ 4&3 \\ \end{bmatrix}$$A\cdot B=\begin{bmatrix}2 &3 \\ 1&5 \\ 4&6 \end{bmatrix} \bullet \begin{bmatrix}1 &2 \\ 4&3 \\ \end{bmatrix}=\begin{bmatrix}14 &13 \\ 21&17 \\ 28 &26 \end{bmatrix}$

矩阵的乘法满足结合律但是不满足交换律

1. $A·(B+C)=A·B+A·C$
2. $(A·B)·C=A·(B·C)$
3. $A·B≠B·A$

代码示例

const int maxn=101;
const int maxm=101;

struct Matrix{
    int n,m;//行 列
    int a[maxn][maxm];
    void clear(){
        n=m=0;
        memset(a,0,sizeof(a));
    }

    Matrix operator +(const Matrix &b) const{
        Matrix tmp;
        tmp.n=n; tmp.m=m;
        for(int i=0;i<n;++i)
            for(int j=0;j<m;++j)
                tmp.a[i][j]=a[i][j]+b.a[i][j];
        return tmp;
    }

    Matrix operator -(const Matrix &b) const{
        Matrix tmp;
        tmp.n=n; tmp.m=m;
        for(int i=0;i<n;++i)
            for(int j=0;j<m;++j)
                tmp.a[i][j]=a[i][j]-b.a[i][j];
        return tmp;
    }

    Matrix operator *(const Matrix &b) const{
        Matrix tmp;
        tmp.clear();
        tmp.n=n; tmp.m=b.m;
        for(int i=0;i<n;++i)
            for(int j=0;j<b.m;++j)
                for(int k=0;k<m;++k)
                    tmp.a[i][j]+=a[i][k]*b.a[k][j];
        return tmp;
    }
};

---

矩阵的简单应用

矩阵经常被应用到递推、动态规划的转移当中，有一个二维状态数组$dp[n][m]$，有如下转移方程$dp[i][j]=a_{j1}dp[i-1][1]+a_{j2}dp[i-1][2]+...+a_{jm}dp[i-1][m]$。

实际上我们正好可以用一个矩阵乘法来表示这个过程，A矩阵是一个m*m的矩阵，我们一般叫做状态转移矩阵

$$\begin{bmatrix}dp[i][1]\\ dp[i][2]\\ ......\\ dp[i][m]\end{bmatrix}=\begin{bmatrix}a[1][1] &a[1][2] & ...... &a[1][m] \\ a[2][1]&a[2][2] &...... & a[2][m]\\ ......& ...... &...... &...... \\ a[m][1]& a[m][2] & ...... &a[m][m] \end{bmatrix}\begin{bmatrix}dp[i-1][1]\\ dp[i-1][2]\\ ......\\ dp[i-1][m]\end{bmatrix}$$

第二维都是1,2,…,m，于是可以简记成$dp[i]=A·dp[i-1]$
通过递推可以得到$dp[n]=A^{n-1}dp[1]$，即我们通过矩阵的乘法表示了dp的转移。
而对于乘法可以采用二分快速幂的方法，求得$A^{n-1}$,从而计算$dp[n]$。

---

矩阵快速幂

形如$dp[n]=A^{n-1}dp[1]$,如果按照一般递推式来正常得到，时间复杂度是$O(nm^2)$的，n是转移的次数，$O(m^2)$是一次转移的复杂度。
我们可以先求出$A^{n-1}$，同普通数的快速幂一样，矩阵的快速幂也是二分的思想，并注意取模操作，直接看代码

时间复杂度为$O(m^3log(n))$

$$A=\begin{bmatrix} 0& 1 & 0&...&0\\ 0& 0 &1 &...&0\\ .. & .. &..&...&..\\ 0& 0&0 &... &1\\ a_{n-1} & a_{n-2} &a_{n-3}&...&a_{0} \end{bmatrix} ,B=\begin{bmatrix} f_{x-n}\\ f_{x-n+1}\\ ..\\ ..\\ f_{x-2}\\ f_{x-1} \end{bmatrix}$$

代码示例

int solve(int a[],int b[],int m,int t){
    //a是常系数数组 b是初值数组 m是数组大小 t是要求解的项数 从第0项开始 所以共有t+1项
    //m为已知递推式的阶数
    //输出函数在第t项的值f(t)
    //调用矩阵类
    Matrix M,F,res;//M是辅助常数矩阵 F是转移矩阵
    M.clear(),F.clear(),res.clear();
    M.n=M.m=m;
    res.n=res.m=m;
    F.n=m,F.m=1;
    for(int i=0;i<m-1;++i)
        M.a[i][i+1]=1;
    for(int i=0;i<m;++i){
        M.a[m-1][i]=a[i];//a[i]先存小项的系数 即递推式最靠末尾的系数
        F.a[i][0]=b[i];//b[i]先存小项的值 即f0(通常)
        res.a[i][i]=1;//初始化为单位矩阵
    }
    if(t<m) return F.a[t][0];
    for(t-=m-1;t;t/=2){//t-=m-1为还要转移的次数
        if(t&1) res=M*res;
        M=M*M;
    }
    F=res*F;
    return F.a[m-1][0];
}

---

例题

用 fib(n) 表示斐波那契数列的第 n 项，现在要求你求 fib(n)mod m，fib(1)=1，fib(2)=1。

斐波那契的递推式是$f[n]=f[n-1]+f[n-2](n>2)$
我们构建转移矩阵A $A=\begin{bmatrix}0 &1 \\ 1&1 \end{bmatrix}$

于是有$\begin{bmatrix}f[i-1]\\ f[i]\end{bmatrix}=A\cdot \begin{bmatrix}f[i-2]\\ f[i-1]\end{bmatrix}=\begin{bmatrix}0 &1 \\1&1 \end{bmatrix}\begin{bmatrix}f[i-2]\\ f[i-1]\end{bmatrix}$

当然矩阵的数的位置可以调换（不唯一），编码时应注意

这样就能将复杂度优化为$O(2^3lgn)$

代码示例

#include<iostream>
#include<string.h>
using namespace std;

const int maxn=101;
const int maxm=101;
int mod;//需要时使用

struct Matrix{
    int n,m;//行 列
    int a[maxn][maxm];//注意函数内数组大小不能超过500000int
    void Clear(){
        n=m=0;
        memset(a,0,sizeof(a));
    }

    Matrix operator +(const Matrix &b) const{
        Matrix tmp;
        tmp.n=n; tmp.m=m;
        for(int i=0;i<n;++i)
            for(int j=0;j<m;++j)
                tmp.a[i][j]=a[i][j]+b.a[i][j];
        return tmp;
    }

    Matrix operator -(const Matrix &b) const{
        Matrix tmp;
        tmp.n=n; tmp.m=m;
        for(int i=0;i<n;++i)
            for(int j=0;j<m;++j)
                tmp.a[i][j]=a[i][j]-b.a[i][j];
        return tmp;
    }

    Matrix operator *(const Matrix &b) const{
        Matrix tmp;
        tmp.Clear();
        tmp.n=n; tmp.m=b.m;
        for(int i=0;i<n;++i)
            for(int j=0;j<b.m;++j)
                for(int k=0;k<m;++k){
                    tmp.a[i][j]+=a[i][k]*b.a[k][j]%mod;
                    tmp.a[i][j]%=mod;
                }
        return tmp;
    }
};

int solve(int a[],int b[],int m,int t){
    //a是常系数数组 b是初值数组 m是数组大小 t是要求解的项数 从第0项开始 所以共有t+1项
    //m为已知递推式的阶数
    //输出函数在第t项的值f(t)
    //调用矩阵类
    Matrix M,F,res;//M是辅助常数矩阵 F是转移矩阵
    M.Clear(),F.Clear(),res.Clear();
    M.n=M.m=m;
    res.n=res.m=m;
    F.n=m,F.m=1;
    for(int i=0;i<m-1;++i)
        M.a[i][i+1]=1;
    for(int i=0;i<m;++i){
        M.a[m-1][i]=a[i];//a[i]先存小项的系数 即递推式最靠末尾的系数
        F.a[i][0]=b[i];//b[i]先存小项的值 即f0(通常)
        res.a[i][i]=1;//初始化为单位矩阵
    }
    if(t<m) return F.a[t][0];
    for(t-=m-1;t;t/=2){//t-=m-1为还要转移的次数
        if(t&1) res=M*res;
        M=M*M;
    }
    F=res*F;//res为最后的常数矩阵
    return F.a[m-1][0];
}


int main()
{
    mod=1000000007;
    int t;
    int a[2],b[2];
    a[0]=a[1]=1;
    b[0]=0;
    b[1]=1;
    while(cin>>t&&t!=-1)
    {
        cout<<solve(a,b,2,t)<<endl;
    }
    return 0;
}