BZOJ3168. 【HEOI2013】钙铁锌硒维生素

题目大意

给定两个$n*n$的矩阵$A,B$，求一个字典序最小的排列$a$满足将任意一个$A_i$换成$B_{a_i}$后，矩阵$A$仍然满秩。

Data Constraint
$n\leq300$

题解

假如我们知道了$A_i$可以被哪一个$B$替换，就可以构出一个二分图，然后匈牙利算法跑出方案即可。
现在的问题是如何构图。
设存在矩阵$C$，$C_{i,j}\neq0$表示$B_i$需要有$A_j$来表出。那么$C$满足$C*A=B$，两边同时左乘$A^{-1}$，$C=B*A^{-1}$。

---

矩阵求逆

如何对矩阵求逆？
考虑构造一个伴随矩阵，初始时伴随矩阵是单位矩阵。对原矩阵进行高斯消元，同时对伴随矩阵做同样的变换。将原矩阵消成单位矩阵后，伴随矩阵就是原矩阵的逆矩阵。

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

---

如何求字典序最小的方案？先跑一遍匈牙利，求出一个初始解，然后跑第二遍匈牙利，贪心地匹配尽可能小的，强制每次只能调整比自己大的匹配。

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

SRC

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std ;

#define N 300 + 10
#define M 90000 + 10
const double eps = 1e-6 ;

double a[N][N] , inv[N][N] , b[N][N] , c[N][N] ;
int from[M] , flag[M] , f[N][N] ;
int n , ans ;

bool GetInv() {
    for (int i = 1 ; i <= n ; i ++ ) {
        int p = -1 ;
        for (int j = 1 ; j <= n ; j ++ )
            if ( fabs( a[i][j] - 0 ) > eps ) { p = i ; break ; }
        if ( p < 0 ) continue ;
        double d = a[i][p] ;
        for (int j = 1 ; j <= n ; j ++ ) {
            a[i][j] /= d ;
            inv[i][j] /= d ;
        }
        for (int j = 1 ; j <= n ; j ++ ) {
            if ( i == j || fabs( a[j][p] - 0 ) <= eps ) continue ;
            double d = a[j][p] ;
            for (int k = 1 ; k <= n ; k ++ ) {
                a[j][k] -= d * a[i][k] ;
                inv[j][k] -= d * inv[i][k] ;
            }
        }
    }
    for (int i = 1 ; i <= n ; i ++ ) if ( !a[i][i] ) return 0 ;
    return 1 ;
}

int DFS( int x , int k ) {
    for (int i = 1 ; i <= n ; i ++ ) {
        if ( !f[x][i] || flag[i] == k ) continue ;
        flag[i] = k ;
        if ( !from[i] || DFS( from[i] , k ) ) {
            from[i] = x ;
            return 1 ;
        }
    }
    return 0 ;
}

int DFS2( int x , int k ) {
    for (int i = 1 ; i <= n ; i ++ ) {
        if ( !f[x][i] || flag[i] == k ) continue ;
        flag[i] = k ;
        if ( from[i] == k || ( from[i] > k && DFS2( from[i] , k ) ) ) {
            from[i] = x ;
            return 1 ;
        }
    }
    return 0 ;
}

int main() {
    scanf( "%d" , &n ) ;
    for (int i = 1 ; i <= n ; i ++ ) {
        for (int j = 1 ; j <= n ; j ++ ) scanf( "%lf" , &a[i][j] ) ;
        inv[i][i] = 1 ;
    }
    if ( !GetInv() ) { printf( "NIE\n" ) ; return 0 ; }
    for (int i = 1 ; i <= n ; i ++ ) {
        for (int j = 1 ; j <= n ; j ++ ) {
            scanf( "%lf" , &b[i][j] ) ;
        }
    }
    for (int i = 1 ; i <= n ; i ++ ) {
        for (int j = 1 ; j <= n ; j ++ ) {
            for (int k = 1 ; k <= n ; k ++ ) {
                c[i][j] += b[i][k] * inv[k][j] ;
            }
        }
    }
    for (int i = 1 ; i <= n ; i ++ ) {
        for (int j = 1 ; j <= n ; j ++ ) {
            if ( fabs( c[j][i] - 0 ) > eps ) f[i][j] = 1 ;
        }
    }
    for (int i = 1 ; i <= n ; i ++ ) ans += DFS( i , i ) ;
    if ( ans < n ) { printf( "NIE\n" ) ; return 0 ; }
    memset( flag , 0 , sizeof(flag) ) ;
    for (int i = 1 ; i <= n ; i ++ )
        ans = DFS2( i , i ) ;
    printf( "TAK\n" ) ;
    for (int i = 1 ; i <= n ; i ++ ) {
        for (int j = 1 ; j <= n ; j ++ ) {
            if ( from[j] == i ) { printf( "%d\n" , j ) ; break ; }
        }
    }
    return 0 ;
}

以上.