Good Bye 2017 D.New Year and Arbitrary Arrangement - 数学

题目大意：你有一个随机数生成器，有p的概率生成字符a，(1-p)的概率生成b。给定p=u/v，k，当生成的字符串中有至少k个形如’ab’的子序列（注意不是子串，例如abaab有4个形如’ab’的子序列）的时候就停止。不难注意到，停止的时候显然可能有多于k个形如’ab’的子序列。现在求期望情况下停止的时候有多少个形如’ab’的子序列。
k<=1000。答案对大质数取模输出。

题解：显然考虑dp……
如果令dp[i][j][t]表示前i个字符有j个a现在’ab’子序列的数量（下文记作序列数量）达到了t的概率。显然是能做的，但是i和t都是无穷大……
不难发现这样一个性质：如果加入一个a，序列数量不变；否则序列数量增加量等价于之前a的数量，显然。
因此你只要纪录a的数量即可，即令dp[i][j]表示当前有i个a，序列数量是j的概率。其中k>j。
转移是显然的，转移过程中如果转移到了j>=k的部分就统计入答案。
这样就完了么？显然没有，因为还有一种情况，就是积攒了很多个a，也就是i>=k，但是现在k>j，这样你再放1个b，序列数量就>=k了，但是在dp过程中没有统计入答案（注意，k>i，并且i+j>=k的时候会在dp过程中直接计入答案）。

这样，不妨考虑dp[i][j]在i>=k的值，不难发现这时候有：（记p=pa/(pa+pb)，q=1-p）

$$dp[i][j]=p*dp[i-1][j]=...=p^{i-k+1}*dp[k][j]$$
这样，这部分对答案的贡献就是：
$$\sum_{j=0}^{k-1}\sum_{i=k}^{+oo}dp[k-1][j]p^{i-k+1}q(i+j)$$
$$\sum_{j=0}^{k-1}q*dp[k-1][j]\sum_{i=k}^{+oo}p^{i-k+1}(i+j)$$
对后面那个数列求和，可以得到上式为：
$$\sum_{j=0}^{k-1}q*dp[i-1][j]*[\frac{p(k+j)}{1-p}+\frac{p^2}{(1-p)^2}]$$
注意到q=1-p，因此有：
$$Answer+=\sum_{j=0}^{k-1}p*dp[i-1][j]*(k+j+\frac{p}{q})$$
然后……就做完了（注意这时dp[0][0]=1/p事后的情况）。
一开始我是直接dp[0][0]=1.这样会发现dp[0][0]可以转移到dp[0][0]，我当时没有处理，导致事后统计答案的事后需要乘上1/p。（如果你令dp[0][0]=1的话以为着什么都不放，也就是这样做出来的答案第一个字符都是a，但是你也可以放若干个b，产生 $\sum_{i=0}^{+oo}q^i=\frac{1}{1-q}=\frac{1}{p}$的概率贡献）。但事实上dp[0][0]应该设为1/p最合适，即考虑先放若干p的情况。
然后……转移时dp[i][j]=p*dp[i-1][j]+q*dp[i][j-i]。
其实写成dp[i+1][j]+=p*dp[i][j],dp[i][i+j]+=q*dp[i][j]，这样i+j>=k的事后直接ans+=dp[i][j] q(i+j)即可。
然后……粘代码（写丑了成了O(n^2lgn)）。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#define N 1010
#define mod 1000000007
#define inv(x) fast_pow(x,mod-2)
#define lint long long
using namespace std;
int dp[N][N];
int fast_pow(int x,int k)
{
    if(k==0) return 1;
    if(k==1) return x%mod;
    int ans=fast_pow(x,k>>1);
    ans=(lint)ans*ans%mod;
    if(k&1) ans=(lint)ans*x%mod;
    return ans;
}
int main()
{
    int k,a,b;scanf("%d%d%d",&k,&a,&b);
    int p=(lint)a*inv(a+b)%mod,q=(lint)b*inv(a+b)%mod;
    dp[0][0]=1;int ans=0;
    for(int i=0,x;i<k;i++)
        for(int j=0;j<=k;j++)
            if(dp[i][j])
            {
                x=dp[i][j],(dp[i+1][j]+=(lint)p*x%mod)%=mod;
                if(i+j<k) (dp[i][i+j]+=(lint)q*x%mod)%=mod;
                else (ans+=(lint)inv(p)*q%mod*x%mod*(i+j)%mod)%=mod;
            }
    dp[0][0]=1;int x=((lint)p*inv(q)%mod+k)%mod;
    for(int i=0;i<k;i++)
        (ans+=(lint)dp[k-1][i]%mod*((i+x)%mod)%mod)%=mod;
    printf("%d\n",ans);return 0;
}