@gym - 100591D@ Fox Rocks

Posted 2021-03-17 tiw-air-oao

目录

• @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 的形式即可。

@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 它高斯消元竟然把一行给消没了？！难道。。。”

然后就发现了问题所在。