P4504 题解

· · 题解

基础知识

基尔霍夫二定律

  1. 节点的流入电流等于流出电流

  2. 每个点有一个电势，两点间的电压等于电势差

欧姆定律：对于一个电阻来说，U=IR。

因此，一个一般性的方法是以 \varphi_i 为变量，构建一个 n 变量的关系。对内部节点有 I=0 即 \sum_{e}(\varphi_x-\varphi_u)/R=0，而对其中一个边界节点可以令 I=1 并求出 \varphi 相除，也可以直接令 \varphi=1，最后随意赋值另一个边界节点即可，因为电势是可以平移的。

但这样解方程是每次 O(n^3) 的，完全过不去ww

一个 O(nq) 做法的尝试

这说明，我们一定需要根据题述电阻网络的特性进行求解。一个很容易想出的结论是，如果一个电阻网络只有两个端点与外界相连，则可以等效为两点之间恰连一个电阻，该电阻和两段测得电阻大小相同。

进一步地，我们有如果只有三个端点与外界相连，则一定等效为两两之间恰连一个电阻，或者三个点分别连向一个额外点某个电阻，使得两两端点间测得电阻均相等，在这两种图之间的变换也被成为星-三变换，其中三取自三角形。

然而，当四端点情况出现时，星形图不再可行，不严谨地来说，这是由于星形图仅有四个电阻自由度，而四端点需要六个自由度，才能满足完美刻画的需求，否则一定存在某种电势关系无法满足，如果想办法塞进六个电阻，即使是简单的完全图情况，图变得丑陋且难以直接计算端对端电阻，更别说用方程反解出各个电阻。

回到本题，一个自然的思路是，将链拎起来，那么链以外的部分可以看做若干子树，只以三个点与这条链相连，因此完全可以看做一个等效“三”或“星”。逐级向上递推。将一个儿子合并到父亲，只需要先忽视父亲的所有结构只保留三个接线柱，再在此条件下分别求解父亲的接线柱间两两的电阻，最后以这样的等效“三”直接并联到原本被忽视的结构上即可，这里最大的难点是需要解出一个 s<->1,2,3 ; 1,2,3<->t ; 1<->2,2<->3,1<->3 之间均有电阻的网络在 s 与 t 之间的电阻值，此时设元求解即可，只是稍微复杂。

然而，更大的困难在于链上。从一端 x_A 开始的合并需要维护 x_A,u_1,u_2,u_3 的四端点结构，每次合并三个节点进入，而这样不仅需要执行 6 次上述大常数的解方程步骤（事实上新的三个节点中间的互相连边是无关紧要的，可以视为最后并联进四端点结构），还需要进行反解参数，恐怕难以通过本题，处理各种父子关系的代码量也比较让人望而却步。

一个刻画多端网络的改进方案

