[BZOJ3168][Heoi2013]钙铁锌硒维生素 高斯消元+匈牙利


原题链接
好像很长时间没更博了...
thuwc之前更一波证明我还活着 当复习了.

这个题就是说,给你两个矩阵A,B,其中矩阵A是满秩的,要求对于矩阵A的每个行向量,在B中找出一个能替代它(使得矩阵仍然满秩)的行向量,要求字典序最小.
注意到矩阵规模非常小,我们只需要处理向量Bi能否替代Aj,然后跑二分图匹配即可.
由于要求替换后整个A矩阵依然是线性无关的,因此Bi能替代Aj,当且仅当Bi和除Aj以外的其他向量线性无关.
换句话说,令,则Bi能否替代Aj取决于Cj是否为0.
因此只需要求出一个C矩阵,使得即可.
也就是,因此求A的逆矩阵,然后用B和它相乘,就得到了C矩阵,它也是二分图的邻接矩阵.然后二分图匹配就好了.
字典序最小怎么办呢?先跑一组解,然后在此基础上看看更靠前的A向量是否可以匹配更靠前的B向量.

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define N 305
using namespace std;
const double eps=1e-6;
inline bool eq(double a,double b)
{
    return max(a-b,b-a)<eps;
}
int n;
double a[N][N<<1],b[N][N],inv[N][N],c[N][N];
int map[N][N];
inline void swap_line(int x,int y)
{
    for(int i=1;i<=n*2;++i)swap(a[x][i],a[y][i]);
}
bool gauss()
{
    int i,j,k;
    double tmp;
    for(i=1;i<=n;++i)
    {
        for(j=i;j<=n&&eq(a[j][i],0);++j);
        if(j>n)return 0;
        if(j!=i)swap_line(i,j);
        tmp=1/a[i][i],a[i][i]=1;
        for(j=i+1;j<=2*n;++j)a[i][j]*=tmp;
        for(j=1;j<=n;++j)if(j!=i)
        {
            tmp=a[j][i],a[j][i]=0;
            for(k=i+1;k<=2*n;++k)a[j][k]-=a[i][k]*tmp;
        }
    }
    for(i=1;i<=n;++i)
        for(j=1;j<=n;++j)inv[i][j]=a[i][j+n];
    return 1;
}
int from[N],to[N],vis[N],tim;
bool xon(int x)
{
    for(int i=1;i<=n;++i)if(map[x][i]&&vis[i]!=tim)
    {
        vis[i]=tim;
        if(!from[i]||xon(from[i]))
        {
            from[i]=x,to[x]=i;
            return 1;
        }
    }
    return 0;
}
bool xon(int x,int y)
{
    for(int i=1;i<=n;++i)if(map[x][i]&&vis[i]!=tim)
    {
        vis[i]=tim;
        if(from[i]==y||(from[i]>=y&&xon(from[i],y)))
        {
            from[i]=x,to[x]=i;
            return 1;
        }
    }
    return 0;
}
int main()
{
    cin>>n;
    int i,j,k;
    for(i=1;i<=n;++i)
        for(j=1;j<=n;++j)scanf("%lf",&a[i][j]);
    for(i=1;i<=n;++i)
        for(j=1;j<=n;++j)scanf("%lf",&b[i][j]);
    for(i=1;i<=n;++i)a[i][i+n]=1;
    if(!gauss())return puts("NIE"),0;
    for(i=1;i<=n;++i)
    {
        for(j=1;j<=n;++j)
        {
            for(k=1;k<=n;++k)c[i][j]+=b[i][k]*inv[k][j];
        }
    }
    for(i=1;i<=n;++i)
        for(j=1;j<=n;++j)if(!eq(c[i][j],0))map[j][i]=1;
    for(i=1;i<=n;++i)
    {
        tim++;
        if(!xon(i))return puts("NIE"),0;
    }
    puts("TAK");
    for(i=1;i<=n;++i)tim++,xon(i,i);
    for(i=1;i<=n;++i)printf("%d\n",to[i]);
    return 0;
}

发表评论

电子邮件地址不会被公开。 必填项已用*标注