分析：矩阵相乘O(n^3), 有k次，则复杂度为n^3*k。
使用矩阵技巧，构造：
B= | A A |
| O I |
则B的乘方的结果其右上会是S，其他三个不变。此时化成了矩阵乘方问题，此时可以使用反复平方法，这样复杂度为(2n)^3*logk
#include<stdio.h>
#define nnum 61
int res[nnum][nnum], te[nnum][nnum], n, m;
void mulmatrix(int a[nnum][nnum], int b[nnum][nnum], int mark) {
int i, j, k, temp[nnum][nnum];
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
temp[i][j] = 0;
for (k = 0; k < n; k++) {
temp[i][j] += a[i][k] * b[k][j];
temp[i][j] %= m;
}
}
}
if (mark) {
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
te[i][j] = temp[i][j];
}
}
} else {
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
res[i][j] = temp[i][j];
}
}
}
}
void cal(int times) {
while (times) {
if (times & 1) {
mulmatrix(te, res, 1);
}
mulmatrix(res, res, 0);
times >>= 1;
}
}
int main() {
#ifndef ONLINE_JUDGE
freopen("t.txt", "r", stdin);
#endif
int k, i, j, re;
while (scanf("%d %d %d", &n, &k, &m) != EOF) {
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
scanf("%d", &res[i][j]);
res[i][j + n] = res[i][j];
}
}
for (i = n; i < n + n; i++) {
for (j = n; j < n; j++) {
res[i][j] = 0;
}
res[i][i] = 1;
}
re = n;
n = n + n;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
te[i][j] = 0;
}
te[i][i] = 1;
}
cal(k);
for (i = 0; i < re; i++) {
for (j = re; j < n; j++) {
printf("%d ", te[i][j]);
}
printf("\n");
}
}
return 0;
}
/*
* hdu1757.c
*
* Created on: 2011-9-12
* Author: bjfuwangzhu
*/
//
//#include<stdio.h>
//#define nmax 10
//typedef struct matrix {
// int num[nmax][nmax];
//} matrix;
//int m;
//matrix res, temp;
//matrix matrix_mul(matrix numa, matrix numb) {
// matrix numc;
// int i, j, k;
// for (i = 0; i < nmax; i++) {
// for (j = 0; j < nmax; j++) {
// numc.num[i][j] = 0;
// for (k = 0; k < nmax; k++) {
// numc.num[i][j] += numa.num[i][k] * numb.num[k][j];
// numc.num[i][j] %= m;
// }
// }
// }
// return numc;
//}
//void solve(int k) {
// while (k) {
// if (k & 1) {
// res = matrix_mul(res, temp);
// }
// temp = matrix_mul(temp, temp);
// k >>= 1;
// }
//}
//int main() {
//#ifndef ONLINE_JUDGE
// freopen("data.in", "r", stdin);
//#endif
// int k, i, j, ans;
// while (scanf("%d %d", &k, &m) != EOF) {
// for (i = 0; i < nmax; i++) {
// for (j = 0; j < nmax; j++) {
// temp.num[i][j] = 0;
// res.num[i][j] = 0;
// }
// res.num[i][i] = 1;
// }
// for (i = 0; i < nmax; i++) {
// scanf("%d", &temp.num[0][i]);
// }
// for (i = 0; i < nmax - 1; i++) {
// temp.num[i + 1][i] = 1;
// }
// if (k < nmax) {
// printf("%d\n", k % m);
// continue;
// }
// solve(k - 9);
// for (i = 0, ans = 0; i < nmax; i++) {
// ans += res.num[0][i] * (nmax - 1 - i);
// ans %= m;
// }
// printf("%d\n", ans);
// }
// return 0;
//}
#include<stdio.h>
#define nmax 10
int matrix_res[nmax][nmax], matrix_temp[nmax][nmax];
void matrix_mul(int a[nmax][nmax], int b[nmax][nmax], int m, int mark) {
int temp[nmax][nmax];
int i, j, k;
for (i = 0; i < nmax; i++) {
for (j = 0; j < nmax; j++) {
temp[i][j] = 0;
for (k = 0; k < nmax; k++) {
temp[i][j] += a[i][k] * b[k][j];
temp[i][j] %= m;
}
}
}
if (mark) {
for (i = 0; i < nmax; i++) {
for (j = 0; j < nmax; j++) {
matrix_res[i][j] = temp[i][j];
}
}
} else {
for (i = 0; i < nmax; i++) {
for (j = 0; j < nmax; j++) {
matrix_temp[i][j] = temp[i][j];
}
}
}
}
void solve(int k, int m) {
if (k < 0) {
return;
}
while (k) {
if (k & 1) {
matrix_mul(matrix_res, matrix_temp, m, 1);
}
matrix_mul(matrix_temp, matrix_temp, m, 0);
k >>= 1;
}
}
int main() {
#ifndef ONLINE_JUDGE
freopen("data.in", "r", stdin);
#endif
int k, m, i, j, res;
while (scanf("%d %d", &k, &m) != EOF) {
for (i = 0; i < nmax; i++) {
for (j = 0; j < nmax; j++) {
matrix_temp[i][j] = 0;
matrix_res[i][j] = 0;
}
matrix_res[i][i] = 1;
}
for (i = 0; i < nmax; i++) {
scanf("%d", &matrix_temp[0][i]);
}
for (i = 0; i < nmax - 1; i++) {
matrix_temp[i + 1][i] = 1;
}
if (k < nmax) {
printf("%d\n", k % m);
continue;
}
solve(k - 9, m);
for (i = 0, res = 0; i < nmax; i++) {
res = res + matrix_res[0][i] * (nmax - i - 1);
res %= m;
}
printf("%d\n", res);
}
return 0;
}
/*
* fzu1683.c
*
* Created on: 2011-9-12
* Author: bjfuwangzhu
*/
#include<stdio.h>
#define nnum 2009
#define nmax 4
typedef struct matrix {
int num[nmax][nmax];
} matrix;
matrix res, temp;
int sum[nmax];
matrix matrix_mul(matrix numa, matrix numb) {
matrix numc;
int i, j, k;
for (i = 0; i < nmax; i++) {
for (j = 0; j < nmax; j++) {
numc.num[i][j] = 0;
for (k = 0; k < nmax; k++) {
numc.num[i][j] += numa.num[i][k] * numb.num[k][j];
numc.num[i][j] %= nnum;
}
}
}
return numc;
}
void solve(int k) {
while (k) {
if (k & 1) {
res = matrix_mul(res, temp);
}
temp = matrix_mul(temp, temp);
k >>= 1;
}
}
int main() {
#ifndef ONLINE_JUDGE
freopen("data.in", "r", stdin);
#endif
int m, n, i, j, k, ant;
while (scanf("%d", &m) != EOF) {
for (i = 1; i <= m; i++) {
scanf("%d", &n);
for (j = 0; j < nmax; j++) {
for (k = 0; k < nmax; k++) {
temp.num[j][k] = 0, res.num[j][k] = 0;
}
res.num[j][j] = 1;
}
temp.num[0][0] = 1, temp.num[0][1] = 3, temp.num[0][2] = 2, temp.num[0][3] =
7;
temp.num[1][1] = 3, temp.num[1][2] = 2, temp.num[1][3] = 7;
temp.num[2][1] = 1, temp.num[3][2] = 1;
sum[3] = 1, sum[2] = 3, sum[1] = 5, sum[0] = 9;
printf("Case %d: ", i);
if (n < nmax - 1) {
for (j = 0, ant = 0; j <= n; j++) {
ant += sum[nmax - 1 - j];
ant %= nnum;
}
printf("%d\n", ant);
continue;
}
solve(n - 2);
for (j = 0, ant = 0; j < nmax; j++) {
ant += res.num[0][j] * sum[j];
ant %= nnum;
}
printf("%d\n", ant);
}
}
return 0;
}
/*
C语言版
*/
#include<stdio.h>
#define nmax 4
#define nnum 10007
typedef struct matrix {
int num[nmax][nmax];
} matrix;
matrix matrix_mul(matrix numa, matrix numb) {
matrix numc;
int i, j, k;
for (i = 0; i < nmax; i++) {
for (j = 0; j < nmax; j++) {
numc.num[i][j] = 0;
for (k = 0; k < nmax; k++) {
numc.num[i][j] = (numc.num[i][j]
+ numa.num[i][k] * numb.num[k][j]) % nnum;
}
}
}
return numc;
}
matrix matrix_exp(matrix numc, matrix numa, int x) {
for (; x; x >>= 1) {
if (x & 1) {
numc = matrix_mul(numc, numa);
}
numa = matrix_mul(numa, numa);
}
return numc;
}
int main() {
#ifndef ONLINE_JUDGE
freopen("data.in", "r", stdin);
#endif
int i, j, N, X, Y, sum[nmax] = { 2, 1, 1, 1 }, ant;
matrix res, temp;
while (scanf("%d %d %d", &N, &X, &Y) != EOF) {
if (N == 0) {
printf("1\n");
continue;
}
if (N == 1) {
printf("2\n");
continue;
}
for (i = 0; i < nmax; i++) {
for (j = 0; j < nmax; j++) {
res.num[i][j] = 0, temp.num[i][j] = 0;
}
res.num[i][i] = 1;
}
X %= nnum,Y %= nnum;
temp.num[0][0] = 1, temp.num[0][1] = X * X % nnum, temp.num[0][2] = Y
* Y % nnum, temp.num[0][3] = 2 * X * Y % nnum;
temp.num[1][1] = X * X % nnum, temp.num[1][2] = Y * Y % nnum, temp.num[1][3] =
2 * X * Y % nnum;
temp.num[2][1] = 1, temp.num[3][1] = X, temp.num[3][3] = Y;
res = matrix_exp(res, temp, (N - 1));
for (i = 0, ant = 0; i < nmax; i++) {
ant += res.num[0][i] * sum[i];
ant %= nnum;
}
printf("%d\n", ant);
}
return 0;
}
/*
* hdu3306.cpp
*
* Created on: 2011-9-13
* Author: bjfuwangzhu
*/
/*
C++语言版 重载运算符
*/
#include<stdio.h>
#define nmax 4
#define nnum 10007
struct matrix {
int num[nmax][nmax];
matrix() {
for (int i = 0; i < nmax; i++) {
for (int j = 0; j < nmax; j++) {
num[i][j] = 0;
}
}
}
};
matrix operator*(matrix numa, matrix numb) {
matrix numc;
int i, j, k;
for (i = 0; i < nmax; i++) {
for (j = 0; j < nmax; j++) {
for (k = 0; k < nmax; k++) {
numc.num[i][j] = (numc.num[i][j]
+ numa.num[i][k] * numb.num[k][j]) % nnum;
}
}
}
return numc;
}
matrix operator^(matrix numa, int x) {
matrix numc;
int i;
for (i = 0; i < nmax; i++) {
numc.num[i][i] = 1;
}
for (; x; x >>= 1) {
if (x & 1) {
numc = numc * numa;
}
numa = numa * numa;
}
return numc;
}
int main() {
#ifndef ONLINE_JUDGE
freopen("data.in", "r", stdin);
#endif
int i, N;
int X, Y, sum[nmax] = { 2, 1, 1, 1 }, ant;
matrix res, temp;
while (scanf("%d %d %d", &N, &X, &Y) != EOF) {
if (N == 0) {
printf("1\n");
continue;
}
if (N == 1) {
printf("2\n");
continue;
}
X %= nnum,Y %= nnum;
temp.num[0][0] = 1, temp.num[0][1] = X * X % nnum, temp.num[0][2] = Y
* Y % nnum, temp.num[0][3] = 2 * X * Y % nnum;
temp.num[1][1] = X * X % nnum, temp.num[1][2] = Y * Y % nnum, temp.num[1][3] =
2 * X * Y % nnum;
temp.num[2][1] = 1, temp.num[3][1] = X, temp.num[3][3] = Y;
res = temp ^ (N - 1);
for (i = 0, ant = 0; i < nmax; i++) {
ant += res.num[0][i] * sum[i];
ant %= nnum;
}
printf("%d\n", ant);
}
return 0;
}