[氧化镍]阴阳两隔

题目

题目描述
你正在听 $\mathsf{T}\mathsf{i}\mathsf{w}$ 讲题。祂有时候讲阴间生成函数，有时候讲阳间（相对而言）数学推导，当然也可能讲了些完全超出你能力范围的东西（显然此类应为绝大多数）。

你认为 $\mathsf{T}\mathsf{i}\mathsf{w}$ 在每讲完一个东西之后，有 $p$ 的概率立刻停止。你认为某个时刻是 “勉强知道自己在干嘛” 的时刻，当且仅当目前 $\mathsf{T}\mathsf{i}\mathsf{w}$ 讲过的阳间数学推导严格多于阴间生成函数。

现在，请问在 $\mathsf{T}\mathsf{i}\mathsf{w}$ 一次讲题的过程中，“勉强知道自己在干嘛” 的时刻的期望数量。可以认为 $\mathsf{T}\mathsf{i}\mathsf{w}$ 每次讲一个东西的时候，有 $p_1$ 的概率是阳间的，有 $p_2$ 的概率是阴间的。（显然有 $1-p_1-p_2$ 的概率讲了一个完全听不懂的东西。）

数据范围与提示
$1\ge p>0$ 且 $p_1,p_2,1-p_1-p_2$ 均为非负。输入的数字不超过 $5$ 位小数。

答案的绝对误差或者相对误差需要不超过 $10^{-8}$ 。

思路

前言：经此一题，我终于意识到，神明是存在的，祂就在我身边。我是 $\mathsf{T}\mathsf{L}\mathsf{Y}$ 太阳神教 的忠实信徒。

那么我们直接上数学推导吧。我希望这不会让你们感到头大。

第一步，随意的考虑某个 “知道自己在跪着听” 的时刻，三种情况（阳间、阴间、听不懂）的个数，然后算概率（因为贡献是 $1$，概率就是期望）。
$$ans=\sum _{a>b\ge 0,c\ge 0}{p}_1^a{p}_2^b(1-p_1-p_2{)}^c(1-p{)}^{a+b+c-1}\left(\genfrac{}{}{0ex}{}{a+b+c}{a+b}\right)\left(\genfrac{}{}{0ex}{}{a+b}{a}\right)$$

神谕：“显然我们会想到把这个东西放进去。” 即，记 $p_1^{\prime }=p_1(1-p)$，然后 $p_2^{\prime },p_3^{\prime }$ 同理，那就可以直接写成 $p_1^{\prime }$ 等的乘积。为了书写方便，一律去掉 ' 符号，还请留心。这里 $p_3$ 对应 $1-p_1-p_2$ 。
$$ans=\frac{1}{1-p}\sum _{a>b\ge 0,c\ge 0}{p}_1^a{p}_2^b\left(\genfrac{}{}{0ex}{}{a+b}{a}\right)p_3^c\left(\genfrac{}{}{0ex}{}{a+b+c}{a+b}\right)$$

这一步需要小心的是，$p=1$ 是有可能的！所以需要特判掉。而从这里开始，所有的 $p_{?}$ 都在 $[0,1)$ 范围内，下面的定义域几乎都满足。

神谕：“发现 $c$ 与另外两个独立，想办法去掉它。如果你对它比较熟悉，你可以直接指出：它就是 $\frac{1}{(1-x{)}^{a+b+1}}$ 的生成函数开形式。” 即，$\frac{1}{(1-x{)}^{a+b+1}}={\sum }_{i=0}^{+\infty }\left(\genfrac{}{}{0ex}{}{i+a+b}{a+b}\right)x^i(x<1)$，令 $x=p_3$，枚举变量 $i$ 和 $c$ 是等价的。
$$ans=\frac{1}{1-p}\sum _{a>b\ge 0}{p}_1^a{p}_2^b\left(\genfrac{}{}{0ex}{}{a+b}{a}\right)\frac{1}{(1-p_3{)}^{a+b+1}}$$

又到了愉快的 “塞进去” 时间！为了书写方便，仍然去掉 ' 符号。
$$ans=\frac{1}{(1-p)(1-p_3)}\sum _{a>b\ge 0}{p}_1^a{p}_2^b\left(\genfrac{}{}{0ex}{}{a+b}{a}\right)$$

做到这里，其实已经有了一个完美的解决方案。设输出的概率是 $P_1,P_2$，令 $P_3=1-P_1-P_2$，式中 $p_1=\frac{P_1(1-p)}{1-P_3(1-p)}$，令 $q=\frac{1}{1-p}>1$，运用一下糖水不等式则有
$$p_1=\frac{P_1}{q-P_3}\le \frac{1-P_3}{q-P_3}\stackrel{\stackrel{sugar}{\mathrm{water}}}{\le }\frac{1}{q}=1-p$$

