@gym - 100591D@ Fox Rocks
2024-08-29 21:00:05
@description@
给定 N 个点以及 P 条单向道路 Ai -> Bi，每条道路有 Ri 块石头（保证\(0 \le \lfloor\frac{Bi}{4}\rfloor-\lfloor\frac{Ai}{4}\rfloor\le 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 它高斯消元竟然把一行给消没了？！难道。。。”
然后就发现了问题所在。
最新文章
- freeCodeCamp:Chunky Monkey
- （转）JS模块化编程之AMD规范
- centos+nginx从零开始配置负载均衡
- Delphi的属性Property
- JS Date.Format
- poj 2312 Battle City
- java.lang.NoClassDefFoundError: JspException
- java集合总结
- Scut：Redis 资源管理器
- C++编写ATM(2)
- windows下配置caffe(环境：win7+vs2013+opencv3.0)
- JspContext对象与PageContext对象
- python的numpy库的学习
- zabbix_agentd客户端安装与配置(windows操作系统)
- 不输入密码执行SUDO命令
- 程序员调 Bug 的样子，非常真实
- 【模板/经典题型】FWT
- 民生银行十五年的数据体系建设，深入解读阿拉丁大数据生态圈、人人BI 是如何养成的？【转】
- Update Node.js Package.json
- Oracle中已知字段名查询所在的表名