UVA 11149 Power of Matrix 构造矩阵

题目大意：意思就是让求A（A是矩阵）+A²⁺A³⁺A⁴⁺A⁵⁺A^{6+······+}A^K，其中矩阵范围n<=40,k<=1000000。

解题思路：由于k的取值范围很大，所以很自然地想到了二分法，用递归逐步将k二分（公式：A+A²+A³+A⁴+A⁵+A⁶ = A+A²+A³ + A³（A+A²+A³））,

　　这种方法只需要注意k是奇数的情况就可以了。

　　最坑的是第二种方法，根据矩阵的性质可以构造出来一个子矩阵，假如有矩阵B=|A  E| ，那么B^K =|A^K   E+ A+A²⁺A³⁺A⁴⁺A⁵⁺A^{6+······+}A^K|

|0  E|                |0          E                                         |

　　呵呵········，这种方法wa了好多次，我曾经开始怀疑线性代数老师是不是讲错了。最后在T巨的提醒下发现 然后还有结束标志，还有每个实例后面都有一个换行。

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <algorithm>
 4 #include <iostream>
 5 using namespace std;
 6 const int maxn = 85;
 7 int  n;
 8 struct mat
 9 {
10     int  p[maxn][maxn];
11 };
12 mat mul (mat a, mat b)
13 {
14     int i, j, k, m;
15     m = n * 2;
16     mat c;
17     memset (c.p, 0, sizeof(c.p));
18     for (i=0; i<m; i++)
19         for (j=0; j<m; j++)
20         {
21             for (k=0; k<m; k++)
22                 c.p[i][j] += a.p[i][k] * b.p[k][j];
23             c.p[i][j] = c.p[i][j] % 10;
24         }
25     return c;
26 }
27
28 mat pow (int n, mat a, mat b)
29 {
30     while (n)
31     {
32         if (n % 2)
33         {
34             b = mul (b, a);
35         }
36         n /= 2;
37         a = mul (a, a);
38     }
39     return b;
40 }
41
42 int main ()
43 {
44     int k;
45     while (scanf ("%d %d", &n, &k))
46     {
47         if (!n)
48             break;
49         mat a, b;
50         memset (b.p, 0, sizeof(b.p));
51         memset (a.p, 0, sizeof(a.p));
52
53         for (int i=0; i<n; i++)
54             for (int j=0; j<n; j++)
55             {
56                 scanf ("%d", &a.p[i][j]);
57                 a.p[i][j] = a.p[i][j] % 10;
58             }
59
60         for (int i=0; i<n; i++)//构造矩阵，使a矩阵的右上，右下成为单位矩阵，把b也初始化为单位矩阵
61             a.p[i][i+n] = a.p[i+n][i+n] = b.p[i][i] = b.p[i+n][i+n] = 1;
62
63         b = pow (k+1, a, b);
64         for (int i=0; i<n; i++)
65             for (int j=0; j<n; j++)
66             {
67                 if (i == j)//在b右上角的那个矩阵减去一个单位矩阵
68                     {
69                         b.p[i][j+n] --;
70                         if (b.p[i][j+n] < 0)//防止出现末尾是零，减去单位矩阵是-1的情况。
71                             b.p[i][j+n] = 9;
72                     }
73                 if (j == n-1)
74                     printf ("%d\n", b.p[i][j+n]);
75                 else
76                     printf ("%d ", b.p[i][j+n]);
77             }
78             printf ("\n");
79     }
80     return 0;
81 }