题意：给定 $n\times n$ 的矩阵 $A$，问有多少 $n\times n$ 的矩阵 $B$ 满足 $AB=BA$。这里运算都在 $\bmod p$ 意义下进行，其中 $p=998244353$。$1\le n\le 500$

• 如果 $A^n=0$？

记 $F=\mathbb F_p$，将 $A$ 视同 $V=F^n$ 上的线性映射。$\ker(T)$ 定义为所有满足 $Tv=0$ 的向量 $v$ 的集合。

我们知道有如下的一系列子空间：$\ker(A)\subseteq \ker(A^2)\subseteq \cdots\subseteq \ker(A^n)=V$。

取 $\ker(A)$ 的一组基 $v_{1,1},\cdots,v_{1,s_1}$，接下来扩充为 $\ker(A^2)$ 的基 $v_{1,1},\cdots,v_{1,s_1},v_{2,1},\cdots,v_{2,s_2}$，依次类推，最后得到 $K$ 的一组基 $v_{i,j}$，其中 $1\le i\le n,1\le j\le s_i$，$s_i=\dim \ker(A^i)-\dim \ker(A^{i-1})$。

事实上，我们可以把基取的「更好看」一点：有 $v_{i,j}\in \ker(A^i),v_{i,j}\not\in \ker(A^{i-1})$，且有 $Av_{i,j}\in\ker(A^{i-1})$。于是可以将 $A$ 看作 $\ker(A^i)/\ker(A^{i-1})\to \ker(A^{i-1})/\ker(A^{i-2})$ 的线性映射（其中 $i\ge 2,A^0=I$），进一步，这个线性映射还是一个单射。

于是我们得到一些推论：必有 $s_{i-1}\ge s_{i}$，且我们可以取 $v_{i,1},\cdots,v_{i,s_i}$ 使得 $Av_{i,j}=v_{i-1,j}$ 对 $1\le j\le s_{i}$ 成立。

• 当然，我们有更简单的方式取这样一组基，也就是先做分解 $V=F[A]v_1\oplus F[A]v_2\oplus\cdots\oplus F[A]v_k$，然后在每个循环子空间 $F[A]v$ 内，假设 $A$ 限制在该子空间内的最小多项式为 $p^e$，其中 $\deg p=r$ 且 $p(x)$ 是一个不可约多项式，且 $v$ 生成整个子空间。设 $\deg p=r$，我们取 $v_{i,j}=p(A)^{i-1}A^{j-1}v$ 即可，其中 $1\le i\le e,1\le j\le r$。最后把这些子空间的基合并到一起就好了。

$B$ 由它在所有 $v_{i,j}$ 上的作用唯一确定。那么 $B$ 只需要满足对任意 $i,j$，有 $ABv_{i,j}=BAv_{i,j}=Bv_{i-1,j}$。对于 $i\ge 2$ 的情况，一旦确定 $Bv_{i,j}$ 是什么，那么 $1\le j\le s_i$ 的 $Bv_{i-1,j}$ 就唯一确定了，而对于 $s_i< j \le s_{i-1}$ 的 $Bv_{i-1,j}$ 则可以任意选取；而在 $i=1$ 的情况，这相当于 $ABv_{1,j}=BAv_{1,j}=0$，即 $Bv_{1,j}\in\ker(A)$。

进而有 $Bv_{1,j}=ABv_{2,j}\in \ker(A)$，那么必须要有 $Bv_{2,j}\in\ker(A^2)$。综上我们发现，$B$ 可以按照如下的步骤确定：对每个 $1\le i\le n$，以及 $s_{i+1}< j \le s_i$，我们确定一个 $Bv_{i,j}\in\ker(A^i)$。那么其自由度是多少呢？就是 $$ \sum_{i=1}^n(s_i-s_{i+1})\times \sum_{j\le i}s_j=\sum_{i=1}^ns_i^2\\ $$ 于是，我们只需要求出 $s_i=\dim\ker(A^i)-\dim\ker(A^{i-1})$，答案就是 $\sum_{i=1}^ns_i^2$。

• 一般情况的解法

