bzoj3168-钙铁锌硒维生素

题目

这道题的题意理解很重要，直接写原题了。
小林把人体需要的营养分成了\(n\)种，他准备了2套厨师机器人，一套厨师机器人有\(n\)个，每个厨师机器人只会做一道菜，这道菜一斤能提供第\(i\)种营养\(x_i\)微克。想要吃这道菜的时候，只要输入一个数，就能吃到对应数量的这道菜了。为防止摄入过量对身体造成的伤害，每个机器人还有防过量摄入药，只要输入一个数，就能生成一定剂量的药，吃了这些药，就能减少相当于食用对应数目的这道菜提供的营养。
小林之所以准备2套厨师机器人，是因为某一个厨师机器人可能坏掉，要是影响了银河队选手的身体，就不好了。因此，第2套厨师机器人被用来做第1套的备用。小林需要为每一个第1套厨师机器人选一个第2套厨师机器人作备份，使得当这个机器人坏掉时，用备份顶替，整套厨师机器人仍然能搭配出任何营养需求，而且，每个第2套厨师机器人只能当一个第1套厨师机器人的备份。求一种备份方案，使得字典序最小（按所备份的机器人编号从小到大输出）。

分析

首先，题中输入的“数”可以是小数，所以可以任意搭配。那“搭配出任何营养需求”是什么意思呢？也就是说，我可以通过改变每种菜的数量来获取任意组合的营养，即：

\[\begin{aligned} y_i=\sum _ {j=1}^n x_j *A _{ji} \\\ \end{aligned} \]

我们把他们都写成矩阵的形式：

\[\begin{aligned} x*A=y \end{aligned} \]

所以判断一个矩阵能否组合出任意营养搭配，就是判断一个矩阵是否能通过消元变成对角阵。这个转化十分关键，其实就是看懂题意。
这是一个一一对应的问题，每个备用机器人只能换一个第一套机器人，然后就想到，这和二分图的匹配模型很相似。以后见到两种物品，一一对应或在一个集合中选择给另一个集合匹配（例如gdkoi2017的演员那题）类似的题目，二分图应该是一种可行的解法。
现在的问题是，我们要求出每一个\(B\_j\)有那些\(A\_i\)可以替换。这里用到了一个很显然的结论。
\(A\)可以通过消元变成对角阵。若存在行向量\(\lambda\)使得\(\lambda\*A=b\)，那么如果\(\lambda \_ i\ne 0\)，那么把行\(A\_i\)替换为\(b\)后，新的矩阵也能消元得出对角阵。这是因为，一行不能被替换，当且仅当替换后这一行会被消成全0，即替换后这一行可以用别的行线性表示（每行乘以一个不为0的系数，行相加），于是就无法消元得出满对角阵。所以我们只能选择\(\lambda\)不为0的行替换。
由于\(\lambda\)有\(n\)个，不如我们把它们全部一起求出来。矩阵\(C\)表示所有的\(\lambda\)，即\(C\_j\)表示\(B\_j\)的替换方案：

\[C*A=B \\ C*A*{A^{-1}}=B*{A^{-1}} \\ C=B*A^{-1} \\ \]

这样我们就求出了所有的\(\lambda\)，即求出了二分图的邻接矩阵，\(C_{ij}\ne 0\)表示左边的\(j\)和右边的\(i\)有连边。这样就可以求二分图匹配了。现在的问题在于，如何求一个矩阵的逆矩阵。
例如原矩阵为：

\[\begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix} \]

我们给出增广炬阵，并把左边消元成单位阵：

\[\left[ \begin{array}{cc|cc} 1 & 2 & 1 & 0 \\ 3 & 4 & 0 & 1 \\ \end{array} \right] \\ \left[ \begin{array}{cc|cc} 1 & 2 & 1 & 0 \\ 0 & -2 & -3 & 1 \\ \end{array} \right] \\ \left[ \begin{array}{cc|cc} 1 & 0 & 1 & 1 \\ 0 & -2 & -3 & 1 \\ \end{array} \right] \\ \left[ \begin{array}{cc|cc} 1 & 0 & -2 & 1 \\ 0 & 1 & \frac{3}{2} & -\frac{1}{2} \\ \end{array} \right] \\ \]

这时，右边得到的矩阵就是原矩阵的逆矩阵。所以矩阵求逆非常简单，只要对左边进行消元，变成对角阵即可。如果一个矩阵不可逆，那么左边无法形成对角阵。
最后一个问题，如何解决字典序最小。通过匈牙利算法的过程可以看出，我们只能保证得到的最后一个匹配的字典序最小，前面的无法保证。因此，我们可以通过一种调整来在匹配数量不变的情况下从前往后保证最小字典序。从小到大，对于左边的一个点，先把原来的匹配边记录下来并断开，再从小到大扫右边的点，跟之前一样找增广路。由于之前断开了匹配边，匹配数量减一；找到增广路后，匹配数量加一，故答案不变。只要找到了增广路，我们就把当前的左边点和右边点删掉（标记），把它们锁定，这样以后不会再更改到它们。这样既可求出最大匹配的最小字典序。其实最小字典序就是一个贪心的东西。

