矩阵快速幂の模板

矩阵快速幂：

当某些转移关系可以用矩阵乘法表示的时候，就可以用矩阵快速幂优化时间复杂度 ， 可以把线性的递推用矩阵快速幂优化成 log

模板：

typedef vector<vector<ll>> M;
// y 相当于 res ， 最后存储答案
// x 相当于底数 a

M cal(M a,M b){
	int l1 = a.size() , l2 = a[0].size() , l3 = b[0].size();
	M c(l1,vector<ll>(l3,0));
	for(int i=0;i<l1;i++){
		for(int j=0;j<l3;j++){
			for(int k=0;k<l2;k++){
				c[i][j] = (c[i][j] + a[i][k] * b[k][j] % p) % p;
			}
		}
	}
	return c;
}

void qp(ll b){
	while(b){
		if(b % 2 != 0) y = cal(y , x);
		x = cal(x , x);
		b /= 2;
	}	
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

例一：

又见斐波那契

状态转移式构造矩阵乘法：

$\left[\begin{array}{cccccc}{F}_{n-2}& F_{n-1}& n^3& n^2& n& 1\end{array}\right]\ast \left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 1& 1& 0& 0& 0& 0\\ 0& 1& 1& 0& 0& 0\\ 0& 1& 3& 1& 0& 0\\ 0& 1& 3& 2& 1& 0\\ 0& 1& 1& 1& 1& 1\end{array}\right]=\left[\begin{array}{cccccc}{F}_{n-1}& F_n& (n+1{)}^3& (n+1{)}^2& (n+1)& 1\end{array}\right]$

1. 首先根据递推式写出先后的矩阵项，由此可以确定出需要构造的矩阵大小。前后的矩阵项中一定要包含递推式中的所有项
2. 然后构造出矩阵，进行矩阵快速幂
3. 矩阵快速幂的难点在于构造快速幂矩阵

#include
using namespace std;
#define fi first
#define se second
#define IOS std::ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
typedef long long ll;
const int N = 2e6+10;
const int p = 1e9 + 7;
typedef pair<int,int>PII;
const int inf = 1e9 + 10;
const double eps = 1e-9;
typedef vector<vector<ll>> M;
M x(6,vector<ll>(6,0));
M y(1,vector<ll>(6,0));


void init(){
	x[0][0] = 0;x[0][1] = 1;x[0][2] = 0;x[0][3] = 0;x[0][4] = 0;x[0][5] = 0; 
	x[1][0] = 1;x[1][1] = 1;x[1][2] = 0;x[1][3] = 0;x[1][4] = 0;x[1][5] = 0; 
	x[2][0] = 0;x[2][1] = 1;x[2][2] = 1;x[2][3] = 0;x[2][4] = 0;x[2][5] = 0; 
	x[3][0] = 0;x[3][1] = 1;x[3][2] = 3;x[3][3] = 1;x[3][4] = 0;x[3][5] = 0; 
	x[4][0] = 0;x[4][1] = 1;x[4][2] = 3;x[4][3] = 2;x[4][4] = 1;x[4][5] = 0; 
	x[5][0] = 0;x[5][1] = 1;x[5][2] = 1;x[5][3] = 1;x[5][4] = 1;x[5][5] = 1; 
	
	y[0][0] = 0;y[0][1] = 1;y[0][2] = 8;y[0][3] = 4;y[0][4] = 2;y[0][5] = 1; 
}

// y 相当于 res ， 最后存储答案
// x 相当于底数 a

M cal(M a,M b){
	int l1 = a.size() , l2 = a[0].size() , l3 = b[0].size();
	M c(l1,vector<ll>(l3,0));
	for(int i=0;i<l1;i++){
		for(int j=0;j<l3;j++){
			for(int k=0;k<l2;k++){
				c[i][j] = (c[i][j] + a[i][k] * b[k][j] % p) % p;
			}
		}
	}
	return c;
}

void qp(ll b){
	while(b){
		if(b % 2 != 0) y = cal(y , x);
		x = cal(x , x);
		b /= 2;
	}	
}

ll t , n;

int main(){
//	IOS

	cin >> t;
	
	while(t--){
		
		cin >> n;
		
		init();
		qp(n);
		
		cout << y[0][0] << "\n" ;
		
	}

	return 0;
}
//freopen("文件名.in","r",stdin);
//freopen("文件名.out","w",stdout);

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74