注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

vfleaking的博客

My name is VFlea King

 
 
 

日志

 
 

3556: [Ctsc2014] 插线板 题解  

2014-05-21 16:22:12|  分类: CTSC |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
【题目大意】
说从前有一棵树,树上每个结点里都含有三个小节点,每条边对应着3*3根电阻不可忽略的导线,连接两端小节点。
每次可以进行如下两个操作之一:
1. 修改一条边上某根导线的电阻值。
2. 询问给定两个小节点问之间的等效电阻。

n <= 10^4, q<=10^4。
电阻均不为0。

【题解】
orz pty……
这题还是很炫酷的233。

首先暴力高斯消元是显然的。
用phi(v)表示v的电势,i(v)表示流出v的电流。
用r(v, u)表示v、u间的电阻。
我们可以列出一个比较显然的方程:
-----
 \     phi(v) - phi(u)
  |   ----------------- = i(v)
 /         r(v, u)
-----
  u
因为欧姆定律嘛啦啦啦。
由于电流平衡,所以有了这个玩意儿我们就能测电阻了。
假设要测电阻的两端为qv, qu。
我们就往他们两端通入1A的电流,测电压就行了。
直接令i(qv) = 1A, i(qu) = -1A,对于其它的v,i(v) = 0A,然后高斯消元,得:
(phi(qv) - phi(qu)) / 1A
即电阻。
当然你用脚丫子想都知道你需要把其中某个点的电势置为0……因为只有电势差有意义。

刚刚那个方程可以写成:
G phi = i
其中G(v, u) =
| v != u:        1
|          - ---------
|             r(v, u)
| v == u:  -----
|           \        1
|            |   ---------
|           /     r(v, u)
|          -----
|            u
然后G居然莫名其妙还有个名字……叫基尔霍夫矩阵。
等等,基尔霍夫矩阵不是生成树计数的那玩意儿吗?

我们设g(v, u) = 1 / r(v, u)   (即v、u间的电导)
那么我们可以把G写得更漂亮:
其中G(v, u) =
| v != u: -g(v, u)
| v == u: sum(g(v, u))
没错!我们熟悉的生成树计数又回来了!
最美妙的时刻莫过于发现一直以来在研究的两种事物其实是一种事物。

不过我很好奇基尔霍夫到底是玩电路的时候发现的这个还是玩生成树计数的时候发现的这个。

不过这么暴力显然会TLE。

我会Δ-Y变换我自豪!
考试时去玩等效电阻去了……
但是没玩出来……

听到题解后吓傻了——直接合并基尔霍夫矩阵就行了囧。

我们来看这玩意儿怎么玩……

3556: [Ctsc2014] 插线板 - vfleaking - vfleaking的博客
 
如图所示……
所以我们就可以一拍脑袋说:只要维护了子树的电路信息然后再在链上搞搞就好了!
= =……所以下面都是怎么合并的故事……

我们定义电路类型Cir1和Cir2。
Cir1有3个端口,Cir2有6个端口。为了方便起见,所有下面表示Cir1的变量为小写,Cir2的变量首字母大写。
我们要记录的是端口电流对端口电势的相应。也就是说端口电流与电势间的函数关系。(假设内部的所有节点的电流平衡)
显然保留一个基尔霍夫矩阵A表示 A phi = i即可。
3556: [Ctsc2014] 插线板 - vfleaking - vfleaking的博客
然后我们手算一下这样几种合并的公式:(下图中的线表示电阻为0的导线)
3556: [Ctsc2014] 插线板 题解 - vfleaking - vfleaking的博客
注意我们可以类似地定义c +a C表示在左边接。当然我们也可以用C +b c表示在右边接。

