【学习笔记】[PA 2022] Drzewa rozpinające

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

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

非常直白的定理。证明基于一个恒等式: [ I 0 v T 1 ] [ I + u v T u 0 1 ] [ I 0 − v T 1 ] = [ I u 0 1 + v T 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} [IvT01][I+uvT0u1][IvT01]=[I0u1+vTu]。这个不难通过手算验证,顺便说一下,把分块矩阵的每个部分看成一个整体来算就好,这和普通的矩阵乘法没有什么区别。两边同时取行列式,就能得到 det ( I + u v T ) = 1 + v T u \text{det}(I+uv^T)=1+v^Tu det(I+uvT)=1+vTu。那么 L H S = det ( A ( I + u v T A − 1 ) ) = det ( A ) det ( I + u ( v T A − 1 ) ) = det ( A ) ( 1 + v T A − 1 u ) 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) LHS=det(A(I+uvTA1))=det(A)det(I+u(vTA1))=det(A)(1+vTA1u),这样就证完了,其实也并不复杂。

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

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

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

细节调了半天,还是太菜了。

#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]); } } //fixed //joker 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; }

__EOF__

本文作者仰望星空的蚂蚁
本文链接https://www.cnblogs.com/cqbzly/p/17529951.html
关于博主:评论和私信会在第一时间回复。或者直接私信我。
版权声明:本博客所有文章除特别声明外,均采用 BY-NC-SA 许可协议。转载请注明出处!
声援博主:如果您觉得文章对您有帮助,可以点击文章右下角推荐一下。您的鼓励是博主的最大动力!
posted @   仰望星空的蚂蚁  阅读(52)  评论(1编辑  收藏  举报  
点击右上角即可分享
微信分享提示