100道动态规划——42 CodeForces 908D New Year and Original Order 概率DP

Good Bye 2017！

概率DP真心不会做......

题意是给出一个k，一个pa和一个pb

一开始当前串是一个空串，每一次操作都有pa/(pa+pb)的概率向当前串附加一个字符'a'，pb/(oa+pb)的概率附加一个字符b，当考虑当前串的子序列，当出现大于等于k个子序列'ab'的时候，就停止，询问这时候子序列中'ab'的期望出现次数。

定义状态dp[i][j]表示当前串中出现了 i个子序列'a'，j个子序列'ab'时，再继续做下去，最后停止时的子序列'ab'期望出现次数

首先要考虑的问题就是，这里i和j是不存在关系的，假若这一点想清楚了，才能继续往下看

按照期望的定义，我们向dp[i][j]后附加一个a时，转移到的状态应该是dp[i+1][j]，附加一个b时，转移到的状态应该是dp[i][i+j]

因此我们的状态转移方程为：

dp[i][j]=pa/(pa+pb)*dp[i+1][j]+pb/(pa+pb)*dp[i][i+j]

我们的初始状态是dp[i][j] = j 这里i随意，j>=k

我们的末状态是dp[0][0]，即对应着一开始的空串

但是！

上面的状态定义其实有两个小问题

1. 我们刚刚说了，初始状态dp[i][j]=j 这里i随意，也就事实上也就意味着字符a可以无限的附加下去。因此我们做一个转化，变成dp[i][j]= j + i +pa/pb  这里i+j>=k

我们来看这个式子，按照定义dp[i][j]停止后ab出现的期望，因为i+j>=k，因此只要出现一个b即会停止，因此我们这里采用一个几何分布的期望，pa/pb代表着出现第一次出现b的时候，a的期望出现次数，然后再加上本来就有i个a，因此最后的ab的数量就是 i+j+pa/pb

2. 由于b可以在第一个a出现之前无穷的出现，因此我们把最终求解的状态设置为 dp[1][0]

3. 一开始就把 pa处理成 pa*(pa+pb)^-1 , 把pb处理成 pb*(pa+pb)^-1 以后直接这样计算即可

实际上代码很短

#include <iostream>
#include <cstring>

using namespace std;
typedef long long ll;
const ll mod=1000000007,maxm=1010;

ll inv(ll a,ll b=mod),dp[maxm][maxm],pa,pb,k;
inline getdp(int i,int j){return j>=k?j:dp[i][j];}
void egcd(ll a,ll b,ll& d,ll& x,ll& y);

int main(){
    ios_base::sync_with_stdio(0);
    cin>>k>>pa>>pb;
    const ll pb_inv=inv(pb),_inv=inv(pa+pb),pa_new=(pa*_inv)%mod,pb_new=(pb*_inv)%mod;
    for(int j=0;j<k;++j)
        dp[k][j]=(j+k+pa*pb_inv)%mod;
    for(int i=k-1;i;--i)
    for(int j=k-1;j>=0;--j)
        dp[i][j]=(pa_new*getdp(i+1,j)+pb_new*getdp(i,i+j))%mod;
    cout<<dp[1][0]<<endl;
    return 0;
}

ll inv(ll a, ll b){
	ll x,y,d=1;
    egcd(a,b,d,x,y);
    return (x+mod)%mod;
}

void egcd(ll a, ll b, ll& d, ll& x, ll& y){
    b?(egcd(b,a%b,d,y,x),y-=x*(a/b)):(d=a,x=1,y=0);
}