代码

用double直接消元比他们的大数取模快几倍，不过大数取模也是很好的方法。通过给定一个永远不会超过的大质数，可以方便地求乘法逆元，解决除法问题，不过每次多一个log。

#include<cstdio>
#include<algorithm>
#include<cctype>
#include<cstring>
#include<cmath>
using namespace std;
int read() {
	int x=0,f=1;
	char c=getchar();
	for (;!isdigit(c);c=getchar()) if (c=='-') f=-1;
	for (;isdigit(c);c=getchar()) x=x*10+c-'0';
	return x*f;
}
const int maxn=305;
const double eps=1e-9;
int b[maxn][maxn],match[maxn];
double a[maxn][maxn],c[maxn][maxn];
double h[maxn][maxn<<1];
int n;
bool alr[maxn];
int gcd(int x,int y) {
	if (x<y) swap(x,y);
	return y?gcd(y,x%y):x;
}
bool elm() {
	for (int i=1;i<n;++i) {
		if (!h[i][i]) {
			for (int j=i+1;j<=n;++j) if (fabs(h[j][i])>eps) {
				for (int k=1;k<=n+n;++k) swap(h[i][k],h[j][k]);
				break;
			}
		}
		for (int j=i+1;j<=n;++j) if (h[j][i]) {
			double tmp=h[j][i]/h[i][i];
			for (int k=1;k<=n+n;++k) h[j][k]-=h[i][k]*tmp,h[j][k]=(fabs(h[j][k])<eps?0:h[j][k]);
		}		
	}
	for (int i=n;i>1;--i) {
		if (!h[i][i]) {
			for (int j=i-1;j>0;--j) if (fabs(h[j][i])<eps) {
				for (int k=1;k<=n+n;++k) swap(h[i][k],h[j][k]);
				break;
			}	
		}	
		for (int j=i-1;j>0;--j) if (h[j][i]) {
			double tmp=h[j][i]/h[i][i];
			for (int k=1;k<=n+n;++k) h[j][k]-=h[i][k]*tmp,h[j][k]=(fabs(h[j][k])<eps?0:h[j][k]);
		}		
	}
	for (int i=1;i<=n;++i) if (fabs(h[i][i])<eps) return false;
	return true;
}
bool abx[maxn],aby[maxn];
bool dfs(int x) {
	for (int i=1;i<=n;++i) if (fabs(c[i][x])>eps && !alr[i] && !aby[i])  {
		alr[i]=true;
		if (!match[i] || dfs(match[i])) {
			match[i]=x;
			return true;
		}
	}
	return false;
}
void change(int x) {
	int k;
	for (int i=1;i<=n;++i) if (match[i]==x) {
		k=i;
		break;
	}
	match[k]=0;
	for (int i=1;i<=n;++i) if (fabs(c[i][x])>eps && !alr[i] && !aby[i]) {
		alr[i]=true;
		if (!match[i] || (!abx[match[i]] && dfs(match[i]))) {
			match[i]=x;
			abx[x]=true;
			aby[i]=true;
			return;
		}
	}
	abx[x]=aby[k]=true;
	match[k]=x;
}
int Ans[maxn];
int main() {
	#ifndef ONLINE_JUDGE
		freopen("test.in","r",stdin);
    		freopen("my.out","w",stdout);
	#endif
	n=read();
	for (int i=1;i<=n;++i) for (int j=1;j<=n;++j) h[i][j]=read();
	for (int i=1;i<=n;++i) for (int j=1;j<=n;++j) b[i][j]=read();
	for (int i=1;i<=n;++i) h[i][i+n]=1;
	bool flag=elm();
	if (!flag) {
		puts("NIE");
		return 0;
	}
	for (int i=1;i<=n;++i) for (int j=1;j<=n;++j) a[i][j]=(double)h[i][j+n]/h[i][i];
	for (int i=1;i<=n;++i) for (int j=1;j<=n;++j) for (int k=1;k<=n;++k) c[i][j]+=b[i][k]*a[k][j];
	int ans=0;
	for (int i=1;i<=n;++i) memset(alr,0,sizeof alr),ans+=dfs(i);
	if (ans!=n) {
		puts("NIE");
		return 0;
	} else puts("TAK");
	for (int i=1;i<=n;++i) {
		memset(alr,0,sizeof alr);
		change(i);
	}
	for (int i=1;i<=n;++i) Ans[match[i]]=i;
	for (int i=1;i<=n;++i) printf("%d\n",Ans[i]);
}