bzoj3168

二分图+矩阵求逆

既然我们考虑b能替换哪些a，那么我们自然要得出b被哪些a表示，这里我们设一个矩阵C,那么C*A = B

为什么呢？直接A*C = B是不可行的，因为都是行向量，不能直接乘，那么我们转置一下，得出At*C=Bt,这样就很科学了，那么再转回来，A*Ct=B,于是Ct=B*A^-1那么矩阵求逆就能得出Ct,于是我们再转置回来就能得出B是由哪些A表示的，然后跑二分图匹配就行了。

二分图匹配理解的不是很清楚，大概是不能用之前已经用过较小的编号来更新自己，还是看CQzhangyu的吧。

矩阵求逆用取模就行了

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 305;
const ll P = 999911657;
int rd()
{
    int x = 0, f = 1; char c = getchar();
    while(c < '0' || c > '9') { if(c == '-') f = -1; c = getchar(); }
    while(c >= '0' && c <= '9') { x = x * 10 + c - '0'; c = getchar(); }
    return x * f;
} 
ll power(ll x, ll t) 
{
    ll ret = 1;
    for(; t; t >>= 1, x = x * x % P) 
        if(t & 1) 
            ret = ret * x % P;
    return ret;
}
void up(ll &x, const ll &t)
{
    x = ((x + t) % P + P) % P;
}
int n;
int vis[N], Map[N][N], match[N];
struct Matrix {
    ll a[N][N];
    Matrix() { memset(a, 0, sizeof(a)); }
    void set() { for(int i = 1; i <= n; ++i) a[i][i] = 1; }
    Matrix friend operator * (const Matrix &a, const Matrix &b) {
        Matrix ret;
        for(int i = 1; i <= n; ++i)
            for(int k = 1; k <= n; ++k) if(a.a[i][k])
                for(int j = 1; j <= n; ++j)
                    up(ret.a[i][j], a.a[i][k] * b.a[k][j] % P);
        return ret;
    }
    Matrix Inverse() 
    {
        Matrix c;
        c.set();
        for(int i = 1; i <= n; ++i)
        {
            int now = i;
            while(now <= n && !a[now][i]) ++now;
            for(int j = 1; j <= n; ++j) 
            {
                swap(a[now][j], a[i][j]);
                swap(c.a[now][j], c.a[i][j]);
            }
            ll Inv = power(a[i][i], P - 2);
            for(int j = 1; j <= n; ++j) 
            {
                a[i][j] = a[i][j] * Inv % P;
                c.a[i][j] = c.a[i][j] * Inv % P;
            }
            for(int k = 1; k <= n; ++k) if(k != i)
            {
                ll t = ((P - a[k][i] % P) % P + P) % P;
                for(int j = 1; j <= n; ++j)
                {
                    a[k][j] = ((a[k][j] + t * a[i][j] % P) % P + P) % P;
                    c.a[k][j] = ((c.a[k][j] + t * c.a[i][j] % P) % P + P) % P;  
                }
            }
        }
        return c;
    }
    void print()
    {
        for(int i = 1; i <= n; ++i)  
        {
            for(int j = 1; j <= n; ++j) printf("%d ", a[i][j]);
            puts("");
        }
    }
} a, b, c, d;
bool dfs(int u) 
{
    for(int i = 1; i <= n; ++i) if(!vis[i] && Map[u][i]) 
    {
        vis[i] = 1;
        if(!match[i] || dfs(match[i])) 
        {
            match[i] = u;
            return true;
        }
    }
    return false;
}
bool dfs(int u, int lev) 
{
    for(int i = 1; i <= n; ++i) if(!vis[i] && Map[u][i]) 
    {
        vis[i] = 1;
        if(match[i] == lev || (match[i] > lev && dfs(match[i], lev)))
        {
            match[i] = u;
            return true;
        }
    }
    return false;
}
int main()
{
//  freopen("ferrous.in", "r", stdin);
//  freopen("ferrous.out", "w", stdout);
    n = rd();
    for(int i = 1; i <= n; ++i)
        for(int j = 1; j <= n; ++j) 
            a.a[i][j] = rd();
    for(int i = 1; i <= n; ++i)
        for(int j = 1; j <= n; ++j)
            b.a[i][j] = rd();
    c = b * a.Inverse();
    for(int i = 1; i <= n; ++i)
        for(int j = 1; j <= n; ++j) if(c.a[i][j]) 
            Map[j][i] = 1;
    int ans = 0;
    for(int i = 1; i <= n; ++i)
    {
        memset(vis, 0, sizeof(vis));
        ans += dfs(i);  
    } 
    if(ans < n) 
    {
        puts("NIE");
        return 0;
    }
    for(int i = 1; i <= n; ++i) 
    {
        memset(vis, 0, sizeof(vis));
        dfs(i, i);
    }
    puts("TAK");
    for(int i = 1; i <= n; ++i) 
        for(int j = 1; j <= n; ++j) if(match[j] == i)
            printf("%d
", j);
//    fclose(stdin);
//    fclose(stdout);
    return 0;
}