矩阵加速

· · 个人记录

题单

突然发现矩阵还不会，昨天学了些，觉得网上资料太碎高深，于是决定自己写一篇罢

矩阵

在数学中，矩阵（Matrix）是一个按照长方阵列排列的复数或实数集合

复数实数什么的我们先不管，总之矩阵就是一堆数，像个二维数组

也就是说，我们要记录矩阵的长，宽，内部的元素

特别的，我们把一个长款相等的矩阵叫做方阵

矩阵加法

首先，两个矩阵能进行一般的加法运算是有前提的

那就是规定两个矩阵是同形矩阵（即矩阵的行数和列数都相同）

我们只要把对应位置的值相加即可

\begin{bmatrix}1&9\\0&8\\7&1\end{bmatrix}+\begin{bmatrix}1&1\\4&5\\1&4\end{bmatrix}=\begin{bmatrix}2&10\\4&13\\8&5\end{bmatrix}

显然矩阵加法满足交换律和结合律，即

A+B=B+A A+(B+C)=(A+B)+C

更详细的定义或者想要学直和还是看这里吧

矩阵减法

把矩阵加法中相加改为相减即可

性质基本和加法相同

\begin{bmatrix}1&9\\0&8\\7&1\end{bmatrix}-\begin{bmatrix}1&1\\4&5\\1&4\end{bmatrix}=\begin{bmatrix}0&8\\-4&3\\6&-3\end{bmatrix}

数乘

即一个数乘一个矩阵

在实数运算中我们并没有数乘这种运算（毕竟本身就是数，直接叫乘法了）

在矩阵中只要只需要将矩阵中的每个元素乘上那个数就可以了

\begin{bmatrix}1&9\\0&8\\7&1\end{bmatrix}\times 4=\begin{bmatrix}4&36\\0&32\\28&4\end{bmatrix}

性质(k是数字，A是矩阵)：

(k_1 \times k_2) \times A=k_1 \times (k_2 \times A) (k_1 + k_2) \times A=k_1 \times A + k_2 \times A k \times (A+B) = k \times A + k \times B

矩阵乘法

即矩阵乘矩阵

要注意

矩阵乘法的前提是前一个矩阵的列数等于后一个矩阵的行数

？

比如A(n * m) B(m * k) C(k * n)

那么只有A\times B,B\times C,C\times A可以，其他都不行~

(也就是说A\times A也是不行哒)

所以我们也可以知道矩阵乘法不满足交换律

（有时候一交换乘都乘不了）

好，那么矩阵怎么乘？

令A为一个n*k的矩阵，B为一个k*m的矩阵

那么它们的乘积C就是一个n*m的矩阵

C_{i,j}=\sum\limits_{p=1}^{k}{A_{i,p}*B_{p,j}}

呃呃诶？

\begin{bmatrix}1*1+1*0+4*7&1*9+1*8+4*1\\5*1+1 *0+4*7&5*9+1*8+4*1\end{bmatrix} =\begin{bmatrix}29&21\\33&57\end{bmatrix}

(建议不太懂的多看几遍，最好手写几遍)

矩阵乘法满足的性质：

(AB)C=A(BC)

(结合律)

(A+B)C=AB+AC C(A+B)=CA+CB

单位矩阵

在普通乘法中，一个数乘 1 还是它本身，而在矩阵乘法中，也有一个这样的矩阵叫单位矩阵

不过值得注意的是对于不同矩阵单位矩阵大小不同

单位矩阵大小为列数 * 列数

即 n*m 的矩阵，对应的单位矩阵大小为 m*m，m*n 的矩阵单位矩阵大小为 n*n

那么单位矩阵肯定是正方形，因为只有这样才能保证矩阵大小不变

那单位矩阵的具体元素是什么？

单位矩阵元素非0即1，当且仅当行列数相同（在左上-右下对角线上）时元素为1，其余为0

差不多是这个样子：

\begin{bmatrix}1&0&0&0&0\\0&1&0&0&0\\0&0&1&0&0\\0&0&0&1&0\\0&0&0&0&1\end{bmatrix}

矩阵没有除法

但是有一个东东叫矩阵的逆

下面这部分和本文主要所讲关系不大，不感兴趣的可以略过直接空降至矩阵快速幂

前置芝士：高斯消元

逆矩阵的定义：

