一個比較慢的做法
首先你要知道矩陣的特征多項式是什么。
直接消元就可以了。
時間復雜度：\(O(n^5)\)或\(O(n^4)\)。
一個稍微快一點的做法
觀察到特征多項式的次數是\(n\)。
我們就可以插值。
具體來說，先求出當\(x=0\ldots n\)時特征多項式對應的點值，然后直接用拉格朗日插值插出來。
時間復雜度：\(O(n^4)\)
一個更快的做法
有一個性質：相似矩陣的特征多項式相同。
所以可以把這個矩陣相似到一個可以快速求特征多項式的矩陣再求。
怎么相似呢？
\(A\)和\(B\)相似意味着\(A=PBP^{-1}\)
假設當前矩陣是\(A\)。在消元過程中會左乘一個初等矩陣\(P\)，這時在\(A\)的右邊乘上這個初等矩陣的逆矩陣\(P^{-1}\)，得到\(PAP^{-1}\)。
顯然這個矩陣和矩陣\(A\)相似。
但是我們只能把所有\(i\geq j+2\)的部分消成\(0\)，也就是說會消成一個上海森堡矩陣。
接下來就要快速求一個上海森堡矩陣的特征多項式。
上海森堡矩陣長這樣：
\[\begin{pmatrix} a_{1,1}&a_{1,2}&a_{1,3}&\cdots&a_{1,n}\\ a_{2,1}&a_{2,2}&a_{2,3}&\cdots&a_{2,n}\\ 0&a_{3,2}&a_{3,3}&\cdots&a_{3,n}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&a_{n,n} \end{pmatrix} \]
設右下角\(i\times i\)的矩陣的特征多項式為\(f_i\)。
考慮對第一列展開。
那么肯定是連續的消掉若干個第二行后再消掉第一行。
即
\[f_i=a_{i,i}f_{i+1}-a_{i+1,i}a_{i,i+1}f_{i+2}+a_{i+1,i}a_{i+2,i+1}a_{i,i+2}f_{i+3}+\cdots \]
注意到只有主對角線上的元素一次的，其他都是常數的，所以前面的系數的次數最多是\(1\)。
這樣可以從后往前求出所有\(f_i\)。時間復雜度是\(O(n^3)\)。
這樣我們就得到了一個時間復雜度是\(O(n^3)\)的優秀做法。
這個做法有什么用？
可以去卡別人。比如用\(O(n^3+n^2\log m)\)的做法卡\(O(n^3\log m)\)的矩陣快速冪。
代碼
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<utility>
using namespace std;
typedef long long ll;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
const ll p=998244353;
struct poly
{
ll a[510];
int n;
poly()
{
memset(a,0,sizeof a);
n=0;
}
ll &operator [](int x)
{
return a[x];
}
};
int n;
ll a[510][510];
void add(poly &a,poly &b,pll c)
{
a.n=max(a.n,b.n+bool(c.first));
int i;
for(i=0;i<=a.n;i++)
{
a[i]=(a[i]+b[i]*c.second)%p;
if(i)
a[i]=(a[i]+b[i-1]*c.first)%p;
}
}
ll fp(ll a,ll b)
{
ll s=1;
for(;b;b>>=1,a=a*a%p)
if(b&1)
s=s*a%p;
return s;
}
void gao1()
{
int i,j,k;
for(i=1;i<=n;i++)
{
for(j=i+1;j<=n;j++)
if(a[i][j])
break;
if(j>n)
continue;
if(j!=i+1)
{
for(k=i;k<=n;k++)
swap(a[i+1][k],a[j][k]);
for(k=1;k<=n;k++)
swap(a[k][i+1],a[k][j]);
}
ll e=fp(a[i+1][i],p-2);
for(j=i+2;j<=n;j++)
if(a[j][i])
{
ll v=e*a[j][i]%p;
for(k=i;k<=n;k++)
a[j][k]=(a[j][k]-a[i+1][k]*v)%p;
for(k=1;k<=n;k++)
a[k][i+1]=(a[k][i+1]+a[k][j]*v)%p;
}
}
}
pll c[510][510];
poly f[510];
pll operator *(pll a,pll b)
{
return pll((a.first*b.second+a.second*b.first)%p,a.second*b.second%p);
}
pll operator *(pll a,ll b)
{
return pll(a.first*b%p,a.second*b%p);
}
void gao2()
{
f[n+1][0]=1;
int i,j;
for(i=n;i>=1;i--)
{
pll v(0,1);
for(j=i+1;j<=n+1;j++)
{
add(f[i],f[j],v*c[i][j-1]*((j-i)&1?1:-1));
v=v*c[j][j-1];
}
}
}
int main()
{
scanf("%d",&n);
int i,j;
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
scanf("%lld",&a[i][j]);
gao1();
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
{
c[i][j].second=-a[i][j];
if(i==j)
c[i][j].first=1;
}
gao2();
for(i=0;i<=n;i++)
printf("%lld ",(f[1][i]+p)%p);
printf("\n");
return 0;
}