【学习笔记】[PA 2022] Drzewa rozpinające-Toy模板网

这篇具有很好参考价值的文章主要介绍了【学习笔记】[PA 2022] Drzewa rozpinające。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

单纯只是为了记录一下这个 $\text{trick}$ 。而且涉及到一点分块矩阵的思想，感觉还挺有意思的。

矩阵行列式引理：设 $A$ 为可逆矩阵， $u, v$ 为列向量，则有： $\text{det}(A+uv^T)=\text{det}(A)(1+v^TA^{-1}u)$ 。

非常直白的定理。证明基于一个恒等式： $\begin{bmatrix}I&0\\ v^T&1\end{bmatrix}\begin{bmatrix}I+uv^T&u\\0&1\end{bmatrix}\begin{bmatrix}I&0\\-v^T&1\end{bmatrix}=\begin{bmatrix}I&u\\0&1+v^Tu\end{bmatrix}$ 。这个不难通过手算验证，顺便说一下，把分块矩阵的每个部分看成一个整体来算就好，这和普通的矩阵乘法没有什么区别。两边同时取行列式，就能得到 $\text{det}(I+uv^T)=1+v^Tu$ 。那么 $LHS=\text{det}(A(I+uv^TA^{-1}))=\text{det}(A)\text{det}(I+u(v^TA^{-1}))=\text{det}(A)(1+v^TA^{-1}u)$ ，这样就证完了，其实也并不复杂。

值得一提的是，如果 $u, v$ 不是向量而是两个矩阵，那么根据分块矩阵的思想上述式子仍然成立，只不过要稍微改动一下： $\text{det}(A+uv^T)=\text{det}(A)\text{det}(I+v^TA^{-1}u)$ 。同样只需要 $\text{Laplace}$ 展开一下即可。

注意看这个公式把 $u,v^T$ 的顺序调换了，这样乘出来的矩阵会发生非常大的变化。比如说这道题，非常熟悉的 $\text{gcd}$ 矩阵和矩阵树定理，考虑构造 $n\times m$ 的矩阵 $U, V$ ，其中 $U_{i,d}=\phi(d)[d|a_i]$ ， $V_{i,d}=[d|a_i]$ ，度数矩阵 $D_{i,i}=\sum \gcd(a_i,a_j)$ ，邻接矩阵 $G_{i,j}=\gcd(a_i,a_j)$ ，那么有 $G=UV^T$ ，根据矩阵树定理我们要求的是 $(n-1)\times (n-1)$ 的代数余子式 $\text{det}(D-G)$ 。

注意到 $LHS=\text{det}(D-UV^T)=\text{det(D)}\text{det}(I_m-V^TD^{-1}U)$ ，记矩阵 $Q=I_m-V^TD^{-1}U$ ，这个矩阵只在 $\text{lcm}(i,j)\le m$ 的地方有值，因此是稀疏的，但是这个稀疏是一个很模糊的概念，而且也没有一个固定的算法。本题当中倒着高消会跑的比较快，大概是因为行的编号越大矩阵越稀疏吧。

细节调了半天，还是太菜了。文章来源地址https://www.toymoban.com/news/detail-527157.html

#include<bits/stdc++.h>
#define ll long long
#define fi first
#define se second
#define pb push_back
#define db double
#define inf 0x3f3f3f3f3f3f3f3f
using namespace std;
const int mod=1e9+7;
int n,m,a[5005],phi[5005],cnt,prime[5005],vis[5005],D[5005],invD[5005],gcd[5005][5005];
ll mat[5005][5005],sum[5005],res=1;
ll fpow(ll x,ll y=mod-2){
    ll z(1);
    for(;y;y>>=1){
        if(y&1)z=z*x%mod;
        x=x*x%mod;
    }
    return z;
}
void init(int n){
    phi[1]=1;
    for(int i=2;i<=n;i++){
        if(!vis[i]){
            prime[++cnt]=i,phi[i]=i-1;
        }
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
            vis[i*prime[j]]=1;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            else{
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
            }
        }
    }
}
void add(ll &x,ll y){x=(x+y)%mod;}
ll det(){
    int rev=1;
    for(int i=m;i;i--){
        if(!mat[i][i]){
            for(int j=i-1;j;j--){
                if(mat[j][i]){
                    swap(mat[i],mat[j]);
                    rev*=-1;
                    break;
                }
            }
        }
        vector<int>pos;for(int j=i;j;j--)if(mat[i][j])pos.pb(j);
        ll inv=fpow(mat[i][i]);
        for(int j=i-1;j;j--){
            if(mat[j][i]){
                ll tmp=mat[j][i]*inv%mod;
                for(auto p:pos){
                    mat[j][p]=(mat[j][p]-mat[i][p]*tmp)%mod;
                }
            }
            
        }
    }
    ll res=1;for(int i=1;i<=m;i++)res=res*mat[i][i]%mod;
    if(rev==1)return res;
    return mod-res;
}
int getgcd(int x,int y){
    if(x<y)swap(x,y);
    return gcd[x][y];
}
signed main(){
    ios::sync_with_stdio(false);
    cin.tie(0),cout.tie(0);
    cin>>n;for(int i=1;i<=n;i++)cin>>a[i],m=max(m,a[i]);
    init(m);
    for(int i=0;i<=m;i++){
        for(int j=0;j<=i;j++){
            if(i==0||j==0)gcd[i][j]=i+j;
            else gcd[i][j]=gcd[j][i%j];
        }
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            D[i]+=getgcd(a[i],a[j]);
        }
    }
    n--;
    for(int i=1;i<=n;i++)invD[i]=fpow(D[i]);
    for(int i=1;i<=n;i++)add(sum[a[i]],invD[i]);
    for(int i=1;i<=m;i++){
        for(int j=i+i;j<=m;j+=i){
            add(sum[i],sum[j]);
        }
    }
    for(int i=1;i<=m;i++)mat[i][i]=1;
    for(int i=1;i<=m;i++){
        for(int j=1;j<=m;j++){
            if(i/getgcd(i,j)*j<=m)add(mat[i][j],-phi[j]*sum[i/getgcd(i,j)*j]);
        }
    }
    for(int i=1;i<=n;i++)res=res*D[i]%mod;
    res=res*det()%mod;
    cout<<(res+mod)%mod;
}