POJ 3233

原题链接

描述

Given a n × n matrix A and a positive integer k, find the sum \(S = A + A^{2} + A^{3} + … + A^{k}.\)

输入

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10e9) and m (m < 104). Then follow n lines each containing n nonnegative integers below 32768, giving A’s elements in row-major order.

输出

Output the elements of S modulo m in the same way as A is given.

样例输入

2 2 4
0 1
1 1

样例输出

1 2
2 3

思路

起初我直接矩阵快速幂，然后把结果相加，可是TLE了，然后去百度看看别人的思路，找到一个神奇的公式。
\[ B = \left[ \begin{matrix} A & E \\ 0 & E \end{matrix} \right] \]

\[ B^{k+1} = \left[ \begin{matrix} A^{k+1} & E + A + A^{2} + A^{3} + … + A^{k} \\ 0 & E \end{matrix} \right] \]
有了这个公式之后，写了个分块矩阵就好，顺便重构了一下我那惨不忍睹的矩阵快速幂模版。
虽然被大佬嫌弃代码啰嗦，但还是放上来吧。

代码

#include <cstdio>
#include<cstring>
#define ll long long
#define maxn 32
using namespace std;

int n, k, mod;

struct Mat
{
    ll f[maxn][maxn];
    void cls(){memset(f, 0, sizeof(f));}//全部置为0 
    Mat() {cls();}
    void myprint(int n)//输出n阶顺序主子式 
    {
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
                printf("%d ", f[i][j] % mod);
            printf("\n");
        }
    }
    friend Mat operator + (Mat a, Mat b)
    {
        Mat res;
        for(int i = 0; i < maxn; i++)
            for(int j = 0; j < maxn; j++)
                res.f[i][j] = a.f[i][j] + b.f[i][j];
        return res;
    }
    friend Mat operator - (Mat a, Mat b)
    {
        Mat res;
        for(int i = 0; i < maxn; i++)
            for(int j = 0; j < maxn; j++)
            {
                res.f[i][j] = a.f[i][j] - b.f[i][j];
                while(res.f[i][j] < 0) res.f[i][j] += mod;
            }
        return res;
    }
    friend Mat operator * (Mat a, Mat b)
    {
        Mat res;
        for(int i = 0; i < maxn; i++)for(int j = 0; j < maxn; j++)
            for(int k = 0; k < maxn; k++)
                (res.f[i][j] += a.f[i][k] * b.f[k][j]) %= mod;
        return res;
    }
} E, I;

struct MatDiv
{
    Mat s[2][2];
    void set(Mat x1, Mat x2, Mat x3, Mat x4)
    {
        s[0][0] = x1; s[0][1] = x2;
        s[1][0] = x3; s[1][1] = x4;
    }
    friend MatDiv operator * (MatDiv a, MatDiv b)
    {
        MatDiv res;
        for(int i = 0; i < 2; i++) for(int j = 0; j < 2;j++)
            for(int k = 0; k < 2; k++)
                res.s[i][j] = a.s[i][k] * b.s[k][j] + res.s[i][j];
        return res;
    }
}; 

MatDiv quick_pow(MatDiv a)  
{  
    MatDiv ans;
    ans.set(E, I, E, I);
    ll b = k;
    while(b != 0)  
    {
        if(b & 1) ans = ans * a;
        b >>= 1;
        a = a * a;
    }
    return ans;  
}

int main()
{
    for(int i = 0; i < maxn; i++)
        E.f[i][i] = 1;
    while(~scanf("%d %d %d", &n, &k, &mod))
    {
        Mat A; k++;
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++)
                scanf("%d", &A.f[i][j]);
        MatDiv B; B.set(A, E, I, E);
        B = quick_pow(B);
        A = B.s[0][1] - E;
        A.myprint(n);
    }
    return 0;
}

转载于:https://www.cnblogs.com/HackHarry/p/8391677.html