LA 3704 (矩阵快速幂 循环矩阵) Cellular Automaton

将这n个格子看做一个向量，每次操作都是一次线性组合，即v_n+1 = Av_n，所求答案为A^kv₀

A是一个n*n的矩阵，比如当n=5，d=1的时候：

不难发现，A是个循环矩阵，也就是将某一行所有元素统一向右移动一位便得到下一行。

而且循环矩阵相乘仍然是循环矩阵，所以只要求出A^k的第一行就行了。

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;

const int maxn = 500 + 10;
typedef long long Vector[maxn];

int n, m, d, k;

inline int value(Vector A, int i, int j)
{//循环矩阵A(i, j)的值
    return A[(((j-i)%n)+n)%n];
}

void matrix_mul(Vector A, Vector B, Vector res)
{
    Vector C;
    memset(C, 0, sizeof(C));
    for(int j = 0; j < n; j++)
        for(int k = 0; k < n; k++)
            C[j] = (C[j] + A[k] * value(B, k, j)) % m;
    memcpy(res, C, sizeof(C));
}

void matrix_pow(Vector A, int n, Vector res)
{
    Vector a, r;
    memcpy(a, A, sizeof(a));
    memset(r, 0, sizeof(r));
    r[0] = 1;
    while(n)
    {
        if(n&1) matrix_mul(r, a, r);
        n >>= 1;
        matrix_mul(a, a, a);
    }
    memcpy(res, r, sizeof(r));
}

void solve(Vector A, Vector v, Vector res)
{
    Vector B;
    memset(B, 0, sizeof(B));
    for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
            B[i] = (B[i] + value(A, i, j) * v[j]) % m;
    memcpy(res, B, sizeof(B));
}

int main()
{
    //freopen("in.txt", "r", stdin);

    while(cin >> n >> m >> d >> k)
    {
        Vector A, v;
        for(int i = 0; i < n; i++) cin >> v[i];
        memset(A, 0, sizeof(A));
        for(int p = -d; p <= d; p++)
        {
            int x = ((p%n)+n)%n;
            A[x] = 1;
        }
        /*for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
                printf("%d ", value(A, i, j));
            printf("
");
        }*/
        matrix_pow(A, k, A);
        solve(A, v, v);
        for(int i = 0; i < n; i++)
        {
            if(i) printf(" ");
            cout << v[i];
        }
        printf("
");
    }

    return 0;
}