「CTSC2017」游戏

题目

点这里看题目。

---

按照如下方式生成一个长度为 \(n\) 的 \(01\) 串 \(s\)：

1. \(s_1\) 由一个参数 \(p_1\) 决定，表示 \(s_1\) 有 \(p_1\) 的概率为 \(\texttt 1\)，有 \(1-p_1\) 的概率为 \(\texttt 0\)。

2. 对于 \(1<k\le n\)，\(s_k\) 将由 \(s_{k-1}\) 和参数 \(p_k,q_k\) 共同决定。

   具体地，我们有下表：

   概率	\(s_{k}=\texttt 0\)	\(s_k=\texttt 1\)
   \(s_{k-1}=\texttt 0\)	\(1-q_k\)	\(q_k\)
   \(s_{k-1}=\texttt 1\)	\(1-p_k\)	\(p_k\)

现在会动态地加入删除一些条件，修改次数总计为 \(m\) 次。一条条件形如 \((x,c)\)，表示 \(s_x=c\)。

在每次条件变化之后，输出 \(s\) 中 \(\texttt 1\) 的个数的期望。

所有数据满足 \(1\le n,m\le 2\times 10^5,0<p_k,q_k<1\)，且输入的所有实数保留至多 \(4\) 位小数。

分析

这道题比较令人困惑的一点是：明明字符串是从前往后生成的，为什么靠后的字符确定后，靠前位置出现某种字符的概率会变化？

我们这样类理解“生成方式”：生成方式没有带来字符之间的必然联系，它仅仅给出了计算字符串的概率测度的方式。所以给出的“条件”实际上就影响了字符串的形态，相应地影响了样本空间，相应地影响了概率。

---

首先，本题的 \(\Omega=\{\texttt 0,\texttt 1\}^n\)，概率测度不方便描述，但可以直接通过“生成方式”计算。我们设 \(A_k\) 表示事件 \(s_k=\texttt 0\) 或者 \(s_k=\texttt 1\)，\(B_k\) 表示事件 \(s_k=\texttt 1\)，\(\overline B_k\) 表示事件 \(s_k=\texttt 0\)。假设所有条件对应的积为 \(R\)，我们相当于要求 \(\sum_{k=1}^nP(B_k|R)\)。

分析一下条件概率：由于“生成方式”为“线性一推一”，我们容易发现：

Observation.

如果 \(j_1<j_2<k\)，则 \(P(A_k|A_{j_1}A_{j_2})=P(A_{k}|A_{j_2})\)。

类似地，如果 \(j_1>j_2>k\)，则 \(P(A_{k}|A_{j_1}A_{j_2})=P(A_{k}|A_{j_2})\)。

此时，对于 \(P(B_{k}|R)\)，我们就可以找出限制中 \(x\) 最靠近 \(k\) 的两条，记它们的 \(x\) 为 \(l,r\)，满足 \(l<k<r\)（先假设这样的 \(l,r\) 均存在），那么答案为 \(P(B_k|R)=P(B_k|A_lA_r)\)。

那么，结合样例解释和 Bayes 公式，我们有：

\[\begin{aligned} P(B_k|A_lA_r) &=\frac{P(B_kA_lA_r)}{P(A_r|A_l)P(A_l)}\\ &=\frac{P(A_r|B_kA_l)P(B_k|A_l)P(A_l)}{P(A_r|A_l)P(A_l)}\\ &=\frac{P(A_r|B_k)P(B_k|A_l)}{P(A_r|A_l)} \end{aligned} \]

最后的结果很友好，所有需要的概率都变成“两两”之间的了。而 \(P(A_r|A_l)\) 这种概率可以用全概率公式慢慢推。例如：

\[P(B_r|B_l)= \begin{cases} p_rP(B_{r-1}|B_l)+q_rP(\overline B_{r-1}|B_l)&r>l+1\\ p_r&r=l+1 \end{cases} \]

---

快进到计算部分。当我们添加一条限制或者删除一条限制时，我们可以计算变化的区间的贡献。因此，我们需要计算的就是——给定 \(l,r\)，求：

\[\frac{1}{P(A_r|A_l)}\sum_{k=l+1}^{r-1}P(A_r|B_k)P(B_k|A_l) \]

结合两两条件概率的计算方式，我们容易想到用矩阵表示递推。考虑已有向量 \(\begin{bmatrix}U_{s,t}&V_{s,t}\end{bmatrix}\)，\(U_{s,t}\) 为 \(2\times 2\) 的概率矩阵，\(V_{s,t}\) 为 \(2\times 2\) 的“概率和”矩阵。举个例子：

\[\begin{aligned} (U_{s,t})_{01}&=P(B_t|\overline B_s)\\ (V_{s,t})_{01}&=\sum_{k=s}^{t-1}P(B_t|B_k)P(B_k|\overline B_s) \end{aligned} \]

Note.

这里的构造比较 tricky。理论上来说把和式定义为左闭右闭会比较常见，但由于我们最终需要乘上 \(r\) 对应的转移矩阵，所以把求和上界定义为 \(t-1\) 会便于最终的计算。

此外，递推向量的结构本身限制了它在一次乘法中无法涉及过多元素，因此如果要算到 \(t\)，那么涉及元素太多，不容易构造。

良好的定义带来了简单的构造过程。最后可以构造出来转移矩阵形如：

