【学习笔记】求矩阵的特征多项式

先膜拜一波神仙yww

给定一个矩阵(没有任何特殊性质)，如何求它的特征多项式？

算法一

直接把(lambda)代入((n+1))个点值，求完行列式之后插值即可。
时间复杂度(O(n^4))

算法二

下面介绍一个更快的做法。
定义 对于矩阵(m A,m B), 若存在可逆矩阵(mPhi)满足(m A=mPhi^{-1}m BmPhi), 则称(m A)和(m B)为相似矩阵，(mPhi)称为桥接矩阵。
定理 相似矩阵的特征多项式相同。
证明: 设(m A)和(m B)为相似矩阵，(f(lambda))和(g(lambda))分别为其特征多项式，则(f(lambda)=|lambda m I-m A|=|lambda m Phi^{-1}mPhi-mPhi^{-1}m BmPhi|=|m Phi^{-1}(lambdam I-m B)mPhi|=|mPhi^{-1}||lambda m I-m B| |mPhi|=|mPhi^{-1}mPhi|g(lambda)=g(lambda)). 证毕。
于是我们可以考虑把没有特殊性质的矩阵化成有特殊性质的与之相似的矩阵来加速运算。
我们知道，对矩阵进行初等行(列)变换相当于将其左(右)乘一个可逆的初等矩阵，那么我们根据相似矩阵的定义，对(m A)进行初等行变换左乘初等矩阵的同时给它右乘同一初等矩阵的逆，不就构造出相似矩阵了吗？
简单推导可知，若左乘某一初等矩阵等价于把原矩阵第(i)行乘以(x)加到第(j)行上，则右乘该初等矩阵的逆矩阵等价于把原矩阵第(j)行乘以(-x)加到第(i)行上。

现在的问题是，我们对这个没有特殊性质的矩阵左乘初等矩阵的同时右乘它的逆，能把它变成什么有特殊性质的矩阵呢？
答案是，上海森堡矩阵。
定义 若一(n imes n)矩阵满足对于任意(1le j<j+1<ile n)的((i,j)), 矩阵第(i)行第(j)列值均为(0)，则其为上海森堡矩阵。
至于为什么是这个答案，首先变成上三角矩阵显然不现实，因为用第(i)行来消第(j)行 ((j>i))的第(i)列时对第(i)行和第(j)行进行了初等行变换，那这对应的初等列变换势必会影响第(j)列和第(i)列。那么我们可以考虑在消第(i)列的时候，用第((i+1))行来消第(j)行 ((j>i+1))的第(i)列，这样进行的初等行变换是第((i+1))行和第(j)行之间的，影响到的是第((i+1))列和第(j)列，而不会影响第(i)列了。

现在问题转化为，如何快速求上海森堡矩阵的特征多项式？
这就非常简单了，最后一行只有两列有值，那么枚举选的是哪一列，如果选第(n)列，那么直接就是((lambda-a_{n,n}))乘以左上角((n-1) imes (n-1))的子矩阵的特征多项式；若选第((n-1))列，那么第((n-1))行选((n-2))列，((n-2))行选((n-3))列……一直到某一行为止。假设这一行是第((j+1))行选第(j)列 ((0<j<n))，那么第(j)行一定选第(n)列，前((j-1))行选的就是前((j-1))列了，而(a_{n-1,n-2},a_{n-2,n-3},...,a_{j+1,j},a_{j,n})都是常数(都不在对角线上)，所以直接递推即可。

时间复杂度(O(n^3)).

代码

(每一项对(998244353)取模, (nle 500))

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cassert>
#include<iostream>
#define llong long long
using namespace std;

inline int read()
{
	int x=0; bool f=1; char c=getchar();
	for(;!isdigit(c);c=getchar()) if(c=='-') f=0;
	for(; isdigit(c);c=getchar()) x=(x<<3)+(x<<1)+(c^'0');
	if(f) return x;
	return -x;
}

const int N = 500;
const int P = 998244353;
llong a[N+3][N+3];
llong c[N+3][N+3];
int n;

llong quickpow(llong x,llong y)
{
	llong cur = x,ret = 1ll;
	for(int i=0; y; i++)
	{
		if(y&(1ll<<i)) {y-=(1ll<<i); ret = ret*cur%P;}
		cur = cur*cur%P;
	}
	return ret;
}
llong mulinv(llong x) {return quickpow(x,P-2);}

void gauss()
{
	for(int i=1; i<=n; i++)
	{
		if(a[i+1][i]==0)
		{
			bool found = false;
			for(int j=i+2; j<=n; j++)
			{
				if(a[j][i]!=0)
				{
					for(int k=i; k<=n; k++) swap(a[i+1][k],a[j][k]);
					for(int k=1; k<=n; k++) swap(a[k][i+1],a[k][j]);
					found = false; break;
				}
			}
			if(found) {continue;}
		}
		for(int j=i+2; j<=n; j++)
		{
			llong coe = P-a[j][i]*mulinv(a[i+1][i])%P;
			for(int k=i; k<=n; k++) a[j][k] = (a[j][k]+coe*a[i+1][k])%P;
			for(int k=1; k<=n; k++) a[k][i+1] = (a[k][i+1]-coe*a[k][j]%P+P)%P;
		}
	}
}

void charpoly()
{
	c[0][0] = 1ll;
	for(int i=1; i<=n; i++)
	{
		for(int j=0; j<=i; j++)
		{
			c[i][j] = (c[i-1][j-1]-a[i][i]*c[i-1][j]%P+P)%P;
		}
		llong coe = P-1,cur = P-a[i][i-1];
		for(int j=i-2; j>=0; j--)
		{
			llong tmp = cur*(P-a[j+1][i])%P;
			tmp = coe*tmp%P;
			for(int k=0; k<=j; k++)
			{
				c[i][k] = (c[i][k]+c[j][k]*tmp)%P;
			}
			cur = cur*(P-a[j+1][j])%P;
			coe = P-coe;
		}
		for(int k=0; k<=i; k++) c[i][k] %= P;
//		printf("%d: ",i); for(int j=0; j<=i; j++) printf("%lld ",c[i][j]); puts("");
	}
}

int main()
{
	scanf("%lld",&n);
	for(int i=1; i<=n; i++) for(int j=1; j<=n; j++) scanf("%lld",&a[i][j]);
	gauss();
//	for(int i=1; i<=n; i++) {for(int j=1; j<=n; j++) printf("%lld ",a[i][j]); puts("");}
	charpoly();
	for(int i=0; i<=n; i++) printf("%lld ",c[n][i]); puts("");
	return 0;
}