HDU 4465繰返しとdoubleの精度

問題は大まかに言わない
 
ここでdp[i][0]は、1つ目のボックスを取った後に2つ目のボックスがi個残る確率を表し、対応する期待はdp[i][0]*iである
dp[i][1]は、2つ目の箱を取った後、1つ目の箱がi個残る確率を表す
 
dp[i][0]  =  p^(n+1) * (1-p)^(n-i) * C(2*n-i , n-i) = p^(n+1) * (1-p)^(n-i) * (2*n-i)!/(n-i)!/n!
dp[i+1][0]  = p^(n+1) * (1-p)^(n-i-1) * C(2*n-i-1 , n-i-1) = p^(n+1) * (1-p)^(n-i-1) * (2*n-i-1)!/(n-i-1)!/n!
 
dp[i][0] = dp[i+1][0] * (1-p) * (2*n-i)/(n-i)
 
dp[i][1]も同じ理屈
dp[n][0]に初期値pow(p,n+1)を最初に付与すると、nが大きすぎると、精度の問題で0が得られる
したがってn+1個のpは計算中に答えが総数を超えたときに1個ずつ乗算される
 

 1 #include <iostream>

 2 #include <cstdio>

 3 #include <cmath>

 4 using namespace std;

 5 const int N = 200005;

 6 double dp[N][2];

 7 

 8 int main()

 9 {

10    // freopen("test.in","rb",stdin);

11 

12     int n,cas = 0;

13     double p;

14     while(scanf("%d%lf",&n,&p)!=EOF){

15         double ans1 = 0;

16         double ans2 = 0;

17         int last[2];

18         last[0] = last[1] = n+1;

19         double q = 1-p;

20         dp[n][0] = 1;

21         dp[n][1] = 1;

22         ans1 += n;

23         ans2 += n;

24         for(int i=n-1;i>=1;i--){

25             double tmp1 =  q * (2*n-i) / (n-i);

26             double tmp2 =  p * (2*n-i) / (n-i);

27             

28             //    ，      ，p    ，                 0

29             //                 p 

30             dp[i][0] = tmp1 * dp[i+1][0];

31             dp[i][1] = tmp2 * dp[i+1][1];

32             ans1 += (dp[i][0] * i);

33             while(ans1>n){

34                 dp[i][0] *= p;

35                 ans1 *= p;

36                 last[0]--;

37             }

38             ans2 += (dp[i][1] * i);

39             while(ans2>n){

40                 dp[i][1]*=q;

41                 ans2 *= q;

42                 last[1]--;

43             }

44         }

45        // cout<<"  dp "<<dp[n]<<endl;

46         ans1 *= pow(p,last[0]);

47         ans2 *= pow(q,last[1]);

48         cas++;

49 

50         printf("Case %d: %.6f
",cas,ans1+ans2);

51     }

52     return 0;

53 }