代数余子式的求法

发表于 分类于

题意

给定一个 $n$ 阶方阵 $A$:

$A=\begin{bmatrix}
a_{1,1} & a_{1,2} & \cdots & a_{1,n} \
a_{2,1} & a_{2,2} & \cdots & a_{2,n}\
\vdots & \vdots & \ddots & \vdots \
a_{n,1} & a_{n,2} & \cdots & a_{n,n}
\end{bmatrix}$

记 $M_{i,j}$ 为矩阵 $A$ 去掉第 $i$ 行和第 $j$ 列后的矩阵或其行列式，$A_{i,j}=(-1)^{i+j}M_{i,j}$，即代数余子式。

求所有的代数余子式 $A_{i,j}$。

$n\le 500$。

做法

记 $r(A)$ 表示矩阵 $A$ 的秩。

一些思路

一个显然的做法，暴力计算每个代数余子式，时间复杂度 $O(n^5)$。

考虑线性代数里的知识：若 $r(A)=n$，则有 $A ^ {-1} = \frac { A^* }{|A|} $，其中： $A^*=\begin{bmatrix}
A_{1,1} & A_{2,1} & \cdots & A_{n,1} \
A_{1,2} & A_{2,2} & \cdots & A_{n,2}\
\vdots & \vdots & \ddots & \vdots \
A_{1,n} & A_{2,n} & \cdots & A_{n,n}
\end{bmatrix}$ ，即 $A$ 的伴随矩阵。

那么当矩阵行列式不为 $0$ 时，通过一次求逆和求行列式即可求出所有的代数余子式。

若 $r(A)\neq n$ 呢？

如果 $r(A)\le n-2$，则说明，任何一个 $n-1$ 阶余子式均为 $0$，于是 $A^*$ 为 $0$ 矩阵。

因此，剩下的只需考虑 $r(A)=n-1$ 时的情况。

$r(A)=n-1$ 时

因 $A$ 不为满秩矩阵，因此行向量是线性相关的。

即对于向量 $v_i=(a_{i,1},a_{i,2},\cdots,a_{i,n})$，存在不全为 $0$ 的数 $p_i$ ，使得 $\sum_{i=1}^mp_iv_i=0$。找到一个不为 $0$ 的 $p_r$。

考虑矩阵：$M_{i,c}$ 和 $M_{r,c}$。其中 $i<r$。

我们将 $M_{i,c}$ 中第 $r-1$ 行插到第 $i-1$ 行的后面去后的矩阵称为 $M ‘ _{i,c}$，那么 $|M ‘ {i,c}|=(-1)^{r-i-1}|M{i,c}|$。

显然，这个 $M’_ {i,c }$ 与 $M_{r,c}$ 只有在第 $i$ 行上是不同的。

设矩阵 $M’$ 为 $M ‘ {i,c}$ 和 $M{r,c}$ 的矩阵中第 $i$ 行的数分别乘上 $p_r$ 或 $p_i$ 相加后，其他行不变的结果。

由行列式性质：$|M ‘ |=p_r(-1)^{r-i-1}|M_{i,c}|+p_i|M_{r,c}|$。

将 $M’$ 中其他行 $j$ 乘上 $p_j$ 后，加到第 $i$ 行中去，则根据 $\sum p_iv_i=0$，第 $i$ 行的结果全为 $0$。因此，$|M’|=0$。

因此：$p_r(-1)^{r-i}|M_{i,c}|=p_i|M_{r,c}|$。

化简可得：$(-1)^{i+c}|M_{i,c}|=\frac{p_i}{p_r}(-1)^{r+c}|M_{r,c}|$。

即：$A_{i,c}=\frac{p_i}{p_r}A_{r,c}$。

当 $i>r$ 时，同理可得上述结果。

当为列向量时，设 $\sum q_{i}v_i=0$。那么同理有：$A_{r,j}=\frac{q_j}{q_c}A_{r,c}$。