如果仅仅追求理论上的解法，我们总可以将 $\mathbb F_p$ 进行扩域（扩域自然不改变解空间维数），使得最小多项式分裂，然后在每个广义特征子空间 $\ker((A-\lambda I)^n)$ 里面做。在这里面，我们可以先把 $A$ 化为 $A-\lambda I$（注意 $AB=BA\iff (A-\lambda I)B=B(A-\lambda I)$），化为幂零的情形。此外，由于每个子空间都是 $A$-不变的，那么它也是 $B$-不变的。最后，我们把每个子空间内的答案加起来就可以了。但是扩域的代价实在太大了。

首先我们考虑如何求解 $A$ 的有理标准型，我们设 $A$ 的最小多项式是 $m_A(x)=\prod p_i^{e_i}$，可以发现这样一个事实：随机选取 $v\in V$，满足 $v$ 相对与 $A$ 的最小多项式 $m_{v}(x)$ 大概率就是 $m_A(x)$，实际上这个概率至少是 $1-\frac{n}{p}$。这是因为首先必有 $m_v(x)\mid m_A(x)$，如果二者不相等那么必须存在 $p_i$ 使得 $m_v(x)\mid m_A(x)/p_i$，也就是说 $v\in \ker(m_A(x)/p_i)$，那么 $v$ 的取值至多是 $\deg m_A(x)$ 个 $n-1$ 维子空间的并，于是取到这其中的概率就 $\le \frac{\deg m_A(x)}{p}$。

那么我们每次随机选取一个向量 $v$，注意向量乘矩阵是 $O(n^2)$ 的，我们不断求 $v,Av,A^2v,\cdots$，直到有一个 $A^kv\in \text{span}(v,Av,\cdots,A^{k-1}v)$ 为止；此外到这里我们还可以顺便求出 $A$ 的最小多项式。

接下来在商空间 $V/F[A]v$ 中接着随机选取向量做下去即可。我们容易在 $O(\deg m_v(x)\times n^2)$ 的时间复杂度内完成一次分解，那么总的时间复杂度是不超过 $O(n^3)$ 的。

这个过程至多进行 $n$ 次，那么我们就有 $(1-\frac{n}{p})^n\ge 1-\frac{n^2}{p}$ 的概率得到 $V$ 的一个循环分解：$V=F[A]v_1\oplus F[A]v_2\oplus\cdots\oplus F[A]v_s$，且我们还顺便求出了所有 $v_i$ 的最小多项式。

• 实际上由于每次进行分解的时候 $\deg m_A(x)$ 的总和其实只有 $n$，因此 union bound 一下发现就有 $\prod (1-\deg m_A(x)/p)\le 1-n/p$，成功概率实际上至少有 $1-n/p$。

这实际上就给出了 $A$ 的不变因子组 $d_k(x)\mid d_{k-1}(x)\mid \cdots\mid d_1(x)$。考虑如何据此求出 $s_i$。

先假装扩域，使所有 $d_i(x)$ 分裂，我们还是只需要对每个特征值 $\lambda$ 对应的广义特征子空间去考虑。假设 $d_i(x)$ 中 $(x-\lambda)$ 的系数为 $c_i$，即 $(x-\lambda)^{c_i}\|d_i(x)$，那么可以发现限制在 $\ker((A-\lambda I)^n)$ 内时，有 $$ s_i=\dim\ker(A^i)-\dim\ker(A^{i-1})=\sum_{j=1}^k[c_j\ge i] $$ 那么有（注意必有 $c_i\ge c_{i+1}$） $$ \sum s_i^2=\sum_i\sum_{x,y}[c_x\ge i][c_y\ge i]=\sum_{x,y}\min(c_x,c_y)=\sum_{j=1}^k(2j-1)c_j $$ 这就是 $\lambda$ 对应的解空间维数。

这时候我们发现，并不需要真的求出分裂后的所有多项式，因为上面这个式子是满足可加性的，所以答案其实就是 $\sum_{j=1}^k(2j-1)\deg d_j(x)$。综上我们就解决了本题，时间复杂度为 $O(n^3)$。

碎碎念：EI 的题解里说要对 $d_i(x)$ 做无平方因子分解，但是好像不需要做这个？

下面的代码写的时候没想清楚，实际上线性基正常写就好了，不用像我这样写。

