POJ 3233 Matrix Power Series(矩阵快速幂)
题意:
求 S = A + A2 + A3 + … + Ak
思路:
和求矩阵幂方法类似,采取二分的思想:
A+A^2+A...+A^(2k+1)= A+A^2+...+A^k + A^(k+1) + A^(k+1)*(A+A^2+...+A^k).
A+A^2+...+A^2k = A+A^2+...+A^k + A^k*(A+A^2+...+A^k).
#include <cstdio>
#include <cstring>
#include <cstring>
#include <iostream>
using namespace std;
const int MAXN = 32;
struct Matrix {
int v[MAXN][MAXN];
} ;
int n, M;
Matrix matrixadd(Matrix a, Matrix b)
{
Matrix c;
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
c.v[i][j] = (a.v[i][j] + b.v[i][j]) % M;
return c;
}
Matrix matrixmul(Matrix a, Matrix b)
{
Matrix c;
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j) {
c.v[i][j] = 0;
for (int k = 0; k < n; ++k)
c.v[i][j] = (a.v[i][k] * b.v[k][j] + c.v[i][j]) % M;
}
return c;
}
Matrix matrixpow(Matrix m, int k)
{
if (k == 1)
return m;
Matrix c;
c = matrixpow(m, k >> 1);
c = matrixmul(c, c);
if (k % 2)
c = matrixmul(c, m);
return c;
}
Matrix matrixcalc(Matrix m, int k)
{
if (k == 1)
return m;
Matrix a, b, c;
a = matrixcalc(m, k >> 1);
if (k % 2) {
b = matrixpow(m, (k+1)>>1);
c = matrixadd(matrixadd(a, b), matrixmul(a, b));
} else {
b = matrixpow(m, k >> 1);
c = matrixadd(a, matrixmul(a, b));
}
return c;
}
int main()
{
Matrix m;
int k;
scanf("%d %d %d", &n, &k, &M);
memset(m.v, 0, sizeof(m.v));
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
scanf("%d", &m.v[i][j]);
m = matrixcalc(m, k);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j)
printf("%d ", m.v[i][j]);
printf("\n");
}
return 0;
}
-------------------------------------------------------
kedebug
Department of Computer Science and Engineering,
Shanghai Jiao Tong University
E-mail: kedebug0@gmail.com
GitHub: http://github.com/kedebug
-------------------------------------------------------