bzoj3168 [Heoi2013]钙铁锌硒维生素

题目链接：bzoj3168
题目大意：
小林准备了2套厨师机器人，一套厨师机器人有n个，每个厨师机器人只会做一道菜，这道菜一斤能提供第i种营养xi微克。第2套厨师机器人被用来做第1套的备用。小林需要为每一个第1套厨师机器人选一个第2套厨师机器人作备份，使得当这个机器人坏掉时，用备份顶替，整套厨师机器人仍然能搭配出任何营养需求，而且，每个第2套厨师机器人只能当一个第1套厨师机器人的备份。如果无法完成任务，输出“NIE”，否则输出“TAK”，并跟着n行，第i行表示第i个第1套机器人的备份是哪一个第2套机器人。如果有多种可能的答案，请给出字典序最小的那一组。

题解：
矩阵的逆+二分图匹配
看后半段是要我们求字典序最小的完备匹配的方案。
那么首先要解决的就是匹配问题。
讲真我看题看了很久，废话太多了。。
首先知道，任意一组线性无关的向量都可以当作基底来表示任意向量。
从上面的加粗部分可知，第一套机器人是线性无关的，同时若第2套中的机器人能当第一套中某个机器人备用，那么也要求替换后的机器人组是线性无关的。
参照ATP的题解
我们设一个参数矩阵$C$，第一套用$A$表示，第二套用$B$表示。
因为$A$是线性无关的向量组，所以每个$B$中的向量都可以用$A$表示。
即有$B_i=\sum_{i=1}^nC_iA_i$，即$B=C\times A$。
要看$B_i$能不能替换$A_j$的话，只要看那个系数为不为0就好了。因为如果系数为0，就是说剩下来的$A$仍能表示出$B_i$，就不满足线性无关了。
而要求$C$，只要求$A^{-1}$，两边同乘$A^{-1}$就好了。
然后这个就是用高斯消元求。不要直接用double会除到0的。

诶各种打错+求字典序的时候随便贪心秒WA了4发。。
二分图匹配要求字典序最小不能一遍做啊。。我还以为跟bzoj1562一样。。从大到小就可以了呢。
应该这样做：先做一遍完备匹配，然后从1到n去从小到大贪心，看看能不能找到不影响前面的交错环来换。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
#define N 310

