P5027 - Barracuda 题解

一道有趣的线性代数题。首先可以枚举错误方程，每次高消跑一遍，\(\mathrm O\!\left(n^4\right)\) 其中 \(n\leq 100\)，卡着时限过就没有意思了（但好像大多数人都是这么做的？这可能是难度只有蓝的原因）。这篇题解讲的是 \(\mathrm O\!\left(n^3\right)\) 的做法。

考虑解线性方程组的另一个解法：克莱姆法则。对于系数矩阵为方阵的线性方程组 \(A\pmb x=\pmb b\)，我们有 \(x_i=\dfrac{|A_i(\pmb b)|}{|A|}\)。直接按照这个做需要求 \(\mathrm O(n)\) 次行列式，这就四方了。再者，系数矩阵可能没满秩，也就是 \(|A|=0\)，那就做不了了。而且行列式的值可能会很大很大，爆 double（虽说消元求行列式仅包含乘法，可以取对数，但那还是麻烦啊）。不难发现普通情况下克莱姆法则解方程组漏洞百出，但是我们还是能修补复杂度的漏洞。

数学书上讲的东西总是很形式化的，那 \(|A_i(\pmb b)|\) 这东西看起来就没有什么快速求的好思路。但你把 \(\pmb b\) 拖拽到最后一列（也就是若干次列对换）就有灵感了：要求的 \(|A_i(\pmb b)|\) 以及 \(|A|\) 这 \(n+1\) 个行列式刚好是增广矩阵分别去掉每一列所得矩阵的行列式乘上 \(\pm1\)（\(\pm1\) 是 trivial 的，不去管它）。去掉一列的行列式，想到了什么？余子式！但是余子式要求去掉一行一列，现在一列有了，对于一行容易想到在增广矩阵下面随便补一行（补完之后恰好又是方阵！），这样要求的就是 \((n+1,1\sim n+1)\) 余子式。

要计算一个方阵每处的余子式，只要计算方阵的伴随矩阵，那么有 \(\mathrm{adj}A=|A|A^{-1}\)，其中行列式和逆矩阵都是好求的。但同时有个条件，那就是方阵行列式非零。但最后一行是我们自己补的，在增广矩阵行满秩的前提下可以刻意构造一行使得 \(n+1\) 阶大方阵满秩。然后发现，克莱姆法则中两个余子式之比，恰好把伴随矩阵公式中的 \(|A|\) 约掉了，只剩逆矩阵元素之比！那就不用担心爆 double 了。总结一下，克莱姆法则解线性方程组有几个漏洞：只能解系数矩阵为方阵的方程组，如果系数矩阵不满秩还判断不了无解还是无限解，无限解还没法求通解，还得动脑子在增广矩阵下面新加一行使得 \(n+1\) 阶方阵满秩。

虽然克莱姆法则很不好，但是我们还是能发现一个优点：它把 \(\mathrm O\!\left(n^2\right)\) 处的余子式都求出来了，但我们只用到了 \(\mathrm O(n)\) 处，这说明这个算法很可扩展。回来看到该题，一开始列出一个 \(n+1\) 阶增广矩阵，每次要去掉一行得到要求的 \(n\times(n+1)\) 增广矩阵。等等，去掉一行？那刚好可以让这一行充当上述算法的新加的一行呀！这样一来不难发现，要解 \(\mathrm O(n)\) 次线性方程组，每次 \(\mathrm O(n)\) 个余子式，一共要 \(\mathrm O\!\left(n^2\right)\) 个余子式，而它们恰好不重不漏地分布在伴随矩阵里！

再看看该题有没有撞上克莱姆法则的几个缺点。系数矩阵为方阵没错了，如果系数矩阵不满秩也根本不用去管它。但是并没有保证 \(n+1\) 阶增广矩阵满秩！事实上，我们可以证明若不满秩，答案一定是 \(\texttt{illegal}\)。如果有自由变量，那把 \(n+1\) 个方程合起来了都有自由变量了，再去掉一个必定没有唯一解。如果没有自由变量且秩等于 \(n\)，那么 \(n+1\) 个方程有唯一解，又分成两类：如果解不合法，那去掉一个方程要么无限解，要么是这个不合法的唯一解，肯定不行；如果解合法，我们努力找到去掉两个方程的时候都得到唯一解，满足这个事情当且仅当去掉的是「冗余的」，即可以被其它行向量线性表出的行，由于 \(n+1\) 行没有满秩，必定有至少一个，而「说明 / 提示」最后一行说 \(m\geq1\)，也就是不可能有零行，那线性表出式的 RHS 中必定有至少一项，那移个项即可得到另一个满足要求的行。

code

#include<bits/stdc++.h>
using namespace std;
#define pb push_back
#define mp make_pair
#define X first
#define Y second
const double eps=1e-6;
const int N=110;
int n;
double a[N][2*N];
void times(int x,double v){
	for(int i=1;i<=2*n+2;i++)a[x][i]*=v;
}
void swp(int x,int y){
	for(int i=1;i<=2*n+2;i++)swap(a[x][i],a[y][i]);
}
void tadd(int x,double v,int y){
	for(int i=1;i<=2*n+2;i++)a[y][i]+=a[x][i]*v;
}
void prt(){
	for(int i=1;i<=n+1;i++)for(int j=1;j<=2*n+2;j++)cout<<a[i][j]<<" \n"[j==2*n+2];
	puts("--------------------");
}
void gauss(){
	for(int i=1;;i++){
		int row=0,col=0;
		for(int j=1;j<=n+1;j++){
			for(int k=i;k<=n+1;k++)if(abs(a[k][j])>eps){row=k,col=j;break;}
			if(row)break;
		}
		if(!row)break;
		swp(i,row);
		for(int j=1;j<=n+1;j++)if(i!=j)tadd(i,-a[j][col]/a[i][col],j);
		times(i,1/a[i][col]);
//		prt();
	}
}
double M(int x,int y){return a[y][x+n+1]*(x+y+(y!=n+1)*(n-y&1)&1?-1:1);}
int main(){
	cin>>n;
	for(int i=1;i<=n+1;i++){
		int m;
		cin>>m;
		while(m--){
			int x;
			scanf("%d",&x);
			a[i][x]=1;
		}
		int x;
		scanf("%d",&x);
		a[i][n+1]=x;
	}
	for(int i=1;i<=n+1;i++)a[i][i+n+1]=1;
//	prt();
	gauss();
	for(int i=1;i<=n+1;i++)if(abs(a[i][i]-1)>eps)return puts("illegal"),0;
	vector<int> legal;
	for(int i=1;i<=n+1;i++){
		if(abs(M(i,n+1))<=eps)continue;
		bool flg=true;multiset<int> st;
		for(int j=1;j<=n;j++){
			double x=M(i,j)/M(i,n+1);
			flg&=x>.5&&abs(round(x)-x)<=eps;
			st.insert(round(x)+.5);
		}
		flg&=st.count(*--st.end())==1;
		if(flg)legal.pb(i);
	}
//	cout<<legal.size()<<"!\n";
	if(legal.size()!=1)return puts("illegal"),0;
	pair<int,int> mx(0,0);
	for(int i=legal[0];i<=legal[0];i++){
		for(int j=1;j<=n;j++){
			double x=M(i,j)/M(i,n+1);
			mx=max(mx,mp(int(round(x)+.5),j));
		}
	}
	cout<<mx.Y;
	return 0;
}