题目描述

给定 $n \times n$ 矩阵 $A$ 和正整数 $k$，求和 $S = A + A ^ 2 + A ^ 3 + \cdots + A ^ k$。

输入格式

输入只包含一个测试用例。

第一行输入包含三个正整数 $n$，$k$ 和 $m$。

接下来 $n$ 行，每行包含 $n$ 个非负整数（均不超过 $32,768$），用以描绘矩阵 $A$。

输出格式

按与描述矩阵 $A$ 相同的方式，输出将 $S$ 中所有元素对 $m$ 取模后得到的矩阵。

数据范围

$1 \leq n \leq 30,$
$1 \leq k \leq 10 ^ 9,$
$1 \leq m < 10 ^ 4$

输入样例:

2 2 4
0 1
1 1

输出样例：

1 2
2 3

算法

(递归) $\mathcal O(n ^ 3 \log ^ 2 k)$

首先假设我们所求和的 $A$ 是一个数，而不是一个矩阵。
那么我们要快速求出 $S = A + A ^ 2 + A ^ 3 + \cdots + A ^ k$。

暴力肯定不可行，$k$ 的最大值为 $10 ^ 9$。

其次想到通项公式，$\displaystyle S = A + A ^ 2 + A ^ 3 + \cdots + A ^ k = A \times \frac{A ^ k - 1}{A - 1}$。
这里用到了除法，转到矩阵上，就是矩阵求逆，可我不会矩阵求逆鸭(wtcl)，所以不可行。

我们观察一下，这玩意好像可以拆开a
设 $S _ k = A + A ^ 2 + A ^ 3 + \cdots + A ^ k$
$\displaystyle S _ k = A + A ^ 2 + \cdots + A ^ \frac k 2 + A ^ {\frac k 2 + 1} + \cdots + A ^ k$
$\displaystyle \quad = A + A ^ 2 + \cdots + A ^ \frac k 2 + A ^ \frac k 2 (A + A ^ 2 + \cdots + A ^ \frac k 2)$
$\displaystyle \quad = (1 + A ^ \frac k 2) (A + A ^ 2 + \cdots + A ^ \frac k 2)$
$\displaystyle \quad = (1 + A ^ \frac k 2) S _ \frac k 2$

$1 + A ^ \frac k 2$ 可以用快速幂算，而 $S _ \frac k 2$ 可以递归算。

这样我们就可以快速求出 $S _ k$ 了。
那么在我们计算的对象是矩阵的时候呢？

同理，
$\displaystyle S _ k = A + A ^ 2 + \cdots + A ^ \frac k 2 + A ^ {\frac k 2 + 1} + \cdots + A ^ k$
$\displaystyle \quad = A + A ^ 2 + \cdots + A ^ \frac k 2 + A ^ \frac k 2 (A + A ^ 2 + \cdots + A ^ \frac k 2)$
$\displaystyle \quad = (E + A ^ \frac k 2) (A + A ^ 2 + \cdots + A ^ \frac k 2)$
$\displaystyle \quad = (E + A ^ \frac k 2) S _ \frac k 2$

那么我们就能快速求出对于矩阵的 $S _ k$ 了，
但这是当 $k$ 为偶数时的理想情况，$k$ 是奇数的时候呢？
当然，可以用 $S _ k = S _ {k - 1} + A ^ k$，但是这样我们就又需要计算一次快速幂。

我们可以用另一种拆法。
设 $k = 2n + 1$
那么 $S _ k = S _ {2n + 1} = A + A ^ 2 + \cdots + A ^ {2n + 1}$
$= A + A(A + \cdots + A ^ {2n})$
$= A + A \times S _ {2n}$
$= A + A \times S _ {k - 1}$

这样我们就可以在不多求一次快速幂的情况下，求出 $S _ k$ 啦~

时间复杂度

emmm这个窝也不是很会算qwq

每次矩阵乘法是 $\mathcal O(n ^ 3)$，
递归一共 $\log k$ 层，
每层要求一次快速幂，
所以总递归的复杂度应该是 $\mathcal O(\log ^ 2 k)$，
那总的时间复杂度就是 $\mathcal O(n ^ 3 \log ^ 2 k)$ 叭。

C++ 代码

这里提供两种写法：
1.结构体重载运算符写法

#include <stdio.h>
#include <string.h>

const int N=32;

int n, k, mod;
struct matrix    // 矩阵结构体
{
    int w[N][N]; // 存矩阵

    void init(bool type) // 矩阵初始化。type 为 true 则初始化为 E，type 为 false 则初始化为 O。
    {
        for (int i = 1; i <= n; i ++ )
            for (int j = 1; j <= n; j ++ )
                if (i != j) w[i][j] = 0;   // 两种矩阵的共同点：不在 i = j 对角线上的数皆为 0
                else    w[i][j] = type;    // 不同点：在 i = j 对角线上，E 为 1，O 为 0
    }

    matrix operator * (const matrix &t) const // 重载矩阵乘法
    {
        matrix ans;
        ans.init(false);
        for (int i = 1; i <= n; i ++ )
            for (int j = 1; j <= n; j ++ )
                for (int k = 1; k <= n; k ++ )
                    ans.w[i][j] = (ans.w[i][j] + t.w[i][k] * w[k][j]) % mod;
        return ans;
    }

