POJ 3150 Cellular Automaton --矩阵快速幂及优化

题意：给一个环，环上有n块，每块有个值，每一次操作是对每个点，他的值变为原来与他距离不超过d的位置的和，问k(10^7)次操作后每块的值。

解法：一看就要化为矩阵来做，矩阵很好建立，大白书P157页有讲，大概为：

[1 1 0 .. 0 1]

[1 1 1 .. .. 0]

...

[1 1 .. .. .. 1]  的循环矩阵，可以证明，循环矩阵的乘积还是循环矩阵，且循环矩阵的性质： a[i][j] = a[i-1][j-1] (循环的) ，所以，我们每次矩阵相乘只需要算出第一行，余下的不需要通过矩阵乘法来算出，直接根据规律推出，这样，矩阵乘法的复杂度就降到了O(n^2)，加上快速幂，总复杂度O(n^2log(k))。

注意：中间相乘的时候a[i][k]*b[k][j]可能会超过int范围，要加一个long long，否则会WA.

代码：(6000+ ms  也是醉了。。)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#define SMod m
#define ll long long
using namespace std;

int n,m,k,d;
struct Matrix
{
    int m[501][501];
    Matrix()
    {
        memset(m,0,sizeof(m));
        for(int i=1;i<=n;i++)
            m[i][i] = 1;
    }
};

Matrix Mul(Matrix a,Matrix b)
{
    Matrix res;
    int i,j,k;
    for(j=1;j<=n;j++)
    {
        res.m[1][j] = 0;
        for(k=1;k<=n;k++)
            res.m[1][j] = (res.m[1][j]+(ll)a.m[1][k]*b.m[k][j]%SMod + SMod)%SMod;
    }
    for(i=2;i<=n;i++)
    {
        for(j=2;j<=n;j++)
            res.m[i][j] = res.m[i-1][j-1];
        res.m[i][1] = res.m[i-1][n];
    }
    return res;
}

Matrix fastm(Matrix a,int b)
{
    Matrix res;
    while(b)
    {
        if(b&1)
            res = Mul(res,a);
        a = Mul(a,a);
        b >>= 1;
    }
    return res;
}

Matrix Muti(Matrix a,Matrix b)
{
    Matrix res;
    int i,j,k;
    for(i=1;i<=n;i++)
    {
        res.m[i][1] = 0;
        for(k=1;k<=n;k++)
            res.m[i][1] = (res.m[i][1]+(ll)a.m[i][k]*b.m[k][1]%SMod + SMod)%SMod;
    }
    return res;
}

int main()
{
    int i,j;
    while(scanf("%d%d%d%d",&n,&m,&d,&k)!=EOF)
    {
        Matrix R;
        memset(R.m,0,sizeof(R.m));
        for(i=1;i<=n;i++)
            scanf("%d",&R.m[i][1]),R.m[i][1]%=SMod;
        Matrix A;
        for(i=2;i<=d+1;i++)
            A.m[1][i] = 1;
        for(i=n;i>=n-d+1;i--)
            A.m[1][i] = 1;
        for(i=2;i<=n;i++)
        {
            for(j=2;j<=n;j++)
                A.m[i][j] = A.m[i-1][j-1];
            A.m[i][1] = A.m[i-1][n];
        }
        /*for(i=1;i<=n;i++)
        {
            for(j=1;j<=n;j++)
            {
                if(min(abs(i-j),n-abs(i-j)) <= d)
                    A.m[j][i] = 1;
                else
                    A.m[j][i] = 0;
            }
        }*/
        Matrix ans = fastm(A,k);
        ans = Muti(ans,R);
        for(i=1;i<=n;i++)
            printf("%d%c",ans.m[i][1]%m,i==n?'\n':' ');
    }
    return 0;
}