POJ3233 – Matrix Power Series(二分+矩阵快速幂)

题目链接：POJ3233

【题意】给出N*N矩阵A，求A^1+A^2+A^3+…+A^K；k(<=10^9)

【分析】对于，A^k次，用快速幂即可，关键是求出Sk = A^1+A^2+A^3+…+A^K;

网上抄来的公式：

然后就是一个n×n的矩阵。

我做这题时用的是类似二分的方法做的，仔细观察可以发现

S_k = A + A² + A³ + … + A^k   

    =(1+A^k/2)*(A + A² + A³ + … + A^k/2  )+{A^k}

    =(1+A^k/2)*(S_k/2 )+{A^k}

当k为偶数时不要大括号里面的数

【AC CODE】1688ms

#include <cstdio>
#include <cstring>
#include <cctype>
#include <cmath>
#include <map>
//#include <unordered_map>
#include <queue>
#include <stack>
#include <vector>
#include <string>
#include <algorithm>
using namespace std;
typedef long long LL;
#define rep(i,a,n) for(int i = a; i < n; i++)
#define repe(i,a,n) for(int i = a; i <= n; i++)
#define per(i,n,a) for(int i = n; i >= a; i--)
#define clc(a,b) memset(a,b,sizeof(a))
const int INF = 0x3f3f3f3f, MAXN = 31;
int n,MOD;
struct MATRIX {
	int num[MAXN][MAXN];
	void init(bool e){
		clc(num,0);
		if(e) rep(i,0,n) num[i][i] = 1;
	}
}a;

MATRIX mul(const MATRIX& a, const MATRIX& b)
{
	MATRIX ans;
	ans.init(0);
	rep(i,0,n)
	{
		rep(j,0,n)
		{
			rep(k,0,n) ans.num[i][j] = (ans.num[i][j]+a.num[i][k]*b.num[k][j])%MOD;
		}
	}
	return ans;
}
MATRIX pow_mod(MATRIX x, int c)
{
	MATRIX ans;
	ans.init(1);
	while(c)
	{
		if(c&1) ans = mul(ans,x);
		x = mul(x,x);
		c >>= 1;
	}
	return ans;
}
MATRIX add(const MATRIX &a, const MATRIX &b)
{
	MATRIX ans;
	rep(i,0,n)
	{
		rep(j,0,n) ans.num[i][j] = (a.num[i][j]+b.num[i][j])%MOD;
	}
	return ans;
}
MATRIX sum(int k)
{
	if(1 == k) return a;
	MATRIX ans;
	ans.init(1);
	ans = add(ans,pow_mod(a,k>>1));//1+A^(k/2)
	ans = mul(ans,sum(k>>1));//(1+A^(k/2))*S(k/2)
	if(k&1) ans = add(ans,pow_mod(a,k));
	return ans;
}
int main()
{
#ifdef SHY
	freopen("e:\\1.txt", "r", stdin);
#endif
	int k;
	scanf("%d %d %d%*c", &n, &k, &MOD);
	rep(i,0,n)
	{
		rep(j,0,n) scanf("%d%*c", &a.num[i][j]);
	}
	MATRIX ans = sum(k);
	rep(i,0,n)
	{
		printf("%d",ans.num[i][0]);
		rep(j,1,n) printf(" %d", ans.num[i][j]);
		putchar('\n');
	}
	return 0;
}