这样子的话,如果要查询v到u之间的电阻,就是要查询v到u之间的Cir2信息,记为Q(v, u)。那么我们就要统计v到u之间的路径p1, p2, p3, ..., pm上的一些信息。
记相邻两点v, u之间的那个3*3的网络的Cir2为E(v, u)。
但是这些是不够的,我们还需要记down(v)为以v为根时,v上的三个端口电流与电势的关系,用一个Cir1表示。
再来点辅助的:若v、f相邻,记down(v, f)为以v为根时,且无视掉f这棵子树,v上的三个端口电流与电势的关系,用一个Cir1表示。
若v、f1相邻,v、f2相邻,记down(v, f1, f2)为以v为根时,且无视掉f1和f2这两棵子树,v上的三个端口电流与电势的关系,用一个Cir1表示。

这样就能写出Q(v, u)的表达式:(注意这里省略括号,例如a + b + c + d表示((a + b) + c) + d)
Q(v, u) = down(p1, p2) +A E(p1, p2) +A down(p2, p1, p3) +B E(p2, p3) +A down(p3, p2, p4) +B E(p3, p4) ... +A down(p[i], p[i - 1], p[i + 1]) +B E(p[i], p[i + 1]) ... +A down(p[m - 1], p[m - 2], p[m]) +B E(p[m - 1], p[m]) +A down(p[m], p[m - 1])
(好一个眼花缭乱……简单来说就是每次接上一个E,再在右端多接个down)
对于任意一个v,怎么求down(v)呢?假设v的儿子分别为u1, u2, ..., um,则:
down(v) = (E(v, u1) +b down(u1, v)) +a (E(v, u2) +b down(u2, v)) ... +a (E(v, uk) +b down(uk, v)) ... +a (E(v, um) +b down(um, v))
我们注意到+a是有逆运算的,我们记为-a
那么:
down(v, f) = down(v) -a (E(v, f) +b down(f, v))
down(v, f1, f2) = down(v) -a (E(v, f1) +b down(f1, v)) -a (E(v, f2) +b down(f2, v))

插播一个搞笑算法:我们就可以强制把1作为根,算出每个v的父亲fa[v],算出down(v, fa[v]),然后每次更新一条边时暴力更新此处到根的所有down。每次查询时暴力找路径,通过down(p[k], fa[p[k]])算出down(v),进而算出所需要的down(p[k], p[k - 1], p[k + 1])。这样居然能跑得很快……最慢的点大概5s……然后貌似已经能AC了……

还是来说树链剖分怎么做吧……
其实也没啥好说的……首先以1为根,把大家剖成链……
算出v的父亲fa[v]以及喜爱的儿子prefer[v]。
原子信息一个奇葩…… down(v, fa[v], prefer[v]) +A E(v, fa[v])
特别地,链底端的原子信息为down(v, fa[v]) +A E(v, fa[v])。
用线段树维护一条链上区间的所有原子信息按顺序 +B 起来的结果。然后就能通过一系列的线段树查询搞出Q(v, u)。
修改的时候也只会修改O(log n)个结点的原子信息。

细节随便搞吧……感觉地球人都知道怎么维护。
详细见代码吧。
mergeA是+A, mergeB是+a,mergeC是+B。
没有去定义+b。但是定义了一个closeL表示把Cir2的左端电势强制设为0从而变为一个Cir1。
faCir[v]即E(v, fa[v])

#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <vector>
using namespace std;

const int MaxN = 10000;

const int ND = 3;

inline int getint()
{
char c;
while (c = getchar(), ('0' > c || c > '9') && c != '-');

if (c != '-')
{
int res = c - '0';
while (c = getchar(), '0' <= c && c <= '9')
res = res * 10 + c - '0';
return res;
}
else
{
int res = 0;
while (c = getchar(), '0' <= c && c <= '9')
res = res * 10 + c - '0';
return -res;
}
}

int sgn(double a)
{
const double Eps = 1e-7;
if (a > Eps)
return 1;
else if (a < -Eps)
return -1;
else
return 0;
}

