$\sum_{i = 1}^k A^i = \sum_{i = 1}^{k / 2 + 1} A^i + A^{k / 2 + 1}(\sum_{i = 1}^{k / 2} A^i)$
- #include<cstdio>
- #include<cstring>
- #include<iostream>
- #include<map>
- #define LL long long
- using namespace std;
- int N, K, mod;
- int mul(int x, int y) {
- if(1ll * x * y > mod) return 1ll * x * y % mod;
- else return 1ll * x * y;
- }
- int add(int x, int y) {
- if(x + y > mod) return x + y - mod;
- else return x + y;
- }
- struct Matrix {
- int m[31][31];
- Matrix() {
- memset(m, 0, sizeof(m));
- }
- bool operator < (const Matrix &rhs) const {
- for(int i = 1; i <= N; i++)
- for(int j = 1; j <= N; j++)
- if(m[i][j] != rhs.m[i][j])
- return m[i][j] < rhs.m[i][j];
- return 1;
- }
- Matrix operator * (const Matrix &rhs) const {
- Matrix ans;
- for(int k = 1; k <= N; k++)
- for(int i = 1; i <= N; i++)
- for(int j = 1; j <= N; j++)
- ans.m[i][j] = add(ans.m[i][j], mul(m[i][k], rhs.m[k][j]));
- return ans;
- }
- Matrix operator + (const Matrix &rhs) const {
- Matrix ans;
- for(int i = 1; i <= N; i++)
- for(int j = 1; j <= N; j++)
- ans.m[i][j] = add(m[i][j], rhs.m[i][j]);
- return ans;
- }
- }a;
- Matrix getbase() {
- Matrix base;
- for(int i = 1; i <= N; i++) base.m[i][i] = 1;
- return base;
- }
- Matrix fp(Matrix a, int p) {
- Matrix base = getbase();
- while(p) {
- if(p & 1) base = base * a;
- a = a * a; p >>= 1;
- }
- return base;
- }
- Matrix solve(int k) {
- if(k == 1) return a;
- Matrix res = solve(k / 2);
- if(k & 1) {
- Matrix po = fp(a, k / 2 + 1);
- return res + po + po * res;
- }
- else return res + fp(a, k / 2) * res;
- }
- main() {
- // freopen("a.in", "r", stdin);
- cin >> N >> K >> mod;
- for(int i = 1; i <= N; i++)
- for(int j = 1; j <= N; j++)
- cin >> a.m[i][j];
- Matrix ans = solve(K);
- for(int i = 1; i <= N; i++, puts(""))
- for(int j = 1; j <= N; j++)
- printf("%d ", ans.m[i][j] % mod);
- }