HDU 7253

题意

$n$ 个点的有向图，求出以每一个点为根的外向生成树的数量，答案对 $10^9+7$ 取模。

$$n\le 500$$

做法

我们先求出外向生成树的 Laplacian 矩阵 $L^{\text{in}} =D^{\text{in}}(G)-A(G)$

所以对于任意的根 $k$ ，根据矩阵树定理可得：

$$t^{\text{leaf}}(G,k)=\text{det} M_{i,i}$$

其中 $M_{i,i}$ 表示 $L^{\text{in}}$ 去掉 $i$ 行 $i$ 列所得的余子式。

发现暴力做这个事情复杂度是 $O(n^4)$ 的，并不能接受，考虑如何快速处理。

先引入 伴随矩阵 的概念:

$A$ 关于第 $i$ 行第 $j$ 列的余子式（记作 $M_{i,j}$ ）是去掉 $A$ 的第 $i$ 行第 $j$ 列之后得到的 $(n − 1)\times(n − 1)$ 矩阵的行列式。

定义：$A$ 关于第 $i$ 行第 $j$ 列的代数余子式是：

$${\displaystyle \mathbf {C} _{ij}=(-1)^{i+j}\mathbf {M} _{ij}}$$

定义：$A$ 的余子矩阵是一个 $n\times n$ 的矩阵 $C$ ，使得其第 $i$ 行第 $j$ 列的元素是 $A$ 关于第 $i$ 行第 $j$ 列的代数余子式。

引入以上的概念后，可以定义：矩阵 $A$ 的伴随矩阵是 $A$ 的余子矩阵的转置矩阵：

$${\displaystyle \mathrm {adj} (\mathbf {A} )=\mathbf {C} ^{T}}$$

对于 $\text{rank}(L^{\text{in}})$ 进行讨论：

• $\text{rank}(L^{\text{in}})=n$ ，那么 $L^{\text{in}}$ 可逆，由伴随矩阵的性质可得 $\mathrm{adj} (A)=\text{det}(A) A^{-1}$，那么只需要求一次逆再求一个行列式即可

• $\text{rank}(L^{\text{in}})\lt n-1$ ，那么由定义可以得到 $\forall i,j$ , $\text{adj} (A)_{i,j}=0$

• $\text{rank}(L^{\text{in}})=n-1$，可以得到 $\text{rank}( \text{adj}(L^{\text{in}}))=1$。实际上在 Laplacian 矩阵中有更好的性质，在 $L^{\text{in}}$ 中每一行都为 $0$ ，可以得到在 $\text{adj} (L^{\text{in}})$ 中每一列都是相同的。

所以只需要在对 $L^{\text{in}}$ 消元的过程中找到出现自由变元的行，然后删去这一行，枚举删去某一列，注意到删去行列后的矩阵为上海森堡矩阵，所以可以 $O(n^2)$ 求其行列式，最后的复杂度是 $O(n^3)$

#include <bits/stdc++.h>

using namespace std;


using int64 = long long;

template<int MOD>
struct ModInt {
    int x;
 
    ModInt() : x(0) {}
    ModInt(int64 y) : x(y >= 0 ? y % MOD : (MOD - (-y) % MOD) % MOD) {}

    inline int inc(const int &x) {
        return x >= MOD ? x - MOD : x;
    }
    inline int dec(const int &x) {
        return x < 0 ? x + MOD : x;
    }
 
    ModInt &operator+= (const ModInt &p) {
        x = inc(x + p.x);
        return *this;
    } 
    ModInt &operator-= (const ModInt &p) {
        x = dec(x - p.x);
        return *this;
    }
 
    ModInt &operator*= (const ModInt &p) {
        x = (int)(1ll * x * p.x % MOD);
        return *this;
    }
    ModInt &operator/= (const ModInt &p) {
        *this *= p.inverse();
        return *this;
    }
 
    ModInt operator-() const { return ModInt(-x); } 

    friend ModInt operator + (const ModInt& lhs, const ModInt& rhs) {
        return ModInt(lhs) += rhs;
    }
    friend ModInt operator - (const ModInt& lhs, const ModInt& rhs) {
        return ModInt(lhs) -= rhs;
    }
    friend ModInt operator * (const ModInt& lhs, const ModInt& rhs) {
        return ModInt(lhs) *= rhs;
    }
    friend ModInt operator / (const ModInt& lhs, const ModInt& rhs) {
        return ModInt(lhs) /= rhs;
    }
 
    bool operator == (const ModInt &p) const { return x == p.x; } 
    bool operator != (const ModInt &p) const { return x != p.x; }
 
