POJ 3233 Matrix Power Series

大意不再赘述。

思路：矩阵快速幂，需要用到两次二分，一次是求A^k，还有一次是求和A^1+A^2+A^3+......A^k时用到。

求和时：

当k为偶数：令t = A^1+A^2+A^3+....A^k/2，b = A^k/2，则S = t + t*b。

当k为奇数：令t = A^1+A^2+A^3+....A^k/2，b = A^k/2+1，则S = t + b + t*b。

递归求和即可。

#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cstring>
using namespace std;

const int MAXN = 101;

typedef struct
{
	int A[MAXN][MAXN];
}M;

int n, k, MOD;

M a, per;

void init()
{
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < n; j++)
		{
			scanf("%d", &a.A[i][j]);
			a.A[i][j] %= MOD;
			per.A[i][j] = (i==j);
		}
	}
}

M add(M a, M b)
{
	M c;
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < n; j++)
		{
			c.A[i][j] = (a.A[i][j] + b.A[i][j]) % MOD;
		}
	}
	return c;
}

M Mul(M a, M b)
{
	M c;
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < n; j++)
		{
			c.A[i][j] = 0;
			for(int k = 0; k < n; k++)
			{
				c.A[i][j] += a.A[i][k] * b.A[k][j];
				if(c.A[i][j] >= MOD) c.A[i][j] %= MOD;
			}
		}
	}
	return c;
}

M power(int k)
{
	M ans = per, p = a;
	while(k)
	{
		if(k & 1) ans = Mul(ans, p);
		p = Mul(p, p);
		k /= 2;
	}
	return ans;
}

M MatraxSum(int k)
{
	if(k == 1) return a;
	M t, b;
	t = MatraxSum(k/2);
	if(k & 1)
	{
		b = power(k/2+1);
		t = add(t, Mul(t, b));
		t = add(t, b);
	}
	else
	{
		b = power(k/2);
		t = add(t, Mul(t, b));
	}
	return t;
}

void print(M a)
{
	for(int i = 0; i < n; i++)
	{
		int first = 1;
		for(int j = 0; j < n; j++)
		{
			if(first)
			{
				printf("%d", a.A[i][j]);
				first = 0;
			}
			else printf(" %d", a.A[i][j]);
		}
		printf("\n");
	}
}

int main()
{
	while(~scanf("%d%d%d", &n, &k, &MOD))
	{
		init();
		M ans = MatraxSum(k);
		print(ans);
	}
	return 0;
}