假设 A 是一个矩阵，I 是单位矩阵，如果存在一个矩阵 A^{-1}，使得 A^{-1}A=I

并且 AA^{-1}=I

那么，矩阵 A 就是可逆的，A^{-1} 称为 A 的逆矩阵

若 A 可逆，则逆唯一

求逆思路：

  1. 把 A 和它等大的单位矩阵 I 放在一个矩阵里

  2. 对 A 进行高斯消元使 A 化成单位矩阵

  3. 此时原来单位矩阵转化成逆矩阵

原理：

[AI] => [I A^{-1}]

由于是在同一矩阵中，因此 A 和 I 得到的操作都是一样的。我们只需要让前面的 A 变换为 I ，那么对于同一过程，I 就会变成矩阵 A 的逆

P4783 【模板】矩阵求逆

矩阵快速幂

P3390 【模板】矩阵快速幂

题目描述

给定 n\times n 的矩阵 A，求 A^k，每个元素对 10^9+7 取模。

对于 100\% 的数据，1\le n \le 100，0 \le k \le 10^{12}，|A_{i,j}| \le 1000。

矩阵乘矩阵我们会，只要把幂拆成乘法就...

k \le 10^{12}

肯定超时！所以我们肯定要优化，如果是求一个数的幂我们可以用快速幂，那矩阵可不可以呢？

题目都告诉你了

快速幂本质上运用了乘法的结合律，但矩阵乘法同样满足结合律

所以没问题！

为了方便，我封装了几个函数而且重载了运算符

还有，别忘开long long

自认为码风不错的我献上代码：

#include<bits/stdc++.h>
using namespace std;
const int N=105;
const int P=1e9+7;
typedef long long ll;

int n;
ll k;
struct node{
    ll a[N][N];
    void _0(){
        memset(a,0,sizeof(a));
    }
    void _1(){
        _0();
        for(int i=1;i<=n;i++)
            a[i][i]=1;
    }
    void read(){
        _0();
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++){
                scanf("%lld",&a[i][j]);
            }
        }
    }
    void print(){
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++){
                printf("%lld ",a[i][j]);
            }
            printf("\n");
        }
    }
}A;
node operator *(const node x,const node y){
    node z;
    z._0();
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++){
            for(int k=1;k<=n;k++)
                z.a[i][j]=(z.a[i][j]+x.a[i][k]*y.a[k][j])%P;
        }
    return z;
}
node ksm(node a,ll b){
    node res;
    res._1();
    for(;b;b>>=1,a=a*a)
        if(b&1) res=res*a;
    return res;
}
int main(){
    scanf("%d%lld",&n,&k);
    A.read();
    A=ksm(A,k);
    A.print();
}

矩阵加速

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 取余的值。

显然不能暴力推，那么肯定要优化，怎么优化呢？

呃呃呃想不出来

想想我们刚刚做的题，它是怎么优化的？它利用矩阵和数字都能用的快速幂优化

但在这题数字显然不能直接快速幂

那这题可不可以把数列转化成矩阵再用矩阵快速幂优化呢？

考虑构造一个矩阵记录a_{x-1},a_{x-2},a_{x-3}

那么我们下一步应该推出a_x,a_{x-1},a_{x-2}

\begin{cases}a_x=a_{x-1}+a_{x-3}\\a_{x-1}=a_{x-1}\\a_{x-2}=a_{x-2}\end{cases}

转化成矩阵就是：

\begin{bmatrix}1&0&1\\1&0&0\\0&1&0\end{bmatrix}

初始矩阵就是：

\begin{bmatrix}a_3\\a_2\\a_1\end{bmatrix}

这样就可以用矩阵快速幂优化啦~

#include<bits/stdc++.h>
using namespace std;
const int N=105;
const int P=1e9+7;
typedef long long ll;