    ModInt inverse() const {
        int a = x, b = MOD, u = 1, v = 0, t;
        while(b > 0) {
            t = a / b;
            swap(a -= t * b, b);
            swap(u -= t * v, v);
        }
        return ModInt(u);
    }

    ModInt pow(int64 n) const {
        ModInt ret(1), mul(x);
        while(n > 0) {
            if(n & 1) ret *= mul;
            mul *= mul;
            n >>= 1;
        }
        return ret;
    }
 
    friend ostream &operator<<(ostream &os, const ModInt &p) {
        return os << p.x;
    }
 
    friend istream &operator>>(istream &is, ModInt &a) {
        int64 t;
        is >> t;
        a = ModInt<MOD>(t);
        return (is);
    }
    static int get_mod() { return MOD; }
};

const int MOD = 1e9 + 7;

using MODint = ModInt<MOD>;
using Varray = vector<vector<MODint>>;

struct Matrix {
  int n , m;
  Varray v;

  void print() const {
    cout << "DEBUG:" << endl;
    for (int i = 1 ; i <= n ; i++) {
      for (int j = 1 ; j <= m ; j++) {
        cout << v[i][j] << ' ';
      }
      cout << endl;
    }
    cout << endl;
  }

  Matrix(int _n , int _m) {
    n = _n , m = _m;
    v = Varray(n + 1 , vector<MODint>(m + 1 , 0));
  }

  Matrix getT() const {
    Matrix ret(m , n);

    for (int i = 1 ; i <= n ; i++) {
      for (int j = 1 ; j <= m ; j++) {
        ret.v[j][i] = v[i][j];
      }
    }
    return ret;
  }

  Matrix operator *(const MODint &k) const {
    Matrix ret = *this;
    for (int i = 1 ; i <= n ; i++) {
      for (int j = 1 ; j <= m ; j++) {
        ret.v[i][j] *= k;
      }
    }
    return ret;
  }

  Matrix operator *(const Matrix &rhs) const {
    Matrix ret(n , rhs.m);
    for (int i = 1 ; i <= n ; i++) {
      for (int k = 1 ; k <= m ; k++) {
        if (v[i][k] == 0) continue;
        for (int j = 1 ; j <= rhs.m ; j++) {
          ret.v[i][j] += v[i][k] * rhs.v[k][j];
        }
      } 
    }
    return ret;
  }

  pair<pair<int , int> , Matrix> gauss(int &swap_cnt) const {
    int r = 0 , p = 0;
    Matrix cur = *this;
    for (int i = 1 ; i <= n ; i++) {
      for (int j = i ; j <= n ; j++) {
        if (cur.v[j][i] != 0) {
          if (i != j) {
            swap(cur.v[i] , cur.v[j]);
            swap_cnt++;
          }
          break;
        }
      }
      if (cur.v[i][i] == 0) {
        p = i;
        continue;
      }
      r++;
      for (int j = i + 1 ; j <= n ; j++) {
        if (cur.v[j][i] == 0) {
          continue;
        }
        MODint d = cur.v[j][i] / cur.v[i][i];
        for (int k = i ; k <= m ; k++) {
          cur.v[j][k] -= d * cur.v[i][k];
        }
      }
    }
    return {{r , p} , cur};
  }

  MODint getval() const {
    MODint ans = 1;
    Matrix cur(n , m);
    for (int i = 1 ; i <= n ; i++) {
      for (int j = 1 ; j <= m ; j++) {
        cur.v[i][j] = v[i][j];
      }
    }
    for (int i = 1 ; i <= n ; i++) {
      for (int j = i ; j <= n ; j++) {
        if (cur.v[j][i] != 0) {
          swap(cur.v[i] , cur.v[j]);
          ans *= -1;
          break;
        }
      }
      for (int j = i + 1 ; j <= n ; j++) {
        if (cur.v[j][i] == 0) {
          continue;
        }
        MODint d = cur.v[j][i] / cur.v[i][i];
        for (int k = i ; k <= n ; k++) {
          cur.v[j][k] -= d * cur.v[i][k];
        }
      }
    }
    for (int i = 1 ; i <= n ; i++) {
      ans *= cur.v[i][i];
    }
    return ans;
  }

