首页 > 其他分享 >Pfaffian and Determinant of $\sum A_iz^i$

Pfaffian and Determinant of $\sum A_iz^i$

时间:2023-12-11 16:44:59浏览次数:35  
标签:Determinant iz int sum Pfaffian lep mul prod buf

// created on 23.12.11

目录

Pfaffian and Determinant of \(\sum A_iz^i\)

求斜对称矩阵的 Pfaffian 即 \(\mathrm{Pf}(A)\),主要用到了性质 \(\mathrm{Pf}(A)|B|=\mathrm{Pf}(BAB^{T})\) 。

我们对矩阵进行消元。若第一列全 \(0\),则矩阵不满秩,\(\mathrm{Pf}(A)=0\) 。否则找到任一第一列元素非 \(0\) 的行,交换到第二行。同时执行 \(B^T\),将列交换。

我们用第二行第一列的非 \(0\) 元素将后面行的第一列全部消掉。在执行 \(B^T\) 时,由于斜对称性,也可以恰好消掉第一列所有元素。

接着用第一行第二列的非 \(0\) 元素将后面行的第二列全部消掉。在执行 \(B^T\) 时,由于斜对称性,也可以恰好消掉第二行所有元素。

此时通过 Pfaffian 的展开式 \(\mathrm{Pf}(A)=\sum\limits_{i\not=1}(-1)^{i}A_{1,i}\mathrm{Pf}(A\backslash\{1,i\})\) 有 \(\mathrm{Pf}(A)=A_{1,2}\mathrm{Pf}(A\backslash\{1,2\})\) 。于是我们得到了 Pfaffian 的 \(O(n^3)\) 解法。

// Pf (a), O (n ^ 3)
template <class T> inline int pfaffian (T a, int n) {
    int prod = 1;
    for (int h = 1; h < n; h += 2) {
        lep (i, h + 1, n) if (a[i][h]) {
            if (i != h + 1) {
                prod = mod - prod, swap (a[i], a[h + 1]);
                lep (j, h, n) swap ( a[j][i], a[j][h + 1] );
            }
            break ;
        }
        if (! a[h][h + 1]) return 0;
        prod = mul (prod, a[h][h + 1]);
        int buf = qpow (a[h][h + 1], mod - 2);
        lep (j, h, n) {
            a[h][j] = mul (a[h][j], buf);
            a[j][h] = mul (a[j][h], buf);
        }
        lep (j, h + 2, n) if (buf = a[j][h]) {
            lep (k, h, n) pls (a[j][k], mul (a[h + 1][k], buf));
            lep (k, h, n) pls (a[k][j], mul (a[k][h + 1], buf));
        }
        lep (j, h + 2, n) if (buf = a[j][h + 1]) {
            lep (k, h, n) sub (a[j][k], mul (a[h][k], buf));
            lep (k, h, n) sub (a[k][j], mul (a[k][h], buf));
        }
    }
    return prod;
}

接着考虑 \(\mathrm{Pf}(\sum A_iz^i)\) 的问题。在 网格图最大流计数 中,\([z^0]\mathrm{Pf}(\sum A_iz^i)=1\) 时必然存在的。于是问题实际转化为求解 \(\det(\sum A_iz_i)\) 。

实际上,只需要通过消元找到 \(MA_m=I\),就可以使得 \(\det(\sum\limits_{i=0}^{m}A_iz^i)\) 转化为 \(\det(M^{-1})\det(M\sum\limits_{i=0}^{m-1}A_iz^i+Iz^m)\) 。在这个过程中,如果找不到某个 \(A_{m,i,i}\not=0\),则将这一列乘 \(z\) 继续考虑。如果最后成为全 \(0\) 列意味着不满秩,原式为 \(0\) 。

最后问题就转化为 \(\det(\sum\limits_{i=0}^{m}A_iz^i+Iz^m)\) 形式的问题。我们设计分块矩阵 \(\begin{bmatrix}Iz& -I& 0&\cdots& 0&0\\0& Iz& -I&\cdots& 0&0\\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&\cdots&Iz&-I\\A_0&A_1&A_2&\cdots& A_{m-2}&A_{m-1}+Iz\end{bmatrix}\),其中 \(I\) 和 \(A_i\) 均为 \(n\) 阶方阵。由于 \(Iz\) 有逆,简单归纳可以发现以上矩阵的行列式就是原问题的答案。而做法是 trivial 的:先消元得到上海森堡矩阵,然后递推即可。这个过程和求解特征多项式是一模一样的。

