ICPC2016 青岛 D - Lucky Coins（概率推公式）

题意：

给出k类硬币和每类硬币的数量和正面朝上的概率，问每类硬币成为幸运硬币的概率是多少。成为幸运硬币的条件是每投一次将所有背面朝上的硬币去掉，继续抛掷，直至剩下一种或者一个都不剩下，那最后一种留下的硬币就是幸运硬币。

分析：

这题公式推导比较新颖。
一般的划分状态在此题因为状态过多开销巨大而不适用，所以考虑以每次投掷硬币为依据划分阶段。

设 $die[i][j]$ 表示第 $i$ 类硬币在第 $j$ 次投掷（包括第 $j$ 次）前就已经被抛弃的概率；
易得：$die[i][j]=(1-p[i{]}^j{)}^{num[i]}$

方便起见，再定一个 $alive[i][j]$，表示第 $i$ 类硬币在第 $j$ 次投掷（包括第 $j$ 次）后至少有一个硬币还存在的概率；
显然：$alive[i][j]=1-die[i][j]$

第 $i$ 类硬币在第 $k$ 轮投掷后成为幸运硬币的概率就是：
$alive[i][k]\times ({\prod }_{j=1,j\ne i}^{n}die[j][k]-{\prod }_{j=1,j\ne i}^{n}die[j][k-1])$
其中 $({\prod }_{j=1,j\ne i}^{n}die[j][k]-{\prod }_{j=1,j\ne i}^{n}die[j][k-1])$ 表示的是其他硬币在第 $k-1$ 轮还至少存在一个，到了第 $k$ 轮全部被抛弃的概率。

因此第 $i$ 类硬币成为幸运硬币的概率就是：
${\sum }_{k=1}^{\infty }alive[i][k]\times ({\prod }_{j=1,j\ne i}^{n}die[j][k]-{\prod }_{j=1,j\ne i}^{n}die[j][k-1])$

这里需要注意特判一下 $n=1$（$n$ 种硬币） 时即只有 $1$ 类硬币的情况。

值得推敲的就是求和公式的 $\infty$ 符号，真的要取到很大很大，大到无穷项吗？ 其实是不需要的，每次投掷硬币有 0.4~0.6 的概率是正面向上，那么当进行几十轮后，这个硬币仍然向上的概率是很小，所以这个求和式很快就收敛了。

因此最终答案就是：
$ans[i]={\sum }_{k=1}^{100}alive[i][k]\times ({\prod }_{j=1,j\ne i}^{n}die[j][k]-{\prod }_{j=1,j\ne i}^{n}die[j][k-1])$

---

其实看很多大佬的博客最后答案写的是：

这个式子表示第 $step$ 步时,其他硬币在 $1$ 到 $step$ 步内死光了,只剩下 $i$ 硬币存活,继续扔, $i$ 硬币在 $step+1$ 步时死光了。我个人觉得是比较难以理解其中的含义。但是实际上我之前写的式子经过变形后也可以得到这个式子，各位读者可以自行选择理解。

代码：

#include<bits/stdc++.h>
using namespace std;

const int N = 15, M = 105;

double p[N], alive[N][M], die[N][M];
int num[N];

double qmi(double a, int b)
{
    double ans = 1;
    for (; b; b >>= 1)
    {
        if (b & 1) ans = ans * a;
        a = a * a;
    }
    return ans;
}

int main()
{
    int t;
    scanf("%d", &t);
    while (t--)
    {
        int n;
        scanf("%d", &n);
        for (int i = 1; i <= n; i++)
            scanf("%d%lf", &num[i], &p[i]);
        
        if (n == 1) puts("1.000000");
        else
        {
            for (int i = 1; i <= n; i++) // 预处理die, alive
            {
                double P = 1;
                for (int j = 0; j <= 100; j++)
                {
                    die[i][j] = qmi(1 - P, num[i]);
                    alive[i][j] = 1 - die[i][j];
                    P *= p[i];
                } 
            }

            for (int i = 1; i <= n; i++)
            {
                double ans = 0;
                for (int k = 1; k <= 100; k++)
                {
                    double x = 1, y = 1;
                    for (int j = 1; j <= n; j++)
                    {
                        if (j == i) continue;
                        x *= die[j][k];
                        y *= die[j][k - 1];
                    }
                    ans += alive[i][k] * (x - y);
                }
                printf("%.6f", ans);
                if (i == n) printf("\n");
                else printf(" ");
            }
        }
    }
    return 0;
}