LA 3704 (矩阵快速幂 循环矩阵) Cellular Automaton
将这n个格子看做一个向量,每次操作都是一次线性组合,即vn+1 = Avn,所求答案为Akv0
A是一个n*n的矩阵,比如当n=5,d=1的时候:
不难发现,A是个循环矩阵,也就是将某一行所有元素统一向右移动一位便得到下一行。
而且循环矩阵相乘仍然是循环矩阵,所以只要求出Ak的第一行就行了。
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
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 }