struct Matrix
{
double a[ND][ND];

friend inline Matrix operator+(const Matrix &lhs, const Matrix &rhs)
{
Matrix res;
for (int i = 0; i < ND; i++)
for (int j = 0; j < ND; j++)
res.a[i][j] = lhs.a[i][j] + rhs.a[i][j];
return res;
}
friend inline Matrix operator-(const Matrix &lhs, const Matrix &rhs)
{
Matrix res;
for (int i = 0; i < ND; i++)
for (int j = 0; j < ND; j++)
res.a[i][j] = lhs.a[i][j] - rhs.a[i][j];
return res;
}

inline Matrix operator-() const
{
Matrix res;
for (int i = 0; i < ND; i++)
for (int j = 0; j < ND; j++)
res.a[i][j] = -a[i][j];
return res;
}

friend inline Matrix operator*(const Matrix &lhs, const Matrix &rhs)
{
Matrix res;
for (int i = 0; i < ND; i++)
for (int j = 0; j < ND; j++)
{
res.a[i][j] = 0;
for (int k = 0; k < ND; k++)
res.a[i][j] += lhs.a[i][k] * rhs.a[k][j];
}
return res;
}

inline Matrix getRev() const
{
double ta[ND][ND];
for (int i = 0; i < ND; i++)
for (int j = 0; j < ND; j++)
ta[i][j] = a[i][j];
Matrix res;
for (int i = 0; i < ND; i++)
for (int j = 0; j < ND; j++)
res.a[i][j] = i == j ? 1 : 0;
for (int i = 0; i < ND; i++)
for (int j = i + 1;j < ND; j++)
if (sgn(ta[j][i]) != 0)
{
double z = ta[j][i] / ta[i][i];
for (int k = i; k < ND; k++)
ta[j][k] -= ta[i][k] * z;
for (int k = 0; k < ND; k++)
res.a[j][k] -= res.a[i][k] * z;
}
for (int i = ND - 1; i >= 0; i--)
{
for (int j = i + 1; j < ND; j++)
for (int k = 0; k < ND; k++)
res.a[i][k] -= ta[i][j] * res.a[j][k];
for (int k = 0; k < ND; k++)
res.a[i][k] /= ta[i][i];
}
return res;
}

friend inline ostream& operator<<(ostream &out, const Matrix &a)
{
for (int i = 0; i < ND; i++)
{
for (int j = 0; j < ND; j++)
out << a.a[i][j] << " ";
out << endl;
}
return out;
}
};

const Matrix MatI = (Matrix){
{
{1, 0, 0},
{0, 1, 0},
{0, 0, 1}
}
};
const Matrix Mat0 = (Matrix){
{
{0, 0, 0},
{0, 0, 0},
{0, 0, 0}
}
};

struct Circuit_1
{
Matrix matA;

Circuit_1(){}
Circuit_1(const Matrix &_matA)
: matA(_matA){}
};
struct Circuit_2
{
bool isEqual;
Matrix matA, matB;
Matrix matC, matD;

Circuit_2(){}
Circuit_2(const Matrix &_matA)
: isEqual(true), matA(_matA){}
Circuit_2(
const Matrix &_matA, const Matrix &_matB,
const Matrix &_matC, const Matrix &_matD)
: isEqual(false), matA(_matA), matB(_matB), matC(_matC), matD(_matD){}

Circuit_2 getRev() const
{
if (isEqual)
return Circuit_2(matA);
else
return Circuit_2(
matD, matC,
matB, matA);
}
};

