poj 3233 Matrix Power Series 矩阵构造+快速幂

点击打开链接

要求A+A²+...+A^k，而不是单个矩阵的幂

即求A(E+A+..A^k-1) 可以不断乘以K次A

那么我们可以构造一个分块的辅助矩阵 S，其中 A 为原矩阵，E 为单位矩阵，O 为0矩阵

即 S=[A,E]   容易用归纳得到S^k=[A^k,A^k-1+A^k-2+...A+E]

 [0 ,E] [0,E]

#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long ll;
const int N=80;
int a[N][N],b[N][N],n,k,mod;
void mul(int x[N][N],int y[N][N])
{
	int i,j,k;
	int c[N][N];
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			c[i][j]=0;
			for(int k=1;k<=n;k++)
			{
				c[i][j]=(c[i][j]+x[i][k]*y[k][j])%mod;
			}
		}
	}
	 for (i=1;i<=n;i++)
        for (j=1;j<=n;j++) 
		x[i][j]=c[i][j];
}
void quickpow(int n)
{
	while(n)
	{
		if(n&1)
		mul(b,a);
		
		mul(a,a);
		n>>=1;
	}
}
int main()
{
	cin>>n>>k>>mod;
	memset(a,0,sizeof(a));
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			cin>>a[i][j];
			a[i][i+n]=a[i+n][i+n]=b[i][i]=b[i+n][i+n]=1;
			//b为单位阵 
		}
	}//构造分块矩阵
	n*=2;
	quickpow(k+1); 
	
	n/=2;
	for(int i=1;i<=n;i++)
	b[i][i+n]--;
	
	
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			printf("%d",b[i][j+n]);
			if(j==n)
			cout<<endl;
			else
			cout<<' ';
		}	
	}
	return 0;
}