矩阵快速幂各类题型总结（一般，共轭，1 * n，矩阵简化）

Reading comprehension

这个不难找出递推式f[n]=f[n−1]+2f[n−2]+1f[n] = f[n - 1] + 2f[n - 2] + 1f[n]=f[n−1]+2f[n−2]+1。

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f;
const double eps = 1e-7;ll n, mod;struct matrix {ll a[3][3];
};matrix operator * (matrix a, matrix b) {matrix ans;for(int i = 0; i < 3; i++) {for(int j = 0; j < 3; j++) {ans.a[i][j] = 0;for(int k = 0; k < 3; k++) {ans.a[i][j] = (ans.a[i][j] + 1ll * a.a[i][k] * b.a[k][j] % mod) % mod;}}}return ans;
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);while(scanf("%lld %lld", &n, &mod) != EOF && (n + mod)) {if(n <= 2) {if(n == 1) printf("%d\n", 1 % mod);else printf("%d\n", 2 % mod);continue;}matrix ans = {2, 1, 1,0, 0, 0,0, 0, 0};matrix a = {1, 1, 0,2, 0, 0, 1, 0, 1};n -= 2;while(n) {if(n & 1) ans = ans * a;a = a * a;n >>= 1;}printf("%lld\n", ans.a[0][0]);}return 0;
}

这里直接写出构造的矩阵了。
A=[a00a10a20…a(n−1)0an010…0…]A =\left[\begin{matrix} a_{00}&a_{10}&a_{20}\dots a_{(n - 1)0}&a_{n0}&1\\ 0& \dots\\ 0& \dots\\ \end{matrix}\right] A=⎣⎡a0000a10……a20…a(n−1)0an01⎦⎤
整体来说就是一个(n+2)×(n+2)(n + 2) \times(n + 2)(n+2)×(n+2)的矩阵吧。
B=[101010…10100011…110001…110…333…331]B =\left[ \begin{matrix} 10&10&10 \dots 10&10&0\\ 0&1&1 \dots 1&1&0\\ 0&0&1\dots1&1&0\\ \dots\\ 3&3&3\dots3&3&1\\ \end{matrix} \right] B=⎣⎢⎢⎢⎢⎡1000…31010310…101…11…13…3101130001⎦⎥⎥⎥⎥⎤
同样的也是一个(n+2)×(n+2)(n + 2) \times(n + 2)(n+2)×(n+2)的矩阵。

我们个a00a_{00}a00初始化为232323，这个时候A×BA \times BA×B发现AAA中的元素会变成第二列的，所以这里就形成了一个递推式了，只要跑一跑矩阵快速幂即可。

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f;
const double eps = 1e-7;const int mod = 10000007;int n, m;struct matrix {ll a[15][15];void init() {memset(a, 0, sizeof a);}
};matrix operator * (matrix a, matrix b) {matrix ans;for(int i = 0; i < n; i++) {for(int j = 0; j < n; j++) {ans.a[i][j] = 0;for(int k = 0; k < n; k++) {ans.a[i][j] = (ans.a[i][j] + 1ll * a.a[i][k] * b.a[k][j] % mod) % mod;}}}return ans;
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);while(scanf("%d %d", &n, &m) != EOF) {matrix ans, a;ans.init(), a.init();for(int i = 1; i <= n; i++) scanf("%lld", &ans.a[0][i]);ans.a[0][0] = 23, ans.a[0][n + 1] = 1;for(int i = 0; i <= n; i++) {a.a[0][i] = 10;a.a[n + 1][i] = 3;}a.a[n + 1][n + 1] = 1;for(int i = 1; i <= n; i++) {for(int j = i; j <= n; j++) {a.a[i][j] = 1;}}n += 2;while(m) {if(m & 1) ans = ans * a;a = a * a;m >>= 1;}printf("%lld\n", ans.a[0][n - 2]);}return 0;
}

B. Jzzhu and Sequences

裸题fn=fn−1−fn−2f_n = f_{n - 1} - f_{n - 2}fn=fn−1−fn−2。

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f;
const double eps = 1e-7;const int mod = 1e9 + 7;int n = 2, m;struct matrix {ll a[2][2];void init() {memset(a, 0, sizeof a);}
};matrix operator * (matrix a, matrix b) {matrix ans;for(int i = 0; i < n; i++) {for(int j = 0; j < n; j++) {ans.a[i][j] = 0;for(int k = 0; k < n; k++) {ans.a[i][j] = (ans.a[i][j] + 1ll * a.a[i][k] * b.a[k][j] % mod) % mod;}}}return ans;
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);ll x, y, n;scanf("%lld %lld %lld", &x, &y, &n);if(n <= 2) {if(n == 1) printf("%d\n", (x % mod + mod) % mod);else printf("%d\n", (y % mod + mod) % mod);return 0;}matrix ans = {y, x, 0, 0};matrix a = {1, 1,-1, 0};n -= 2;while(n) {if(n & 1) ans = ans * a;a = a * a;n >>= 1;}printf("%lld", (ans.a[0][0] % mod + mod) % mod);return 0;
}