inline Circuit_2 mergeA(const Circuit_2 &cirL, const Circuit_1 &cirR)
{
if (cirL.isEqual)
return Circuit_2(cirL.matA + cirR.matA);
else
return Circuit_2(
cirL.matA, cirL.matB,
cirL.matC, cirL.matD + cirR.matA);
}
inline Circuit_2 mergeA(const Circuit_1 &cirL, const Circuit_2 &cirR)
{
if (cirR.isEqual)
return Circuit_2(cirL.matA + cirR.matA);
else
return Circuit_2(
cirR.matA + cirL.matA, cirR.matB,
cirR.matC, cirR.matD);
}
inline Circuit_2 mergeB(const Circuit_2 &cirL, const Circuit_2 &cirR)
{
if (cirL.isEqual)
return Circuit_2(
cirL.matA + cirR.matA, cirR.matB,
cirR.matC , cirR.matD);
else if (cirR.isEqual)
return Circuit_2(
cirL.matA, cirL.matB,
cirL.matC, cirL.matD + cirR.matA);
else
{
Matrix matMR = (cirL.matD + cirR.matA).getRev();
return Circuit_2(
cirL.matA - cirL.matB * matMR * cirL.matC, - cirL.matB * matMR * cirR.matB,
- cirR.matC * matMR * cirL.matC, cirR.matD - cirR.matC * matMR * cirR.matB);
}
}
inline Circuit_1 mergeC(const Circuit_1 &cirL, const Circuit_1 &cirR)
{
return Circuit_1(cirL.matA + cirR.matA);
}
inline Circuit_1 subC(const Circuit_1 &cirL, const Circuit_1 &cirR)
{
return Circuit_1(cirL.matA - cirR.matA);
}
inline Circuit_1 closeL(const Circuit_2 &cir)
{
Matrix matAR = cir.matA.getRev();
return Circuit_1(cir.matD - cir.matC * matAR * cir.matB);
}

struct halfEdge
{
int u;
halfEdge *next;
};
halfEdge adj_pool[(MaxN - 1) * 2], *adj_tail = adj_pool;

int n;
int fa[MaxN + 1];
int faR[MaxN + 1][ND][ND];
halfEdge *adj[MaxN + 1];
Circuit_2 faCir[MaxN + 1];

void addEdge(int v, int u)
{
adj_tail->u = u, adj_tail->next = adj[v], adj[v] = adj_tail++;
}

int bfsSeq[MaxN];

int weight[MaxN + 1];
int depth[MaxN + 1];
int preferredChild[MaxN + 1];
bool isPreferred[MaxN + 1];
int trSeq[MaxN];
int trPos[MaxN + 1];
int chainL[MaxN + 1], chainR[MaxN + 1];

Circuit_1 down[MaxN + 1];
Circuit_1 downD[MaxN + 1];

Circuit_2 seg[MaxN * 4 + 1];

void seg_build(int p, int pL, int pR)
{
if (pL + 1 == pR)
seg[p] = mergeA(down[trSeq[pL]], faCir[trSeq[pL]]);
else
{
int pMid = (pL + pR) / 2;
seg_build(p << 1 | 0, pL, pMid);
seg_build(p << 1 | 1, pMid, pR);
seg[p] = mergeB(seg[p << 1 | 0], seg[p << 1 | 1]);
}
}
void seg_change(int p, int pL, int pR, int q)
{
if (pL + 1 == pR)
{
seg[p] = mergeA(down[trSeq[pL]], faCir[trSeq[pL]]);
return;
}
int pMid = (pL + pR) / 2;
if (q < pMid)
seg_change(p << 1 | 0, pL, pMid, q);
else
seg_change(p << 1 | 1, pMid, pR, q);
seg[p] = mergeB(seg[p << 1 | 0], seg[p << 1 | 1]);
}
Circuit_2 seg_query(int p, int pL, int pR, int qL, int qR)
{
if (qL <= pL && pR <= qR)
return seg[p];
int pMid = (pL + pR) / 2;
Circuit_2 res(Mat0);
if (qL < pMid)
res = mergeB(res, seg_query(p << 1 | 0, pL, pMid, qL, qR));
if (pMid < qR)
res = mergeB(res, seg_query(p << 1 | 1, pMid, pR, qL, qR));
return res;
}

