poj3233Matrix Power Series矩阵快速幂

/*
 * 思路：这题数据k<10^9，很显然需要O(log n)的算法，自然就是矩阵
 * 快速幂．矩阵快速幂和快速幂差不多就是求矩阵A的k次方时，每次折半
 * 类似于二分的思想，比如要求A的１００次方，我们先求出A的５０次方，
 * 然后用得到的结果平方就好了．试想如果要算16次方，每次折半，我们
 * 只需要做4次矩阵乘法，但是如果采用朴素的连乘则需要15次，效率对比
 * 可想而知．
 * 但是还有一个问题，我们用矩阵快速幂可以求得一个矩阵的n次方．但是
 * 这题要求的是一系列的矩阵幂之和，就算每个矩阵幂都用上上面的快速
 * 幂也会超时啊！表示俺线性代数学得太糟糕了，看了网友的构造方法，
 * 才学会转换．如下:
 * 构造一个矩阵
 * B = |A En|
 *     |0 En|
 * 那么B^(k+1) = |A^(k+1) En+A+A^2+A^3+...+A^k|
 *              |  0            En           |
 * 其中En表示n阶单位矩阵．
 * 所以要求A+A^2+...+A^k，只需要求出B^(k+1)，然后右上角那个分
 * 块矩阵减单位矩阵就可以了．方法真是看似简单，但是炒鸡厉害．．．
 * 计算思想乳此感人，代码就算了，各有各的风格．如果你看到了本文，
 * 程序千万要自己实现，别直接模仿或是抄写别人的代码．
 */
#include <cstdio>
#include <cstring>
#include <vector>
#include <iostream>
using namespace std;

const int MAX = 65;
typedef long long ll;
int n, k, m;
int SIZE;

struct mat {
	int a[MAX][MAX];
	mat(){}
	mat(int b) {
		memset(a, 0, sizeof(a));
		for (int i = 0; i < MAX; ++i) {
			a[i][i] = b;
		}
	}//init

	mat operator*(mat b) {
		mat tmp(0);
		int sum;
		/************矩阵乘法***************/
		/*for (int i = 0; i < MAX; ++i) {
			for (int j = 0; j < MAX; ++j) {
				sum = 0;
				for (int k = 0; k < MAX; ++k) {
					sum += a[i][k]*b.a[k][j]; sum %= m;
				}
				tmp.a[i][j] = sum;
			}
		}*/
		/**********优化版，原理以后学了计组再说*********/

		for (int i = 0; i < MAX; ++i) {
			for (int k = 0; k < MAX; ++k) {
				for (int j = 0; j < MAX; ++j) {
					tmp.a[i][j] += a[i][k]*b.a[k][j]; tmp.a[i][j] %= m;
				}
			}
		}

		return tmp;
	}

	mat operator=(mat b) {
		memcpy(a, b.a, sizeof(b.a));
		return *this;
	}
};
mat YY(1);

mat fast_pow(mat a, int x) {
	mat res(1);
    while (x) {
		if (x&1) res = res*a;
		a = a*a;
		x >>= 1;
    }
    return res;
}

int main() {
	mat M(MAX);
	while (~scanf(" %d %d %d", &n, &k, &m)) {
		SIZE = n << 1;
		for (int i = 0; i < n; ++i) {
			for (int j = 0; j < n; ++j) {
				scanf(" %d", &M.a[i][j]);
			}
		}
		/* 我们要的是构造后的矩阵B，上面输入存入的只是B左上角的
		 * 分块矩阵，另外三块还要手动搞
		 */
		/* 右上角n阶单位矩阵 */
		for (int i = 0; i < n; ++i) {
			for (int j = n; j < SIZE; ++j) {
				if (j-n == i) M.a[i][j] = 1;
				else M.a[i][j] = 0;
			}
		}
		/* 左下角０矩阵 */
		for (int i = n; i < SIZE; ++i) {
			for (int j = 0; j < n; ++j) {
				M.a[i][j] = 0;
			}
		}
		/* 右下角n阶单位矩阵 */
		for (int i = n; i < SIZE; ++i) {
			for (int j = n; j < SIZE; ++j) {
				M.a[i][j] = (i == j ? 1 : 0);
			}
		}

		mat res = fast_pow(M, k+1);

		/* 减去单位矩阵 */
		for (int i = 0; i < n; ++i) {
			res.a[i][i+n]--;
			if (res.a[i][i+n] < 0) res.a[i][i+n] += m;
		}

		/* 打印 */
		for (int i = 0; i < n; ++i) {
			for (int j = n; j < SIZE; ++j) {
				printf("%d%c",res.a[i][j], (j^(SIZE-1) ? ' ' : '\n'));
			}
		}
	}
	return 0;
}