hdu 4089 Activation

/*
dp[i][j] 表示当前总共有i人， 番茄站在第j位， 到达目标状态的概率。
dp[1][1] = p1*dp[1][1] + p2*dp[1][1] + p4 = p4 / (1 - p1 - p2)
j == 1时，
dp[i][1] = dp[i][1]*p1 + p2*dp[i][i] + p4;
2<=j<=k时,
dp[i][j] = dp[i][j]*p1 + p2*dp[i][j-1] + p3*dp[i-1][j-1] + p4;
k < j <=i
dp[i][j] = dp[i][j]*p1 + p2*dp[i][j-1] + p3*dp[i-1][j-1];

令
p21 = (p2 / (1 - p1));
p31 = (p3 / (1 - p1));
p41 = (p4 / (1 - p1));
先化简
dp[i][1] =  p21*dp[i][i] + p41;
dp[i][j] =  p21*dp[i][j-1] + p31*dp[i-1][j-1] + p41;
dp[i][j] =  p21*dp[i][j-1] + p31*dp[i-1][j-1];

再令
j <= k,
c[j] = p31*dp[i-1][j-1] + p41;
j > k
c[j] = p31*dp[i-1][j-1];
则 j>1, dp[i][j] = p21*dp[i-1][j-1] + c[j];

发现这个递推式是循环的， 采用迭代先将dp[i][i]先求出来
dp[i][i] = dp[i-1][i-1]*p21 + c[i];
dp[i][i-1]*p21 = dp[i-1][i-2]*p21^2 + c[i-1]*p21;  //左右两边同乘p21
...
dp[i][2]*p21^(i-2) = dp[i][1]*p21^(i-1) + c[2]*p21^(i-2);  // 左右两边同乘p21^(i-2)
dp[i][1]*p21^(i-1) = dp[i][i]*p21^(i) + c[1]*p21^(i-1);  // 左右两边同乘p21^(i-1)

将上式左右两边同时相加，
dp[i][i] = dp[i][i]*p21^(i) + c[i] + c[i-1]*p21 + ... + c[1]*p21^(i-1).

求出dp[i][i]后， 根据递推式求出dp[i][j].
最后dp[n][m]就是答案。
*/

#include <cstdio>
#include <cstring>
#include <iostream>
#include <queue>
#include <vector>
#include <map>
#include <set>
#include <cmath>
#include <algorithm>
#include <functional>
#include <cmath>
#include <bitset>
using namespace std;
#define sqr(x) ((x)*(x))

const double eps = 1e-10;
const int maxn = 2222;
double dp[maxn][maxn], c[maxn];

int n,m,k;
double p1,p2,p3,p4;
double p21,p31,p41;

int sgn(double x)
{
    return (x > eps) - (x < -eps);
}

int main()
{
    while(~scanf("%d%d%d",&n,&m,&k))
    {
        scanf("%lf%lf%lf%lf", &p1,&p2,&p3,&p4);
        if(sgn(p4) == 0)
        {
            puts("0.00000");
            continue;
        }
        p21 = (p2 / (1 - p1));
        p31 = (p3 / (1 - p1));
        p41 = (p4 / (1 - p1));
        dp[1][1] = p4 / (1 - p1 - p2);
        c[1] = p41;
        for(int i=2; i <= n; i++)
        {
            for(int j=2; j<=i; j++)
            {
                if( j <= k) c[j] = dp[i-1][j-1]*p31 + p41;
                else c[j] = dp[i-1][j-1]*p31;
            }

            double tmp = 0, p = 1;
            for(int j=i; j >= 1; j--)
            {
                tmp += c[j]*p;
                p *= p21;
            }
            dp[i][i] = tmp / (1 - p);
            dp[i][1] = dp[i][i]*p21 + c[1];
            for(int j=2; j<i; j++)
            {
                dp[i][j] = dp[i][j-1]*p21 + c[j];
            }
        }

        printf("%.5f\n", dp[n][m]);
    }
    return 0;
}