int n=3;
ll k;
struct node{
    ll a[N][N];
    void _0(){
        memset(a,0,sizeof(a));
    }
    void _1(){
        _0();
        for(int i=1;i<=n;i++)
            a[i][i]=1;
    }
    void init(){
        _0();
        a[1][1]=a[1][3]=a[2][1]=a[3][2]=1;
    }
    void print(){
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++){
                printf("%lld ",a[i][j]);
            }
            printf("\n");
        }
    }
}A;
node operator *(const node x,const node y){
    node z;
    z._0();
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++){
            for(int k=1;k<=n;k++)
                z.a[i][j]=(z.a[i][j]+x.a[i][k]*y.a[k][j])%P;
        }
    return z;
}
node ksm(node a,ll b){
    node res;
    res._1();
    for(;b;b>>=1,a=a*a)
        if(b&1) res=res*a;
    return res;
}
int main(){
    int T,k;
    scanf("%d",&T);
    while(T--){
        scanf("%d",&k);
        node A,B;
        A._0();
        A.a[3][1]=1;//a3
        A.a[2][1]=1;//a2
        A.a[1][1]=1;//a1
        B.init();
        A=A*ksm(B,k);
//      A.print();
//其中A.a[1][3]是a_k A.a[1][2]是a_{k+1} A.a[1][1]是a_{k+2} 
        printf("%lld\n",A.a[1][3]);
    }
}

相关习题

P2044 [NOI2012] 随机数生成器

题目描述

四个非负整数参数 m,a,c,X_0，按照下面的公式生成出一系列随机数 \{X_n\}：

X_{n+1}=(aX_n +c)\bmod m

其中 \bmod m 表示前面的数除以 m 的余数。求 X_n \bmod g

对于 100\% 的数据，a,c,X_0 是非负整数，m,n,g 是正整数。n,m,a,c,X_0\leq 10^{18}，1\leq g\leq 10^8，n,m\geq 1，a,c,X_0\geq 0。

题目就是说形如 x_{i+1}=ax_i+c 的式子

显然也要推柿子

考虑怎么把这东东变成矩阵

在查阅了资料后可以得知

构造如下初始矩阵(x_0可能没有意义)

\begin{bmatrix}x_1&1\\x_0&1\end{bmatrix}

构造以下转移矩阵

\begin{bmatrix}a&0\\c&1\end{bmatrix}

因为题目要求第 n 项，那么只要将转移矩阵快速幂 n 次方，再将初始矩阵乘以快速幂得到的结果矩阵就好啦~

注意一下细节，要求变换 n 次，那要输出结果矩阵的 [2][1]

但是 n,m,a,c,X_0\leq 10^{18}

怎么了？我记得long long存的下 10^{18} 的呀

但你有乘法啊，一乘不就炸了

高精警告

为了避免代码变得过于冗长，我们也可以在这种情况下避免高精

那就是龟速乘！

利用快速幂的思想把乘法拆成一步一步加法，这样就可以一遍算一遍取模

利用这种方法可以在模数不炸long long时不用高精

完了完了我已经不会高精了考试考了怎么办

#include<bits/stdc++.h>
using namespace std;
const int N=105;
typedef long long ll;

int n=2;
ll k,m,X0,_a,_c;
struct node{
    ll a[N][N];
    void _0(){
        memset(a,0,sizeof(a));
    }
    void _1(){
        _0();
        for(int i=1;i<=n;i++)
            a[i][i]=1;
    }
    void init(){
        _0();
        a[1][1]=_a;
        a[2][1]=_c;
        a[2][2]=1;
    }
    void print(){
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++){
                printf("%lld ",a[i][j]);
            }
            printf("\n");
        }
    }
}A;
ll mul(ll x,ll y){
    ll res=0;
    while(y){
        if(y&1) res=(res+x)%m;
        x=(x+x)%m;
        y>>=1;
    }
    return res;
}
node operator *(const node x,const node y){
    node z;
    z._0();
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++){
            for(int k=1;k<=n;k++)
                z.a[i][j]=(z.a[i][j]+mul(x.a[i][k],y.a[k][j]))%m;
        }
    return z;
}
node ksm(node a,ll b){
    node res;
    res._1();
    for(;b;b>>=1,a=a*a)
        if(b&1) res=res*a;
    return res;
}

int main(){
    ll n,g;
    scanf("%lld%lld%lld%lld%lld%lld",&m,&_a,&_c,&X0,&n,&g);
    node A,B;
    A._0();
    A.a[1][1]=(mul(_a,X0)+_c)%m;
    A.a[2][1]=X0%m;
    A.a[1][2]=A.a[2][2]=1;
    B.init();
    A=A*ksm(B,n);
//  A.print();
    printf("%lld\n",A.a[2][1]%g);//变换n次从0->n 
}

相关习题

更多进阶详见 OI-wiki