void init()
{
int bfsSeq_n = 0;
bfsSeq[bfsSeq_n++] = 1;
for (int i = 0; i < bfsSeq_n; i++)
{
int v = bfsSeq[i];
for (halfEdge *e = adj[v]; e; e = e->next)
bfsSeq[bfsSeq_n++] = e->u;
}

depth[0] = 0;
for (int i = 0; i < n; i++)
{
int v = bfsSeq[i];
depth[v] = depth[fa[v]] + 1;
}

for (int i = n - 1; i >= 0; i--)
{
int v = bfsSeq[i];
weight[v] = 1;
for (halfEdge *e = adj[v]; e; e = e->next)
weight[v] += weight[e->u] + 1;

int u = 0;
for (halfEdge *e = adj[v]; e; e = e->next)
if (!u || weight[e->u] > weight[u])
u = e->u;
preferredChild[v] = u;
if (u)
isPreferred[u] = true;
}

int trSeq_n = 0;
for (int v = 1; v <= n; v++)
if (!adj[v])
{
int u = v;
for (u = v; isPreferred[u]; u = fa[u])
trSeq[trPos[u] = trSeq_n++] = u;
trSeq[trPos[u] = trSeq_n++] = u;

for (int i = trPos[v]; i < trSeq_n; i++)
chainL[trSeq[i]] = trPos[v], chainR[trSeq[i]] = trSeq_n;
}

downD[0] = Mat0;
for (int i = n - 1; i >= 0; i--)
{
int v = bfsSeq[i];
down[v] = Mat0;
for (halfEdge *e = adj[v]; e; e = e->next)
if (e->u != preferredChild[v])
down[v] = mergeC(down[v], downD[e->u]);
if (fa[v] && trPos[v] == chainR[v] - 1)
{
downD[v] = Mat0;
Circuit_2 cur(Mat0);
for (int i = chainL[v]; i < chainR[v]; i++)
cur = mergeB(cur, mergeA(down[trSeq[i]], faCir[trSeq[i]]));
downD[v] = closeL(cur);
}
}
seg_build(1, 0, n);
}

Circuit_2 tr_query(int qL, int qR)
{
if (qL == qR)
return Mat0;
return seg_query(1, 0, n, qL, qR);
}

void calcFaCir(int v)
{
faCir[v] = Circuit_2(
Mat0, Mat0,
Mat0, Mat0);

for (int i = 0; i < ND; i++)
for (int j = 0; j < ND; j++)
{
faCir[v].matA.a[i][i] += faR[v][i][j];
faCir[v].matB.a[i][j] -= faR[v][i][j];
faCir[v].matC.a[j][i] -= faR[v][i][j];
faCir[v].matD.a[j][j] += faR[v][i][j];
}
}

inline int getLca(int v, int u)
{
while (chainL[v] != chainL[u])
{
if (depth[trSeq[chainR[v] - 1]] < depth[trSeq[chainR[u] - 1]])
swap(v, u);
v = fa[trSeq[chainR[v] - 1]];
}
return depth[v] < depth[u] ? v : u;
}

Circuit_1 getDownD(int v)
{
if (!v)
return Mat0;
return closeL(tr_query(chainL[v], trPos[v] + 1));
}
Circuit_1 getDown(int v)
{
return mergeC(down[v], getDownD(preferredChild[v]));
}
Circuit_1 getDown(int v, int forbV)
{
Circuit_1 res = getDown(v);
if (forbV)
res = subC(res, getDownD(forbV));
return res;
}
Circuit_1 getDown(int v, int forbV1, int forbV2)
{
Circuit_1 res = getDown(v, forbV1);
if (forbV2)
res = subC(res, getDownD(forbV2));
return res;
}
Circuit_2 goUp(int v, int target, int &lastV)
{
Circuit_2 res(Mat0);
while (chainL[v] != chainL[target])
{
res = mergeA(res, getDown(v, lastV));
res = mergeB(res, faCir[v]);
res = mergeB(res, tr_query(trPos[v] + 1, chainR[v]));
lastV = trSeq[chainR[v] - 1], v = fa[lastV];
}
if (v != target)
{
res = mergeA(res, getDown(v, lastV));
res = mergeB(res, faCir[v]);
res = mergeB(res, tr_query(trPos[v] + 1, trPos[target]));
lastV = trSeq[trPos[target] - 1];
}
return res;
}
Circuit_1 getUp(int v)
{
if (!fa[v])
return Circuit_1(Mat0);
else
{
int lastV = v;
Circuit_2 res = goUp(fa[v], 1, lastV);
res = mergeA(res, getDown(1, lastV));
return closeL(mergeB(faCir[v], res).getRev());
}
}

