题解 AT4690 e

题目传输门

更好的阅读体验

题意简述

• 一张有 \(M\) 个位置的长椅，初始为空。
• 有 \(M\) 个人依次落座，每个人会在长椅中满足 自己与相邻共三个位置均未被占据 的位置中随机选择一个坐下。若不存在满足条件的位置则不落座。
• 随机选择 \(N\) 个连续的位置，求当 \(M\) 趋近于无穷大时，这 \(N\) 个位置的落座位置与给出序列 \(s\) 一致的概率的极限（\(s\) 中 X 表示有人，- 表示无人）。
• 可以证明答案可以表示成 \(p+\frac{q}{e}+\frac{r}{e^2}\)，其中 \(p,q,r\) 均为有理数，对 \(10^9+7\) 取模。
• \(1 \leq |S|=N \leq 1000\)。

题解

Part 1

先考虑 \(s=X\)（样例 \(1\)），即求落座比例的期望。

题目中的过程较为繁琐，每个人会随机选择自己的位置。我们不妨令每个人只对应一个位置，这样最后的情况只与落座顺序有关。

具体而言，随机排列 \(P_1,\cdots,P_M\)，然后让 \(M\) 个人依次落座 \(P_1,\cdots,P_M\)（如果符合题意的话）。

但是由于 \(M\) 趋近于无穷大，枚举排列不好处理。考虑给每个位置赋一个 \([0,1]\) 中的值 \(X_i\) 来代替排列 \(P\)。让 \(M\) 个人按 \(X_i\) 从小到大的顺序依次落座，仍是等价的。

考虑一个位置 \(i\) 最终是否被占据。

设其向左与向右的极长连续下降长度分别为 \(L\) 与 \(R\)（即 \(X_{i-L-1}<X_{i-L}<\cdots<X_i<\cdots<X_{i+R}>x_{i+R+1}\)），那么 \(i-L\) 与 \(i+R\) 一定会被占据，同理 \(i-L+2\) 与 \(i+R-2\) 也会被占据……

故 \(i\) 被占据当且仅当 \(L\) 与 \(R\) 都是偶数。

设 \(X_i=x\)，那么 \(L\) 为偶数的概率为：

\[\begin{aligned} P(2 \mid L)=P(2 \mid R) & =1-P(X_{i-1}<X_i)+P(X_{i-2}<X_{i-1}<X_i)-P(X_{i-3}<X_{i-2}<X_{i-1}<X_i)+\cdots\\ & =1-x+\dfrac{x^2}{2}-\dfrac{x^3}{6}+\cdots\\ & =e^{-x} \end{aligned} \]

故答案为 \(\int_0^1 e^{-2x}=\dfrac{1}{2}-\dfrac{1}{2e^2}\)。

Part 2

考虑原问题，可以考虑 DP。

假设选出的 \(N\) 个位置就是 \(1,2,\cdots,N\)，我们来从左向右填 \(X_i\)（\(-\infty\) 到 \(0\) 均已填好了），同时计算与 \(s\) 匹配的概率。

这样会有一个问题：在 DP 过程中，我们只知道 \(L\) 的奇偶性，但 \(R\) 的奇偶性还未确定。

但是 \(R\) 的奇偶性只有两种，不妨倒推。

设 \(dp_{i,id_1,id_2,id_3}\) 为关于 \(x\) 的概率密度函数，四个变量分别表示到第 \(i\) 个位置，\(L\) 的奇偶性，\(R\) 奇偶性为 \(0\) 时是否合法，\(R\) 奇偶性为 \(1\) 时是否合法。

转移时仍然设 \(X_i=x\)，考虑 \(X_{i-1}\) 与 \(x\) 的大小关系。

