正确的求解代数余子式的方法

前言

在联合省选day2t3中，存在一种使用行列式求导来计算生成树边权和的方法。需要计算出每个位置的代数余子式
常见的做法均是套用$A\times A^{\ast }=|A|\times I$的等式求解，这种做法在模意义下矩阵$rank(A)=n-1$时不能正确得出伴随矩阵
以下给出正确的求解伴随矩阵的方法

代数余子式

给$n$阶方阵$A=(a_{i,j})$，定义$a_{i,j}$的余子式$M_{i,j}$为将$i$行$j$列划去后的行列式，$a_{i,j}$的代数余子式$A_{i,j}$为$(-1{)}^{i+j}{M}_{i,j}$

伴随矩阵

若$A_{i,j}$为代数余子式矩阵，我们定义
$$A^{\ast }=(A_{ij}{)}^{\mathrm{T}}$$
$A^{\ast }$为矩阵$A$的伴随矩阵

我们存在一个结论，即
$$A\times A^{\ast }=A^{\ast }\times A=|A|\times I_n$$
其中$I_n$为$n$阶单位矩阵，$|A|$为矩阵$A$的行列式

计算

若$rank(A)=n$，则$A^{-1}$是良定义的，可以直接用$A^{\ast }=|A|\times A^{-1}$来进行计算。

若$rank(A)\le n-2$，不难发现任意$n-1$阶子式行列式均为$0$，故$A^{\ast }=O$，其中$O$为空矩阵

我们只需要特殊考虑$rank(A)=n-1$的情况

不妨设$A$的列向量组为$c_1,c_2\cdots c_n$，由其线性相关，可以得到存在一组不全为$0$的系数$q_1,q_2\cdots q_n$，满足$\sum q_i{c}_i=0$。不妨假设$q_c\ne 0$

考虑同行两余子式$M_{r,i}$与$M_{r,c}$的关系（在此后的表述中，$M_{r,i}$可能为矩阵，也可能为该矩阵的行列式。视语境而定）。不失一般性地令$i<c$，对于$i>c$的情况是同理推导的

若我们将$M_{r,i}$的第$c-1$列向前置换至第$i$列，同时原本第$[i,c-2]$列向后循环移位。我们能够得到与$M_{r,c}$十分相似的一个矩阵。唯一的不同在于$M_{r,i}$的第$i$列为原本的第$c$列去掉行$r$，而$M_{r,c}$的第$i$列为原本的第$i$列去掉行$r$

我们构造新矩阵$M^{\prime }$，将$M_{r,i}$（循环移位后）的第$i$列乘上$q_c$，将$M_{r,c}$的第$i$列乘上$q_i$，令$M^{\prime }$为这两个矩阵的和。这里两个矩阵的和的计算方式是这样的：将唯一不同的两列求和，其他所有的列照写。
由行列式的性质，不难发现有
$$|M^{\prime }|=(-1{)}^{c-i-1}{q}_c|M_{r,i}|+q_i|M_{r,c}|$$
若我们将$M^{\prime }$的其他列乘上原本对应的列$c$的$q_c$，再加到$M^{\prime }$的第$i$列上，行列式不变。而第$i$列全为$0$，此时可以证明$|M^{\prime }|=0$，因而有
$$\begin{gathered} (-1{)}^{c-i}{q}_c|M_{r,i}|=q_i|M_{r,c}| \\ |M_{r,i}|=\frac{q_i}{q_c}(-1{)}^{c-i}|M_{r,c}| \\ (-1{)}^{r+i}|M_{r,i}|=\frac{q_i}{q_c}(-1{)}^{c-i+r+i}|M_{r,c}| \\ {A}_{r,i}=\frac{q_i}{q_c}{A}_{r,c} \end{gathered}$$
通过对$i>c$的情况加以分析，不难发现，上述结论仍然成立。故有
$$A_{r,i}=\frac{q_i}{q_c}{A}_{r,c}$$
对于行向量组的分析，也是相似的。