在第一个式子中取 $c=j$，联立两式可得：$A_{i,j}=\frac{p_iq_j}{p_rq_c} A_{r,c}$。

于是，只需求出 $p,q$ 和一个 $A_{r,c}$ ，即可求出所有的代数余子式了。

求 $p,q$ 相当于求解齐次线性方程组，$A_{r,c}$ 相当于求一代数余子式。时间复杂度均为 $O(n^3)$。

程序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
#include <bits/stdc++.h>
using namespace std;
#define FO(x) freopen(#x".in","r",stdin),freopen(#x".out","w",stdout)
#define fo(i,j,k) for(int i=(j),end_i=(k);i<=end_i;i++)
#define ff(i,j,k) for(int i=(j),end_i=(k);i< end_i;i++)
#define fd(i,j,k) for(int i=(j),end_i=(k);i>=end_i;i--)
#define DEBUG(x) cerr<<#x<<"="<<x<<endl
#define all(x) (x).begin(),(x).end()
#define cle(x) memset(x,0,sizeof(x))
#define lowbit(x) ((x)&-(x))
#define ll long long
#define ull unsigned ll
#define db double
#define lb long db
#define pb push_back
#define mp make_pair
#define fi first
#define se second
inline int read()
{
	int x=0; char ch=getchar(); bool f=0;
	for(;ch<'0'||ch>'9';ch=getchar()) if(ch=='-') f=1;
	for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<3)+(x<<1)+(ch^48);
	return f?-x:x;
}
#define CASET fo(___,1,read())
const ll mod=998244353;
inline ll Add(ll x,ll y){x+=y; return (x<mod)?x:x-mod;}
inline ll Dec(ll x,ll y){x-=y; return (x<0)?x+mod:x;}
inline ll Mul(ll x,ll y){return x*y%mod;}
inline ll Pow(ll x,ll y)
{
	y%=(mod-1);ll ans=1;for(;y;y>>=1,x=x*x%mod)if(y&1) ans=ans*x%mod;
	return ans;
}
const int N=502;
struct matrix{
	ll a[N][N];
	matrix(){memset(a,0,sizeof(a));}
	void clear(){memset(a,0,sizeof(a));}
};
namespace Mat{
	static ll a[N][N],b[N][N+N],c[N][N],d[N][N];
	static matrix B,C;
	inline ll det(const matrix &A,int n)
	{
		fo(i,1,n) fo(j,1,n) a[i][j]=A.a[i][j];
		ll d=1,iv,tmp;
		fo(i,1,n)
		{
			int k=i;
			fo(j,i+1,n) if(a[j][i]) {k=j; break;}
			if(k!=i) {fo(j,i,n) swap(a[k][j],a[i][j]); d=(mod-d)%mod;}
			if(!a[i][i]) return 0;
			iv=Pow(a[i][i],mod-2);
			fo(j,i+1,n)
			{
				tmp=a[j][i]*iv%mod;
				fo(k,i,n) a[j][k]=Dec(a[j][k],a[i][k]*tmp%mod);
			}
			d=d*a[i][i]%mod;
		}
		return d;
	}
	inline matrix inv(const matrix &A,int n)
	{
		fo(i,1,n) fo(j,1,n) b[i][j+n]=0,b[i][j]=A.a[i][j];
		fo(i,1,n) b[i][i+n]=1;
		ll iv,tmp;
		fo(i,1,n)
		{
			int k=n+1;
			fo(j,i,n) if(b[j][i]) {k=j; break;}
			if(k==n+1) continue;
			if(k!=i) fo(j,1,n+n) swap(b[i][j],b[k][j]);
			iv=Pow(b[i][i],mod-2);
			fo(j,i+1,n)
			{
				tmp=b[j][i]*iv%mod;
				fo(k,i,n+n) b[j][k]=Dec(b[j][k],b[i][k]*tmp%mod);
			}
		}
		fd(i,n,1)
		{
			iv=Pow(b[i][i],mod-2);
			fo(j,i,n+n) b[i][j]=b[i][j]*iv%mod;
			fd(j,i-1,1)
				if(b[j][i])
				{
					tmp=b[j][i];
					fo(k,i,n+n) b[j][k]=Dec(b[j][k],b[i][k]*tmp%mod);
				}
		}
		B.clear();
		fo(i,1,n) fo(j,1,n) B.a[i][j]=b[i][j+n];
		return B;
	}
	int r(const matrix A,int n)
	{
		fo(i,1,n) fo(j,1,n) a[i][j]=A.a[i][j];
		int d=0;
		ll iv,tmp;
		fo(i,1,n)
		{
			int k=n+1;
			fo(j,i,n) if(a[j][i]) {k=j; break;}
			if(k==n+1) continue;
			d++;
			if(k!=i) fo(j,i,n) swap(a[i][j],a[k][j]);
			iv=Pow(a[i][i],mod-2);
			fo(j,i+1,n)
			{
				tmp=a[j][i]*iv%mod;
				fo(k,i,n) a[j][k]=Dec(a[j][k],tmp*a[i][k]%mod);
			}
		}
		return d;
	}
	static ll v[N],w[N];
	inline bool ins(ll *v,int n,int id,ll *ans)
	{
		fo(i,1,n) w[i]=0;
		w[id]=1;
		ll tmp;
		fd(i,n,1)
			if(v[i])
			{
				if(!c[i][i])
				{
					fd(j,i,1) c[i][j]=v[j];
					fo(j,1,n) d[i][j]=w[j];
					return 0;
				}
				tmp=Pow(c[i][i],mod-2)*v[i]%mod;
				fd(j,i,1) v[j]=Dec(v[j],c[i][j]*tmp%mod);
				fo(j,1,n) w[j]=Dec(w[j],d[i][j]*tmp%mod);
			}
		fo(i,1,n) ans[i]=w[i];
		return 1;
	}
	inline void get_G(const matrix &A,int n,ll *p)
	{
		fo(i,1,n) fo(j,1,n) a[i][j]=A.a[i][j];
		memset(c,0,sizeof(c)); memset(d,0,sizeof(d));
		fo(i,1,n)
		{
			fo(j,1,n) v[j]=a[j][i];
			if(ins(v,n,i,p)) return;
		}
	}
	inline matrix solve(const matrix &A,int n)
	{
		int rank=r(A,n);
		if(rank == n)
		{
			ll d=det(A,n);
			B=inv(A,n);
			fo(i,1,n) fo(j,1,n) C.a[i][j]=B.a[j][i]*d%mod;
			return C;
		}
		else if(rank <= n-2)
		{
			C.clear();
			return C;
		}
		else
		{
			static ll p[N],q[N];
			get_G(A,n,q);
			fo(i,1,n) fo(j,1,n) B.a[j][i]=A.a[i][j];
			get_G(B,n,p);
			int c=0,r=0;
			fo(i,1,n) if(q[i]) {c=i; break;}
			fo(i,1,n) if(p[i]) {r=i; break;}
			fo(i,1,n)
				if(i!=r)
					fo(j,1,n)
						if(j!=c)
							B.a[i-(i>r)][j-(j>c)]=A.a[i][j];
			ll d=det(B,n-1);
			C.a[r][c]=((r+c)%2==1)?(mod-d)%mod:d;
			ll iv=Pow(q[c]*p[r]%mod,mod-2);
			fo(i,1,n)
				fo(j,1,n)
					C.a[i][j]=C.a[r][c]*iv%mod*p[i]%mod*q[j]%mod;
			return C;
		}
	}
}
using Mat::solve;

matrix A,B;
int n;
int main()
{
	CASET
	{
		n=read();
		fo(i,1,n) fo(j,1,n) A.a[i][j]=read();
		B=solve(A,n);
		fo(i,1,n)
		{
			fo(j,1,n) printf("%lld ",((i+j)%2==1)?(mod-B.a[i][j])%mod:B.a[i][j]);
			printf("\n");
		}
		printf("\n");
	}
	return 0;
}