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

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

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

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

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

 1 #include <iostream>
 2 #include <cstdio>
 3 #include <cstring>
 4 using namespace std;
 5
 6 const int maxn = 500 + 10;
 7 typedef long long Vector[maxn];
 8
 9 int n, m, d, k;
10
11 inline int value(Vector A, int i, int j)
12 {//循环矩阵A(i, j)的值
13     return A[(((j-i)%n)+n)%n];
14 }
15
16 void matrix_mul(Vector A, Vector B, Vector res)
17 {
18     Vector C;
19     memset(C, 0, sizeof(C));
20     for(int j = 0; j < n; j++)
21         for(int k = 0; k < n; k++)
22             C[j] = (C[j] + A[k] * value(B, k, j)) % m;
23     memcpy(res, C, sizeof(C));
24 }
25
26 void matrix_pow(Vector A, int n, Vector res)
27 {
28     Vector a, r;
29     memcpy(a, A, sizeof(a));
30     memset(r, 0, sizeof(r));
31     r[0] = 1;
32     while(n)
33     {
34         if(n&1) matrix_mul(r, a, r);
35         n >>= 1;
36         matrix_mul(a, a, a);
37     }
38     memcpy(res, r, sizeof(r));
39 }
40
41 void solve(Vector A, Vector v, Vector res)
42 {
43     Vector B;
44     memset(B, 0, sizeof(B));
45     for(int i = 0; i < n; i++)
46         for(int j = 0; j < n; j++)
47             B[i] = (B[i] + value(A, i, j) * v[j]) % m;
48     memcpy(res, B, sizeof(B));
49 }
50
51 int main()
52 {
53     //freopen("in.txt", "r", stdin);
54
55     while(cin >> n >> m >> d >> k)
56     {
57         Vector A, v;
58         for(int i = 0; i < n; i++) cin >> v[i];
59         memset(A, 0, sizeof(A));
60         for(int p = -d; p <= d; p++)
61         {
62             int x = ((p%n)+n)%n;
63             A[x] = 1;
64         }
65         /*for(int i = 0; i < n; i++)
66         {
67             for(int j = 0; j < n; j++)
68                 printf("%d ", value(A, i, j));
69             printf("\n");
70         }*/
71         matrix_pow(A, k, A);
72         solve(A, v, v);
73         for(int i = 0; i < n; i++)
74         {
75             if(i) printf(" ");
76             cout << v[i];
77         }
78         printf("\n");
79     }
80
81     return 0;
82 }

代码君