    matrix operator + (const matrix &t) const // 重载矩阵加法
    {
        matrix ans;
        for (int i = 1; i <= n; i ++ )
            for (int j = 1; j <= n; j ++ )
                ans.w[i][j] = (w[i][j] + t.w[i][j]) % mod;
        return ans;
    }

    void print() // 将该矩阵输出
    {
        for (int i = 1; i <= n; i ++ ) // 不写奇怪的 for 循环了 2333
        {
            for (int j = 1; j <= n; j ++ )
                printf("%d ", w[i][j]);
            puts("");
        }
    }

    void read() // 读入该矩阵
    {
        for (int i = 1; i <= n; i ++ )
            for (int j = 1; j <= n; j ++ )
            {
                scanf("%d", &w[i][j]);
                w[i][j] %= mod;
            }
    }
};

matrix E, v; // E 即 E 矩阵，v 为读入矩阵

matrix operator ^ (matrix a, int k) // 重载矩阵快速幂。由于要用到乘法，所以在结构体外重载。
{
    matrix ans;
    ans.init(true);
    while (k)
    {
        if (k & 1) ans = ans * a;
        a = a * a, k >>= 1;
    }
    return ans;
}

matrix S(int k)        // 计算 S[k]
{
    if (k == 1) return v;        // 如果 k 为 1，那么返回 v，终止递归。
    if (k & 1) return v + v * S(k - 1);    // 如果 k 是奇数，那么返回上述 A + A * S(k - 1)
    return (E + (v ^ k >> 1)) * S(k >> 1); // 否则返回上述 (E + A ^ (k / 2)) * S(k / 2)
}

int main()
{
    scanf("%d%d%d", &n, &k, &mod); // 读入三个数值
    E.init(true); // 初始化矩阵 E
    v.read();     // 读入矩阵 v
    S(k).print(); // 计算结果 S(k) 并输出
    return 0;
}

2.纯数组模拟写法

#include <stdio.h>
#include <string.h>

const int N = 31;

int n, k;
int mod;
int E[N][N];
int w[N][N];

void init(int v[][N], bool type)   // 初始化，同上。
{
    for (int i = 1; i <= n; i ++ )
        for (int j = 1; j <= n; j ++ )
            v[i][j] = i == j ? type : 0;
}

void add(int A[][N], int B[][N])   // 计算矩阵 A + B 的结果，并存储在 A 中。
{
    for (int i = 1; i <= n; i ++ )
        for (int j = 1; j <= n; j ++ )
            A[i][j] = (A[i][j] + B[i][j]) % mod;
}

void multi(int A[][N], int B[][N]) // 计算矩阵 A * B 的结果，并存储在 A 中。
{
    static int ans[N][N];
    init(ans, false);
    for (int i = 1; i <= n; i ++ )
        for (int j = 1; j <= n; j ++ )
            for (int k = 1; k <= n; k ++ )
                ans[i][j] = (ans[i][j] + A[i][k] * B[k][j]) % mod;
    memcpy(A, ans, sizeof ans);
}

void qmi(int A[][N], int k)        // 计算 A 的 k 次方的结果，并存储在 A 中。
{
    static int ans[N][N];
    init(ans, true);
    while (k)
    {
        if (k & 1) multi(ans, A);
        multi(A, A);
        k >>= 1;
    }
    memcpy(A, ans, sizeof ans);
}

void S(int k)                 // 计算 S[k]，并存储在矩阵 w 中。 
{
    if (k == 1) return ;
    int tmp[N][N];
    memcpy(tmp, w, sizeof w); // 在计算前要先将 w 复制到 tmp 中
    if (k & 1)                // 如果 k 是奇数
    {
        S(k - 1);             // 先计算 S(k - 1)，并存储在矩阵 w 中。
        multi(w, tmp);        // 计算 w * 原w 的结果，并存储在矩阵 w 中。
        add(w, tmp);          // 计算 w + 原w 的结果，并存储在矩阵 w 中。
    }
    else
    {
        S(k >> 1);            // 计算 S(k / 2)，并存储在矩阵 w 中。
        qmi(tmp, k >> 1);     // 计算 原w 的 k / 2 次方，并存储在矩阵 tmp 中。
        add(tmp, E);          // 计算求完快速幂后的 原w + E 的结果，并存储在 tmp 中。
        multi(w, tmp);        // 计算 w * (E + 原w ^ (k / 2)) 的结果，并存储在 w 中。
    }
}

int main()
{
    scanf("%d%d%d", &n, &k, &mod);
    init(E, true);            // 初始化矩阵 E
    for (int i = 1; i <= n; i ++ )
        for (int j = 1; j <= n; j ++ )
        {
            scanf("%d", &w[i][j]);
            w[i][j] %= mod;
        }
    S(k);
    for (int i = 1; i <= n; i ++ )
    {
        for (int j = 1; j <= n; j ++ )
            printf("%d ", w[i][j]);
        puts("");
    }
    return 0;
}