M斐波那契数列

这是个乘法的递推式，显然乘法可以变成指数相加，所以还是一个矩阵快速幂，统计最后一项有多少个A0,A1A_0,A_1A0,A1即可。

当然矩阵中的取模得用费马小定理。

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f;
const double eps = 1e-7;const int mod = 1000000006;int n = 4, m;struct matrix {ll a[4][4];void init() {memset(a, 0, sizeof a);}
};matrix operator * (matrix a, matrix b) {matrix ans;for(int i = 0; i < n; i++) {for(int j = 0; j < n; j++) {ans.a[i][j] = 0;for(int k = 0; k < n; k++) {ans.a[i][j] = (ans.a[i][j] + 1ll * a.a[i][k] * b.a[k][j] % mod) % mod;}}}return ans;
}ll quick_pow(ll a, int n) {ll ans = 1;const int mod = 1e9 + 7;while(n) {if(n & 1) ans = ans * a % mod;a = a * a % mod;n >>= 1;}return ans;
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);ll x, y, n;while(scanf("%lld %lld %lld", &x, &y, &n) != EOF) {const int mod = 1e9 + 7;if(n <= 1) {if(n == 0) {printf("%lld\n", x % mod);}else {printf("%lld\n", y % mod);}continue;}matrix ans = {1, 0, 0, 1,0, 0, 0, 0,0, 0, 0, 0,0, 0, 0, 0};matrix a = {1, 0, 1, 0,0, 1, 0, 1,1, 0, 0, 0,0, 1, 0, 0};n -= 1;while(n) {if(n & 1) ans = ans * a;a = a * a;n >>= 1;}ll res = 1ll * quick_pow(x, ans.a[0][1]) * quick_pow(y, ans.a[0][0]) % mod;printf("%lld\n", res);}return 0;
}

So Easy!

共轭矩阵构造的经典题了。

假设An=(a+b)n,Bn=(a−b)nA_n = (a + \sqrt b) ^n,B_n = (a - \sqrt b) ^nAn=(a+b)n,Bn=(a−b)n，一定有An+BnA_n + B_nAn+Bn是个有理数

并且有Bn<1B_n < 1Bn<1，所以有⌈Sn⌉=An+Bn\lceil S_n \rceil = A_n + B_n⌈Sn⌉=An+Bn

2aSn=((a+b)n+(a−b)n)((a+b)+(a−b))=Sn+1+(a2−b)Sn−12aS_n = ((a + \sqrt b) ^n + (a - \sqrt b) ^n)((a + \sqrt b) + (a - \sqrt b)) = S_{n + 1} +(a ^2 - b)S_{n - 1}2aSn=((a+b)n+(a−b)n)((a+b)+(a−b))=Sn+1+(a2−b)Sn−1

得到Sn=2aSn−1+(b−a2)Sn−2S_n = 2aS_{n - 1} +(b - a ^ 2)S_{n - 2}Sn=2aSn−1+(b−a2)Sn−2，于是递推式就得到了进行矩阵快速幂即可。

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f;
const double eps = 1e-7;// const int mod = 1000000006;ll n = 2, mod;struct matrix {ll a[2][2];void init() {memset(a, 0, sizeof a);}
};matrix operator * (matrix a, matrix b) {matrix ans;for(int i = 0; i < n; i++) {for(int j = 0; j < n; j++) {ans.a[i][j] = 0;for(int k = 0; k < n; k++) {ans.a[i][j] = (ans.a[i][j] + 1ll * a.a[i][k] * b.a[k][j] % mod) % mod;}}}return ans;
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);ll a, b, n;while(scanf("%lld %lld %lld %lld", &a, &b, &n, &mod) != EOF) {ll a1 = 2ll * a % mod, a2 = 2ll * (a * a % mod + b) % mod;if(n <= 2) {if(n == 1) printf("%lld\n", a1);else printf("%lld\n", a2);continue;}n -= 2;matrix ans = {a2, a1, 0, 0};matrix fat = {2ll * a % mod, 1,((b - a * a % mod) % mod + mod) % mod, 0};while(n) {if(n & 1) ans = ans * fat;fat = fat * fat;n >>= 1;}printf("%lld\n", ans.a[0][0]);}return 0;
}

Problem of Precision