$p_2$ 同理。而 $1-p\le 0.99999$，在 $a+b>3\times 10^6$ 时，幂已经小于 $10^{-12}$ 了，不会引起误差。而 $a+b\le 3\times 10^6$ 时，可以很简单的 $\mathtt{d}\mathtt{p}$ 求概率。这种方法的代码我附在最后，好奇的可以自行研究。

可是这里就体现出神与人的差距了。神无意吊打我，祂只是不自知地吊打了出题人。更不幸的是：$\mathsf{T}\mathsf{L}\mathsf{Y}$ 就是神。

记 $S_a={\sum }_{b=0}^{a-1}\left(\genfrac{}{}{0ex}{}{a+b}{a}\right)p_2^b$，虽然这里并不是一个直接的生成函数展开式，但是我们可以 “错位相减”。这可能就是神的直觉吧！
$$\begin{gathered} S_a-p_2\cdot S_a=\sum _{b=0}^{a-1}\left[\left(\genfrac{}{}{0ex}{}{a+b}{a}\right)-\left(\genfrac{}{}{0ex}{}{a+b-1}{a}\right)\right]\cdot p_2^b-\left(\genfrac{}{}{0ex}{}{2a-1}{a}\right)p_2^a \\ =\sum _{b=0}^{a-1}\left(\genfrac{}{}{0ex}{}{a-1+b}{a-1}\right)\cdot p_2^b-(\cdots )=S_{a-1}+\left(\genfrac{}{}{0ex}{}{2a-2}{a-1}\right){p_2}^{a-1}-(\cdots ) \end{gathered}$$

为了防止过长，有所省略。令
$$r_a=\left(\genfrac{}{}{0ex}{}{2a-2}{a-1}\right)p_2^{a-1}-\left(\genfrac{}{}{0ex}{}{2a-1}{a}\right)p_2^a$$

递推式简写为
$$S_a=\frac{S_{a-1}+r_a}{1-p_2}$$

一个很不错的 $\mathrm{o}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$ 是 $1-p_2$ 是一个定值。再考虑到它的递归出口 $\{\begin{array}{l}{S}_0=0\\ r_1=1-p_2\end{array}$，显然我们可以直接考虑每个 $r_a$ 最后得到的系数。所以
$$\begin{gathered} (1-p)(1-p_3)\cdot \mathrm{a}\mathrm{n}\mathrm{s}=\sum _{a=1}^{+\infty }{p}_1^a{S}_a \\ =\sum _{a=1}^{+\infty }{p}_1^a\sum _{i=1}^a\frac{r_i}{(1-p_2{)}^{a-i+1}} \\ =\sum _{i=1}^{+\infty }{r}_i{p}_1^i\sum _{a=i}^{+\infty }\frac{p_1^{a-i}}{(1-p_2{)}^{a-i+1}} \\ =\frac{1}{1-p_2}\sum _{i=1}^{+\infty }{r}_i{p}_1^i\frac{1-p_2}{1-p_2-p_1} \\ =\frac{1}{1-p_1-p_2}\sum _{i=1}^{+\infty }{r}_i{p}_1^i \end{gathered}$$

那个地方，公比 $\frac{p_1}{1-p_2}$ 可靠吗？当然。因为 $p_1+p_2<1$，具体原因就去 $\mathrm{s}\mathrm{u}\mathrm{g}\mathrm{a}\mathrm{r}\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}$ 那一节比对吧。

于是我们竟然只需要求 ${\sum }_{i=1}^{+\infty }{r}_i{p}_1^i$ 了耶！两个组合数分开考虑。第一项：
$$\sum _{i=1}^{+\infty }\left(\genfrac{}{}{0ex}{}{2i-2}{i-1}\right)p_2^{i-1}{p}_1^i$$

这东西似乎和卡特兰数很像。考虑用生成函数化简，记
$$C_1(x)=\sum _{i=0}^{+\infty }\left(\genfrac{}{}{0ex}{}{2i}{i}\right)\cdot x^i$$

同样试着写递推式。现实意义仍然是 $n$ 个 $0$、$n$ 个 $1$ 的序列。还是可以考虑前缀 $0,1$ 数量相同之处。不妨设第一位是 $0$，将其反转就是第一位为 $1$ 的。那么前半部分是卡特兰数，后半部分还是 $C_1$ 。记卡特兰数的生成函数 $\frac{1-\sqrt{1-4x}}{2x}=C_0(x)$，则有
$$\begin{gathered} C_1(x)=2x\cdot C_0(x)\cdot C_1(x)+1 \\ \Rightarrow C_1(x)=\frac{1}{\sqrt{1-4x}} \end{gathered}$$

第二项是要减去的，绝对值为
$$\begin{gathered} \sum _{i=1}^{+\infty }\left(\genfrac{}{}{0ex}{}{2i-1}{i}\right)\cdot p_2^i\cdot p_1^i \\ {C}_2(x)=\sum _{i=1}^{+\infty }\left(\genfrac{}{}{0ex}{}{2i-1}{i}\right)\cdot x^i \end{gathered}$$

