[POJ 3233] Matrix Power Series （矩阵快速幂）

链接

POJ 3233

---

题意

A是一个矩阵，求A + A^2 + A^3 + … A^k。

---

题解

有一个经典解法，构造矩阵：
|A E |
|0 E |
该矩阵自乘k次，得到的矩阵为：
|A^k E + A + A^1 + … + A^(k-1)|
|0 E |
很容易验证。

感觉它的巧妙之处在于向量(A^k, sum(k-1))和向量T(E, E)相乘后能得到sum(k)， 和向量T(A, 0)相乘后能得到A^(k+1)，这样向后迭代就好。
效率上用快速幂优化一下就好。

另外还有一个经典的矩阵二分求解，复杂度是一样的，感觉效率上会快一些，毕竟这种构造的做法需要将矩阵规模扩大一倍。
po个链接吧：
二分做法

---

代码

#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;
int mod;
struct Mat
{
    #define maxk (100)
    int a[maxk][maxk], n;
    void clear() { memset(a, 0, sizeof(a)); }
    Mat(int k, int type) {
        n = k;
        clear();
        if(type) for(int i = 0; i < n; i++)
            a[i][i] = 1;
    }
    Mat operator*(const Mat& b) const
    {
        Mat o = Mat(n, 0);
        for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
        for(int k = 0; k < n; k++)
        {
            o.a[i][j] += (a[i][k] * b.a[k][j]) % mod;
            o.a[i][j] %= mod;
        }
        return o;
    }
    friend Mat operator^(Mat tmp, int k)
    {
        Mat o = Mat(tmp.n, 1);
        while(k)
        {
            if(k & 1) o = o * tmp;
            tmp = tmp * tmp;
            k >>= 1;
        }
        return o;
    }
};

int main()
{
    int n, k, m;
    while(cin >> n >> k >> m)
    {
        mod = m;
        Mat mat = Mat(n*2, 0);
        for(int i = 0; i < n; i++)
        for(int j = 0; j < n; j++)
        scanf("%d", &mat.a[i][j]);
        for(int i = 0; i < n; i++)
        {
            mat.a[i][n + i] = 1;
            mat.a[n + i][n + i] = 1;
        }
        Mat ans = mat^(k + 1);
        for(int i = 0; i < n; i++)
        ans.a[i][n + i] = (ans.a[i][n + i] + mod - 1) % mod;
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                if(j) putchar(' ');
                printf("%d", ans.a[i][n + j]);
            }
            printf("\n");
        }
    }
    return 0;
}