@gym

目录

• @description@
• @solution@
• @accepted code@
• @details@

---

@description@

给定 N 个点以及 P 条单向道路 Ai -> Bi，每条道路有 Ri 块石头（保证(0 le lfloorfrac{Bi}{4} floor-lfloorfrac{Ai}{4} floorle 1)）。

接下来 D 天每天有一个询问，分为三类：
（1）给定 X, Y, Z，修建一条 X -> Y 且有 Z 块石头的道路。
（2）给定 X, Y，拆除 X -> Y 的原有道路。
（3）查询从 X 开始随机游走（从它的出边中等概率选一块石头，并往那条边走）可以到达 Y 的概率。

链接。

@solution@

一杯茶，一包烟，一道破题调一天。

如果我们把连续 4 个点看成一个块，则连边情况应是块内连边 + 块向下一块连边。
如果暴力着做概率 dp，可以从起点 X 出发，向后递推到 Y 所在块，每一块的块内进行高斯消元。

显然可以用线段树加速这一过程：每个结点 [l, r] 维护从第 l 个块的第 i 个点随机游走，第一次到达第 r 个块是第 r 个块的第 j 个点的概率 p[i][j]。
线段树上 merge 的过程就是个高斯消元。
最后维护一下对 Y 所在块做一个高斯消元，从 Y 所在块的第 i 个点到第 j 个点的概率 q[i][j]。

然而。。。有一个小问题。
概率 dp 的一般形式是：给定若干终止状态，保证每个点都可以到达一种终止状态，且每条转移边有一定概率。
定义 dp[i] 表示到达想要的终止状态的概率，则 dp[i] 要么游走到它可以到达的点，要么直接到达终止状态。

而这道题终止状态有 “没有出度”，“到达 Y”，“到达 Y 的下一个块”，还有一个 “进入一个无法到达 Y 的环”。
前三种可以通过 dp 的转移式子表达出来，然而最后一个并不能直接表示出来。。。
方法要么是通过简单 dfs 判断，要么还有一种特殊的判断：如果高斯消元中主对角线 = 0，则一定出现了最后一种情况。此时直接把这个方程修改成 xi = 0 的形式即可。

注意还有一个情况：询问中可能 Y 所在块 < X 所在块。此时直接输出 0。

@accepted code@

#include <cmath>
#include <cstdio>
#include <algorithm>
using namespace std;

const int MAXN = 50000;
const double EPS = 1E-9;

double dcmp(double x) {return fabs(x) <= EPS ? 0 : (x > 0 ? 1 : -1);}

struct node{double a[4][4];};

int K;

double d[MAXN + 5];
node f[MAXN + 5], g[MAXN + 5], h[MAXN + 5];

double M[100][100];
void gauss(int n, int m) {
	int r = 0, c = 0;
	while( r < n && c < m ) {
		if( dcmp(M[r][c]) ) {
			double k = M[r][c];
			for(int j=c;j<m;j++)
				M[r][j] /= k;
			for(int i=0;i<n;i++) {
				if( i == r ) continue;
				k = M[i][c];
				for(int j=c;j<m;j++)
					M[i][j] -= k*M[r][j];
			}
			r++;
		}
		else {
			for(int j=c;j<m;j++)
				M[r][j] = 0;
			M[r][c] = 1;
			r++;
		}
		c++;
	}
}

void get(int x, node &p) {
	for(int t=0;t<4;t++) {
		for(int r=0;r<4;r++)
			for(int c=0;c<4;c++)
				M[r][c] = 0;
		for(int r=0;r<4;r++) {
			if( r == t ) {
				for(int c=0;c<4;c++)
					M[r][c] = (r == c);
				M[r][4] = 1;
			}
			else {
				int k = 4*x + r;
				for(int c=0;c<=4;c++)
					M[r][c] = (r == c);
				if( d[k] == 0 ) continue;
				for(int c=0;c<4;c++)
					M[r][c] -= f[x].a[r][c] / d[k];
			}
		}
		gauss(4, 4+1);
		for(int s=0;s<4;s++) 
			p.a[s][t] = M[s][4];
	}
}

struct segtree{
	#define lch (x << 1)
	#define rch (x << 1 | 1)
	
