[校内模拟]光影交错

Description

有一个盒子，每个时刻有pl的概率往里面放一个白球，有pd的概率往里面放一个黑球，1-pl-pd的概率什么都不干，然后每个时刻末尾有p的概率直接结束过程
问所有满足“白球数量大于黑球数量”的时刻的数量的期望
所有读入的实数均保留5位小数

Solution

其实是JS的省队集训
考虑设f[i]表示，所有满足有i个球的时刻的期望，g[i]表示，i个球中，白球>黑球的概率
那么$Ans={\sum }_{i\ge 1}f[i]\ast g[i]$
设F(x)为f的生成函数，我们有$$F(x)=\sum _{i\ge 1}(1-p{)}^{i-1}[(pl+pd)x+(1-pl-pd){]}^i$$
设$a=1-p,b=(pl+pd)x+1-pl-pd$，那么$F(x)={\sum }_{i\ge 1}{a}^{i-1}{b}^i=\frac{b}{1-ab}$
考虑将F(x)写成$\frac{b}{c-dx}=\frac{b/c}{1-d/cx}$，可以发现F(x)存在一个一阶的线性递推
只需要求出$f[1]={\sum }_{i\ge 1}(1-p{)}^{i-1}(1-pl-pd{)}^{i-1}i(pl+pd)=\frac{pl+pd}{(1-(1-p)(1-pl-pd){)}^2}$，即可递推计算
注意到当x大于某个界的时候f[x]已经很小了，所以只用保留前面的项
再考虑g[i]，注意这里计算的是条件概率需要/总概率，即$p{l}^{\prime }=pl/(pl+pd),p{d}^{\prime }=pd/(pl+pd)$
当n为奇数时，$g[n]=g[n-1]\ast (p{l}^{\prime }+p{r}^{\prime })+p{l}^{\prime }\ast \left(\genfrac{}{}{0ex}{}{n-1}{\frac{n-1}{2}}\right)\ast p{l}^{\prime \frac{n-1}{2}}\ast p{r}^{\prime \frac{n-1}{2}}$，表示加上前n-1白=黑，且第n个为白
当n为偶数时，$g[n]=g[n-1]\ast (p{l}^{\prime }+p{r}^{\prime })+p{r}^{\prime }\ast \left(\genfrac{}{}{0ex}{}{n-1}{\frac{n}{2}}\right)p{l}^{\prime \frac{n}{2}}\ast p{r}^{\prime \frac{n}{2}-1}$，表示减去前n-1里白=黑+1，且第n个为黑
后面的组合数部分也可以直接递推
然后就做完了

Code

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
using namespace std;

typedef double db;

const int N=1e7+5;
const db eps=1e-15;

int id,ty;
db pl,pd,p,f[N],g[N],c[N];

int main() {
	freopen("augury.in","r",stdin);
	freopen("augury.out","w",stdout);
	scanf("%d",&id);
	for(scanf("%d",&ty);ty;ty--) {
		scanf("%lf%lf%lf",&pl,&pd,&p);
		db P=(1-p)*(1-pl-pd);
		db a=1-(1-p)*(1-pl-pd),b=(1-p)*(pl+pd);a=b/a;
		f[1]=1/(1-P)/(1-P)*(pl+pd);
		int n=1;
		while (1) {
			n++;
			f[n]=f[n-1]*a;
			if (f[n]<eps) break;
		}
		a=pl/(pl+pd);b=pd/(pl+pd);
		db ans=0;
		fo(i,0,n) {
			if (i&1) {
				if (i==1) c[i]=a;
				else c[i]=c[i-2]*a*b*(i-1)*i/(i>>1)/((i+1)>>1);
				g[i]=g[i-1]+c[i-1]*a;
			} else {
				if (!i) c[i]=1;
				else c[i]=c[i-2]*4*a*b*(i-1)/i;
				g[i]=g[i-1]-c[i-1]*b;
			}
			ans+=f[i]*g[i];
		}
		printf("%.10lf\n",ans);
	}
	return 0;
}