HDU4686 Arc of dream

矩阵快速幂

ai = ai-1 * Ax + Ay;

bi = bi-1 * Bx + By;

ai * bi = ai-1 * bi-1 * (Ax * Bx) + ai-1 * (Ax * By) + bi-1 * (Bx * Ay) + AyBy;

Si = Si - 1 + ai * bi = Si - 1  + ai-1 * bi-1 * (Ax * Bx) + ai-1 * (Ax * By) + bi-1 * (Bx * Ay) + AyBy;

根据SI ， 一定需要  Si - 1，  ai - 1,  bi - 1,  ai - 1 * bi - 1,  常数。

矩阵大小为5，

构造 初始矩阵A 为 记  S0

[1， A0,  B0,  A0 * B0,  A0 * B0]  即 [1, a0, b0, a0 * b0, S0]

[0 …… 0]

转换矩阵 B.   使得  A * B = S1 = [1,  a1,  b1, a1 * b1 , S1]

根据上面4个公式得。

第一列为  1， 0， 0， 0， 0；

第二列为  Ay,  Ax,   0,   0,   0;

第三列为  By,  0,   Bx,  0,   0;

第四列为  AyBy, Ax * By,   Bx * Ay,  Ax * Bx, 0;

第五列为  AyBy, Ax * By,   Bx * Ay,  Ax * Bx, 1;

即

N次递推得

现在要求 Sn - 1 ,   所以  init * Pow ^ n - 1.   然后取第一行， 第5个元素。

[cpp]  view plain copy print ?

1. #include <cstdio>  
2. #include <cstring>  
3. #include <algorithm>  
4. using namespace std;  
5. const int V = 100000 + 50;  
6. const int mod = 1000000000 + 7;  
7. const int MaxN = 5;  
8. struct Matrix{  
9.     __int64 mat[MaxN][MaxN];  
10. };  
11. Matrix init, Pow;  
12. __int64 n, A0, Ax, Ay, B0, Bx, By;  
13. Matrix multi(Matrix a, Matrix b) {  
14.     Matrix ans;  
15.     for(int i = 0; i < MaxN; ++i)  
16.     for(int j = 0; j < MaxN; ++j) {  
17.         __int64 sum = 0;  
18.         for(int k = 0; k < MaxN; ++k)  
19.             sum = (sum + a.mat[i][k] * b.mat[k][j] % mod) % mod;  
20.         ans.mat[i][j] = sum;  
21.     }  
22.     return ans;  
23. }  
24. Matrix MatrixQuickPow(Matrix a, __int64 b) {  
25.     Matrix ans = init;  
26.     while(b) {  
27.         if(b & 1)  
28.             ans = multi(ans, a);  
29.         b /= 2;  
30.         a = multi(a, a);  
31.     }  
32.     return ans;  
33. }  
34. int main() {  
35.     while(~scanf("%I64d", &n)) {  
36.         scanf("%I64d%I64d%I64d%I64d%I64d%I64d", &A0, &Ax, &Ay, &B0, &Bx, &By);  
37.         if(!n) {  
38.             printf("0\n");  
39.             continue;  
40.         }  
41.         init.mat[0][0] = 1;  
42.         init.mat[0][1] = A0;  
43.         init.mat[0][2] = B0;  
44.         init.mat[0][3] = init.mat[0][4] = A0 * B0 % mod;  
45.         Pow.mat[0][0] = Pow.mat[4][4] = 1;  
46.         Pow.mat[0][1] = Ay;  
47.         Pow.mat[0][2] = By;  
48.         Pow.mat[0][3] = Pow.mat[0][4] = Ay * By % mod;  
49.         Pow.mat[1][1] = Ax;  
50.         Pow.mat[2][2] = Bx;  
51.         Pow.mat[1][3] = Pow.mat[1][4] = Ax * By % mod;  
52.         Pow.mat[2][3] = Pow.mat[2][4] = Ay * Bx % mod;  
53.         Pow.mat[3][3] = Pow.mat[3][4] = Ax * Bx % mod;  
54.         /* 
55. for(int i = 0; i < MaxN; ++i) { 
56.             for(int j = 0; j < MaxN; ++j) 
57.                 printf("%I64d ", Pow.mat[i][j]); 
58.             printf("\n"); 
59.         }*/  
60.         Matrix ans = MatrixQuickPow(Pow, n - 1);  
61.         printf("%I64d\n", ans.mat[0][4]);  
62.     }  
63. }