  Matrix getinv() const {
    Matrix cur(n , n + n);
    for (int i = 1 ; i <= n ; i++) {
      for (int j = 1 ; j <= n ; j++) {
        cur.v[i][j] = v[i][j];
      }
      cur.v[i][i + n] = 1;
    }
    for (int i = 1 ; i <= n ; i++) {
      for (int j = i ; j <= n ; j++) {
        if (cur.v[j][i] != 0) {
          swap(cur.v[j] , cur.v[i]);
          break;
        }
      }
      MODint d = 1 / cur.v[i][i];
      for (int j = i ; j <= n + n ; j++) {
        cur.v[i][j] *= d;
      } 
      for (int j = 1 ; j <= n ; j++) {
        if (i == j) continue;
        MODint d = cur.v[j][i];
        for (int k = i ; k <= n + n ; k++) {
          cur.v[j][k] -= d * cur.v[i][k];
        }
      }
    }
    Matrix ret(n , n);
    for (int i = 1 ; i <= n ; i++) {
      for (int j = 1 ; j <= n ; j++) {
        ret.v[i][j] = cur.v[i][j + n];
      }
    }
    return ret;
  }

  Matrix getR(int r) const {
    Matrix ret(n - 1 , m);
    int R = 0 , C = 0;
    for (int i = 1 ; i <= n ; i++) {
      if (i == r) continue;
      R++ , C = 0;
      for (int j = 1 ; j <= n ; j++) {
        C++;
        ret.v[R][C] = v[i][j];
      }
    }
    return ret;
  }

  Matrix getRC(int r , int c) const {
    Matrix ret(n - 1 , n - 1);
    int R = 0 , C = 0;
    for (int i = 1 ; i <= n ; i++) {
      if (i == r) continue;
      R++ , C = 0;
      for (int j = 1 ; j <= n ; j++) {
        if (j == c) continue;
        C++;
        ret.v[R][C] = v[i][j];
      }
    }
    return ret;
  }

  Matrix getC(int c) const {
    Matrix ret(n , m - 1);
    for (int i = 1 ; i <= n ; i++) {
      int C = 0;
      for (int j = 1 ; j <= m ; j++) {
        if (j == c) continue;
        C++;
        ret.v[i][C] = v[i][j];
      }
    }
    return ret;
  }

  MODint getval_H() const { // Hessenberg
    MODint ans = 1;
    Matrix cur = *this;
    for (int i = 1 ; i + 1 <= n ; i++) {
      if (cur.v[i + 1][i] == 0) continue;
      if (cur.v[i][i] == 0) {
        swap(cur.v[i] , cur.v[i + 1]);
        ans *= -1;
      }
      MODint d = cur.v[i + 1][i] / cur.v[i][i];
      for (int j = i ; j <= n ; j++) {
        cur.v[i + 1][j] -=  d * cur.v[i][j];
      }
    }
    for (int i = 1 ; i <= n ; i++) {
      ans *= cur.v[i][i];
    }
    return ans;
  }

  Matrix getadj(int type = 0) {
    int swap_cnt = 0;
    auto P = gauss(swap_cnt);

    int r = P.first.first , p = P.first.second;
    auto mat = P.second;

    if (r < n - 1) {
      return Matrix(n , n);
    } else if (r == n - 1) {
      Matrix ans(n , n);
      Matrix t = getR(p);

      int swap_cnt2 = 0;
      auto PI = t.gauss(swap_cnt2);
      auto M = PI.second;
      for (int i = 1 ; i <= p ; i++) {
        auto M1 = M.getC(i);
        ans.v[p][i] = M1.getval_H();
        if ((p + i) % 2 == 1) {
          ans.v[p][i] = -ans.v[p][i];
        }
        if (swap_cnt2 % 2 == 1) {
          ans.v[p][i] = -ans.v[p][i];
        }
      }

      if (type == 0) { // Laplacian
        for (int i = p - 1 ; i >= 1 ; i--) {
          for (int j = 1 ; j <= n ; j++) {
            ans.v[i][j] = ans.v[i + 1][j];
          }
        }
        ans = ans.getT();
      } else {
      }
      return ans;
    } else { // r == n
      return getinv() * getval();
    }
  } 
};


int main() {		
  cin.tie(nullptr)->sync_with_stdio(0);

  int T;
  cin >> T;

  while (T--) {
    int n , m;
    cin >> n >> m;

    vector<vector<int>> a(n + 1 , std::vector<int>(n + 1));
    vector<vector<int>> d(n + 1 , std::vector<int>(n + 1));

    for (int i = 1 ; i <= m ; i++) {
      int u , v;
      cin >> u >> v;

      a[u][v]++;
      d[v][v]++;
    }

    Matrix cur(n , n);
    for (int i = 1 ; i <= n ; i++) {
      for (int j = 1 ; j <= n ; j++) {
        cur.v[i][j] = d[i][j] - a[i][j];
      }
    }

    if (n == 1) {
      cout << "1\n";
      continue;
    }

    auto ans = cur.getadj(0);
    for (int i = 1 ; i <= n ; i++) {
      cout << ans.v[i][i] << " \n"[i == n];
    }
  }
  return 0;
}