【XSY3154】入门多项式 高斯消元

题目大意

　　给你一个 \(n\times n\)的矩阵 \(A\)，求次数最小且最高次项为 \(1\) 的多项式 \(F(x)\)，满足 \(F(A)=0\)。

　　所有操作都对 \(p\) 取模。

　　\(n\leq 70,n<p\leq 998244353\)

题解

　　显然特征多项式满足条件，但不一定是最优的。

　　设答案为 \(F(x)=\sum_{i\geq 0}f_ix^i\)。

　　那么

\[\begin{cases} f_0{(A^0)}_{1,1}+f_1{(A^1)}_{1,1}+\cdots+f_n{(A^n)}_{1,1}&=0\\ f_0{(A^0)}_{1,2}+f_1{(A^1)}_{1,2}+\cdots+f_n{(A^n)}_{1,2}&=0\\ \vdots\\ f_0{(A^0)}_{n,n}+f_1{(A^1)}_{n,n}+\cdots+f_n{(A^n)}_{n,n}&=0 \end{cases} \]

　　这就是一个方程组，可以通过高斯消元来求解。

　　观察高斯消元的过程。

　　如果在消第 \(i\) 列的时候找不到主元，就说明这个矩阵的前 \(i\) 列不满秩，那么就可以钦定 \(f_{i-1}=1\)，从而得到一组解。

　　否则前 \(i\) 列是满秩的，唯一可能的解为 \(f_0=f_1=\ldots=f_{i-1}=0\)

　　时间复杂度：\(O(n^4)\)

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=80;
int n;
ll 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;
}
struct mat
{
	ll a[N][N];
	mat()
	{
		memset(a,0,sizeof a);
	}
	ll *operator [](int x)
	{
		return a[x];
	}
};
mat operator *(mat a,mat b)
{
	mat c;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
		{
			__int128 s=0;
			for(int k=1;k<=n;k++)
				s+=(ll)a[i][k]*b[k][j];
			c[i][j]=s%p;
		}
	return c;
}
mat a[N];
ll ans[N];
ll c[N*N][N];
int m;
void gao(int x)
{
	ans[x]=1;
	for(int i=1;i<x;i++)
		ans[i]=(-c[i][x]*fp(c[i][i],p-2)%p+p)%p;
	printf("%d\n",x-1);
	for(int i=1;i<=x;i++)
		printf("%lld ",ans[i]);
}
void gao()
{
	for(int i=1;i<=n+1;i++)
	{
		int flag=0;
		for(int j=i;j<=m;j++)
			if(c[j][i])
			{
				flag=j;
				break;
			}
		if(!flag)
		{
			gao(i);
			return;
		}
		if(flag!=i)
		{
			for(int k=i;k<=n+1;k++)
				swap(c[i][k],c[flag][k]);
		}
		ll inv=fp(c[i][i],p-2);
		for(int j=1;j<=m;j++)
			if(j!=i&&c[j][i])
			{
				ll v=c[j][i]*inv%p;
				for(int k=i;k<=n+1;k++)
					c[j][k]=(c[j][k]-v*c[i][k])%p;
			}
	}
}
int main()
{
#ifndef ONLINE_JUDGE
	freopen("a.in","r",stdin);
	freopen("a.out","w",stdout);
#endif
	scanf("%d%lld",&n,&p);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			scanf("%lld",&a[1][i][j]);
	for(int i=1;i<=n;i++)
		a[0][i][i]=1;
	for(int i=2;i<=n;i++)
		a[i]=a[i-1]*a[1];
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
		{
			m++;
			for(int k=0;k<=n;k++)
				c[m][k+1]=a[k][i][j];
		}
	gao();
	return 0;
}