这道题目是向下取整，跟上一道题目一样构造，唯一的不同就是Sn=An+Bn−1S_n = A_n + B_n - 1Sn=An+Bn−1，所以我们还是按照SnS_nSn递推，最后在答案上减一即可。

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f;
const double eps = 1e-7;const int mod = 1024;int n = 2;struct matrix {ll a[2][2];void init() {memset(a, 0, sizeof a);}
};matrix operator * (matrix a, matrix b) {matrix ans;for(int i = 0; i < n; i++) {for(int j = 0; j < n; j++) {ans.a[i][j] = 0;for(int k = 0; k < n; k++) {ans.a[i][j] = (ans.a[i][j] + 1ll * a.a[i][k] * b.a[k][j] % mod) % mod;}}}return ans;
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int T;scanf("%d", &T);while(T--) {ll n;scanf("%lld", &n);ll a1 = 10, a2 = 98;if(n <= 2) {if(n == 1) printf("%lld\n", (a1 - 1) % mod);else printf("%lld\n", (a2 - 1) % mod);continue;}n -= 2;matrix ans = {a2, a1,0, 0};matrix fat = {10, 1,1023, 0};while(n) {if(n & 1) ans = ans * fat;fat = fat * fat;n >>= 1;}printf("%lld\n", (ans.a[0][0] - 1 + mod) % mod);}return 0;
}

C. Partial Sums

由于递推矩阵的特殊关系，所以可以转化为1×n1 \times n1×n的矩阵递推，最后达到n2log(n)n ^ 2log(n)n2log(n)的复杂度。

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f;
const double eps = 1e-7;const int N = 2e3 + 10, mod = 1e9 + 7;int n, k;struct matrix {ll a[N];
};matrix operator * (matrix a, matrix b) {matrix ans;for(int i = 1; i <= n; i++) {ans.a[i] = 0;for(int j = 1; j <= i; j++) {ans.a[i] = (ans.a[i] + a.a[j] * b.a[i - j + 1] % mod) % mod;}}return ans;
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);scanf("%d %d", &n, &k);matrix ans, a;for(int i = 1; i <= n; i++) {scanf("%d", &ans.a[i]);a.a[i] = 1;}while(k) {if(k & 1) ans = ans * a;a = a * a;k >>= 1;}for(int i = 1; i <= n; i++) {printf("%lld%c", ans.a[i], i == n ? '\n' : ' ');}return 0;
}

Fast Matrix Calculation

ABn2=A(BA)n2−1BAB^{n^2} = A(BA)^{n^2 - 1}BABn2=A(BA)n2−1B,转化为k×kk \times kk×k的矩阵快速幂，这样就可以保证不超时了。

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f;
const double eps = 1e-7;const int N = 1e3 + 10, mod = 6;int n, m;struct matrix {ll a[6][6];void init() {memset(a, 0, sizeof a);}
}c;struct Matrix {ll a[N][N];
}a, b, Ans1, Ans2;matrix operator * (matrix a, matrix b) {matrix ans;for(int i = 0; i < m; i++) {for(int j = 0; j < m; j++) {ans.a[i][j] = 0;for(int k = 0; k < m; k++) {ans.a[i][j] = (ans.a[i][j] + a.a[i][k] * b.a[k][j] % mod) % mod;}}}return ans;
}matrix mult() {matrix ans;for(int i = 0; i < m; i++) {for(int j = 0; j < m; j++) {ans.a[i][j] = 0;for(int k = 0; k < n; k++) {ans.a[i][j] = (ans.a[i][j] + b.a[i][k] * a.a[k][j] % mod) % mod;}}}return ans;
}int main() {// freopen("in.txt", "r", stdin);   // freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);while(scanf("%d %d", &n, &m) && (n + m)) {for(int i = 0; i < n; i++) {for(int j = 0; j < m; j++) {scanf("%lld", &a.a[i][j]);}}for(int i = 0; i < m; i++) {for(int j = 0; j < n; j++) {scanf("%lld", &b.a[i][j]);}}c = mult();matrix ans;ans.init();for(int i = 0; i < m; i++) ans.a[i][i] = 1;int n1 = n * n - 1;while(n1) {if(n1 & 1) ans = ans * c;c = c * c;n1 >>= 1;}for(int i = 0; i < n; i++) {for(int j = 0; j < m; j++) {Ans1.a[i][j] = 0;for(int k = 0; k < m; k++) {Ans1.a[i][j] = (Ans1.a[i][j] + a.a[i][k] * ans.a[k][j] % mod) % mod;}}}for(int i = 0; i < n; i++) {for(int j = 0; j < n; j++) {Ans2.a[i][j] = 0;for(int k = 0; k < m; k++) {Ans2.a[i][j] = (Ans2.a[i][j] + Ans1.a[i][k] * b.a[k][j] % mod) % mod;}}}int res = 0;for(int i = 0; i < n; i++) {for(int j = 0; j < n; j++) {res += Ans2.a[i][j];}}printf("%d\n", res);}return 0;
}