题解 P4100 【[HEOI2013]钙铁锌硒维生素 】
cosmicAC
2019-02-07 14:26:54
题意：给定n个线性无关的向量`A[1..n]`(如果不是线性无关直接输出无解即可),另外n个向量`B[1..n]`，求能否给A中的每一个向量选择一个B中的备用向量，使得任意两个备用向量在B中编号不同，且A中的一个向量的备用向量和A中其余向量线性无关。
显然可以先对A中的每一个向量确定哪些向量可以备用，然后跑一个二分图最小字典序完美匹配即可。没有完美匹配则无解。
首先可以考虑一个系数矩阵V,`V[i][j]`表示“B中第i个向量用A的线性组合表示时，A[j]项的系数”。容易证明`A[i]`可以使用`B[j]`作为备用向量当且仅当$V_{ji}\ne 0$。（如果$V_{ji}=0$,`B[j]`是A中其余向量的线性组合；否则可以移一下项把A[i]用B[j]表示）
可以写出一个等式：$$B_{ij}=\sum_{k=1}^n{V_{ik}A_{kj}}$$即$$B=VA;V=BA^{-1}$$
所以矩阵求逆即可。不会矩阵求逆可以看[我的博客](https://www.luogu.org/blog/474D/solution-p4783)
但是这里不会矩阵求逆也没有关系。可以在把A和B拼成一个`2n*n`的矩阵（还是叫矩阵A），高斯消元的同时维护一下当前的`A[i]`由初始的`A[1..n]`线性组合表示的各项系数。复杂度还是$O(n^3)$。（也许这个算法本质上就是伴随矩阵？）
然后需要求二分图最小字典序完美匹配。我一开始觉得先跑一个Dinic，然后退流就行了，但是退流好像不能凭空加入一条新的匹配边。所以我最后写了一个类似匈牙利的方法，维护当前哪些点的匹配已经被固定了，找到一条不经过已固定的点的交错路（这里其实是一个环），然后增广。复杂度仍然是$O(n^3)$，实际跑的飞快。
C++11代码：
```cpp
// luogu-judger-enable-o2
#include<bits/stdc++.h>
using namespace std;
int n,vis[610];
double a[610][310],b[610][310],sm[610];
struct Dinic{
struct edge{int v,cap;};
vector<edge> e,e1;
vector<int> v[610];
int dis[610],cur[610],que[610],f[610],s,t;
void addEdge(int x,int y){
e.insert(e.end(),{{y,1},{x}});
v[x].push_back(e.size()-2);
v[y].push_back(e.size()-1);
}
int bfs(int s,int t){
int l=1,r=1;
memset(dis,0x3f,sizeof dis);
que[1]=s;dis[s]=0;
while(l<=r){
int p=que[l++];
for(int i:v[p])if(e[i].cap && !vis[e[i].v] && dis[e[i].v]>1e9)
f[e[i].v]=i^1,dis[e[i].v]=dis[p]+1,que[++r]=e[i].v;
}
return dis[t]<1e9;
}
void wtf(int x,int y){
for(int p=x;p!=y;p=e[f[p]].v)e[f[p]].cap=1,e[f[p]^1].cap=0;
}
int dfs(int p,int a){
if(p==t || !a)return a;
int f,to,fl=0;
for(int&i=cur[p];i<(int)v[p].size();++i){
edge &E=e[v[p][i]];
if(dis[to=E.v]==dis[p]+1 && (f=dfs(to,min(a,E.cap)))){
fl+=f,a-=f,e[v[p][i]^1].cap+=f,E.cap-=f;
if(!a)break;
}
}
return fl;
}
int dinic(int s,int t){
this->s=s,this->t=t;
int f=0;
while(bfs(s,t))
memset(cur,0,sizeof cur),f+=dfs(s,1e9);
return f;
}
}sol;
int main(){
scanf("%d",&n);
for(int i=0;i<2*n;i++)
for(int j=0;j<n;j++)
scanf("%lf",a[i]+j);
for(int i=0;i<n;i++)
b[i][i]=1,sol.addEdge(2*n,i),sol.addEdge(i+n,2*n+1);
for(int i=0;i<n;i++){
double mx=-1;int mp;
for(int j=0;j<2*n;j++)if(!vis[j]){
sm[j]=0;
for(int k=0;k<n;k++)
sm[j]+=b[j][k]*a[k][i];
if(j<n){if(fabs(sm[j])>mx)mx=fabs(sm[j]),mp=j;}
else sm[j]+=a[j][i];
}
if(mx<1e-10){puts("NIE");return 0;}
vis[mp]=1;
for(int j=0;j<2*n;j++)if(!vis[j] && fabs(sm[j])>1e-10)
for(int k=0;k<n;k++)
b[j][k]-=b[mp][k]*sm[j]/sm[mp];
}
for(int i=n;i<2*n;i++)
for(int j=0;j<n;j++)if(fabs(b[i][j])>1e-10)
sol.addEdge(j,i);
memset(vis,0,sizeof vis);
if(sol.dinic(2*n,2*n+1)<n){puts("NIE");return 0;}
puts("TAK");
vis[2*n]=1;
for(int i=0,to;i<n;i++)
for(int j:sol.v[i])if(!vis[to=sol.e[j].v]&&(!sol.e[j].cap||sol.bfs(to,i))){
printf("%d\n",to-n+1);vis[to]=1;
if(sol.e[j].cap)
sol.wtf(i,to),sol.e[j].cap=0,sol.e[j^1].cap=1;
break;
}
return 0;
}
```