对于一组行向量$p=(p_1{p}_2\cdots p_n)$，列向量$q=(q_1{q}_2\cdots q_n)$，我们求出这样的两个向量使得$pA=0,Aq=0$。不妨设$p_r,q_c\ne 0$，就有
$$A_{i,j}=\frac{p_i{q}_j}{p_r{q}_c}{A}_{r,c}$$
求解向量$p,q$，代数余子式$A_{r,c}$均可以在$\Theta (n^3)$的时间内用高斯消元法求解。综上所述，求解一个矩阵的伴随矩阵的时间复杂度为$\Theta (n^3)$

代码

调用Mat::get_adjoint_matrix(a,n)即可求解$n$阶矩阵$a$的伴随矩阵
可以在这里验证你的代码
https://vjudge.net/problem/OpenJ_POJ-C19J

const int MAXN=35;
const int MAXM=152505;
const int mod=998244353;
inline void ad(int &x,int y){x+=y;if(x>=mod)x-=mod;}
inline void dl(int &x,int y){x-=y;if(x<0)x+=mod;}
inline int pow_mod(int a,int b)
{
	int ret=1;
	for(;b;b>>=1,a=1LL*a*a%mod)if(b&1)ret=1LL*ret*a%mod;
	return ret;
}
struct matrix{int m[MAXN][MAXN];matrix(){memset(m,0,sizeof(m));}};
namespace Mat
{
	int A[MAXN][MAXN],B[MAXN][2*MAXN];matrix mat_ret,mat_lin;
	int dat(const matrix &a,int n)
	{
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)A[i][j]=a.m[i][j];
		int ret=1;
		for(int i=1;i<=n;i++)
		{
			int u=i;for(int j=i+1;j<=n;j++)if(A[j][i]){u=j;break;}
			if(u!=i){for(int j=i;j<=n;j++)swap(A[u][j],A[i][j]);ret=1LL*ret*(mod-1)%mod;}
			if(!A[i][i])return 0;int Iv=pow_mod(A[i][i],mod-2);
			for(int j=i+1;j<=n;j++)
			{
				int temp=1LL*A[j][i]*Iv%mod;
				for(int k=i;k<=n;k++)dl(A[j][k],1LL*A[i][k]*temp%mod);
			}ret=1LL*ret*A[i][i]%mod;
		}return ret;
	}
	matrix mat_inv(const matrix &a,int n)
	{
		memset(B,0,sizeof(B));
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)B[i][j]=a.m[i][j];
		for(int i=1;i<=n;i++)B[i][i+n]=1;
		for(int i=1;i<=n;i++)
		{
			int u=i;for(int j=i+1;j<=n;j++)if(B[j][i]){u=j;break;}
			for(int j=1;j<=2*n;j++)swap(B[i][j],B[u][j]);
			if(!B[i][i])continue;int Iv=pow_mod(B[i][i],mod-2);
			for(int j=i+1;j<=n;j++)
			{
				int temp=1LL*B[j][i]*Iv%mod;
				for(int k=i;k<=2*n;k++)dl(B[j][k],1LL*B[i][k]*temp%mod);
			}
		}
		for(int i=n;i>=1;i--)
		{
			int Iv=pow_mod(B[i][i],mod-2);
			for(int j=i;j<=2*n;j++)if(B[i][j])B[i][j]=1LL*B[i][j]*Iv%mod;
			for(int j=i-1;j>=1;j--)if(B[j][i])
			{
				int Z=B[j][i];
				for(int k=i;k<=2*n;k++)if(B[i][k])dl(B[j][k],1LL*B[i][k]*Z%mod);
			}
		}
		memset(mat_ret.m,0,sizeof(mat_ret.m));
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)
			mat_ret.m[i][j]=B[i][j+n];
		return mat_ret;
	}
	int calculate_r(const matrix &a,int n)
	{
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)A[i][j]=a.m[i][j];int ret=0;
		for(int i=1;i<=n;i++)
		{
			if(!A[i][i])
			{
				int u=-1;
				for(int j=i+1;j<=n;j++)if(A[j][i]){u=j;break;}
				if(u==-1)continue;
				for(int j=i;j<=n;j++)swap(A[i][j],A[u][j]);
			}++ret;int Iv=pow_mod(A[i][i],mod-2);
			for(int j=i+1;j<=n;j++)
			{
				int temp=1LL*A[j][i]*Iv%mod;
				for(int k=i;k<=n;k++)dl(A[j][k],1LL*A[i][k]*temp%mod);
			}
		}return ret;
	}
	matrix get_adjoint_matrix_full(const matrix &a,int n)
	{
		matrix Iv=mat_inv(a,n);int Z=dat(a,n);
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)
			mat_ret.m[j][i]=1LL*Iv.m[i][j]*Z%mod;
		return mat_ret;
	}
	
	
	int us[MAXN][MAXN],Inv[MAXN];
	VI insert(VI Z,int n,int id)
	{
		VI bel(n+1,0);bel[id]=1;bel[0]=-1;
		for(int i=n;i>=1;i--)if(Z[i])
		{
			if(!A[i][i])
			{
				for(int j=i;j>=1;j--)A[i][j]=Z[j];Inv[i]=pow_mod(A[i][i],mod-2);
				for(int j=1;j<=n;j++)us[i][j]=bel[j];
				return bel;
			}
			int temp=1LL*Inv[i]*Z[i]%mod;
			for(int j=i;j>=1;j--)dl(Z[j],1LL*A[i][j]*temp%mod);
			for(int j=1;j<=n;j++)dl(bel[j],1LL*us[i][j]*temp%mod);
		}bel[0]=0;return bel;
	}
	VI get_G(const matrix &a,int n,int opt)
	{
		VI ret;
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)A[i][j]=a.m[i][j];
		if(opt)for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)B[j][i]=A[i][j];
		else for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)B[i][j]=A[i][j];
		
		memset(A,0,sizeof(A));memset(us,0,sizeof(us));
		for(int i=1;i<=n;i++)
		{
			VI Z(n+1,0);
			for(int j=1;j<=n;j++)Z[j]=B[j][i];
			ret=insert(Z,n,i);
			if(ret[0]!=-1)return ret;
		}assert(false);
	}
	matrix get_adjoint_matrix_one(const matrix &a,int n)
	{
		VI Q=get_G(a,n,0);
		VI P=get_G(a,n,1);
		int c=-1,r=-1;
		for(int i=1;i<=n;i++)if(Q[i]){c=i;break;}assert(c!=-1);
		for(int i=1;i<=n;i++)if(P[i]){r=i;break;}assert(r!=-1);
		
		for(int i=1;i<=n;i++)if(i!=r)for(int j=1;j<=n;j++)if(j!=c)
			mat_lin.m[i-(i>r)][j-(j>c)]=a.m[i][j];
		int Ivq=pow_mod(Q[c],mod-2),Ivp=pow_mod(P[r],mod-2);
		int Iv=1LL*Ivq*Ivp%mod;
		mat_ret.m[r][c]=dat(mat_lin,n-1);
		if((r+c)&1)mat_ret.m[r][c]=mod-mat_ret.m[r][c];
		for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)
		{
			int val=1LL*P[i]*Q[j]%mod*Iv%mod*mat_ret.m[r][c]%mod;
			mat_ret.m[i][j]=val;
		}return mat_ret;
	}
	
	matrix get_adjoint_matrix(const matrix &a,int n)
	{
		int Z=calculate_r(a,n);
		if(Z==n)return get_adjoint_matrix_full(a,n);
		else if(Z==n-1)return get_adjoint_matrix_one(a,n);
		else
		{
			memset(mat_ret.m,0,sizeof(mat_ret.m));
			return mat_ret;
		}
	}
}