	int le[4*MAXN + 5], ri[4*MAXN + 5]; node h[4*MAXN + 5];
	void merge(const node &f1, const node &f2, node &f3, int m) {
		for(int i=0;i<4;i++)
			for(int j=0;j<4;j++)
				f3.a[i][j] = 0;
		for(int t=0;t<4;t++) {
			for(int r=0;r<4;r++) {
				int k = 4*m + r;
				for(int c=0;c<4;c++) {
					if( d[k] ) M[r][c] = (r == c) - f[m].a[r][c]/d[k];
					else M[r][c] = (r == c);
				}
				if( r == t ) M[r][4] = 1;
				else M[r][4] = 0;
			}
			gauss(4, 5);
			int k = 4*m + t;
			if( d[k] ) {
				double a[4] = {}, b[4] = {};
				for(int p=0;p<4;p++)
					for(int s=0;s<4;s++)
						a[p] += M[s][4]*f1.a[p][s];
				for(int q=0;q<4;q++)
					for(int x=0;x<4;x++)
						b[q] += g[m].a[t][x]/d[k]*f2.a[x][q];
				for(int p=0;p<4;p++)
					for(int q=0;q<4;q++)
						f3.a[p][q] += a[p]*b[q];
			//f3.a[p][q] += M[s][4]*div(g[m].a[t][x], d[k])*f1.a[p][s]*f2.a[x][q];
			}
		}
	}
	void pushup(int x) {
		int m = (le[x] + ri[x]) >> 1;
		merge(h[lch], h[rch], h[x], m);
	}
	void build(int x, int l, int r) {
		le[x] = l, ri[x] = r;
		if( l == r ) {
			for(int i=0;i<4;i++)
				for(int j=0;j<4;j++)
					h[x].a[i][j] = (i == j);
			return ;
		}
		int m = (le[x] + ri[x]) >> 1;
		build(lch, l, m), build(rch, m + 1, r);
		pushup(x);
	}
	void update(int x, int p) {
		if( le[x] == ri[x] ) return ;
		int m = (le[x] + ri[x]) >> 1;
		if( p <= m ) update(lch, p);
		else update(rch, p);
		pushup(x);
	}
	node query(int x, int ql, int qr) {
		if( ql <= le[x] && ri[x] <= qr )
			return h[x];
		int m = (le[x] + ri[x]) >> 1;
		if( qr <= m ) return query(lch, ql, qr);
		else if( ql > m ) return query(rch, ql, qr);
		else {
			node a = query(lch, ql, qr), b = query(rch, ql, qr), c;
			merge(a, b, c, m); return c;
		}
	}
}T;


void update(int A, int B, double R) {
	d[A] += R;
	if( A / 4 == B / 4 ) f[A/4].a[A%4][B%4] += R;
	else g[A/4].a[A%4][B%4] += R;
	get(A/4, h[A/4]), T.update(1, A/4);
}
void remove(int A, int B) {
	if( A / 4 == B / 4 ) update(A, B, -f[A/4].a[A%4][B%4]);
	else update(A, B, -g[A/4].a[A%4][B%4]);
}

int N, P, D;

void solve() {
	scanf("%d%d%d", &N, &P, &D), K = N, N = (N + 3) / 4;
	for(int i=0;i<4*N;i++) d[i] = 0;
	for(int i=0;i<N;i++)
		for(int j=0;j<4;j++)
			for(int k=0;k<4;k++)
				f[i].a[j][k] = g[i].a[j][k] = 0;
	for(int i=1;i<=P;i++) {
		int A, B, R; scanf("%d%d%d", &A, &B, &R);
		d[A] += R;
		int del = B - (A/4*4);
		if( del < 4 ) f[A/4].a[A%4][del] += R;
		else g[A/4].a[A%4][del-4] += R;
	}
	for(int i=0;i<N;i++)
		get(i, h[i]);
	T.build(1, 0, N - 1);
	for(int i=1;i<=D;i++) {
		int op; scanf("%d", &op);
		if( op == 1 ) {
			int X, Y, Z; scanf("%d%d%d", &X, &Y, &Z);
			update(X, Y, Z);
		}
		else if( op == 2 ) {
			int X, Y; scanf("%d%d", &X, &Y);
			remove(X, Y);
		}
		else {
			int X, Y; scanf("%d%d", &X, &Y);
			if( X / 4 <= Y / 4 ) {
				node p = T.query(1, X / 4, Y / 4);
				double ans = 0;
				for(int i=0;i<4;i++)
					ans += p.a[X%4][i] * h[Y/4].a[i][Y%4];
				printf(" %.6f", ans);
			}
			else printf(" %.6f", 0.0);
		}
	}
}

int main() {
//	freopen("data.in", "r", stdin);
//	freopen("data.out", "w", stdout);
	int t; scanf("%d", &t);
	for(int i=1;i<=t;i++)
		printf("Case #%d:", i), solve(), puts("");
}
/*
5
8 0 10
1 0 1 1
1 1 2 1
1 2 3 1
3 0 3
1 0 4 1
3 0 3
1 0 3 1
3 0 3
2 1 2
3 0 3

4 4 10
0 1 1
1 0 1
1 2 1
0 3 1
3 0 2
2 0 1
1 0 1 2
3 0 2
2 0 1
1 0 1 3
3 0 2
2 0 1
1 0 1 4
3 0 2

8 7 5
0 1 1
1 2 1
2 3 1
0 4 1
1 5 1
2 6 1
3 7 1
3 0 5
3 0 7
1 3 0 1
3 0 5
3 0 7

8 4 10
4 5 1
5 6 1
6 7 1
7 4 1
1 0 4 1
1 0 1 4
3 0 7
1 1 6 1
3 0 7
1 1 2 1
1 2 3 1
3 0 7
1 2 0 1
3 0 7

12 5 7
0 4 1
4 8 1
0 1 1
1 2 1
2 3 1
3 0 8
1 1 4 1
3 0 8
1 2 4 1
3 0 8
1 3 4 1
3 0 8
*/

@details@

做题心路历程：

“woc 这是什么恶心题，算了思路挺简单的，先写再说”
“woc 怎么过不了样例，woc 我题给读错了，它不能往前面那个块跑。没事儿反正也是往简单的方向改”
“woc 怎么卡在第二个数据上了，woc 它竟然 WA 了！怎么办该不会是精度问题吧”
“（去偷了一份正确代码，生成了一点大数据）！怎么会有 nan 这种东西，看来得调试了”
“woc 它高斯消元竟然把一行给消没了？！难道。。。”

然后就发现了问题所在。