POJ 3233 Matrix Power Series（矩阵快速幂前缀和）

题意：给出一个矩阵A求 S = A + A² + A³ + … + A^k.

思路：矩阵快速幂最重要的就是转移过程的实现，由当前状态转移到下一个状态的中间变化即使转移矩阵。

这道题目由两种解法:
1.直接找出S 和Aⁿ关系

先看看前面几个 S2=S1+A² S3=S2+A³ …

Sn                    E        E                 Sn-1
             =                                      
A^(n+1)               0        A                 A^n

这里我们直接把矩阵行列扩大一倍就可以直接做了

#include<stdio.h>
#include<iostream>
#include<cmath>
#include<string.h>
#define fi first
#define se second
#define log2(a) log(n)/log(2)
#define show(a) cout<<a<<endl;
#define show2(a,b) cout<<a<<" "<<b<<endl;
#define show3(a,b,c) cout<<a<<" "<<b<<" "<<c<<endl;
using namespace std;

typedef long long ll;

const ll inf = 1e17 + 10;
const int N = 1e6 + 10;
const ll mod = 1e9+7;
const int base=131;
const double pi=acos(-1);




int n,m,k;

struct mat{
	ll x[70][70];
	mat(){memset(x,0,sizeof x);}
};

mat mul(mat a,mat b)
{
	mat ans;
	for(int i=1;i<=2*n;i++)
	{
		for(int j=1;j<=2*n;j++)
		{
			for(int k=1;k<=2*n;k++)
			{
				ans.x[i][j]=(ans.x[i][j]+a.x[i][k]*b.x[k][j])%m;
			}
		}
	}
	return ans;
}

mat ksm(mat x,ll k)
{
	mat ans;
	for(int i=1;i<=2*n;i++) ans.x[i][i]=1;
	while(k)
	{
		if(k&1) ans=mul(ans,x);
		k>>=1;
		x=mul(x,x);
	}

	return ans;
}
int main( )
{

	ios::sync_with_stdio(false);
	cin.tie(0);

	cin>>n>>k>>m;
	mat ans,tr,a;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			cin>>a.x[i][j];
			tr.x[i+n][j+n]=ans.x[i][j]=a.x[i][j];
		}
	}
	a=mul(a,a);
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			ans.x[i+n][j]=a.x[i][j];
		}
	}

	for(int i=1;i<=n;i++)
	{
		tr.x[i][i]=tr.x[i][i+n]=1;
	}

	ans=mul(ksm(tr,k-1),ans);
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<n;j++)
		{
			cout<<ans.x[i][j]<<" ";
		}
		cout<<ans.x[i][n]<<endl;
	}



}

方法二：

就是通过不断的二分操作，按照如图判断n的奇偶性，进行上述公式的求和，显然这种方法会慢许多。

#include<stdio.h>
#include<iostream>
#include<cmath>
#include<string.h>
#define fi first
#define se second
#define log2(a) log(n)/log(2)
#define show(a) cout<<a<<endl;
#define show2(a,b) cout<<a<<" "<<b<<endl;
#define show3(a,b,c) cout<<a<<" "<<b<<" "<<c<<endl;
using namespace std;

typedef long long ll;

const ll inf = 1e17 + 10;
const int N = 1e6 + 10;
const ll mod = 1e9 + 7;
const int base = 131;
const double pi = acos ( -1 );




int n, m, k;
struct mat
{
	ll x[35][35];
	mat()
	{
		memset ( x, 0, sizeof x );
	}
};

mat ans, tr, a;

mat mul ( mat a, mat b )
{
	mat ans;
	for ( int i = 1; i <= n; i++ )
	{
		for ( int j = 1; j <= n; j++ )
		{
			for ( int k = 1; k <= n; k++ )
			{
				ans.x[i][j] = ( ans.x[i][j] + a.x[i][k] * b.x[k][j] ) % m;
			}
		}
	}
	return ans;
}
mat add ( mat a, mat b )
{
	mat ans;
	for ( int i = 1; i <= n; i++ )
	{
		for ( int j = 1; j <= n; j++ )
		{
			ans.x[i][j] = ( a.x[i][j] + b.x[i][j] ) % m;

		}
	}
	return ans;
}
mat ksm ( mat x, ll k )
{
	mat ans;
	for ( int i = 1; i <= n; i++ ) ans.x[i][i] = 1;
	while ( k )
	{
		if ( k & 1 ) ans = mul ( ans, x );
		k >>= 1;
		x = mul ( x, x );
	}
	return ans;
}

mat solve ( int k )
{
	if ( k == 1 ) return a;
	mat s = solve ( k / 2 );

	mat tmp = ksm ( a, k / 2 );

	mat q = mul ( tmp, s );

	s = add ( s, q );

	if ( k & 1 )
	{
		s = add ( s, ksm ( a, k ) );
	}
	return s;
}
int main( )
{

	ios::sync_with_stdio ( false );
	cin.tie ( 0 );

	cin >> n >> k >> m;
	for ( int i = 1; i <= n; i++ )
	{
		for ( int j = 1; j <= n; j++ )
		{
			cin >> a.x[i][j];

		}
	}
	ans = solve ( k );
	
	for ( int i = 1; i <= n; i++ )
	{
		for ( int j = 1; j < n; j++ )
		{
			cout << ans.x[i][j] << " ";
		}
		cout << ans.x[i][n] << endl;
	}

}