一个比较慢的做法
首先你要知道矩阵的特征多项式是什么。
直接消元就可以了。
时间复杂度:(O(n^5))或(O(n^4))。
一个稍微快一点的做法
观察到特征多项式的次数是(n)。
我们就可以插值。
具体来说,先求出当(x=0ldots n)时特征多项式对应的点值,然后直接用拉格朗日插值插出来。
时间复杂度:(O(n^4))
一个更快的做法
有一个性质:相似矩阵的特征多项式相同。
所以可以把这个矩阵相似到一个可以快速求特征多项式的矩阵再求。
怎么相似呢?
(A)和(B)相似意味着(A=PBP^{-1})
假设当前矩阵是(A)。在消元过程中会左乘一个初等矩阵(P),这时在(A)的右边乘上这个初等矩阵的逆矩阵(P^{-1}),得到(PAP^{-1})。
显然这个矩阵和矩阵(A)相似。
但是我们只能把所有(igeq j+2)的部分消成(0),也就是说会消成一个上海森堡矩阵。
接下来就要快速求一个上海森堡矩阵的特征多项式。
上海森堡矩阵长这样:
[egin{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 imes 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^2log m))的做法卡(O(n^3log 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("
");
return 0;
}