[Medium-Hard] AT_arc190_d [ARC190D] Matrix Pow Sum

简要题意

给定一个 $n\times n$ 的方阵 $A$，若 $A_{ij}=0$，表示 $A_{ij}$ 不确定。

给出整数 $p$，你需要将每个不确定的位置填上一个 $[1,p)$ 的整数得到矩阵 $B$。求全体 $B^p$ 的和。答案对 $p$ 取模。

$1\leq n\leq 100,p\in\mathbb{P},p\leq 10^9$。

思路

对于 $p=2$，求解是平凡的，这里只考虑 $p>2$ 的情况。

不妨将每个不确定的位置 $A_{ij}$ 看成一个自由元 $x_{ij}$。先考察 $A_{ij}^k$，则每个位置都是一个多元至多 $p$ 次的多项式。

考虑每一个单元格的多项式的每一项，应当形如 $C x_1^{k_1} x_2^{k_2} \cdots$，则这一项产生的贡献也应当形如 $C\left(\sum_{i=1}^{p-1} i^{k_1}\right)\left(\sum_{i=1}^{p-1} i^{k_2}\right)\cdots$。这启发我们往 $F_p(k)=\left(\sum_{i=1}^{p-1} i^k\right)\bmod p$ 去思考。

直觉告诉你这个 $F_p(k)$ 一定是好求的，考虑打表：

F[p_] := Table[Mod[Sum[i ^ k, {i, 1, p - 1}], p], {k, 0, p - 1}]
F[17]

得到 {16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16}，于是做出猜想：

$$F_p(k) = \begin{cases} p - 1& p - 1\mid k\\ 0& \mathrm{otherwise} \end{cases}$$

尝试证明这个结论。取模 $p$ 意义下的某一个原根 $g$，若 $p-1\nmid k$ 则得到：

$$F_p(k)=\sum_{i=1}^{p-1}i^k\bmod p=\sum_{i=0}^{p-2} g^{ik} \bmod p=\frac{g^{k(p-1)}-1}{g^k-1} \bmod p=\frac{1-1}{g^k-1}=0$$

否则，$i^k\equiv 1\pmod{p}$，所以 $F_p(k)=p-1$，证毕。

所以只有形如 $Cx_1^{p-1}$ 的项或常数项才会对答案产生贡献。

先来考虑常数项，这是容易的，$A^p$ 就是答案，然后尝试求出形如 $x_1^{p-1}$ 的系数之和。

看起来十分困难，不妨模拟矩阵乘法的过程，分析每一项是如何得到的。

对于 $p>3$，我们只有下面的两种方法：

• $x_{ii}\cdot \cdots\cdot x_{ii}\cdot A_{ij}$。
• $A_{ji}\cdot x_{ii}\cdot\cdots\cdot x_{ii}$。

对于 $p=3$，还有一种方法 $x_{ij}\cdot A_{ji}\cdot x_{ij}$。依照这种方法统计系数即可。

时间复杂度 $O(\mathrm{MatrixMul}(n)\log p)$。

代码

#include <bits/stdc++.h>
using namespace std;

int n, mod;

int Add(int x, int y){ return (x + y) >= mod ? (x + y - mod) : (x + y); }
int Sub(int x, int y){ return (x - y) < 0 ? (x - y + mod) : (x - y); }
int Mul(int x, int y){ return 1ll * x * y % mod; }
int fastpow(int x, int y){
    int ret = 1;
    for(;y;y>>=1,x=Mul(x, x)){ if(y & 1) ret = Mul(ret, x); }
    return ret;
}

const int N = 105;

struct Matrix {
    int a[N][N];
    Matrix(){ memset(a, 0, sizeof(a)); }
    int* operator[](int x){ return a[x]; }
    void idt(){
        for(int i=1;i<=n;i++) a[i][i] = 1;
    }
} A, B;

Matrix operator*(Matrix A, Matrix B){
    Matrix C;
    for(int i=1;i<=n;i++){
        for(int k=1;k<=n;k++){
            for(int j=1;j<=n;j++) C[i][j] = Add(C[i][j], Mul(A[i][k], B[k][j]));
        }
    }
    return C;
}

Matrix operator^(Matrix A, int k){
    Matrix B; B.idt();
    while(k){
        if(k & 1) B = B * A;
        A = A * A, k >>= 1;
    }
    return B;
}

signed main(){
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);
    cin >> n >> mod; int cnt = 0;
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++) cin >> A[i][j], cnt += (A[i][j] == 0);
    }
    if(mod == 2){
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++) A[i][j] = 1;
        }
        B = A ^ mod;
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++) cout << B[i][j] << ' ';
            cout << '\n';
        }
        return 0;
    }
    B = A ^ mod;
    for(int i=1;i<=n;i++){
        if(A[i][i]) continue;
        for(int j=1;j<=n;j++){
            B[i][j] = Add(B[i][j], A[i][j]);
            B[j][i] = Add(B[j][i], A[j][i]);
        }
    }
    if(mod == 3){
        for(int i=1;i<=n;i++){
            for(int j=1;j<=n;j++){
                if(A[i][j]) continue;
                B[i][j] = Add(B[i][j], A[j][i]);
            }
        }
    }
    int tmp = cnt & 1 ? mod - 1 : 1;
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++) cout << Mul(B[i][j], tmp) << ' ';
        cout << '\n';
    }
    return 0;
}

// Written by xiezheyuan

Submission。