const LL mod=1e9+7;
// const double eps=1e-5;
LL a[N][N],b[N][N],inv[N][N];int n;
bool map[N][N];int bf[N],gf[N],vis[N],tim;
LL qpow(LL x,int t)
{
    LL ret=1;
    while (t)
    {
        if (t&1) ret=(ret*x)%mod;
        x=(x*x)%mod;t>>=1;
    }return ret;
}
void gauss()
{
    for (int i=1;i<=n;i++) inv[i][i]=1;
    for (int i=1;i<=n;i++)
    {
        int now=i;
        for (int j=i+1;j<=n;j++)
          if (abs(a[now][i])<abs(a[j][i])) now=j;
        if (now!=i)
        {
            for (int j=1;j<=n;j++)
            {
                swap(a[now][j],a[i][j]);
                swap(inv[now][j],inv[i][j]);
            }
        }
        LL tt=qpow(a[i][i],mod-2)%mod;
        // double tt=1.0/a[i][i];
        for (int j=1;j<=n;j++)
        {
            // a[i][j]*=tt;inv[i][j]*=tt;
            a[i][j]=a[i][j]*tt%mod;
            inv[i][j]=inv[i][j]*tt%mod;
        }
        for (int j=1;j<=n;j++)
         if (j!=i && abs(a[j][i])>0)
         {
             // tt=a[j][i];
             tt=a[j][i]%mod;
             for (int k=1;k<=n;k++)
             {
                 // a[j][k]-=a[i][k]*tt;
                 // inv[j][k]-=inv[i][k]*tt;
                a[j][k]=(a[j][k]-a[i][k]*tt)%mod;
                inv[j][k]=(inv[j][k]-inv[i][k]*tt)%mod;
             }
         }
    }
}
void get_c()
{
    for (int i=1;i<=n;i++)
     for (int j=1;j<=n;j++)
     {
        LL now=0;
        for (int k=1;k<=n;k++)
            // now+=b[i][k]*inv[k][j];
         now=(now+b[i][k]*inv[k][j]%mod)%mod;
        if (now) map[j][i]=true;
        // if (fabs(now)>eps) map[i][j]=true;
        else map[j][i]=false;
     }
}
bool ffind(int x)
{
    for (int i=1;i<=n;i++)
     if (map[x][i] && vis[i]!=tim)
     {
         vis[i]=tim;
         if (bf[i]==-1 || ffind(bf[i]))
         {
             bf[i]=x;gf[x]=i;
             return true;
         }
     }
    return false;
}
bool ffind2(int x,int fr)
{
    for (int i=1;i<=n;i++)
     if (map[x][i] && vis[i]!=tim)
     {
         vis[i]=tim;
         if (bf[i]==fr || bf[i]>fr && ffind2(bf[i],fr))
         //这里的bf[i]>fr就是要保证不影响之前的只做fr之后的
         {
             bf[i]=x;gf[x]=i;
             return true;
         }
     }
    return false;
}
int main()
{
    //freopen("a.in","r",stdin);
    //freopen("a.out","w",stdout);
    int i,j;bool flag=true;
    scanf("%d",&n);
    for (i=1;i<=n;i++)
     for (j=1;j<=n;j++)
      scanf("%lld",&a[i][j]);
    for (i=1;i<=n;i++)
     for (j=1;j<=n;j++)
      scanf("%lld",&b[i][j]);
    gauss();get_c();
    // for (i=1;i<=n;i++)
    // {
       // for (j=1;j<n;j++)
        // printf("%d ",map[i][j]);
       // printf("%d\n",map[i][n]);
    // }
    tim=0;memset(bf,-1,sizeof(bf));
    for (i=n;i>=1;i--)
    {
        tim++;
        if (!ffind(i)) {flag=false;break;}
    }
    if (!flag) printf("NIE\n");
    else
    {
        printf("TAK\n");
        for (i=1;i<=n;i++)
        {
            tim++;
            ffind2(i,i);
        }
        for (i=1;i<=n;i++)
         printf("%d\n",gf[i]);
    }
    return 0;
}

附：高斯消元求矩阵的逆的求法
学习资料：ATP的总结
什么叫矩阵的逆？
设一个矩阵$A$的话，满足$A\times A^{-1}=$单位矩阵 的那个$A^{-1}$就是$A$的逆了吧。
如果我们 把$A$经过一系列的操作变成单位矩阵（大概算是除以$A$）的这个过程用在单位矩阵身上，那么得出来的矩阵就是$A^{-1}$了。
而把 A <script type="math/tex" id="MathJax-Element-23">A</script>变成单位矩阵就是用高斯消元来做啊(单位矩阵就像什么只有一位是1其余为0的线性基)。然后让单位矩阵也做同样的操作就好了。
单独拉出来的代码就是：

void gauss()
{
    for (int i=1;i<=n;i++) inv[i][i]=1;
    for (int i=1;i<=n;i++)
    {
        int now=i;
        for (int j=i+1;j<=n;j++)
          if (abs(a[now][i])<abs(a[j][i])) now=j;//大于小于别写错
        if (now!=i)
        {
            for (int j=1;j<=n;j++)
            {
                swap(a[now][j],a[i][j]);
                swap(inv[now][j],inv[i][j]);
            }
        }
        LL tt=qpow(a[i][i],mod-2)%mod;
        for (int j=1;j<=n;j++)
        {
            a[i][j]=a[i][j]*tt%mod;
            inv[i][j]=inv[i][j]*tt%mod;
        }
        for (int j=1;j<=n;j++)
         if (j!=i && abs(a[j][i])>0)
         {
             tt=a[j][i]%mod;
             for (int k=1;k<=n;k++)
             {
                a[j][k]=(a[j][k]-a[i][k]*tt)%mod;
                inv[j][k]=(inv[j][k]-inv[i][k]*tt)%mod;
             }
         }
    }
}