POJ 3233 Matrix Power Series(矩阵乘法+分治/矩阵加速递推/倍增)

POJ 3233

矩阵乘法+分治/矩阵加速递推/倍增
分治:
$$\begin{gathered} sum=A^1+A^2+...+A^k \\ =A^1+A^2+...+A^{k/2}+A^{k/2+1}+...+A^k \\ =(I+A^{k/2})(A^1+A^2+...+A^{k/2}) \end{gathered}$$(I为单位矩阵)
对于$su{m}_i=A^1+A^2+A^3+...+A^k$,可以将其分为两部分求解;
若k为偶数,则使$sum=(I+A^{k/2})(A^1+A^2+...+A^{k/2})$,用矩阵快速幂求出$A^{k/2}$,其余继续递推下去;
若k为奇数,则计算出$A^k$,然后递推剩余的$sum(1--k-1)$;
(2547MS)

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#define RG register
using namespace std;
const int maxn=35;
int n,m,k;

int read()
{
	int res=0,fh=1;char ch=getchar();
	while(ch>'9'||ch<'0') {if(ch=='-') fh=-1; ch=getchar();}
	while(ch>='0'&&ch<='9') res=(res<<3)+(res<<1)+(ch^'0'),ch=getchar();
	return fh*res;
}

struct matrix{
	int a[maxn][maxn];
	matrix() { memset(a,0,sizeof(a)); }
	matrix operator +(matrix b) {
		matrix c;
		for(RG int i=1;i<=n;++i)
			for(RG int j=1;j<=n;++j)
				c.a[i][j]=(1ll*a[i][j]+b.a[i][j])%m;
		return c;
	}
	matrix operator *(matrix b) {
		matrix c;
		for(RG int i=1;i<=n;++i)
			for(RG int j=1;j<=n;++j)
				for(RG int s=1;s<=n;++s)
					c.a[i][j]=(1ll*c.a[i][j]+1ll*a[i][s]*b.a[s][j])%m;
		return c;
	}
}I,A;

matrix pw(int s)
{
	matrix t=I,B=A;
	for(;s;s>>=1) {
		if(s&1) t=t*B;
		B=B*B;
	}
	return t;
}

matrix work(int s)
{
	if(s==1) return A;
	if(s&1) return work(s-1)+pw(s);
	else return work(s>>1)*(I+pw(s>>1));
}

int main()
{
	n=read(),k=read(),m=read();
	for(RG int i=1;i<=n;++i)
		for(RG int j=1;j<=n;++j)
			A.a[i][j]=read()%m,I.a[i][j]= (i==j);
	matrix S=work(k);
	for(RG int i=1;i<=n;++i)
	{
		for(RG int j=1;j<n;++j)
			printf("%d ",S.a[i][j]);
		printf("%d\n",S.a[i][n]);
	}
	
	return 0;
}

矩阵加速递推:
设$f[i]=A^1+...+A^i$
则$f[i]=A\ast f[i-1]+A$
$\left[\begin{array}{cc}f[i-1]& A\end{array}\right]\ast \left[\begin{array}{cc}A& 0\\ I& I\end{array}\right]=\left[\begin{array}{cc}f[i]& A\end{array}\right]$
用矩阵快速幂即可;
(1344MS)

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
const int maxn=35;
int n,m,k;

int read()
{
	int res=0,fh=1;char ch=getchar();
	while(ch>'9'||ch<'0') {if(ch=='-') fh=-1; ch=getchar();}
	while(ch>='0'&&ch<='9') res=(res<<3)+(res<<1)+(ch^'0'),ch=getchar();
	return fh*res;
}

