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巨的提醒下发现 然后还有结束标志，还有每个实例后面都有一个换行。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <iostream>
using namespace std;
const int maxn = 85;
int  n;
struct mat
{
    int  p[maxn][maxn];
};
mat mul (mat a, mat b)
{
    int i, j, k, m;
    m = n * 2;
    mat c;
    memset (c.p, 0, sizeof(c.p));
    for (i=0; i<m; i++)
        for (j=0; j<m; j++)
        {
            for (k=0; k<m; k++)
                c.p[i][j] += a.p[i][k] * b.p[k][j];
            c.p[i][j] = c.p[i][j] % 10;
        }
    return c;
}

mat pow (int n, mat a, mat b)
{
    while (n)
    {
        if (n % 2)
        {
            b = mul (b, a);
        }
        n /= 2;
        a = mul (a, a);
    }
    return b;
}

int main ()
{
    int k;
    while (scanf ("%d %d", &n, &k))
    {
        if (!n)
            break;
        mat a, b;
        memset (b.p, 0, sizeof(b.p));
        memset (a.p, 0, sizeof(a.p));

        for (int i=0; i<n; i++)
            for (int j=0; j<n; j++)
            {
                scanf ("%d", &a.p[i][j]);
                a.p[i][j] = a.p[i][j] % 10;
            }

        for (int i=0; i<n; i++)//构造矩阵，使a矩阵的右上，右下成为单位矩阵，把b也初始化为单位矩阵
            a.p[i][i+n] = a.p[i+n][i+n] = b.p[i][i] = b.p[i+n][i+n] = 1;

        b = pow (k+1, a, b);
        for (int i=0; i<n; i++)
            for (int j=0; j<n; j++)
            {
                if (i == j)//在b右上角的那个矩阵减去一个单位矩阵
                    {
                        b.p[i][j+n] --;
                        if (b.p[i][j+n] < 0)//防止出现末尾是零，减去单位矩阵是-1的情况。
                            b.p[i][j+n] = 9;
                    }
                if (j == n-1)
                    printf ("%d
", b.p[i][j+n]);
                else
                    printf ("%d ", b.p[i][j+n]);
            }
            printf ("
");
    }
    return 0;
}