#include<bits/stdc++.h>

#define ll long long
#define mk make_pair
#define fi first
#define se second

using namespace std;

#define gc getchar
inline int read(){
    int x=0,f=1;char c=getchar();
    for(;(c<'0'||c>'9');c=getchar()){if(c=='-')f=-1;}
    for(;(c>='0'&&c<='9');c=getchar())x=x*10+(c&15);
    return x*f;
}

const int mod=998244353;
int ksm(int x,ll y,int p=mod){
    int ans=1;y%=(p-1);
    for(int i=y;i;i>>=1,x=1ll*x*x%p)if(i&1)ans=1ll*ans*x%p;
    return ans%p;
}
int inv(int x,int p=mod){return ksm(x,p-2,p)%p;}
mt19937 rnd(20070819);
int randint(int l,int r){
    uniform_int_distribution<>dist(l,r);
    return dist(rnd);
}
void add(int &x,int v){x+=v;if(x>=mod)x-=mod;}
void Mod(int &x){if(x>=mod)x-=mod;}
int cmod(int x){if(x>=mod)x-=mod;return x;}

template<typename T>void cmax(T &x,T v){x=max(x,v);}
template<typename T>void cmin(T &x,T v){x=min(x,v);}

int n;
vector<vector<int> >B;
vector<int>I;
void Ins(vector<int>v){
    for(int i=0;i<n;i++)if(v[i]!=0){B[i]=v,I[i]=inv(v[i]);break;}
}
vector<int>rec(vector<int>v){
    for(int i=0;i<n;i++)if(v[i]!=0){
        if(B[i][i]==0)return v;
        int r=1ll*I[i]*v[i]%mod;
        for(int j=i;j<n;j++)add(v[j],mod-1ll*B[i][j]*r%mod);
    }
    return v;
}

struct basis{
    vector<vector<int> >Vs;
    vector<int>ids,invs,vis;int D;
    bool ins(vector<int>v){
        if(vis.empty())vis.resize(n,0),D=0;
        for(int i=0;i<D;i++){
            int p=ids[i];
            for(int j=0;j<n;j++){
                if(v[j]==0&&Vs[p][j]==0)continue;
                if(v[j]==0&&Vs[p][j]!=0)break;
                if(v[j]!=0&&Vs[p][j]==0){
                    Vs.emplace_back(v),vis[j]=1;D++;
                    ids.insert(ids.begin()+i,D-1);
                    invs.emplace_back(inv(v[j]));
                    Ins(v);
                    return true;
                }

                int r=1ll*invs[p]*v[j]%mod;
                for(int k=j;k<n;k++)add(v[k],mod-1ll*r*Vs[p][k]%mod);
                assert(v[j]==0);break;
            }
        }
        for(int j=0;j<n;j++)if(v[j]!=0){
            Vs.emplace_back(v),ids.emplace_back(D);
            D++,invs.emplace_back(inv(v[j])),vis[j]=1;
            Ins(v);
            return true;
        }
        return false;
    }
};

vector<vector<int> >A;
vector<int>mul(vector<int>v){
    vector<int>res(n,0);
    for(int i=0;i<n;i++)for(int j=0;j<n;j++)add(res[j],1ll*v[i]*A[i][j]%mod);
    return res;
}
vector<int> work(vector<int>v){
    basis W;int cnt=0;
    while(W.ins(v))v=mul(v),v=rec(v);
    return W.vis;
}

signed main(void){

    n=read();A.resize(n,vector<int>(n)),B.resize(n,vector<int>(n)),I.resize(n);
    for(int i=0;i<n;i++)for(int j=0;j<n;j++)A[i][j]=read();

    vector<int>ban(n),degs;int sum=0;
    while(sum<n){
        vector<int>vc(n,0);
        for(int i=0;i<n;i++)if(!ban[i])vc[i]=randint(0,mod-1);
        auto vis=work(vc);int now=0;
        for(int i=0;i<n;i++)if(vis[i])ban[i]=1,now++;
        degs.emplace_back(now),sum+=now;
    }

    int ans=0;
    for(int x:degs)for(int y:degs)ans+=min(x,y);
    cout<<ans<<endl;

    return 0;
}