Atcoder abc323

發表於 更新於 分類於 ， Disqus： 文章字數： 2.5k 所需閱讀時間 ≈ 8 分鐘
G(基尔霍夫定律（矩阵树定理），特征多项式（Hessenberg算法）)

https://atcoder.jp/contests/abc323

G - Inversion of Tree

给定一个 1~n的排列p

1
2
for k = 0..n-1:
	求 对于n点树，满足((边端点 点的idx大小 和 p[点idx]大小 是相反的) 的边的个数 == k) 的树的个数 mod 998244353

n 500

4s

1024mb

我的思路

因为这里要的相反不是 父子和p相反，而是index和p[index]相反， 所以如果要选根，选1就行

然后问题是 如果在树上做dp,不光要记录满足树和子树的点树，还需要记录哪些点已经用了，感觉这状态有点多

注意到$n\le 500$是有点小的，可能有$n^2,n^3$的考虑的方向，但实际上这里还有一个k,所以可能是$nk,n^2k,nk^2$一类的

但直接的思路还是没有

然后就是插入新点的想法

如果有一个n-1点的方案，那么对于 p[1]...p[n-1]需要把所有大于$p[n]$的p进行减1,这样就是规模减1的子问题，

但是问题是，这样的话，并没有看出n-1到n的转移

首先钦定了根是1，所以(n,p(n))不是根，如果它是叶子，那么 一个n-1树上原来叶子p[叶] >= p[n] 都会贡献到对应的k+1, 而原来p[叶] < p[n] 会贡献到k

那还不用说非叶子的，光是这里，就需要记录叶子的p和更大的p之间的关系的个数，那么更小的层数时要记的关联太多了

似乎连一个暴力的计数法都没有

然后 如果 == k难求，可以考虑求 $\le k$ 然后做相邻相减

题解

Kirchhoff’s theorem: 无自环可重边无向图的 基尔霍夫矩阵（拉普拉斯矩阵）的任意余子式的行列式 = 树的个数

而 定理中的证明过程有用到

$B$是行表示点，列表示边的，值+1/-1表示大小点的邻接矩阵

$L=B\cdot B^T$

$L$去掉第i行第i列的主余子式 $= \sum_{选n-1列} \mathrm{det}(B_{去掉i行选n-1列})\mathrm{det}(B_{去掉i行选n-1列}^T)$

而如果我们让 满足本题统计的 逆向边乘上 $\sqrt{x}$, 其中$x$是生成函数里一样的$x$没有具体值的一个单位

那么一个具体的树的贡献 就是

$x^{满足题意统计的个数}$

那么$[x^k] L$的主余子式,就是要求的$\mathrm{ans}_k$

因此 $L_{i,j}$为$1$或$x$（看i和j的关系）求出的$L$的主余子式还是一个 生成函数

因此问题变成了计算一个$M_0+M_1x$的行列式的问题，（因为本题对应的L只包含1(不满足题意统计),n-1(对角线的度),x(满足题意统计)

$O(N^4)$的办法

把$x=0\cdots N$带入 求行列式

再反过来求系数

每次求行列式是$O(N^3)$因此总的是$O(N^4)$

$O(N^3)$的办法

把$\mathrm{det}(xI-M)$称作矩阵$M$的 特征多项式

而特征多项式可以$O(N^3)$ https://judge.yosupo.jp/problem/characteristic_polynomial

所以需要转化

如果$M_0+M_1x$中的$M_1$是regular的(行列式不为0)，可以通过初等行列变换(elementary row and column operations.)成I，即存在$I=AM_1B$

那么有$\det({M_0 + M_1 x })= \frac{1}{\det A \det B}\det(A (M_0 + M_1 x)B) = \frac{1}{\det A \det B}\det(AM_0B + I x),$
那么只要找到$AM_0B$再用 计算特征多项式的算法

如果$M_1$是 not regular

随机取一个$a$

令$y = x-a$有$x=y+a$

$\det({M_0 + M_1 x })= \det({M_0 + M_1 (y+a) })= \det({(M_0+M_1a) + M_1y})$

注意到$\det({(M_0+M_1a) + M_1y})$与$\det({(M_0+M_1a)z + M_1})$ 在对应的特征多项式 的y幂次与z的幂次 正好 逆幂次对应系数

所以 可以求$\det({(M_0+M_1a)z + M_1})$ , 而$M_0+M_1a$因为随机取的关系，有高概率是 regular的