于是我们得到了 \(\det(\sum\limits_{i=0}^{m}A_iz^i)\) 的 \(O((mn)^3)\) 解法。

// det (xI + A), O (n ^ 3)
template <class T> inline vector <int> char_poly (T a, int n) {
    lep (i, 1, n - 1) {
        int t = i + 1; for ( ; t <= n && ! a[t][i]; ++ t);
        if (t > n) continue ;
        if (t != i + 1) {
            lep (j, i, n) swap (a[t][j], a[i + 1][j]);
            lep (j, 1, n) swap (a[j][t], a[j][i + 1]);
        }
        int inv = qpow (a[i + 1][i], mod - 2);
        lep (j, i + 2, n) if (a[j][i]) {
            int buf = mul (a[j][i], inv);
            lep (k, i, n) sub (a[j][k], mul (a[i + 1][k], buf));
            lep (k, 1, n) pls (a[k][i + 1], mul (a[k][j], buf));
        }
    }
    vector <vector <int> > p (n + 2, vector <int> (n + 1, 0) );
    p[n + 1][0] = 1;
    rep (i, n, 1) {
        lep (j, 0, n + 1 - i) p[i][j] = j ? p[i + 1][j - 1] : 0;
		for (int j = i, t = 1; ; t = mul (t, a[j + 1][j]), ++ j) {
			int coef = mul (mul (t, a[i][j]), (((j - i) & 1) ? mod - 1 : 1));
			lep (k, 0, n - j) pls (p[i][k], mul (p[j + 1][k], coef));
            if (j == n) break ;
		}
    }
    return p[1];
}
// det (a0 + a1 x + a2 x ^ 2 + ... ak x ^ k), O ((nk) ^ 3)
template <class T> inline vector <int> det (T a, int k, int n) {
    int prod = 1, offset = 0;
    lep (h, 1, n) {
        lep (i, 1, h - 1) {
            int buf = a[k][i][h];
            lep (d, 0, k) lep (j, 1, n) {
                sub (a[d][j][h], mul (a[d][j][i], buf));
            }
        }
        for ( ; ; ) {
            if (a[k][h][h]) break ;
            lep (i, h + 1, n) if (a[k][i][h]) {
                prod = mod - prod;
                lep (d, 0, k) lep (j, 1, n) swap (a[d][h][j], a[d][i][j]);
            }
            if (a[k][h][h]) break ;
            if ( ++ offset > k * n ) return vector <int> (k * n + 1, 0);

            rep (d, k, 1) lep (j, 1, n) a[d][j][h] = a[d - 1][j][h];
            lep (j, 1, n) a[0][j][h] = 0;

            lep (i, 1, h - 1) {
                int buf = a[k][i][h];
                lep (d, 0, k) lep (j, 1, n) {
                    sub (a[d][j][h], mul (a[d][j][i], buf));
                }
            }
        }
        if (! a[k][h][h]) return vector <int> (k * n + 1, 0);
        prod = mul (prod, a[k][h][h]);
        int inv = qpow (a[k][h][h], mod - 2);
        lep (d, 0, k) lep (j, 1, n) a[d][h][j] = mul (a[d][h][j], inv);
        lep (i, 1, n) if (i != h) {
            int buf = a[k][i][h];
            lep (d, 0, k) lep (j, 1, n) {
                sub (a[d][i][j], mul (a[d][h][j], buf));
            }
        }
    }

    vector <vector <int> > b (k * n + 1, vector <int> (k * n + 1, 0) );
    lep (i, 1, (k - 1) * n) b[i][n + i] = mod - 1;
    lep (d, 0, k - 1) {
        lep (i, 1, n) lep (j, 1, n) {
            b[(k - 1) * n + i][d * n + j] = a[d][i][j];
        }
    }
    auto ret = char_poly (b, k * n);
    lep (i, offset, k * n) ret[i - offset] = mul (ret[i], prod);
    lep (i, k * n - offset + 1, k * n) ret[i] = 0;
    return ret;
}

