BZOJ 3168 Heoi2013 钙铁锌硒维生素 矩阵求逆+匈牙利算法

题目大意：给定一个$n*n$的满秩矩阵$A$和一个$n*n$的矩阵$B$，求一个字典序最小的$1...n$的排列$a$满足将任意一个$A_i$换成$B_{a_i}$后矩阵$A$仍然满秩

我们考虑建立一个二分图，如果$A_i$能换成$B_j$，就在$i->j$之间连接一条边

那么这个图怎么建呢？

考虑一个行向量$B_i$，我们在$A$中找到最小的行向量集合满足$B_i$可以被这些行向量线性表出，那么显然$B_i$只能替换这些行向量

我们可以设矩阵$C$满足$C*A=B$，那么$C=B*A^{-1}$
$C_{i,j}\neq 0$表示$B_i$的线性表出需要$A_j$，因此$C^T$就是这个二分图的邻接矩阵

现在我们有了一个二分图，如何求字典序最小的完备匹配呢？

我们可以枚举每一条边，然后判断剩余的图是否存在一个完备匹配，但是这样做是$O(n^4)$的

我们可以跑两遍匈牙利算法，第一遍求出任意一个完备匹配，第二遍对于每个点贪心选最小的出边判断能否找到不影响前面点的交错环

总时间复杂度$O(n^3)$

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define M 330
#define MOD 999911657
using namespace std;
int n;
bool state[M];
int result[M];
long long Quick_Power(long long x,int y)
{
    long long re=1;
    while(y)
    {
        if(y&1) (re*=x)%=MOD;
        (x*=x)%=MOD; y>>=1;
    }
    return re;
}
struct Matrix{
    int a[M][M];
    Matrix() {}
    Matrix(bool flag)
    {
        int i;
        memset(a,0,sizeof a);
        for(i=1;i<=n;i++)
            a[i][i]=flag;
    }
    int* operator [] (int x)
    {
        return a[x];
    }
    friend istream& operator >> (istream &_,Matrix &a)
    {
        int i,j;
        for(i=1;i<=n;i++)
            for(j=1;j<=n;j++)
                scanf("%d",&a[i][j]);
        return _;
    }
    friend Matrix operator * (Matrix x,Matrix y)
    {
        Matrix z(false);
        int i,j,k;
        for(i=1;i<=n;i++)
            for(j=1;j<=n;j++)
                for(k=1;k<=n;k++)
                    (z[i][j]+=(long long)x[i][k]*y[k][j]%MOD)%=MOD;
        return z;
    }
    friend Matrix Get_Inv(Matrix a)
    {
        Matrix re(true);
        int i,j,k;
        for(i=1;i<=n;i++)
        {
            for(k=i;k<=n;k++)
                if(a[k][i])
                    break;
            for(j=1;j<=n;j++)
            {
                swap(a[i][j],a[k][j]);
                swap(re[i][j],re[k][j]);
            }
            long long inv=Quick_Power(a[i][i],MOD-2);
            for(j=1;j<=n;j++)
            {
                a[i][j]=a[i][j]*inv%MOD;
                re[i][j]=re[i][j]*inv%MOD;
            }
            for(k=1;k<=n;k++)
                if(k!=i)
                {
                    long long temp=(MOD-a[k][i])%MOD;
                    for(j=1;j<=n;j++)
                    {
                        (a[k][j]+=a[i][j]*temp%MOD)%=MOD;
                        (re[k][j]+=re[i][j]*temp%MOD)%=MOD;
                    }
                }
        }
        return re;
    }
}A,B,C,f;
bool DFS1(int x)
{
    int i;
    for(i=1;i<=n;i++)
        if(f[x][i]&&!state[i])
        {
            state[i]=true;
            if( !result[i] || DFS1(result[i]) )
            {
                result[i]=x;
                return true;
            }
        }
    return false;
}
bool DFS2(int x,int from)
{
    int i;
    for(i=1;i<=n;i++)
        if(f[x][i]&&!state[i])
        {
            state[i]=true;
            if( result[i]==from || result[i]>from && DFS2(result[i],from) )
            {
                result[i]=x;
                return true;
            }
        }
    return false;
}
int main()
{
    int i,j;
    cin>>n>>A>>B;
    C=B*Get_Inv(A);
    for(i=1;i<=n;i++)
        for(j=1;j<=n;j++)
            f[i][j]=(bool)C[j][i];
    for(i=1;i<=n;i++)
    {
        memset(state,0,sizeof state);
        if( !DFS1(i) )
            return puts("NIE"),0; 
    }
    puts("TAK");
    for(i=1;i<=n;i++)
    {
        memset(state,0,sizeof state);
        DFS2(i,i);
    }
    for(i=1;i<=n;i++)
        for(j=1;j<=n;j++)
            if(result[j]==i)
                printf("%d\n",j);
    return 0;
}