稍微麻烦点。它是在坐标系内行走到 $(i,i-1)$ 的方案数。如果第一步是向右，那么它就是上面的 $C_1$ 的情形；如果第一步是向上，那么必然经过 $\ell :x=y$ 这条直线，分成两部分，前半部分是卡特兰数，后半部分又是 $C_2$ 。所以
$$\begin{gathered} C_2(x)=x\cdot C_1(x)+x\cdot C_0(x)\cdot C_2(x) \\ \Rightarrow C_2(x)=\frac{2x}{1-4x+\sqrt{1-4x}} \end{gathered}$$

于是终于写出
$$ans=\frac{p_1\cdot C_1(p_1{p}_2)-C_2(p_1{p}_2)}{(1-p)(1-p_3)(1-p_1-p_2)}$$

当然，这里 $p_1{p}_2\le \frac{1}{4}$ 才行欸！其实这是成立的。因为 $\{\begin{array}{l}{p}_1=\frac{P_1(1-p)}{1-(1-P_1-P_2)(1-p)}\\ p_2=\frac{P_2(1-p)}{1-(1-P_1-P_2)(1-p)}\\ p_3=(1-P_1-P_2)(1-p)\end{array}$，那么你很容易发现
$$\begin{gathered} q=\frac{1}{1-p}>1,P_3=1-P_1-P_2\le 1 \\ {p}_1{p}_2=\frac{P_1{P}_2}{(q-P_3{)}^2}\stackrel{A.G.}{\le }\frac{(\frac{P_1+P_2}{2}{)}^2}{(q-P_3{)}^2}\le \frac{(1-P_3{)}^2}{4(q-P_3{)}^2}\stackrel{\stackrel{sugar}{\mathrm{water}}}{\le }\frac{1}{4q^2}<\frac{1}{4} \end{gathered}$$

如果你愿意的话，你也可以把 $ans$ 写成关于输入数据的式子。

就这样，这道题就被 无精度误差（除了 $\mathtt{d}\mathtt{o}\mathtt{u}\mathtt{b}\mathtt{l}\mathtt{e}$ 自动丢失精度）算出来了。$\mathsf{T}\mathsf{L}\mathsf{Y}$ 太阳神 永生！

代码

#include <bits/stdc++.h>
int main(){
	freopen("augury.in","r",stdin);
	freopen("augury.out","w",stdout);
	double p1, p2, p3, p; int T;
	for(scanf("%*d%d",&T); T; --T){
		scanf("%lf %lf %lf",&p1,&p2,&p);
		if(p == 1){ // must end immediately
			printf("%.10f\n",p1); continue;
		} else p3 = 1-p1-p2;
		double ans = 1/(1-p)/(1-p3*(1-p));
		p1 /= 1/(1-p)-p3, p2 /= 1/(1-p)-p3;
		double k = 1-4*p1*p2; // 1-4x
		ans *= p1/sqrt(k)-2*p1*p2/(k+sqrt(k));
		printf("%.10f\n",ans/(1-p1-p2));
	}
	return 0;
}

以及那个卡精度的方法：

#include <cstdio>
#include <iostream>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
inline int readint(){
	int a = 0; char c = getchar(), f = 1;
	for(; c<'0'||c>'9'; c=getchar())
		if(c == '-') f = -f;
	for(; '0'<=c&&c<='9'; c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}
double qkpow(double b,int q){
	double a = 1;
	for(; q; q>>=1,b*=b)
		if(q&1) a *= b;
	return a;
}

const int MaxN = 3000005;
double f[MaxN]; 

int main(){
	readint();
	for(int T=readint(); T; --T){
		double p1, p2, p;
		scanf("%lf %lf %lf",&p1,&p2,&p);
		if(p == 1){ // must end immediately
			printf("%.10f\n",p1);
			continue;
		}
		double p3 = 1-p1-p2;
		p1 *= (1-p), p2 *= (1-p), p3 *= (1-p);
		p1 /= (1-p3), p2 /= (1-p3);
		f[0] = 0, f[1] = p1, f[2] = p1*p1;
		// tmp = C(i-1,i/2)*pow(p1,i/2)*pow(p2,i/2)
		double tmp = p1*p2;
		double ans = f[1]+f[2];
		for(int i=3; i<MaxN; ++i){
			if(!(i&1)) tmp *= p1*p2*
				(i-1-(i>>1))/(i>>1);
			tmp *= (i-1); tmp /= (i-1-(i/2));
			if(i&1) f[i] = f[i-1]*(p1+p2)+tmp*p1;
			else f[i] = f[i-1]*(p1+p2)-tmp;
			ans += f[i]; // add them up
		}
		ans /= (1-p), ans /= (1-p3);
		printf("%.10f\n",ans);
	}
	return 0;
}