[BZOJ3168][Heoi2013]钙铁锌硒维生素（矩阵求逆+二分图匹配）

=== ===

这里放传送门

=== ===

题解

这个题的意思是给出一组线性无关的A向量和一组B向量，要求A里面的每一个向量从B里面找到一个替换，然后把A的某一个向量用它的替换代替掉以后整个A向量组仍然是线性无关的。同一个B向量不能作为两个A向量的替换。最后还要求输出字典序最小的解。

首先可以发现因为这个题的数据范围比较小，那么如果能够知道每一个B向量可不可以作为每一个A向量的替换，然后做二分图匹配就可以了。那么关键就是求向量$B_i$能不能作为向量$A_j$的替换。
可以发现因为A和B都是n维向量组，并且都有n个，并且A还是线性无关的，那么也就是说A向量组实际上形成了这个n维空间的一组基底，任何一个n维向量都可以由A中的向量线性表示。那么对于B里面的某一个向量$B_i$，设$B_i=c_1A_1+c_2A_2+...+c_nA_n$。考虑$B_i$能否替换掉$A_j$，可以发现如果$c_j$不为0，也就是说去掉了$A_j$以后剩下的向量不能线性表示出$B_i$，就说明$B_i$和剩下的向量是线性无关的，即$B_i$可以替换掉$A_j$。但是如果$c_j=0$，也就是剩下那些向量可以线性表示出$A_i$，那么$B_i$不能替换掉$A_j$。

那么问题就变成了怎么求出这个系数。可以发现$A$和$B$都是行向量组，并且B中的向量可以被A中的向量线性表示。这就说明存在可逆矩阵$C$，使得$C \times A=B$。写成矩阵的形式也就是

$$\begin{pmatrix} c_{1,1} & c_{1,2} & c_{1,3} & \cdots & c_{1,n} \\ c_{2,1} & c_{2,2} & c_{2,3} & \cdots & c_{2,n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_{n,1} & c_{n,2} & c_{n,3} & \cdots & c_{n,n} \\ \end{pmatrix} \times \begin{pmatrix} A_1\\A_2\\\vdots \\A_n\\ \end{pmatrix}= \begin{pmatrix} B_1\\B_2\\\vdots \\B_n\\ \end{pmatrix}$$
那么我们只需要求出A矩阵的逆矩阵，然后用它右乘B，就可以求出C矩阵啦。然后根据C矩阵的系数 $C_{i,j}$就能判断 $B_i$能否替换掉 $A_j$啦。

---

然而这样这题还没完。。。
它要求输出字典序最小的一组解，但普通的匈牙利算法做出来的匹配是不能保证字典序最小的。那么需要用别的方法来求这个东西。首先用普通的匈牙利跑一遍求出一组“基准解”，这样的话首先可以判断它是不是无解。然后我们从1到n枚举每一个A部分中的点依次考虑改变这个点当前的匹配，看看它能不能选到一个编号更小的点。

首先因为这实际上是一个贪心的思路，让A中编号小的匹配的B的编号也尽量小。也就是当枚举到点$A_i$的时候$A_1..A_{i-1}$的匹配都已经确定下来不会再改变了。那么当枚举到i的时候首先要把点$A_1..A_{i-1}$以及它们的匹配点都从图中删掉，这样才不会影响前面的结果。然后按照编号从小到大枚举每一个B部分中的点，如果从这个点出发能够找到一条新的增广路，那么就可以把这条路径取反然后把$A_i$的匹配改成它，原理和普通的匈牙利算法是类似的。

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const long long Mod=1e9+7;
int n,cnt,link[610],vis[610],tot;
long long v[310][310],rec[310];
bool ext[310][310];
struct element{
    long long p[310];
    int id;
}A[310],B[310],Inv[310];
long long powww(long long a,int t){
    long long ans=1;
    while (t!=0){
        if (t&1) ans=(ans*a)%Mod;
        a=(a*a)%Mod;t>>=1;
    }
    return ans;
}
long long Abs(long long x){return (x<0)?-x:x;}
void Gauss_Eli(){
    int num;
    long long t;
    for (int i=1;i<=n;i++) Inv[i].p[i]=1;
    for (int i=1;i<=n;i++){
        num=i;
        for (int j=i;j<=n;j++)
          if (Abs(A[j].p[i])>Abs(A[num].p[i])) num=j;
        swap(A[i],A[num]);
        swap(Inv[i],Inv[num]);
        t=powww(A[i].p[i],Mod-2);
        for (int j=1;j<=n;j++){
            A[i].p[j]=(A[i].p[j]*t)%Mod;
            Inv[i].p[j]=(Inv[i].p[j]*t)%Mod;
        }
        for (int j=1;j<=n;j++)
          if (i!=j&&A[j].p[i]!=0){
              t=A[j].p[i];
              for (int k=1;k<=n;k++){
                  A[j].p[k]=(A[j].p[k]-(A[i].p[k]*t)%Mod)%Mod;
                  Inv[j].p[k]=(Inv[j].p[k]-(Inv[i].p[k]*t)%Mod)%Mod;
              }//注意逆矩阵也要减去同样的量
          }
    }
}
void Multi(){
    memset(A,0,sizeof(A));
    for (int i=1;i<=n;i++)
      for (int j=1;j<=n;j++)
        for (int k=1;k<=n;k++)
          A[i].p[j]=(A[i].p[j]+(B[i].p[k]*Inv[k].p[j])%Mod)%Mod;
}
bool find(int u,int k){
    for (int i=1;i<=n;i++)
      if (ext[u][i]&&vis[i]!=k){
          vis[i]=k;
          if (link[i]==-1||find(link[i],k)){
              link[i]=u;return true;
          }
      }
    return false;
}
bool find_again(int u,int st,int k){
    if (u==st) return true;
    for (int i=1;i<=n;i++)
      if (ext[u][i]==true&&vis[i]!=k){
          vis[i]=k;
          if (find_again(link[i],st,k)){
              link[i]=u;return true;
          }
      }
    return false;
}
int main()
{
    scanf("%d",&n);
    for (int i=1;i<=n;i++){
        A[i].id=i;
        for (int j=1;j<=n;j++)
          scanf("%d",&A[i].p[j]);
      }
    for (int i=1;i<=n;i++){
        B[i].id=i;
        for (int j=1;j<=n;j++)
          scanf("%d",&B[i].p[j]);
    }
    Gauss_Eli();Multi();
    for (int i=1;i<=n;i++)
      for (int j=1;j<=n;j++)
        if (A[i].p[j]!=0) ext[j][i]=true;
    memset(link,-1,sizeof(link));
    for (int i=1;i<=n;i++)
      if (find(i,i)) ++cnt;
    if (cnt!=n){printf("NIE\n");return 0;}
    printf("TAK\n");
    memset(vis,0,sizeof(vis));
    for (int i=1;i<=n;i++){
        for (int j=1;j<=n;j++) ext[i-1][j]=0;//注意要删去前面的边防止影响结果
        for (int j=1;j<=n;j++)
          if (ext[i][j]&&vis[j]!=i){
              vis[j]=i;
              if (link[j]>=i&&find_again(link[j],i,i)){
                  link[j]=i;printf("%d\n",j);
                  for (int k=1;k<=n;k++) ext[k][j]=0;
                  break;//找到可行解就及时break
              }
          }
    }
    return 0;
}