光影交錯

題目傳送門

Description

Solution

可以看出，假設我們設定一個臨界點 \(n\) ，當 \(n\) 足夠大的時候，\(n\) 步操作之後對答案的影響就不再精度考慮範圍之類了。

我們設 \(f(i)\) 分別表示 \(i\) 次操作時非中性靈氣的期望出現次數，\(g(i)\) 表示 \(i\) 次非中性靈氣是總體陽間的概率。

那麼答案就是 \(\sum_{i=1}^{n} f(i)\times g(i)\) 。

可以得到：

\[g(x)=\sum_{i=\lfloor\frac{x}{2}\rfloor+1}^{x} \binom{x}{i}qL^iqR^{x-i} \]

我們可以得到遞推式：

\[\left\{\begin{array}{l} g(x)=0,\text{when }x=0\\ g(x)=g(x-1)-pD\binom{x-1}{\lfloor\frac{x}{2}\rfloor} pL^{x/2}pD^{x/2-1},\text{when }x\equiv 0\pmod{2}\\ g(x)=g(x-1)+pL\binom{x-1}{\lfloor\frac{x}{2}\rfloor} pL^{(x-1)/2}pD^{(x-1)/2},\text{when }x\equiv 1\pmod{2} \end{array}\right. \]

考慮如何求解 \(f(x)\)，可以考慮對其構造生成函式 \(F(x)\)

。那麼我們就有：

\[F(x)=\sum_{i=0}^{\infty} f(i)\times x^i=\sum_{i=1}^{\infty} (1-p)^{i-1}((pL+pD)x+(1-pL-pD))^i=\frac{1-p}{1-(1-p)[(pL+pD)x+(1-pL-pD)]} \]

然後你就可以得到：

\[f(x)=\frac{f(x-1)(1-p)(pL+pD)}{1-(1-p)(1-pL-pD)} \]

那麼你就可以 \(\Theta(n)\) 解決這個問題，\(n\) 大概是 \(1e7\) 左右。

Code

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

#define Int register int
#define MAXN 10000005

template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c != ' ' && c != '\n') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}
template <typename T,typename ... Args> inline void read (T &t,Args&... args){read (t);read (args...);}
template <typename T> inline void write (T x){if (x < 0){x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}
template <typename T> void chkmax (T &a,T b){a = max (a,b);}
template <typename T> void chkmin (T &a,T b){a = min (a,b);}

#define eps 1e-10
double p,pl,pd,f[MAXN],g[MAXN],drw[MAXN],tmp[MAXN];

void work (){
	scanf ("%lf%lf%lf",&pl,&pd,&p);
	double now = 1,ans = 0;int n = 0;tmp[0] = 1,tmp[1] = f[0] = f[1] = 0;
	while (now > eps){
		for (Int i = 1;~i;-- i){
			tmp[i] *= (1 - pl - pd);
			if (i) tmp[i] += tmp[i - 1] * (pl + pd);
			f[i] += tmp[i] * now;
		}
		now *= (1 - p),n ++;
	}
	double h = 1 - (1 - p) * (1 - pl - pd);
	for (Int i = 2;i <= n;++ i) f[i] = f[i - 1] * (1 - p) * (pl + pd) / h;
	double sum = pl + pd;pl /= sum,pd /= sum;
	for (Int i = 0;i <= n;++ i){
		if (i & 1){
			g[i] = g[i - 1] + drw[i - 1] * pl;
			if (i == 1) drw[i] = pl;
			else drw[i] = drw[i - 2] * pl * pd * i * (i - 1) / (i / 2) / (i / 2 + 1);
			ans += g[i] * f[i];
		}
		else{
			if (i == 0) drw[i] = 1;
			else drw[i] = 4 * drw[i - 2] * (i - 1) / i * pl * pd;
			if (i) g[i] = g[i - 1] - drw[i - 1] * pd;
			ans += g[i] * f[i];
		}
	}
	printf ("%.10f\n",ans);
}

signed main(){
	freopen ("augury.in","r",stdin);
	freopen ("augury.out","w",stdout);
	int num,T;read (num,T);
	while (T --> 0) work ();
	return 0;
}