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

题目大意：给定一个n?n的满秩矩阵A和一个n?n的矩阵B。求一个字典序最小的1...n的排列a满足将随意一个Ai换成Bai后矩阵A仍然满秩

我们考虑建立一个二分图。假设Ai能换成Bj。就在i?>j之间连接一条边

那么这个图怎么建呢？

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

我们能够设矩阵C满足C?A=B，那么C=B?A?1

Ci,j≠0表示Bi的线性表出须要Aj，因此CT就是这个二分图的邻接矩阵

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

我们能够枚举每一条边，然后推断剩余的图是否存在一个完备匹配。可是这样做是O(n4)的

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

总时间复杂度O(n3)

#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;
}