double solveCircuit_2(const Circuit_2 &cir, int v, int u)
{
const int NVar = ND * 2;
double mat[NVar][NVar + 1];
for (int i = 0; i < ND; i++)
for (int j = 0; j < ND; j++)
mat[i][j] = cir.matA.a[i][j];
for (int i = 0; i < ND; i++)
for (int j = 0; j < ND; j++)
mat[i][ND + j] = cir.matB.a[i][j];
for (int i = 0; i < ND; i++)
for (int j = 0; j < ND; j++)
mat[ND + i][j] = cir.matC.a[i][j];
for (int i = 0; i < ND; i++)
for (int j = 0; j < ND; j++)
mat[ND + i][ND + j] = cir.matD.a[i][j];
for (int i = 0; i < NVar; i++)
mat[i][NVar] = 0;

mat[v][NVar] = 1, mat[u][NVar] = -1;

for (int j = 0; j < NVar; j++)
mat[0][j] = 0;
mat[0][0] = 1, mat[0][NVar] = 0;

for (int i = 0; i < NVar; i++)
for (int j = i + 1; j < NVar; j++)
if (sgn(mat[j][i]) != 0)
{
double z = mat[j][i] / mat[i][i];
for (int k = 0; k <= NVar; k++)
mat[j][k] -= mat[i][k] * z;
}
for (int i = NVar - 1; i >= 0; i--)
{
for (int j = i + 1; j < NVar; j++)
mat[i][NVar] -= mat[i][j] * mat[j][NVar];
mat[i][NVar] /= mat[i][i];
}
return mat[v][NVar] - mat[u][NVar];
}

void updateVer(int v)
{
calcFaCir(v);
seg_change(1, 0, n, trPos[v]);

for (v = trSeq[chainR[v] - 1]; fa[v]; v = trSeq[chainR[fa[v]] - 1])
{
down[fa[v]] = subC(down[fa[v]], downD[v]);
downD[v] = getDownD(v);
down[fa[v]] = mergeC(down[fa[v]], downD[v]);
seg_change(1, 0, n, trPos[fa[v]]);
}
}

double query(int v, int vl, int u, int ul)
{
int lca = getLca(v, u);
int lastV = 0, lastU = 0;
Circuit_2 cirL = goUp(v, lca, lastV);
Circuit_2 cirR = goUp(u, lca, lastU);
Circuit_1 cirM = mergeC(getDown(lca, lastV, lastU), getUp(lca));

Circuit_2 res = mergeB(mergeA(cirL, cirM), cirR.getRev());
return solveCircuit_2(res, vl, ND + ul);
}

int main()
{
freopen("circuit.in", "r", stdin);
freopen("circuit.out", "w", stdout);

cin >> n;

fa[1] = 0;
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
faR[1][i][j] = 100000;
calcFaCir(1);

for (int v = 2; v <= n; v++)
{
fa[v] = getint();
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
faR[v][i][j] = getint();
calcFaCir(v);
addEdge(fa[v], v);
}

init();

int nQ;
cin >> nQ;
while (nQ--)
{
int type;
scanf("%d", &type);
if (type == 1)
{
int v = getint(), vl = getint(), fl = getint(), nr = getint();
vl--, fl--;
faR[v][vl][fl] = nr;
updateVer(v);
}
else
{
int v = getint(), vl = getint(), u = getint(), ul = getint();
vl--, ul--;
printf("%.6f\n", query(v, vl, u, ul));
}
}

return 0;
}


  评论这张
 
阅读(2034)| 评论(1)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017