• [模板][题解][Luogu1939]矩阵乘法加速递推(详解)


    题目传送门

    题目大意:计算数列a的第n项,其中:

    [a[1] = a[2] = a[3] = 1 ]

    [a[i] = a[i-3] + a[i - 1] ]

    [(n ≤ 2 imes 10^9) ]

    一般的递推是O(n)的,显然时间和空间都不能承受。

    由于每一步递推都是相同的。这句话包含了2个层面:首先,递推式是相同的;其次,递推的条件也要是相同的。综合来说,每一步的递推都是相同的。这是应用矩阵加速递推的充分条件。

    那么怎么进行矩阵加速呢?我们首先观察,第(i)项和哪些项有关? 与(i-3)(i-1)优化,也就是和前3项有关。为了能够**仅仅利用矩阵就能推出(a[i]), 我们需要记录前3项,于是,我们构造一个3*3的矩阵:

    [A=egin{bmatrix} 1 & 0 & 1 \ 1 & 0 & 0 \ 0 & 1 & 0 \ end{bmatrix} ]

    有同学会问:这个矩阵是怎么构造出来的呢?

    我们首先要构造出类似于DP的状态来存下所有计算过程中可能会用到的信息,对于这道题,我们需要记录:(假设我们要从(a[i])推到(a[i+1]))

    [B=egin{bmatrix} a[i] \ a[i-1] \ a[i-2] \ end{bmatrix} ]

    这个矩阵要推到:

    [C=egin{bmatrix} a[i+1] \ a[i] \ a[i-1] \ end{bmatrix} ]

    也就是说,我们需要构造一个矩阵(A),使得(A*B = C),根据线性代数的相关定义,A一定是一个(3*3)的矩阵,没错吧。

    那好,我们已经得到:

    [egin{bmatrix} ? & ? & ? \ ? & ? & ? \ ? & ? & ? \ end{bmatrix}*egin{bmatrix} a[i] \ a[i-1] \ a[i-2] \ end{bmatrix} =egin{bmatrix} a[i+1] \ a[i] \ a[i-1] \ end{bmatrix} \A*B=C ]

    我们只需要根据递推式,把矩阵(A)中的数填满就可以了,比如说:

    由于$a[i+1] = a[i-2] +a[i] (,而根据矩阵,)a[i+1] = (0,0) * a[i] + (0,1) * a[i-1] + (0,2) * a[i-2]$,所以,矩阵的第一行是(1,0,1),以此类推,就可以吧矩阵填满了。

    然后,我们可以得到:

    [egin{bmatrix} 1 & 0 & 1 \ 1 & 0 & 0 \ 0 & 1 & 0 \ end{bmatrix}*egin{bmatrix} a[i] \ a[i-1] \ a[i-2] \ end{bmatrix} =egin{bmatrix} a[i+1] \ a[i] \ a[i-1] \ end{bmatrix} \A*B=C ]

    可是,有了这个,怎么从(a[1])推到(a[n])呢?

    我们知道:(a[1] = a[2] = a[3] = 1),如果把它们代入矩阵(B)(就是中间的那个矩阵),我们会得到:

    [egin{bmatrix} 1 & 0 & 1 \ 1 & 0 & 0 \ 0 & 1 & 0 \ end{bmatrix}*egin{bmatrix} a[3]=1 \ a[2]=1 \ a[1]=1 \ end{bmatrix} =egin{bmatrix} a[4]=1*a[3] + 1 * a[1] = 2 \ a[3] = 1 \ a[2] = 1 \ end{bmatrix} \A * B = C ]

    一开始我们只知道(a[1], a[2], a[3]),但是两个矩阵相乘后,我们得到了一个新的值(a[4] = 2),很开心有木有。如果我们用矩阵(A)去乘矩阵(C),会得到一个新的矩阵,暂且叫(C'),你会得到有一个新的值(a[5]),我们有点兴奋起来,这有点想多米诺骨牌,推到第一个,会一直向前倒,知道最后一个。我相信你脑子一定有了这样一个式子:

    […… egin{bmatrix} 1 & 0 & 1 \ 1 & 0 & 0 \ 0 & 1 & 0 \ end{bmatrix}* egin{bmatrix} 1 & 0 & 1 \ 1 & 0 & 0 \ 0 & 1 & 0 \ end{bmatrix}* egin{bmatrix} 1 & 0 & 1 \ 1 & 0 & 0 \ 0 & 1 & 0 \ end{bmatrix}*egin{bmatrix} a[3]=1 \ a[2]=1 \ a[1]=1 \ end{bmatrix} =egin{bmatrix} a[n] \ a[n-1] \ a[n-2] \ end{bmatrix} ]

    矩阵乘法有结合律(但没有交换律,不要问我为什么),所以左边一堆相同的矩阵(不妨叫系数矩阵)可以用一个括号括起来,就像这样:

    [left(…… egin{bmatrix} 1 & 0 & 1 \ 1 & 0 & 0 \ 0 & 1 & 0 \ end{bmatrix}* egin{bmatrix} 1 & 0 & 1 \ 1 & 0 & 0 \ 0 & 1 & 0 \ end{bmatrix}* egin{bmatrix} 1 & 0 & 1 \ 1 & 0 & 0 \ 0 & 1 & 0 \ end{bmatrix} ight)*egin{bmatrix} a[3]=1 \ a[2]=1 \ a[1]=1 \ end{bmatrix} =egin{bmatrix} a[n] \ a[n-1] \ a[n-2] \ end{bmatrix} ]

    想到了什么?

    [egin{bmatrix} 1 & 0 & 1 \ 1 & 0 & 0 \ 0 & 1 & 0 \ end{bmatrix}^k*egin{bmatrix} a[3]=1 \ a[2]=1 \ a[1]=1 \ end{bmatrix} =egin{bmatrix} a[n] \ a[n-1] \ a[n-2] \ end{bmatrix} ]

    我们可以得到(k = n - 3)(想想为什么?),由于n很大,能不能快速求矩阵k次方呢?

    在mod p意义下?矩阵乘法满足结合律?

    快速幂!

    为什么可以用快速幂这个黑科技呢?

    普通的快速幂是用来求(b^k mod p)的,其原理是把(k)二进制拆分成(k=2^{a_1}+2^{a_2}+ ……+2^{a_k}),从而得到(b^k mod p = b^{2^{a_1}} * b^{2^{a_2}} * ……*b^{2^{a_k}} mod p = ((b^{2^{a_1}} mod p) *(b^{2^{a_2}} mod p) * …… * (b^{2^{a_k}} mod p)) mod p) ,只要满足乘法结合律,就可以进行快速幂。

    矩阵快速幂通常用来加速递推。比如说斐波那契数列的第n项mod p的值也可以用矩阵快速幂来求。但是并不是所有的递推都可以用矩阵快速幂,只有那些转移具有周期性的递推才能使用。

    代码模块

    1、矩阵的定义(结构体)

    struct Square{
    	int mat[3][3];
    	void clear(){
    		memset(mat, 0, sizeof(mat));
    	}
    } Base, Result;
    

    2、方阵的乘法

    void Times(Square &A, Square B){
    	Square C; C.clear();
    	for (int i = 0; i <= 2; ++i)
    		for (int j = 0; j <= 2; ++j)
    			for (int k = 0; k <= 2; ++k)
    				(C.mat[i][j] += (LL)A.mat[i][k] * B.mat[k][j] % p) %= p;
    	A = C;
    }
    

    3、矩阵快速幂

    void SquareQpow(Square Base, int k){
    	if (k == 1){
    		Result = Base;
    		return;
    	}
    	Result.clear();
    	SquareQpow(Base, k / 2);
    	Times(Result, Result);
    	if (k % 2 == 1) Times(Result, Base);	
    }
    

    4、矩阵初始化

    void init(){
    	Base.clear();
    	Base.mat[0][0] = 1; Base.mat[0][2] = 1;
    	Base.mat[1][0] = 1; Base.mat[2][1] = 1;
    }
    

    易错点

    1. 乘法时需进行强制类型转换:(C.mat[i][j] += (LL)A.mat[i][k] * B.mat[k][j] % p) %= p;
    2. C++数组从0开始的哦QAQ
    3. 计算答案时注意要加3项,不要只加2项

    完整代码

    #include<iostream>
    #include<cstdio>
    #include<cstring>
    #define LL long long
    using namespace std;
    const int p = 1e9 + 7;
    struct Square{
    	int mat[3][3];
    	void clear(){
    		memset(mat, 0, sizeof(mat));
    	}
    } Base, Result;
    
    void init(){
    	Base.clear();
    	Base.mat[0][0] = 1; Base.mat[0][2] = 1;
    	Base.mat[1][0] = 1; Base.mat[2][1] = 1;
    }
    
    void Times(Square &A, Square B){
    	Square C; C.clear();
    	for (int i = 0; i <= 2; ++i)
    		for (int j = 0; j <= 2; ++j)
    			for (int k = 0; k <= 2; ++k)
    				(C.mat[i][j] += (LL)A.mat[i][k] * B.mat[k][j] % p) %= p;
    	A = C;
    }
    
    void SquareQpow(Square Base, int k){
    	if (k == 1){
    		Result = Base;
    		return;
    	}
    	Result.clear();
    	SquareQpow(Base, k / 2);
    	Times(Result, Result);
    	if (k % 2 == 1) Times(Result, Base);	
    }
    
    int main(){
    	int T; scanf("%d", &T);
    	while (T--){
    		int n;
    		scanf("%d", &n);
    		if (n <= 3) printf("1
    ");
    		else{
    			init();
    			SquareQpow(Base, n-3);
    			printf("%d
    ", ((Result.mat[0][0] + Result.mat[0][1]) % p + Result.mat[0][2]) % p);
    		}
    	}
    	return 0;
    }
    
  • 相关阅读:
    windows设置文件夹显示缩略图
    windows合并文件夹窗口
    crm高速开发之Entity
    在Eclipse中搭建Dagger和Dagger2使用环境
    CCEditBox/CCEditBoxImpl
    failed to sync branch You might need to open a shell and debug the state of this repo.
    五个方法:让站点用心服务于每一位客户
    Mobile first! Wijmo 5 + Ionic Framework之:Hello World!
    ST Nucleo mbed套件开发 一 MBED环境使用 以Nucleo-F401为例 (二)
    关于Windows通过远程桌面訪问Ubuntu
  • 原文地址:https://www.cnblogs.com/YJZoier/p/9442086.html
Copyright © 2020-2023  润新知