题解：P10216 【模板】Pfaffian

qianyuzhe · 2025-04-19 17:40:08 · 题解

upd：修改了几处错误，添加了外代数证法。

称一个长度为 2n 的排列 \pi 是完美匹配，当且仅当其满足

• \forall 1\le i\le n,\pi_{2i-1}<\pi_{2i}
• \forall 1\le i<n,\pi_{2i-1}<\pi_{2i+1}

或者一种更加简单的记法：所有完美排列表示将 1\sim 2n 分成 n 对的所有情况。

记 \textup{inv }\pi 表示 \pi 的逆序对数，\textup{sgn }\pi=(-1)^{\textup{inv }\pi}，M_{2n} 表示全体长度为 2n 的完美匹配构成的集合。

令 \mathbf{A}=(a_{i,j})_{1\le i<j\le 2n} 是一个反对称矩阵，定义 \mathbf{A} 的 \text{Pfaffian} 为

举个例子，二阶反对称矩阵

它的 \text{Pfaffian} 值 \textup{Pf}(\mathbf{A})=a。

四阶反对称矩阵

它的 \text{Pfaffian} 值 \textup{Pf}(\mathbf{B})=af-be+cd。

至此还没有什么发现，我们把它们平方一下试试：

是否对于一切偶数阶反对称矩阵都有如上结论呢？

Proof 1. \textup{Pf}(\mathbf{A})^2=\det(\mathbf{A})

证法一（组合证法）：

首先易得

构造 \tau_i=\begin{cases}\pi_i&2\nmid i\\\sigma_i&2|i\end{cases}，则：

再考虑行列式的组合定义：

对于每个排列 \pi，我们对其作轮换分解：

由 \mathbf{A} 的反对称性，

• 若 \exist|\sigma_x|=1，设 \sigma_x=(y)，则 \prod\limits_{i=1}^{n}a_{i,\pi_i} 中有因数 a_{y,y}=0，消去。

• 若 \exist|\sigma_x|>2，设 \pi^{-1}=\prod\limits\sigma_i^{-1}，则 \sigma_x\not=\sigma_x^{-1}，于是 \pi\not=\pi^{-1}，进而 (\textup{sgn }\pi)\prod\limits_{i=1}^{n}a_{i,\pi_i}+(\textup{sgn }\pi^{-1})\prod\limits_{i=1}^{n}a_{i,\pi^{-1}_i}=0，消去。

于是只有可分解为对换之积的排列（即完美匹配）不会被消去，故

定理得证。

证法二（外代数证法）：

要阅读本部分内容，请确认您了解外代数。

将 n 阶反对称矩阵 A 视作二次外形式 \omega(e_i,e_j)（e_i,e_j 是外代数的标准基），有 \omega^n=\land^nV。这个式子是一维的，又由 \text{Pfaffian} 的定义知 \omega^n=\textup{Pf(A)}\land^n。

考虑 \det(A) 的几何定义，熟知 (\omega^n)^2=\det(A)\land^{2n}，于是有 \textup{Pf}(\mathbf{A})^2=\det(\mathbf{A})。

定理得证。

（对于奇数阶反对称矩阵，其行列式值为 0，\text{Pfaffian} 也为 0，证明可类比偶数阶反对称矩阵的证明。）

考虑一个更强的定理：

Proof 2. \text{Pf}(\mathbf{A})\det(\mathbf{Q})=\text{Pf}(\mathbf{Q^T}\mathbf{A}\mathbf{Q})

证明：

先证 [\text{Pf}(\mathbf{A})\det(\mathbf{Q})]^2=[\text{Pf}(\mathbf{Q^T}\mathbf{A}\mathbf{Q})]^2。

再证 \text{Pf}(\mathbf{A})\det(\mathbf{Q})=\text{Pf}(\mathbf{Q^T}\mathbf{A}\mathbf{Q})。

显然当 \mathbf{Q}=\mathbf{E} 时成立，又因为 \det(\mathbf{Q}) 和 \text{Pf}(\mathbf{Q^T}\mathbf{A}\mathbf{Q}) 均为关于 \mathbf{Q} 的多项式函数，且在 \mathbf{Q} 可逆时连续，因此符号必须恒为正（否则会导致不连续性），于是 \text{Pf}(\mathbf{A})\det(\mathbf{Q})=\text{Pf}(\mathbf{Q^T}\mathbf{A}\mathbf{Q})。

定理得证。

