自动搬运

查看原文

来自洛谷，原作者为

	Edward1002001 这个家伙很菜，什么也没有留下

搬运于2025-08-24 22:00:49，当前版本为作者最后更新于2025-07-21 23:16:58，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

基础知识

基尔霍夫二定律

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 提供的出题人彭天翼的题解