• \(X_{i-1}<x\)。则 \(dp_{i-1,id_1 \oplus 1,1,0/1} \rightarrow dp_{i,id_1,X,X}\)，具体而言：
  • 如果该位为 -：\(dp_{i-1,1,1,0/1} \rightarrow dp_{i,0,0,1}\)；\(dp_{i-1,0,1,0/1} \rightarrow dp_{i,1,1,1}\)。
  • 如果该位为 X：\(dp_{i-1,1,1,0/1} \rightarrow dp_{i,0,1,0}\)。
• \(X_{i-1}>x\)。则 \(dp_{i-1,0/1,id_2,id_3} \rightarrow dp_{i,0,X,X}\)，具体而言：
  • 如果该位为 -：\(dp_{i-1,0/1,1,0/1} \rightarrow dp_{i,0,0,1}\)。
  • 如果该位为 X：\(dp_{i-1,0/1,0/1,1} \rightarrow dp_{i,0,1,0}\)。

初始状态为 \(dp_{0,0,1,1}=e^{-x}\) 及 \(dp_{0,1,1,1}=1-e^{-x}\)。

实现时需要计算 \(\int_0^x dp_{i-1,X,X,X} dy\) 或是 \(\int_x^1 dp_{i-1,X,X,X} dy\)，注意到 \(c \int e^{-x}=-ce^{-x}+C\)，故每个 DP 值一定是 \(P(x)+Q(x)e^{-1}+ce^{-x}\)，其中 \(P(x)\) 与 \(Q(x)\) 均为关于 \(x\) 的不超过 \(i\) 次多项式。

计算答案时，要将 DP 状态乘上 \(e^{-x}\) 或是 \(1\) 或是 \(1-e^{-x}\) 再求在 \([0,1]\) 上的定积分。

只需要计算两个积分：\(\int_0^1 x^k e^{-x} dx\) 及 \(\int_0^1 e^{-2x}dx=\dfrac{1}{2}-\dfrac{1}{2e^2}\)。

分部积分法：

\[\begin{aligned} \int x^k e^{-x} dx & =-\int x^k d e^{-x}\\ & =-x^ke^{-x}+\int e^{-x} d x^k\\ & =-x^ke^{-x}+k \int x^{k-1} e^{-x} d x\\ & =-x^ke^{-x}-k x^{k-1}e^{-x}+k(k-1) \int x^{k-2} e^{-x} d x\\ & =\cdots\\ & =-e^{-x} \sum\limits_{i=0}^k k^{\underline i} x^{k-i} \end{aligned} \]

\[\int_0^1 x^k e^{-x} dx=k!- e^{-1} \sum\limits_{i=0}^k k^{\underline i} \]

总时间复杂度 \(O(N^2)\)。

代码

#include<bits/stdc++.h>
using namespace std;
const long long N=1100,mod=1e9+7,inv2=(mod+1)/2;
long long n;
char S[N];
long long invv[N],ans[3];

struct abc{
    long long P[N],Q[N],c;
    abc(){
        memset(P,0,sizeof(P));
        memset(Q,0,sizeof(Q));
        c=0;
    }
}dp[N][2][2][2];

abc operator +(abc X,abc Y){
	abc T;
	for(long long i=0;i<=n;i++) T.P[i]=(X.P[i]+Y.P[i])%mod,T.Q[i]=(X.Q[i]+Y.Q[i])%mod;
	T.c=(X.c+Y.c)%mod;
	return T;
}
abc operator -(abc X,abc Y){
	abc T;
	for(long long i=0;i<=n;i++) T.P[i]=(X.P[i]+mod-Y.P[i])%mod,T.Q[i]=(X.Q[i]+mod-Y.Q[i])%mod;
	T.c=(X.c+mod-Y.c)%mod;
	return T;
}