标签:Determinant,iz,int,sum,Pfaffian,lep,mul,prod,buf
From: https://www.cnblogs.com/qiulyqwq/p/17894761.html

相关文章

  • FAILED: ParseException line 1:65 cannot recognize input near 'row' 'formatted' &
    hive报FAILED:ParseExceptionline1:65cannotrecognizeinputnear'row''formatted''delimited'intablerowformatspecification错误语句:insertoverwritelocaldirectory'/home/ljpbd/datas/student'rowformatteddel......
  • FAILED: ParseException line 1:17 cannot recognize input near 'student2' 'select'
     hive向表中插入数据时报错:FAILED:ParseExceptionline1:17cannotrecognizeinputnear'student2''select''id'indestinationspecification错误:insertoverwritestudent2selectid,namefromstudent;正确:insertoverwritetablest......
  • ApplicationContextInitializer在Spring容器执行refresh之前执行
    ApplicationContextInitializer用于在刷新Spring容器之前的回调接口。ApplicationContextInitializer是Spring框架原有的概念,这个类的主要目的就是在ConfigurableApplicationContext类型(或者子类型)的ApplicationContext进行刷新refresh之前,允许我们对ConfigurableApplicatio......
  • CF1777C Quiz Master 题解
    题意:思路:由于需要维护极差,因此将$a$排序;由于相同的数对因子的种类和极差的贡献重复,因此将$a$去重。设满足条件且极差最小的方案为:$a_1$,$a_3$,$a_4$,$a_7$,$a_9$,该方案等价于区间$[1,9]$。因此,满足条件且极差最小的方案一定为一个连续的区间$[l,r......
  • [ARC164E] Segment-Tree Optimization 题解
    题目链接题目链接题目解法一个自认为比较自然的解法这种一段序列切成两部分的问题首先考虑区间\(dp\)令\(f_{l,r}\)为\([l,r]\)能构成的最小深度,\(g_{l,r}\)为在\(f_{l,r}\)最小的情况下最少的最大深度的点的个数转移枚举\(k\)即可需要在下面是叶子的时候处理一......
  • Java synchronized
    synchronized是Java中最基本的线程同步机制之一,通过在方法或代码块上添加synchronized关键字,可以确保只有一个线程可以访问该方法或代码块。它是Java中实现线程安全的重要机制之一。synchronized关键字的使用方式有两种:1、修饰实例方法当synchronized关键字修饰一个实例方法时,......
  • Java synchronized 、ReentrantLock和Semaphore
    synchronized在Java中,使用synchronized关键字可以实现对代码块或方法的同步访问,以确保多个线程不会同时访问共享资源。当一个线程获取了对象的锁(即进入了synchronized代码块),其他线程如果也希望获取该对象的锁,它们将被阻塞,直到拥有锁的线程执行完毕并释放锁。因此,在某种意义上,使......
  • Graph regularized non-negative matrix factorization with prior knowledge consist
    Graphregularizednon-negativematrixfactorizationwithpriorknowledgeconsistencyconstraintfordrug-targetinteractionspredictionJunjunZhang 1, MinzhuXie 2 3Affiliations expandPMID: 36581822 PMCID: PMC9798666 DOI: 10.1186/s1285......
  • Predict potential miRNA-disease associations based on bounded nuclear norm regul
    PredictpotentialmiRNA-diseaseassociationsbasedonboundednuclearnormregularizationYidongRao 1, MinzhuXie 1, HaoWang 1Affiliations expandPMID: 36072658 PMCID: PMC9441603 DOI: 10.3389/fgene.2022.978975 SigninFreePMCa......
  • Graph regularized non-negative matrix factorization with [Formula: see text] nor
    Graphregularizednon-negativematrixfactorizationwith[Formula:seetext]normregularizationtermsfordrug-targetinteractionspredictionJunjunZhang 1, MinzhuXie 2 3Affiliations expandPMID: 37789278 PMCID: PMC10548602 DOI: 10.11......