[poj 3233]Matrix Power Series --- 构造 + 矩阵快速幂

传送门:poj3233

---

题目大意

给定一个$n\times n$的矩阵$A$,求矩阵$S$,其中$S=A+A^2+\cdots +A^k$,且对$m$取模.

---

分析

乍一看,是个等比数列求和,然后扯上除法、逆元,然后$\dots$就没有然后了,可以直接弃疗了.
　　对此不得不换个思路,考虑一下构造递推式.
　　由秦九韶算法可知$S_i=A\ast S_{i-1}+A$,于是可以考虑像fibonacci那样构造一个辅助矩阵来优化递推
　　于是就有了如下的式子：
$$\left[\begin{array}{cc}{S}_i& A\end{array}\right]\times \left[\begin{array}{cc}A& O\\ I& I\end{array}\right]=\left[\begin{array}{cc}{S}_{i+1}& A\end{array}\right]$$
　　　　注:$I$为单位矩阵，$O$为零矩阵
　　因此直接快速幂即可
　　PS:网上也有直接自乘的,但其实差不多($S_1=A$)

---

代码

#include <cstdio>
#include <cstdlib>

#define IL inline

using namespace std;

IL int read()
{
    int sum = 0;
    bool k = 1;
    char c = getchar();
    for(;'0' > c || c > '9'; c = getchar())
    if(c == '-')
        k = 0;
    for(;'0' <= c && c <= '9'; c = getchar())
        sum = sum * 10 + (c - '0');
    return k ? sum : -sum;
}

typedef long long ll;

int mo;

struct matrix
{   
    int n, m;
    ll s[65][65];
    
    IL matrix(int n_ = 0, int m_ = 0)
    {
        n = n_; m = m_;
        for(int i = 1; i <= n; ++i)
        for(int j = 1; j <= m; ++j)
            s[i][j] = 0;
    }
    
    IL matrix operator * (const matrix &b)
    {
        matrix c(n, b.m);
        for(int i = 1; i <= n; ++i)
        for(int k = 1; k <= m; ++k)
        if(s[i][k])
        for(int j = 1; j <= b.m; ++j)
        if(b.s[k][j])
        {
            c.s[i][j] += s[i][k] * b.s[k][j] % mo;
            if(c.s[i][j] >= mo) c.s[i][j] -= mo;
        }
        return c;
    }
    
    IL void out()
    {
        for(int i = 1; i <= n; ++i)
        {
            for(int j = 1; j <= n; ++j)
                printf("%lld ", s[i][j]);
            printf("\n");
        }
    }
}a, b;

IL matrix power(matrix x, int t)
{
    matrix sum(x.n, x.m);
    
    for(int i = 1; i <= x.n; ++i) sum.s[i][i] = 1;
    
    for(; t; t >>= 1)
    {
        if(t & 1) sum = sum * x;
        x = x * x;
    }
    return sum;
}


int main()
{
    int n = read(), k = read(); 
    mo = read();
    
    a.n = n; b.n = b.m = a.m = (n << 1);
    
    for(int i = 1; i <= n; ++i)
    {
        for(int j = 1; j <= n; ++j)
            a.s[i][j] = a.s[i][j + n] = b.s[i][j] = read();
        b.s[i + n][i] = b.s[i + n][i + n] = 1;
    }
    
    if(k == 1) a.out(); else (a * power(b, k - 1)).out();
    
    return 0;
}