「CTSC2017」游戏
题目
点这里看题目。
按照如下方式生成一个长度为 \(n\) 的 \(01\) 串 \(s\)：
-
\(s_1\) 由一个参数 \(p_1\) 决定，表示 \(s_1\) 有 \(p_1\) 的概率为 \(\texttt 1\)，有 \(1-p_1\) 的概率为 \(\texttt 0\)。
-
对于 \(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 公式，我们有：
最后的结果很友好，所有需要的概率都变成“两两”之间的了。而 \(P(A_r|A_l)\) 这种概率可以用全概率公式慢慢推。例如：
快进到计算部分。当我们添加一条限制或者删除一条限制时，我们可以计算变化的区间的贡献。因此，我们需要计算的就是——给定 \(l,r\)，求：
结合两两条件概率的计算方式，我们容易想到用矩阵表示递推。考虑已有向量 \(\begin{bmatrix}U_{s,t}&V_{s,t}\end{bmatrix}\)，\(U_{s,t}\) 为 \(2\times 2\) 的概率矩阵，\(V_{s,t}\) 为 \(2\times 2\) 的“概率和”矩阵。举个例子：
Note.
这里的构造比较 tricky。理论上来说把和式定义为左闭右闭会比较常见，但由于我们最终需要乘上 \(r\) 对应的转移矩阵，所以把求和上界定义为 \(t-1\) 会便于最终的计算。
此外，递推向量的结构本身限制了它在一次乘法中无法涉及过多元素，因此如果要算到 \(t\)，那么涉及元素太多，不容易构造。
良好的定义带来了简单的构造过程。最后可以构造出来转移矩阵形如：
需要注意，初始向量需要单独构造。
上述讨论我们保证 \(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;
}
