POJ 3150 Cellular Automaton（矩阵快速幂）

Description

给出一个n个元素组成的环状序列，现在对这个序列进行k次同样的操作，每次操作，从第一个元素开始到最后一个元素结束，将每个元素的值更新为这个元素周围d范围内所有元素(2*d+1个)的和模m后的值，输出操作完后的序列

Input

第一行为四个整数n,m,d,k,第二行为n个整数表示这个环状序列

Output

输出操作完后的序列

Sample Input

sample input #1

5 3 1 1

1 2 2 1 2

sample input #2

5 3 1 10

1 2 2 1 2

Sample Output

sample output #1

2 2 2 2 1

sample output #2

2 0 0 2 2

Solution

显然是要构造转移矩阵，首先将原序列序列放入一个1*n矩阵

B=


下面来看转移矩阵的构造，以n=7,d=2为例，转移矩阵如下



轻易得出这个转移矩阵的每一行都是由上一行得到，即有转移方程A[i][j]=A[(i-1+n)%n][(j-1+n)%n]，这种性质十分有用，对于最初转移矩阵的构造遗迹矩阵快速幂时每次矩阵的乘法，我们都只需算出第一行的数据，之后递推求出后面n-1行即可，这样可以在矩阵乘法部分把时间复杂度由O(n^3)降到O(n^2)，下面来证明，两个具有这种性质的矩阵相乘之后仍然具有这种性质，



所以这种性质是具有积性的，故构造完转移矩阵A后用矩阵快速幂算出A^k然后用B*A^k即为答案，总时间复杂度为O(n^2*log k)

注意这道题要开的二维数组非常大，在函数中没法开，只能开全局，所以函数中要用指针传递数组

Code

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
using namespace std;
#define maxn 520
typedef long long ll;
int n,m,d,k;
ll num[maxn];
ll temp[maxn][maxn],ans[maxn][maxn],f[maxn][maxn];
void mod_mul(ll a[][maxn],ll b[][maxn])//矩阵乘法,结果放在b数组中
{
memset(temp,0,sizeof(temp));//初始化
for(int k=0;k<n;k++)//算出第一行
if(a[0][k])//优化
for(int j=0;j<n;j++)
{
temp[0][j]+=a[0][k]*b[k][j];
temp[0][j]%=m;
}
for(int i=1;i<n;i++)//递推算出剩下n-1行
for(int j=0;j<n;j++)
temp[i][j]=temp[i-1][(j-1+n)%n];
for(int i=0;i<n;i++)//复制结果
for(int j=0;j<n;j++)
b[i][j]=temp[i][j];
}
void mod_pow(ll a[][maxn])//矩阵快速幂
{
//初始化一个单位矩阵
memset(ans,0,sizeof(ans));
for(int i=0;i<n;i++) ans[i][i]=1;
while(k)
{
if(k&1) mod_mul(a,ans);
mod_mul(a,a);
k>>=1;
}
}
int main()
{
scanf("%d%d%d%d",&n,&m,&d,&k);
for(int i=0;i<n;i++)
{
scanf("%lld",&num[i]);
f[0][i]=0;
}
f[0][0]=1;
for(int i=1;i<=d;i++)//构造转移矩阵的第一行
f[0][i]=f[0][n-i]=1;
for(int i=1;i<n;i++)//递推算出剩下n-1行
for(int j=0;j<n;j++)
f[i][j]=f[i-1][(j-1+n)%n];
mod_pow(f);
ll end[maxn];//储存结果的数组
memset(end,0,sizeof(end));//初始化
for(int i=0;i<n;i++)//初始序列组成的1*n矩阵与转移矩阵相乘即为答案
for(int j=0;j<n;j++)
{
end[i]+=ans[i][j]*num[j];
end[i]%=m;
}
for(int i=0;i<n;i++)
printf("%lld%c",end[i],i==n-1?'\n':' ');
return 0;
}