台湾菜

Matrix Power Series

twcai • geek stuff

经小HH大牛指点,重新过了一遍POJ 3233,效率有很大的提高
题目给出一个矩阵A, 要求求出A^1 + A^2 + ... + A^k
这里的做法是构造一个两倍于A的矩阵B[math]matrix{2}{2}{e A 0 A}[/math]
然后求B^K次即可,所要的矩阵位于B的右上角
#include
using namespace std;
//A^1 + A^2 + ... + A^k
int modnum;

class matrix{
public :
int rw, cn, ma[65][65];

matrix() {
rw = 0, cn = 0;
}

matrix(int x, int y) {
int i, j;
rw = x;
cn = y;
for(i = 0; i < rw; i++) {
for(j = 0; j < cn; j++) { ma[i][j] = 0; }
}
}

matrix(const matrix &x) {
int i, j;
rw = x.rw, cn = x.cn;
for(i = 0; i < rw; i++){
for(j = 0; j < cn; j++) { ma[i][j] = x.ma[i][j]; }
}
}

matrix operator*(matrix &x)
{
int i, j, k;
matrix ret(cn, x.rw);
for(i = 0; i < rw ; i++) {
for(j = 0; j < x.cn ; j++) {
for(k = 0; k < cn; k++) {
ret.ma[i][j] += ma[i][k] * x.ma[k][j];
if(ret.ma[i][j] >= modnum)
ret.ma[i][j] %= modnum;
}
}
}
return ret;
}

matrix power(int x) {
matrix ret(rw, cn);
if(x == 1) return *this;
if(x == 2) {
int i, j ,k;
for(i = 0; i < rw; i++) {
for(j = 0; j < cn; j++) {
for(k = 0; k < cn; k++) {
ret.ma[i][j] += ma[i][k] * ma[k][j];
if(ret.ma[i][j] >= modnum)
ret.ma[i][j] %= modnum;
}
}
}
}
else if(!(x & 1)) {
ret = this->power(x >> 1);
ret = ret.power(2);
}
else {
ret = this->power(x >> 1);
ret = ret.power(2);
ret = ret * (*this);
}
return ret;
}

inline matrix operator+= (const matrix &x) {
int i, j;
for(i = 0; i < rw; i++) {
for(j = 0; j < cn; j++)
{
ma[i][j] += x.ma[i][j];
if(ma[i][j] >= modnum)
ma[i][j] %= modnum;
}
}
return *this;
}

inline matrix operator = (const matrix &x) {
int i, j;
rw = x.rw, cn = x.cn;
for(i = 0; i < rw; i++) {
for(j = 0; j < cn; j++) { ma[i][j] = x.ma[i][j]; }
}
return *this;
}
};

int main()
{
int n, m, k;
int i, j;
while(scanf("%d %d %d", &n, &k, &m) != EOF) {
matrix org(2 * n, 2 * n), ret(2 * n, 2 * n);
for(i = 0; i < n; i++) {
org.ma[i][i] = 1;
for(j = 0; j < n; j++) {
scanf("%d", &org.ma[i][n + j]);
org.ma[n + i][n + j] = org.ma[i][n + j];
}
}
modnum = m;
ret = org.power(k);
for(i = 0; i < n; i++) {
for(j = 0; j < n; j++) {
if(j) printf(" ");
printf("%d", ret.ma[i][n + j]);
}
printf("\n");
}
}
return 0;
}

comments powered by Disqus