另一个想法是整体的想法, 当我们在 对$M_0+M_1x$ 进行基础行列变换时，它的整体是一个

1
2
3
4
5
6
7
              i当前要处理的列
x+v00   v01  v02 ...
  v10 x+v11  v12 ...
  v20   v21  v22                       <-- i当前要处理的行
   .
   .
   .

如果对于$M_1$做变换时 在列上 找不到非0的系数，说明该列全为不带x的值，不妨把该列直接乘上x

1
2
3
4
5
6
7
8
              i当前要处理的列
x+v00   v01  v02x ...
  v10 x+v11  v12x ...
  v20   v21  v22x                       <-- i当前要处理的行
   .         v23x
   . 
   .

这样整个行列式 被乘上了x, 最后做偏移量即可，比随机取一个a更确定性，而且对应的操作也就是把$M_0$的对应列搬到$M_1$，如果操作数大于n 直接表达式就为0了

特征多项式O(N^3)

$\mathrm{det}(I_nx-A_{n\times n})$

  1. 利用相似矩阵特征多项式相同的性质将给定矩阵消成 上 海森堡 矩阵。
  2. 使用行列式展开定理递推计算它的特征多项式。

相似矩阵行列式相等

$A\sim B$那么

$|A|=|P^{-1}BP|=|P^{-1}||B||P|=|P^{-1}||P||B|=|P^{-1}P||B|=|I||B|=|B|$

相似矩阵特征多项式相等

$A\sim B$那么

$|\lambda I -A|=|P^{-1}(\lambda I -B)P|=|P^{-1}||\lambda I-B||P|=|P^{-1}||P||\lambda I-B|=|\lambda I -B|$

我们用上面的这个思路把矩阵变成 这样的形式（上 海森堡 矩阵）

左乘行变换，右乘列变换，例如

  • （第3n行减去k倍第2行）A（第2列加上k倍的3n列）就能完成第1列的 下方0
  • （第4n行减去k倍第3行）A（第3列加上k倍的4n列）就能完成第2列的 下方0

$$\begin{pmatrix} a_{1,1}&a_{1,2}&a_{1,3}&\cdots&a_{1,n}\\ a_{2,1}&a_{2,2}&a_{2,3}&\cdots&a_{2,n}\\ 0&a_{3,2}&a_{3,3}&\cdots&a_{3,n}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&a_{n,n} \end{pmatrix}$$

$\displaystyle f_n(x)=(x-a_{n,n})f_{n-1}(x)-\sum_{i=1}^{n-1}f_{n-i-1}(x)\left( \prod_{j=n-i+1}^n a_{j,j-1} \right) a_{n-i,n}$

可以$O(n^3)$算出

代码