struct matrix{
	int a[maxn][maxn];
	matrix() { memset(a,0,sizeof(a)); }
	matrix operator +(matrix b) {
		matrix c;
		for(int i=1;i<=n;++i)
			for(int j=1;j<=n;++j)
				c.a[i][j]=(1ll*a[i][j]+b.a[i][j])%m;
		return c;
	}
	matrix operator *(matrix b) {
		matrix c;
		for(int i=1;i<=n;++i)
			for(int j=1;j<=n;++j)
				for(int s=1;s<=n;++s)
					c.a[i][j]=(1ll*c.a[i][j]+1ll*a[i][s]*b.a[s][j]%m)%m;
		return c;
	}
}A,f[2],S[2][2];

void mul()
{
	matrix m[2];
	for(int j=0;j<2;++j)
		for(int k=0;k<2;++k)
			m[j]=m[j]+f[k]*S[k][j];
	f[0]=m[0];
}

void mulself()
{
	matrix m[2][2];
	for(int i=0;i<2;++i)
		for(int j=0;j<2;++j)
			for(int k=0;k<2;++k)
				m[i][j]=m[i][j]+S[i][k]*S[k][j];
	for(int i=0;i<2;++i)
		for(int j=0;j<2;++j)
			S[i][j]=m[i][j];
}

int main()
{
	n=read(),k=read(),m=read();
	for(int i=1;i<=n;++i)
		for(int j=1;j<=n;++j)
		{
			A.a[i][j]=f[1].a[i][j]=S[0][0].a[i][j]=read();
			S[1][0].a[i][j]=S[1][1].a[i][j]= (i==j);
		}
	for(;k;k>>=1)
	{
		if(k&1) mul();
		mulself();
	}
	for(int i=1;i<=n;++i)
	{
		for(int j=1;j<n;++j)
			printf("%d ",f[0].a[i][j]);
		printf("%d\n",f[0].a[i][n]);
	}
	
	return 0;
}

倍增:
设$F[i]=A^{2^i},sum[i]=A^1+A^2+A^3+...+A^{2^i}$
则对于i有$F[i]=F[i-1]\ast F[i-1],sum[i]=sum[i-1]\ast F[i-1]+sum[i-1]$
对于k,其在二进制下若第i位为1,则答案$t=t\ast F[i]+sum[i]$
(344MS)

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
const int maxn=35;
int n,m,k;

int read()
{
	int res=0,fh=1;char ch=getchar();
	while(ch>'9'||ch<'0') {if(ch=='-') fh=-1; ch=getchar();}
	while(ch>='0'&&ch<='9') res=(res<<3)+(res<<1)+(ch^'0'),ch=getchar();
	return fh*res;
}

struct matrix{
	int a[maxn][maxn];
	matrix() { memset(a,0,sizeof(a)); }
	matrix operator +(matrix b) {
		matrix c;
		for(int i=1;i<=n;++i)
			for(int j=1;j<=n;++j)
				c.a[i][j]=(1ll*a[i][j]+b.a[i][j])%m;
		return c;
	}
	matrix operator *(matrix b) {
		matrix c;
		for(int i=1;i<=n;++i)
			for(int j=1;j<=n;++j)
				for(int s=1;s<=n;++s)
					c.a[i][j]=(1ll*c.a[i][j]+1ll*a[i][s]*b.a[s][j]%m)%m;
		return c;
	}
}A[maxn],sum[maxn],t;//sum_i---A^1+A^2+A^3+...+A^2^i

int main()
{
	n=read(),k=read(),m=read();
	for(int i=1;i<=n;++i)
		for(int j=1;j<=n;++j)
			A[0].a[i][j]=sum[0].a[i][j]=read();
	
	for(int i=1; (1<<i)<=k; ++i)
		A[i]=A[i-1]*A[i-1],sum[i]=sum[i-1]*A[i-1]+sum[i-1];
	for(int i=0;(1<<i)<=k;++i)
		if(k&(1<<i)) t=t*A[i],t=t+sum[i];
	for(int i=1;i<=n;++i)
	{
		for(int j=1;j<n;++j)
			printf("%d ",t.a[i][j]);
		printf("%d\n",t.a[i][n]);
	}
	
	return 0;
}