这启示我们运用合同变换将原矩阵化为特殊矩阵求解。我们知道，对称矩阵可以通过合同变换化为对角矩阵，对于偶数阶反对称矩阵，由于对其施加任意次合同变换后其依然为反对称矩阵，因此我们不能将其化为对角矩阵，但是我们可以将其化为如下分块对角矩阵：

其中 \mathbf{B_1},\dots,\mathbf{B_n} 均为二阶反对称矩阵。

也就是说，我们可以构造一个矩阵 \mathbf{Q}=\begin{bmatrix}\mathbf{E}&\mathbf{P}\\\mathbf{O}&\mathbf{E}\end{bmatrix} 来将 A 化为 \mathbf{Q^T}。

考虑如下矩阵 \mathbf{A}。

假设现在我们要求出它的 \text{Pfaffian} 值。

我们先从 a_{1,3} 消起：

要使 a_{1,3} 变为 0，有两种选择：

1. r_1-\dfrac{2}{3}r_2
2. c_3-2c_2

为了与后面保持一致，这里我们采用第 2 种方法消去。注意由于是合同变换，我们还要把 r_3-2r_2。

再看 a_{1,4}，我们仍然用 c_2 来消，用 c_4-3c_2,r_4-3r_2。

类似地，我们用 c_5-6c_2,r_5-6r_2。

用 c_6-5c_2,r_6-5r_2。

我们忽略 r_2，直接看 r_3，由于 a_{3,4}=0，我们让 c_4\leftrightarrow c_5,r_4\leftrightarrow r_5，得：

注意此处我们对换了矩阵的两行列，相当于在 Q 上实行一次对换，故结果需要变号。不要误因为在原矩阵上实行了两次对换而忘记变号。

2. c_4+4c_1,r_4+4r_1
3. c_5+5c_1,r_5+5r_1
4. c_6+\tfrac{14}{3}c_1,r_6+\tfrac{14}{3}r_1

5. c_5-3c_3,r_5-3r_3
6. c_6+\tfrac{7}{9}c_3,r_6+\tfrac{7}{9}r_3

可以看到，对于我们希望知道的位置的元素，其数值是不变的。

对于这个矩阵，其行列式为 (-1)^1\times[1\times(-1)\times9\times(-9)\times6\times(-6)]=(1\times9\times6)^2，所以 \textup{Pf}(A)=1\times9\times6=54。

我们总结一下算法步骤：初始计数器 cnt=0，对于每一个奇数行 r_{2i-1}，若该行元素全部为 0，直接输出 0 并终止程序。否则，找到一个 j 使 a_{2i-1,j}\not=0 并对换 c_{i} 和 c_j，对换 r_{i} 和 r_j，若 j\not=2i 则令 cnt\leftarrow cnt+1。然后对于 j\in[2i+1,2n] 将 c_j-\tfrac{a_{2i-1,j}}{a_{2i-1,2i}}c_{2i}。最后输出 (-1)^{cnt}\prod\limits_{i=1}^n a_{2i-1,2i}。

时间复杂度 O(n^3)，空间复杂度 O(n^2)。

下面是你们最喜欢的代码：

#include<bits/stdc++.h>
using namespace std;
long long n,w,r,t,q,i,j,k,p=1e9+7,a[505][505];
int main(){
    cin.tie(0)->sync_with_stdio(0);
    cin>>n;
    for(i=1;i<n;i++)
        for(j=i+1;j<=n;j++){
            cin>>a[i][j];
            a[j][i]=-a[i][j];//补全反对称矩阵 
        }
    for(r=i=1;i<n;i+=2){
        for(j=i+1;j<=n&&!a[i][j];j++);//找到非0位置
        if(j>n)return cout<<0,0;//遇到全0行直接输出0
        swap(a[i+1],a[j]);//r[i+1] <-> r[j]
        for(k=1;k<=n;k++)swap(a[k][i+1],a[k][j]);//c[i+1] <-> c[j]
        if(j>i+1)r=-r;//结果取反
        for(j=i+2;j<=n;j++){
            t=a[i][i+1];
            q=a[i][j];
            for(k=p-2;k;k>>=1){//计算乘法逆元 
                if(k&1)q=q*t%p;
                t=t*t%p;
            }
            for(k=1;k<=n;k++){
                (a[j][k]-=a[i+1][k]*q)%=p;//c[j] - (a[i][j]/a[i][i+1])*c[i+1]
                (a[k][j]-=a[k][i+1]*q)%=p;//r[j] - (a[i][j]/a[i][i+1])*r[i+1]
            }
        }
        r=r*a[i][i+1]%p;//ans <- ans*a[i][i+1]
    }
    cout<<(r+p)%p;//为避免负数再加模取模 
}