P5027

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

考虑解线性方程组的另一个解法：克莱姆法则。对于系数矩阵为方阵的线性方程组 (Apmb 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,1sim n+1)) 余子式。

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

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

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