https://atcoder.jp/contests/abc323/submissions/49822377

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
#include <bits/stdc++.h>
#include <atcoder/modint>
const int MOD=998244353;
using namespace std;
typedef long long ll;
#define rep(i,a,n) for (ll i=a;i<(ll)(n);i++)
#define per(i,a,n) for (ll i=n;i-->(ll)(a);)
ll read(){ll r;scanf("%lld",&r);return r;}
template<class T> // det(Ix-A), T is your custom modint
std::vector<T> characteristic_polynominal(std::vector<std::vector<T>> A){ // copy(need modify), A[0-index]
  int n = A.size();
  if(n == 0) return vector<T>{1};
  assert(n == (int)A[0].size());
  // 先转化成 上 海森堡 矩阵
  rep(j,0,n-1) { // 列
    rep(i,j+1,n) if (A[i][j] != 0) { // 取一个j列 行主元 非0的行
      swap(A[j+1], A[i]); // 交换
      rep(k,0,n) swap(A[k][j+1], A[k][i]); // 列交换
      break;
    }
    if (A[j+1][j]!=0){ // 基础行/列变换,[i=j+2..n][j]置为0
      T ajj_inv = A[j+1][j].inv();
      rep(i,j+2,n){
        T c = A[i][j] * ajj_inv;
        rep(k,j,n) A[i][k  ] -= A[j+1][k] * c; // row[i]  -=A[i][j]*A[j+1][j].inv()*row[j+1]
        rep(k,0,n) A[k][j+1] += A[k  ][i] * c; // col[j+1]+=A[i][j]*A[j+1][j].inv()*col[i]
      }
    }
  }
  // 暴力展开最后一列
  std::vector<std::vector<T> > f(n+1); // 1-index
  f[0] = {1}; // |[0x0]| = 1
  // f[i+1] = x*f[i] - sum ( A[j=i..0][i] * (prod A[_=i..j+1][_-1]) * f[j])
  rep(i,0,n) {
    f[i+1].push_back(0);
    rep(j,0,size(f[i])) f[i+1].push_back(f[i][j]); // *x
    T prod=1;
    per(j,0,i+1){ // i..0
      T t = prod*A[j][i];
      rep(k,0,size(f[j])) f[i+1][k] -= f[j][k]*t;
      if(j) prod *= A[j][j-1];
    }
  }
  return f[n];
}
template<class T> // det(Ax+B), T is your custom modint
std::vector<T> matrix_det(std::vector<std::vector<T>> A,std::vector<std::vector<T>>B){ // copy(need modify), A[0-index]x+B
  int n = A.size();
  auto op=[&](auto fn){fn(A);fn(B);}; // 同时操作A,B
  // Ax+B => LARx + LBR => Ix + B' 的形式
  int mul=0; // det(Ix+B') = det(Ix+newB)/x^mul
  T coef=1; // 从矩阵中统一提取出来的系数?
  rep(i,0,n) {
    // 行变换 让A[i][i]!=0
    auto _=[&](){ rep(j,i,n) if(A[j][i] != 0){ // A,B 同时行变换, 让A[i][i] != 0
      op([&](auto &M){ swap(M[i],M[j]);});
      if(i!=j) coef*=-1;
      break;
    } };
    if(A[i][i]==0) { // i列全为不带x, A[0~n][i]==0, 则i列乘上x, 最后det/x
      rep(j,0,n) swap(A[j][i],B[j][i]);
      mul++;
      _();
    }
    if(A[i][i]==0) return vector<T>(n+1,0); // 该列A,B全是0
    T inv=A[i][i].inv();
    coef *= A[i][i];
    // A的[i][i] 行变成 1
    op([&](auto&M){ rep(j,0,n) M[i][j]=M[i][j]*inv; });
    rep(j,0,n) if(j!=i and A[j][i]!=0) { // if(A[j][i]!=0) { 把A中i列上，除了i行的全行变换 变成0
      T x=A[j][i];
      op([&](auto&M){ rep(k,0,n) M[j][k]-=M[i][k]*x; });
    }
    rep(j,i+1,n) if(A[i][j] != 0){ // 列变换 可以省略吗？省略了似乎不能那样搬列了?
      T x=A[i][j];
      op([&](auto&M){ rep(k,0,n) M[k][j]-=M[k][i]*x; });
    }
  }
  // Ix+B => Ix-B
  rep(i,0,n) rep(j,0,n) B[i][j]=-B[i][j];
  auto cp=characteristic_polynominal(B); // 即求B的特征方程
  vector<T> ans(n+1,0);
  rep(i,mul,n+1) ans[i-mul]=cp[i]*coef;
  return ans;
}
// ------------
using mint = atcoder::static_modint<MOD>;
int main(){
  int n=read();
  vector<int> p(n,0);
  rep(i,0,n) p[i]=read();
  vector A(n,vector<mint>(n,0)); // kirchhoff's matrix, B+Ax
  vector B(n,vector<mint>(n,0));
  rep(i,0,n) rep(j,i+1,n) { // i < j
    vector<vector<mint> > *ptr= p[i] > p[j] ? &A: &B;
    (*ptr)[i][i]++;
    (*ptr)[j][j]++;
    (*ptr)[i][j]--;
    (*ptr)[j][i]--;
  }
  auto op=[&](auto fn){fn(A);fn(B);}; // 同时操作A,B
  // 求n-1阶主余子式 A[n-1,n-1] x + B[n-1,n-1]
  op([&](auto&M){M.pop_back();for(auto &row:M)row.pop_back();});
  n--;
  auto ans=matrix_det(A,B); // det(Ax+B)
  rep(i,0,n+1) printf("%d ",ans[i].val());
  return 0;
}

总结

G:

abc253出现过一次矩阵树定理

所以 图的生成树问题就要向 基尔霍夫定理（矩阵树）想，第二次遇到了，回忆是做过但自己没想到

Hessenberg算法第一次学 https://judge.yosupo.jp/problem/characteristic_polynomial

矩阵快速幂加速 https://codeforces.com/problemset/problem/923/E

https://oi-wiki.org/math/linear-algebra/char-poly/