LA 3704细胞自动机——循环矩阵&&矩阵快速幂

题目

一个细胞自动机包含 $n$ 个格子，每个格子的取值为 $0 sim m-1$。给定距离 $d$，则每次操作是将每个格子的值变为到它的距离不超过 $d$ 的所有格子的在操作之前的值的和除以 $m$ 的余数。给出 $n, m, d, k$ 和自动机各个格子的初始值。你的任务是计算 $k$ 次操作以后各格子的值。（$1 leq nleq 500, 1 leq mleq 10^6, 0 leq dleq n/2, 1leq kleq 10^7$）.

分析

如果我们把 $t$ 次操作以后的各格子值写成列向量 $v_t$，不难发现 $v_{t+1}$ 的每一维都是 $v_t$ 中各维的线性组合，其中的加法和乘法都是在模 $m$ 的剩余系中完成。

每次操作相当于乘以一个 $n imes n $ 矩阵，直接使用矩阵快速幂的复杂度为 $O(n^3logk)$，

由于这里的矩阵比较特殊，是循环矩阵（从第二行开始每一行都是上一行循环右移），

可以证明，两个循环矩阵的乘积仍然为循环矩阵。

因此在存储时只需保存第一行，而计算矩阵乘法时也只需算出第一行即可。这样，矩阵乘法的时间复杂度降为 $O(n^2)$。总时间降为 $O(n^2log k)$，可以承受。

用FFT优化的话可做到 $0(nlognlogk)$.　　

$$egin{bmatrix}
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
end{bmatrix} imes
egin{bmatrix}
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
end{bmatrix} =
egin{bmatrix}
3 & 2 & 1 & 1 & 2\
2 & 3 & 2 & 1 & 1\
1 & 2 & 3 & 2 & 1\
1 & 1 & 2 & 3 & 2\
2 & 1 & 1 & 2 & 3
end{bmatrix}$$

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

typedef long long ll;
const int maxn = 500+10;
struct matrix
{
    int n;
    ll mat[maxn];
    matrix(){
        memset(mat, 0, sizeof(mat));
    }
};
ll n, p, d, k;
ll a[maxn];

matrix mul(matrix A, matrix B)   //矩阵相乘,这里A=B，且都是n x n的方阵
{
    matrix ret;
    ret.n = A.n;
    for(int i = 0;i < A.n;i++)
        for(int j = 0;j < B.n;j++)  ret.mat[i] = (ret.mat[i] +  A.mat[j] * B.mat[(j-i+A.n)%A.n]) % p;
    return ret;
}

matrix mpow(matrix A, int n)
{
    matrix ret;
    ret.n = A.n;
    ret.mat[0]=1;
    while(n)
    {
        if(n & 1)  ret = mul(ret, A);
        A = mul(A, A);
        n >>= 1;
    }
    return ret;
}

int  main()
{
    while(scanf("%lld%lld%lld%lld", &n, &p,&d, &k) == 4)
    {
        for(int i = 0;i < n;i++)  scanf("%lld", &a[i]);
        matrix A;
        A.n = n;
        for(int i = 0;i <= d;i++) A.mat[i] = 1;
        for(int i = n-d; i < n;i++)  A.mat[i] = 1;

        A = mpow(A, k);

        for(int i = 0;i < A.n;i++)
        {
            ll tmp = 0;
            for(int j = 0;j < A.n;j++)  tmp = (tmp + A.mat[(j-i+n) % n] * a[j]) % p;
            printf("%lld%c", tmp, i == n-1 ? '
' : ' ');
        }
    }
    return 0;
}

记得开 long long !!

参考链接： https://vjudge.net/status/#un=&OJId=UVALive&probNum=3704&res=0&orderBy=run_id&language=