http://poj.org/problem?id=3233

ps转：

用二分方法求等比数列前n项和：即

原理：

(1)若n==0

(2)若n%2==0 

(3)若n%2==1

代码如下：

LL sum(LL p,LL n)
{
    if(n==0) return 1;
    if(n&1)  return (1+pow(p,(n>>1)+1))*sum(p,n>>1);
    else     return (1+pow(p,(n>>1)+1))*sum(p,(n-1)>>1)+pow(p,n>>1);
}

那么类比等比数列，我们可以二分等比矩阵中的幂K,每次分成一半求。

下面是二分方法，这份代码可以做模板啦。。

#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <set>
#include <map>
#include <vector>
#include <math.h>
#include <string.h>
#include <queue>
#include <string>
#include <stdlib.h>
#define LL long long
#define _LL __int64
#define eps 1e-8
#define PI acos(-1.0)
using namespace std;
const int maxn = 35;

int mod,n,k,m;

struct matrix
{
	int mat[maxn][maxn];
	//初始化为单位矩阵
	void init()
	{
		memset(mat,0,sizeof(mat));
		for(int i = 1; i <= maxn; i++)
			mat[i][i] = 1;
	}
}a,res;
//矩阵相加
matrix matrixAdd(matrix x, matrix y)
{
	matrix tmp;

	for(int i = 1; i <= n; i++)
		for(int j = 1; j <= n; j++)
		{
			tmp.mat[i][j] = x.mat[i][j] + y.mat[i][j];
			if( tmp.mat[i][j] >= mod )
				tmp.mat[i][j] %= mod;
		}
	return tmp;
}
//矩阵相乘
matrix matrixMul(matrix x, matrix y)
{
	matrix tmp;
	memset(tmp.mat,0,sizeof(tmp.mat));

	for(int i = 1; i <= n; i++)
	{
		for(int k = 1; k <= n; k++)
		{
			if(x.mat[i][k] == 0) continue;
			for(int j = 1; j <= n; j++)
			{
				tmp.mat[i][j] += x.mat[i][k] * y.mat[k][j];
				if(tmp.mat[i][j] >= mod)
					tmp.mat[i][j] %= mod;
			}
		}
	}
	return tmp;
}
//矩阵求幂
matrix Mul(matrix x, int k)
{
	matrix tmp;
	tmp.init();

	while(k)
	{
		if(k & 1)
			tmp = matrixMul(tmp,x);
		x = matrixMul(x,x);
		k >>= 1;
	}
	return tmp;
}
// 求A + A^1 + A^2 + ... + A^k
matrix Sum(matrix x, int k)
{
	if(k == 1)
		return x;
	matrix tmp;
	tmp.init();

	tmp = matrixAdd(tmp,Mul(x,k>>1));
	tmp = matrixMul(tmp,Sum(x,k>>1));
	if(k&1)
		tmp = matrixAdd(tmp,Mul(x,k));
	return tmp;
}

void output()
{
	for(int i = 1; i <= n; i++)
	{
		for(int j = 1; j <= n; j++)
		{
			if(j < n)
				printf("%d ",res.mat[i][j]);
			else printf("%d\n",res.mat[i][j]);
		}
	}
}

int main()
{
	scanf("%d %d %d",&n,&k,&m);
	mod = m;
	for(int i = 1; i <= n; i++)
		for(int j = 1; j <= n; j++)
			scanf("%d",&a.mat[i][j]);

	res = Sum(a,k);
	output();
	return 0;
}

还有一种更快的方法，便是构造矩阵。可惜不是我自己想出来的。膜拜，矩阵太强大了。

我们要求的矩阵设为A，先构造这样的矩阵

B = |A A|，可以计算B^2 = |A^2    A+A^2| .... B^k = |A^k   A+A^2+...+A^k| ，

       |0  1|                              |0              1   |                 |0                1              |

因此我们只需求出B ^k，然后取其右上角的n*n的矩阵便是答案。。

#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <set>
#include <map>
#include <vector>
#include <math.h>
#include <string.h>
#include <queue>
#include <string>
#include <stdlib.h>
#define LL long long
#define _LL __int64
#define eps 1e-8
#define PI acos(-1.0)
using namespace std;
const int maxn = 100;

int mod,n,k,m;

struct matrix
{
	int mat[maxn][maxn];
	//初始化为单位矩阵
	void init()
	{
		memset(mat,0,sizeof(mat));
		for(int i = 1; i <= n; i++)
			mat[i][i] = 1;
	}
}a,res;
//矩阵相加
matrix matrixAdd(matrix x, matrix y)
{
	matrix tmp;

	for(int i = 1; i <= n; i++)
		for(int j = 1; j <= n; j++)
		{
			tmp.mat[i][j] = x.mat[i][j] + y.mat[i][j];
			if( tmp.mat[i][j] >= mod )
				tmp.mat[i][j] %= mod;
		}
	return tmp;
}
//矩阵相乘
matrix matrixMul(matrix x, matrix y)
{
	matrix tmp;
	memset(tmp.mat,0,sizeof(tmp.mat));

	for(int i = 1; i <= n; i++)
	{
		for(int k = 1; k <= n; k++)
		{
			if(x.mat[i][k] == 0) continue;
			for(int j = 1; j <= n; j++)
			{
				tmp.mat[i][j] += x.mat[i][k] * y.mat[k][j];
				if(tmp.mat[i][j] >= mod)
					tmp.mat[i][j] %= mod;
			}
		}
	}
	return tmp;
}
//矩阵求幂
matrix Mul(matrix x, int k)
{
	matrix tmp;
	tmp.init();

	while(k)
	{
		if(k & 1)
			tmp = matrixMul(tmp,x);
		x = matrixMul(x,x);
		k >>= 1;
	}

	return tmp;
}


void output()
{
	for(int i = 1; i <= n; i++)
	{
		for(int j = n+1; j <= n*2; j++)
		{
			if(j < n*2)
				printf("%d ",res.mat[i][j]);
			else printf("%d\n",res.mat[i][j]);
		}
	}
}

int main()
{
	while(~scanf("%d %d %d",&n,&k,&m))
	{
		mod = m;
		//扩大A矩阵，看成是4个小矩阵组合而成
		for(int i = 1; i <= n; i++)
			for(int j = 1; j <= n; j++)
				scanf("%d",&a.mat[i][j]);

		for(int i = 1; i <= n; i++)
			for(int j = n+1; j <= n*2; j++)
				a.mat[i][j] = a.mat[i][j-n];

		for(int i = n+1; i <= n*2; i++)
			for(int j = 1; j <= n; j++)
				a.mat[i][j] = 0;

		for(int i = n+1; i <= n*2; i++)
			for(int j = n+1; j <= n*2; j++)
			
		n = n << 1;
		res = Mul(a,k);
		n = n >> 1;
		output();
	}
	return 0;
}