abc INTx(abc X){
    abc T;
    for(long long i=0;i<=n;i++) T.P[i+1]=X.P[i]*invv[i+1]%mod;
    for(long long i=0;i<=n;i++) T.Q[i+1]=X.Q[i]*invv[i+1]%mod;
    T.c=(mod-X.c)%mod;
    return T;
}//求 x 处积分
abc INT0(abc X){
    abc T;
    T.P[0]=(mod-X.c)%mod;
    return T;
}//求 0 处积分
abc INT1(abc X){
    abc T;
    for(long long i=0;i<=n;i++) T.P[0]=(T.P[0]+X.P[i]*invv[i+1]%mod)%mod;
    for(long long i=0;i<=n;i++) T.Q[0]=(T.Q[0]+X.Q[i]*invv[i+1]%mod)%mod;
    T.Q[0]=(T.Q[0]+mod-X.c)%mod;
    return T;
}//求 1 处积分

abc cal1(abc X){
    return INTx(X)-INT0(X);
}
abc cal2(abc X){
    return INT1(X)-INTx(X);
}

void solve1(abc X,long long w){
	abc T=INT1(X)-INT0(X);
	ans[0]=(ans[0]+T.P[0]*w%mod)%mod;
	ans[1]=(ans[1]+T.Q[0]*w%mod)%mod;
}

void solvex(abc X,long long w){
	ans[0]=(ans[0]+X.c*inv2%mod*w%mod)%mod;
	ans[2]=(ans[2]+mod-X.c*inv2%mod*w%mod)%mod;
	
	ans[0]=(ans[0]+X.P[0]*w%mod)%mod;
	ans[1]=(ans[1]+mod-X.P[0]*w%mod)%mod;
	
	ans[1]=(ans[1]+X.Q[0]*w%mod)%mod;
	ans[2]=(ans[2]+mod-X.Q[0]*w%mod)%mod;
	
	for(long long k=1;k<=n;k++){
		long long tot=1;
		long long nw=1;
		for(long long i=k;i>=1;i--){
			nw=nw*i%mod;
			tot=(tot+nw)%mod;
		}
		ans[1]=(ans[1]+mod-X.P[k]*tot%mod*w%mod)%mod;
		ans[0]=(ans[0]+X.P[k]*nw%mod*w%mod)%mod;
		
		ans[2]=(ans[2]+mod-X.Q[k]*tot%mod*w%mod)%mod;
		ans[1]=(ans[1]+X.Q[k]*nw%mod*w%mod)%mod;
	}
}//求答案

int main(){
    invv[1]=1;
    for(long long i=2;i<N;i++) invv[i]=(mod-mod/i)*invv[mod%i]%mod;

    cin>>n;
    scanf("\n%s",S+1);
	
	dp[0][0][1][1].c=1;
	dp[0][1][1][1].P[0]=1,dp[0][1][1][1].c=mod-1;
	
	for(long long i=1;i<=n;i++){
		if(S[i]=='-'){
			dp[i][0][0][1]=dp[i][0][0][1]+cal1(dp[i-1][1][1][0]+dp[i-1][1][1][1]);
			dp[i][1][1][1]=dp[i][1][1][1]+cal1(dp[i-1][0][1][0]+dp[i-1][0][1][1]);
			
			dp[i][0][0][1]=dp[i][0][0][1]+cal2(dp[i-1][0][1][0]+dp[i-1][1][1][0] +dp[i-1][0][1][1]+dp[i-1][1][1][1]);
		}
		else{
			dp[i][0][1][0]=dp[i][0][1][0]+cal1(dp[i-1][1][1][0]+dp[i-1][1][1][1]);
			
			dp[i][0][1][0]=dp[i][0][1][0]+cal2(dp[i-1][0][0][1]+dp[i-1][1][0][1] +dp[i-1][0][1][1]+dp[i-1][1][1][1]);
		}
	}
	
	for(long long i=0;i<=1;i++)
		for(long long j=0;j<=1;j++)
			for(long long k=0;k<=1;k++){
				if(j && k) solve1(dp[n][i][j][k],1);
				if(j && !k) solvex(dp[n][i][j][k],1);
				if(!j && k) solve1(dp[n][i][j][k],1),solvex(dp[n][i][j][k],mod-1);
			}
	cout<<ans[0]<<" "<<ans[1]<<" "<<ans[2];
}