\[\begin{bmatrix} U_{s,t-1}&V_{s,t-1} \end{bmatrix} \begin{bmatrix} 1-q_{t}&q_{t}&0&0\\ 1-p_{t}&p_{t}&1-p_{t}&p_{t}\\ 0&0&1-q_{t}&q_{t}\\ 0&0&1-p_{t}&p_{t} \end{bmatrix}= \begin{bmatrix} U_{s,t}&V_{s,t} \end{bmatrix} \]

需要注意，初始向量需要单独构造。

---

上述讨论我们保证 \(l,r\) 存在。为了延续这一点，我们可以引入 \(s_0\) 和 \(s_{n+1}\)，并且将 \(B_0\) 和 \(B_{n+1}\) 始终作为固定条件，这样 \(l,r\) 就始终存在了。

最后求矩阵区间乘积可以使用 zkw 线段树，不知道可不可以求逆。复杂度为 \(O((n+m)\log n)\)，（我的）常数比较大。

代码

#include <set>
#include <cstdio>

#define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ )
#define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- )

const int MAXN = 2e5 + 5;

template<typename _T>
inline void Read( _T &x ) {
    x = 0; char s = getchar(); bool f = false;
    while( s < '0' || '9' < s ) { f = s == '-', s = getchar(); }
    while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
    if( f ) x = -x;
}

template<typename _T>
inline void Write( _T x ) {
    if( x < 0 ) putchar( '-' ), x = -x;
    if( 9 < x ) Write( x / 10 );
    putchar( x % 10 + '0' );
}

struct Matrix {
    double mat[4][4];

    Matrix(): mat{} {}
};

std :: set<int> s;

Matrix tre[MAXN << 2];
int bas;

double P[MAXN], Q[MAXN];
int fix[MAXN];

int N, M;

inline Matrix operator * ( const Matrix &a, const Matrix &b ) {
    Matrix ret;
    rep( i, 0, 3 ) rep( k, 0, 3 ) rep( j, 0, 3 )    
        ret.mat[i][j] += a.mat[i][k] * b.mat[k][j];
    return ret;
}

void Build( const int &n ) {
    for( bas = 1 ; bas <= n + 2 ; bas <<= 1 );
    rep( i, 1, n + 1 ) {
        int x = i + bas;
        tre[x].mat[0][0] = 1 - Q[i], tre[x].mat[0][1] = Q[i];
        tre[x].mat[1][0] = 1 - P[i], tre[x].mat[1][1] = P[i];
        rep( b, 0, 1 ) rep( c, 0, 1 )
            tre[x].mat[b + 2][c + 2] = tre[x].mat[b][c];
        tre[x].mat[1][2] = 1 - P[i], tre[x].mat[1][3] = P[i];
    }
    per( i, bas - 1, 1 )
        tre[i] = tre[i << 1] * tre[i << 1 | 1];
}

inline Matrix Query( int l, int r ) {
    Matrix lef, rig;
    l += bas - 1, r += bas + 1;
    rep( i, 0, 3 ) lef.mat[i][i] = rig.mat[i][i] = 1;
    while( l ^ r ^ 1 ) {
        if( ! ( l & 1 ) ) lef = lef * tre[l ^ 1];
        if( r & 1 ) rig = tre[r ^ 1] * rig;
        l >>= 1, r >>= 1;
    }
    return lef * rig;
}

inline double Ask( const int &l, const int &r ) {
    if( l > r ) return 0;
    double ans = 0;
    Matrix tmp = Query( l + 1, r + 1 );
    if( fix[l - 1] ) {
        ans = ( 1 - P[l] ) * tmp.mat[0][fix[r + 1] + 2] + P[l] * tmp.mat[1][fix[r + 1] + 2];
        ans /= ( 1 - P[l] ) * tmp.mat[0][fix[r + 1]] + P[l] * tmp.mat[1][fix[r + 1]];
    } else {
        ans = ( 1 - Q[l] ) * tmp.mat[0][fix[r + 1] + 2] + Q[l] * tmp.mat[1][fix[r + 1] + 2];
        ans /= ( 1 - Q[l] ) * tmp.mat[0][fix[r + 1]] + Q[l] * tmp.mat[1][fix[r + 1]];
    }
    return ans;
}

int main() {
    scanf( "%d %d %*c", &N, &M );
    scanf( "%lf", P + 1 );
    rep( i, 2, N ) 
        scanf( "%lf %lf", P + i, Q + i );
    P[N + 1] = Q[N + 1] = 1;
    fix[0] = fix[N + 1] = 1;

    Build( N );
    double ans = Ask( 1, N );
    s.insert( 0 ), s.insert( N + 1 );
    while( M -- ) {
        char buf[10];
        scanf( "%s", buf );
        if( buf[0] == 'a' ) {
            int x, y; 
            Read( x ), Read( y );
            std :: set<int> :: iterator 
                it = s.lower_bound( x );
            int succ = *it, pred = * -- it;
            ans -= Ask( pred + 1, succ - 1 );
            fix[x] = y, ans += y;
            ans += Ask( pred + 1, x - 1 );
            ans += Ask( x + 1, succ - 1 );
            s.insert( x );
        }
        if( buf[0] == 'd' ) {
            int x; 
            Read( x ), s.erase( x );
            std :: set<int> :: iterator
                it = s.lower_bound( x );
            int succ = *it, pred = * -- it;
            ans -= Ask( pred + 1, x - 1 );
            ans -= Ask( x + 1, succ - 1 );
            ans -= fix[x], fix[x] = -1;
            ans += Ask( pred + 1, succ - 1 );
        }
        printf( "%.10lf\n", ans );
    }
    return 0;
}