POJ 3150 Cellular Automaton(矩陣快速冪)

題目大意：給定n（1<=n<=500）個數字和一個數字m，這n個數字組成一個環(a0,a1.....an-1)。如果對ai進行一次d-step操作，那麼ai的值變為與ai的距離小於d的所有數字之和模m。求對此環進行K次d-step（K<=10000000）後這個環的數字會變為多少。

看了一篇部落格：http://www.cppblog.com/varg-vikernes/archive/2011/02/08/139804.html說的很清楚。

拿樣例來說：

a矩陣：

a = 1 2 2 1 2

構造出來一個矩陣：

b = 
1 1 0 0 1
1 1 1 0 0
0 1 1 1 0
0 0 1 1 1
1 0 0 1 1

(a4+a0+a1) (a0+a1+a2) (a1+a2+a3) (a2+a3+a4) (a3+a4+a0)得到的矩陣e = 5 5 5 5 4

所以題意就是：a*(b^k) = e。求出e就是通過解b^k。因為直接求b^k會超時，所以要二分求冪。

這個題的矩陣還有一個性質：

利用矩陣A，B具有A[i][j]=A[i-1][j-1],B[i][j]=B[i-1][j-1]（i-1<0則表示i-1+n,j-1<0則表示j-1+n）
我們可以得出矩陣C=a*b也具有這個性質
C[i][j]=sum(A[i][t]*B[t][j])=sum(A[i-1][t-1],B[t-1][j-1])=sum(A[i-1][t],B[t][j-1])=C[i-1][j-1] 
“”“

所以可以求出C[0][0],以後的遞推得到。

Description

Input

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

#include <algorithm>
#include <iostream>
#include <stdlib.h>
#include <string.h>
#include <iomanip>
#include <stdio.h>
#include <string>
#include <queue>
#include <cmath>
#include <stack>
#include <map>
#include <set>
#define eps 1e-10
#define LL __int64
///#define LL long long
///#define INF 0x7ffffff
#define INF 0x3f3f3f3f
#define PI 3.1415926535898
#define zero(x) ((fabs(x)<eps)?0:x)

const int maxn = 550;

using namespace std;

LL e[maxn][maxn];
LL c[maxn][maxn];
LL xmul[maxn][maxn];
LL num[maxn];
LL ans[maxn];
int n, m, d;

void Mul(LL a[][maxn], LL b[][maxn])
{
    memset(c, 0, sizeof(c));
    for(int i = 0; i < n; i++)
    {
        if(!a[0][i])
            continue;
        for(int j = 0; j < n; j++)
        {
            c[0][j] += a[0][i]*b[i][j];
            c[0][j] %= m;
        }
    }
    for(int i = 1; i < n; i++)
        for(int j = 0; j < n; j++) c[i][j] = c[i-1][(j-1+n)%n];
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++) b[i][j] = c[i][j];
}

void expo(LL a[][maxn], int k)
{
    if(k == 1)
    {
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++) e[i][j] = a[i][j];
        return;
    }
    memset(e, 0, sizeof(e));
    for(int i = 0; i < n; i++) e[i][i] = 1;
    while(k)
    {
        if(k&1) Mul(a, e);
        Mul(a, a);
        k /= 2;
    }
}
int main()
{
    int k;
    while(cin >>n>>m>>d>>k)
    {
        for(int i = 0; i < n; i++) cin >>num[i];
        memset(xmul, 0, sizeof(xmul));
        xmul[0][0] = 1;
        for(int i = 1; i <= d; i++) xmul[0][i] = xmul[0][n-i] = 1;
        for(int i = 1; i < n; i++)
            for(int j = 0; j < n; j++) xmul[i][j] = xmul[i-1][(j-1+n)%n];
        expo(xmul, k);
        memset(ans, 0, sizeof(ans));
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                ans[i] += e[i][j]*num[j];
                ans[i] %= m;
            }
        }
        for(int i = 0; i < n-1; i++) cout<<ans[i]<<" ";
        cout<<ans[n-1]<<endl;
    }
    return 0;
}