求矩阵的特征多项式

一个比较慢的做法

　　首先你要知道矩阵的特征多项式是什么。

　　直接消元就可以了。

　　时间复杂度：\(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;
}