我们不再使用所谓的“端对端电阻”来刻画一个电阻网络。相反，我们令 G_{ij}=\left\{ \begin{aligned} -\sum\limits_k 1/R_{ik} &(i=j) \\ 1/R_{ij}&(i \neq j)\end{aligned} \right.，于是，令 I 列向量为 x 号点额外流出网络的电流，U 列向量为 x 号点的电势，则 GU=I。

此时，我们发现一个 n\times n 的对称矩阵足以刻画整个电阻网络，且这个矩阵的自由度恰好为 \binom n2。相当于此时我们用一种更具结构的工具重新刻画了多端网络。

引入 LCT

以上的工具给予我们能力处理一般的六端网络，于是我们用LCT的一个节点存储包括其所辖链和链内子树中所有点的六端网络矩阵，主要困难在于三端网络合并入六端网络、六端网络转三端网络和六端网络的串接。

不妨令靠“上”（离当前根近）的一个端定为矩阵 0\sim 2 的位置。

先讲三端网络合并入六端网络，只需将六阶矩阵对应位置直接加上三阶矩阵即可，因为加入新部分电流是叠加的。

再说六端网络转三端网络，这里并非直接取左上角，而是需要满足下端 I_2=0 的限制。不妨令矩阵分块为 \left(\begin{array}{c|c} A&B \\ \hline C&D \end{array}\right)，于是有 AU_1+BU_2=I_1,CU_1+DU_2=0，因此 U_2=-D^{-1}CU_1，故 I_2=(A+BD^{-1}C)U_1，令三端网络为 U_0 前系数即可。

最后阐述六端的串接，使 P 矩阵靠上，Q 靠下拼接，于是 P=\left(\begin{array}{c|c} A_p&B_p \\ \hline C_p&D_p \end{array}\right),Q=\left(\begin{array}{c|c} A_q&B_q \\ \hline C_q&D_q \end{array}\right)，且 A_pU_1+B_pU_2=I_1,C_pU_1+D_pU_2=I_2,A_qU_2+B_qU_3=-I_2,C_qU_2+D_qU_3=I_3，你想要用 U_1 和 U_3 表示出 I_1 与 I_3，因此应消去 U_2 与 I_2，故 U_2=-(D_p+A_q)^{-1}(C_pU_1+B_qU_3)，带入即可得到最终有 A=A_p-B_p(D_p+A_q)^{-1}C_p,B=-B_p(D_p+A_q)^{-1}B_q,C=-C_q(D_p+A_q)^{-1}C_p,D=D_q-C_q(D_p+A_q)^{-1}B_q。

于是，LCT只剩最后一些细节尚未处理，将此时的零号点初始化为 1/R_{i,i+3}=10^{10}，而其他 1/R=0 即可，而换根引起的区间翻转操作只需将该位置权值和六端答案的所有 i 与 i+3 交换即可。最后由于我们需要处理的是边权而不是点权，新建 n-1 个点将权值放到虚拟“边点”上即可解决。

最后一个问题：已知矩阵 G，如何求出两点间电阻？

注意到 I=GU，我们想要 U=G_vI，然而 U 并非 I 所能确定，还有上下移动的差值，因此我们从 G 中去掉一行一列（记得不要把 A 或 B 去掉了），此时 G 满秩，故对其求逆，求出 U=G_t^{-1}I，令 I_A=-1,I_B=1，\varphi_A-\varphi_B 即为所求。

代码如下

#include<cstdio>
#include<cmath>
#include<vector>
#include<algorithm>
const int N=10007,M=2*N;
typedef long double ld;
const ld ginf=1e10,eps=1e-9;
struct M33{ld a[3][3];bool Z;inline ld*operator[](int x){return a[x];}};
struct M66
{
    ld a[6][6];bool Z;
    inline ld*operator[](int x){return a[x];}
#define Ma(t,bd) a[t][bd+0],a[t][bd+1],a[t][bd+2]
#define Mb(ad,bd) {Ma(ad+0,bd),Ma(ad+1,bd),Ma(ad+2,bd)}
    inline M33 A(){return Mb(0,0);}
    inline M33 B(){return Mb(0,3);}
    inline M33 C(){return Mb(3,0);}
    inline M33 D(){return Mb(3,3);}
    inline void swp()
    {
        for(int i=0;i<3;++i)for(int j=0;j<3;++j)
        {
            std::swap(a[i][j],a[i+3][j+3]);
            std::swap(a[i+3][j],a[i][j+3]);
        }
    }
};
template<int n,class T>
T ginv(T x)
{
    T res{0};
    for(int i=0;i<n;++i)res[i][i]=1;
    for(int i=0;i<n;++i)
    {
        for(int j=i;j<n;++j)if(std::abs(x[j][i])>eps)
        {
            std::swap(x.a[i],x.a[j]);
            std::swap(res.a[i],res.a[j]);
            break;
        }
        assert(std::abs(x[i][i])>eps);
        for(int j=0;j<n;++j)if(j!=i)
        {
            ld cof=x[j][i]/x[i][i];
            for(int k=0;k<n;++k)x[j][k]-=cof*x[i][k],res[j][k]-=cof*res[i][k];
        }
    }
    for(int i=0;i<n;++i)
        for(int j=0;j<n;++j)
            res[i][j]/=x[i][i];
    return res;
}
M33 add(M33 x,M33 y)
{
#define Mc(a,b) x[a][b]+y[a][b]
#define Md(a) Mc(a,0),Mc(a,1),Mc(a,2)
    return {Md(0),Md(1),Md(2)};
}
M33 sub(M33 x,M33 y)
{
#define Me(a,b) x[a][b]-y[a][b]
#define Mf(a) Me(a,0),Me(a,1),Me(a,2)
    return {Mf(0),Mf(1),Mf(2)};
}
M33 mul(M33 x,M33 y)
{
    M33 res{0};
    for(int i=0;i<3;++i)
        for(int j=0;j<3;++j)if(x[i][j])
            for(int k=0;k<3;++k)
                res[i][k]+=x[i][j]*y[j][k];
    return res;
}
#define ginv3 ginv<3,M33>
M66 gv[M];int ch[M][2],fa[M];
M33 gs[M];M66 fs[M];bool srv[M];
inline bool nrt(int x){return x==ch[fa[x]][0]||x==ch[fa[x]][1];}
inline M33 subt(int x)
{
    if(fs[x].Z)return {0};
    M33 res=sub(fs[x].A(),mul(mul(fs[x].B(),ginv3(fs[x].D())),fs[x].C()));
    return res;
}
M66 meg(M66 p,M66 q)//p up q down
{
    if(p.Z)return q;
    if(q.Z)return p;
    M33 Ap=p.A(),Bp=p.B(),Cp=p.C(),Dp=p.D();
    M33 Aq=q.A(),Bq=q.B(),Cq=q.C(),Dq=q.D();
    M33 Is=ginv3(add(Dp,Aq)),U1s=mul(Is,Cp),U3s=mul(Is,Bq);

    M33 Aa=sub(Ap,mul(Bp,U1s)),Ba=mul(Bp,U3s);
    M33 Ca=mul(Cq,U1s),Da=sub(Dq,mul(Cq,U3s));

    M66 ans;ans.Z=0;
    for(int i=0;i<3;++i)for(int j=0;j<3;++j)ans[i+0][j+0]= Aa[i][j];
    for(int i=0;i<3;++i)for(int j=0;j<3;++j)ans[i+0][j+3]=-Ba[i][j];
    for(int i=0;i<3;++i)for(int j=0;j<3;++j)ans[i+3][j+0]=-Ca[i][j];
    for(int i=0;i<3;++i)for(int j=0;j<3;++j)ans[i+3][j+3]= Da[i][j];

    return ans;
}
bool tg[M];
inline void rev(int x)
{
    if(!x)return;
    std::swap(ch[x][0],ch[x][1]);tg[x]^=1;srv[x]^=1;
    gv[x].swp();fs[x].swp();
}
inline void push(int x){if(tg[x])rev(ch[x][0]),rev(ch[x][1]),tg[x]=0;}
void pull(int x)
{
    M66 t=fs[ch[x][1]];t.Z=0;
    for(int i=0;i<3;++i)for(int j=0;j<3;++j)t[i][j]+=gs[x][i][j];
    fs[x]=meg(fs[ch[x][0]],meg(gv[x],t));
}
void rot(int x)
{
    int y=fa[x],z=fa[y],k=(x==ch[y][1]),w=ch[x][!k];
    fa[x]=z;if(nrt(y))ch[z][ch[z][1]==y]=x;
    ch[y][k]=w;if(w)fa[w]=y;fa[y]=x;ch[x][!k]=y;pull(y);
}
int que[M];
void splay(int x)
{
    int y=x,z=0;que[++z]=y;
    while(nrt(y))que[++z]=y=fa[y];
    while(z)push(que[z--]);
    while(nrt(x))
    {
        y=fa[x];z=fa[y];
        if(nrt(y))rot((ch[z][1]==y)^(ch[y][1]==x)?x:y);
        rot(x);
    }
    pull(x);
}
void access(int x)
{
    for(int y=0;x;x=fa[y=x])
    {
        splay(x);
        gs[x]=add(gs[x],subt(ch[x][1]));
        ch[x][1]=y;
        gs[x]=sub(gs[x],subt(ch[x][1]));
        pull(x);
    }
}
void mkrt(int x){access(x);splay(x);rev(x);}
void upd(int x,int A,int B,int z)
{
    access(x);splay(x);
    if(srv[x])std::swap(A,B);
    ld D=z-gv[x][A+3][B];
    gv[x][A+3][B]+=D;gv[x][B][A+3]+=D;
    gv[x][A+3][A+3]-=D;gv[x][B][B]-=D;
    pull(x);
}
std::vector<int> to[M];
void dfs(int i){for(int v:to[i])dfs(v),gs[i]=add(gs[i],subt(v));pull(i);}
int main()
{
    int n;scanf("%d",&n);
    fs[0].Z=1;
    for(int i=0;i<3;++i)
    {
        fs[0][i][i+3]=fs[0][i+3][i]=ginf;
        fs[0][i][i]-=ginf;fs[0][i+3][i+3]-=ginf;
    }
    for(int i=0;i<=n;++i)gv[i]=*fs;
    for(int i=2;i<=n;++i)
    {
        scanf("%d",fa+i);fa[n+i]=fa[i];fa[i]=n+i;
        for(int A=0,z;A<3;++A)for(int B=0;B<3;++B)
        {
            scanf("%d",&z);
            gv[n+i][A+3][B]=z;gv[n+i][B][A+3]=z;
            gv[n+i][A+3][A+3]-=z;gv[n+i][B][B]-=z;
        }
    }
    for(int i=2;i<=n+n;++i)to[fa[i]].push_back(i);
    dfs(1);

    int m;scanf("%d",&m);
    while(m--)
    {
        int op,x,y,z,A,B;
        scanf("%d",&op);
        if(op==1)scanf("%d%d%d%d",&x,&A,&B,&z),upd(n+x,A-1,B-1,z);
        else
        {
            scanf("%d%d%d%d",&x,&A,&y,&B);--A;B+=2;
            mkrt(x);access(y);splay(y);
            M66 ans=fs[y];
            int del=0;while(del==A||del==B+3)++del;
            for(int i=0;i<6;++i)for(int j=0;j<6;++j)ans[i][j]=ans[i+(i>=del)][j+(j>=del)];
            M66 tans=ginv<5,M66>(ans);A-=(A>=del);B-=(B>=del);
            ld zans=(tans[B][A]-tans[B][B])-(tans[A][A]-tans[A][B]);
            printf("%Lf\n",zans);
        }
    }
    return 0;
}

致谢

感谢 Qingyu 提供的出题人彭天翼的题解