Let
$$Z=\sum_{n\ge 1}\frac{B_n}{n^2},\qquad B_n\in\{0,1\},\qquad \Pr(B_n=0)=\Pr(B_n=1)=\tfrac12.$$
The problem asks for the tail probability
$$p(a)=\Pr(Z\gt a)$$
at \(a=0.5\), to high precision. Since \(0\le Z\le \zeta(2)=\pi^2/6\) and the series contains infinitely many independent bits, direct enumeration is hopeless. The implementation instead rewrites the distribution in a symmetric form, derives its characteristic function, and then recovers the probability by Fourier inversion.
The whole method is a clean chain: center the random series, compute the characteristic function term by term, convert the tail probability into a one-dimensional integral, and then approximate that integral numerically.
Each random bit contributes either \(0\) or \(1/n^2\), so its mean is \(1/(2n^2)\). Therefore
$$\mu=\mathbb E[Z]=\sum_{n\ge 1}\frac{1}{2n^2}=\frac{\zeta(2)}{2}=\frac{\pi^2}{12}.$$
Write \(B_n=\frac{1+\varepsilon_n}{2}\), where \(\varepsilon_n\in\{-1,+1\}\) with equal probability. Then
$$Z=\mu+\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2}.$$
So the centered variable is
$$Y=Z-\mu=\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2},$$
which is symmetric around \(0\). Also every term is nonnegative and the largest possible sum is \(\zeta(2)\), so
$$a\lt 0 \implies p(a)=1,\qquad a\ge \zeta(2)\implies p(a)=0.$$
For a fixed \(n\), the centered summand has characteristic function
$$\mathbb E\left[e^{it\varepsilon_n/(2n^2)}\right]=\frac{e^{it/(2n^2)}+e^{-it/(2n^2)}}{2}=\cos\left(\frac{t}{2n^2}\right).$$
The summands are independent, so the characteristic function of the full series factorizes into an infinite product:
$$\varphi_Y(t)=\prod_{n\ge 1}\cos\left(\frac{t}{2n^2}\right).$$
This function is real and even, reflecting the symmetry \(Y\overset{d}= -Y\). That symmetry is the reason the final inversion formula contains only a sine factor and no imaginary part to evaluate separately.
Set
$$x=\mu-a.$$
Then \(p(a)=\Pr(Y\gt -x)\). Gil-Pelaez inversion gives
$$p(a)=\frac12+\frac{1}{\pi}\int_0^{\infty}\frac{\sin(tx)}{t}\,\varphi_Y(t)\,dt.$$
So an infinite random sum has been reduced to a deterministic one-dimensional integral. The hard combinatorics disappear; the remaining work is numerical analysis.
The implementation cannot keep infinitely many cosine factors, so it replaces \(\varphi_Y(t)\) by
$$\varphi_{Y,N}(t)=\prod_{n=1}^{N}\cos\left(\frac{t}{2n^2}\right),$$
with a large cutoff \(N\). It also replaces the interval \([0,\infty)\) by a finite interval \([0,T]\). This is effective because the omitted cosine factors are extremely close to \(1\) once \(n\) is large, and because the oscillatory integral contributes less and less beyond a sufficiently large upper limit.
If \(T=mh\) with even \(m\), the numerical quadrature uses
$$\int_0^T f(t)\,dt\approx \frac{h}{3}\left(f(0)+f(T)+4\sum_{j=1,3,\dots,m-1}f(jh)+2\sum_{j=2,4,\dots,m-2}f(jh)\right),$$
where
$$f(t)=\frac{\sin(tx)}{t}\,\varphi_{Y,N}(t).$$
At \(t=0\), the expression \(\sin(tx)/t\) is interpreted by continuity:
$$\lim_{t\to 0}\frac{\sin(tx)}{t}=x.$$
That removes the apparent singularity at the origin and makes Simpson's rule stable there.
If every bit \(B_n\) is replaced by \(1-B_n\), then
$$\sum_{n\ge 1}\frac{1-B_n}{n^2}=\zeta(2)-Z.$$
So the distribution of \(Z\) is symmetric around \(\mu=\zeta(2)/2\). It follows that
$$p(\mu)=\Pr(Z\gt \mu)=\frac12,$$
and more generally
$$p(\mu-u)+p(\mu+u)=1.$$
In particular, the thresholds \(1\) and \(\zeta(2)-1\) are symmetric around \(\mu\), so both have probability \(1/2\). With the finer numerical settings used by the implementation, the target quantity comes out as
$$p(0.5)\approx 0.565654540708545,$$
which is the value checked before the final rounded answer is printed.
The C++, Python, and Java implementations all compute the same truncated Fourier integral. The C++ implementation performs the numerical work directly: it precomputes the coefficients \(1/(2n^2)\) up to the cutoff, evaluates the truncated cosine product at each Simpson node, and accumulates the weighted node values over the chosen interval.
The interior Simpson nodes are divided across available threads, so the wall-clock time improves on multi-core hardware while the mathematical result stays unchanged. The endpoint at \(t=0\) is handled through the limit above, and the number of subintervals is adjusted so Simpson's rule always uses an even count.
The Python and Java implementations are thin front ends that delegate to the same compiled numerical core, so all three languages share the same checkpoints, the same truncation strategy, and the same final probability once their numerical parameters agree.
If the integration limit is \(T\), the step size is \(h\), and the cosine-product cutoff is \(N\), then the number of Simpson subintervals is about \(T/h\). Each node evaluation multiplies \(N\) cosine factors, so the total running time is
$$O\left(\frac{T}{h}\,N\right).$$
The memory cost is \(O(N)\) for the precomputed coefficients, plus a small extra buffer for per-thread partial sums. Multithreading changes the constant factor in runtime, not the asymptotic order.
Gegeben ist die Zufallsreihe
$$Z=\sum_{n\ge 1}\frac{B_n}{n^2},\qquad B_n\in\{0,1\},\qquad \Pr(B_n=0)=\Pr(B_n=1)=\tfrac12.$$
Gesucht ist die Schwanzwahrscheinlichkeit
$$p(a)=\Pr(Z\gt a)$$
für \(a=0.5\) mit hoher Genauigkeit. Da \(0\le Z\le \zeta(2)=\pi^2/6\) gilt und unendlich viele unabhängige Binärentscheidungen beteiligt sind, ist vollständiges Durchprobieren unmöglich. Die Lösung bringt die Reihe deshalb in eine symmetrische Form, bestimmt ihre charakteristische Funktion und gewinnt daraus die Wahrscheinlichkeit per Fourier-Inversion.
Die Methode besteht aus vier gedanklichen Schritten: zentrieren, faktorisieren, invertieren und numerisch integrieren.
Jeder Summand trägt im Mittel \(1/(2n^2)\) bei. Also ist
$$\mu=\mathbb E[Z]=\sum_{n\ge 1}\frac{1}{2n^2}=\frac{\zeta(2)}{2}=\frac{\pi^2}{12}.$$
Schreibe \(B_n=\frac{1+\varepsilon_n}{2}\) mit \(\varepsilon_n\in\{-1,+1\}\) gleichverteilt. Dann folgt
$$Z=\mu+\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2}.$$
Damit ist die zentrierte Zufallsvariable
$$Y=Z-\mu=\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2}$$
symmetrisch um \(0\). Zugleich ist \(Z\) immer zwischen \(0\) und \(\zeta(2)\), also gilt sofort
$$a\lt 0 \implies p(a)=1,\qquad a\ge \zeta(2)\implies p(a)=0.$$
Für einen festen Summanden erhält man
$$\mathbb E\left[e^{it\varepsilon_n/(2n^2)}\right]=\frac{e^{it/(2n^2)}+e^{-it/(2n^2)}}{2}=\cos\left(\frac{t}{2n^2}\right).$$
Wegen der Unabhängigkeit multiplizieren sich diese Faktoren, also
$$\varphi_Y(t)=\prod_{n\ge 1}\cos\left(\frac{t}{2n^2}\right).$$
Die Funktion ist reell und gerade. Genau diese Symmetrie reduziert die Inversionsformel später auf ein einziges Sinusintegral.
Setze
$$x=\mu-a.$$
Dann ist \(p(a)=\Pr(Y\gt -x)\). Mit dem Satz von Gil-Pelaez erhält man
$$p(a)=\frac12+\frac{1}{\pi}\int_0^{\infty}\frac{\sin(tx)}{t}\,\varphi_Y(t)\,dt.$$
Damit wird aus einer probabilistischen Frage über eine unendliche Reihe ein eindimensionales bestimmtes Integral.
Numerisch wird \(\varphi_Y(t)\) durch
$$\varphi_{Y,N}(t)=\prod_{n=1}^{N}\cos\left(\frac{t}{2n^2}\right)$$
ersetzt. Außerdem wird nur bis zu einer großen Obergrenze \(T\) integriert. Das ist sinnvoll, weil die ausgelassenen Kosinusfaktoren für große \(n\) fast \(1\) sind und weil der oszillierende Integrand jenseits von \(T\) immer weniger zum Endwert beiträgt.
Ist \(T=mh\) mit geradem \(m\), so verwendet die Implementierung
$$\int_0^T f(t)\,dt\approx \frac{h}{3}\left(f(0)+f(T)+4\sum_{j=1,3,\dots,m-1}f(jh)+2\sum_{j=2,4,\dots,m-2}f(jh)\right),$$
wobei
$$f(t)=\frac{\sin(tx)}{t}\,\varphi_{Y,N}(t).$$
Am Ursprung wird der Grenzwert benutzt:
$$\lim_{t\to 0}\frac{\sin(tx)}{t}=x.$$
Damit verschwindet die scheinbare Singularität bei \(t=0\).
Ersetzt man jedes Bit durch \(1-B_n\), dann wird \(Z\) zu
$$\sum_{n\ge 1}\frac{1-B_n}{n^2}=\zeta(2)-Z.$$
Die Verteilung ist also symmetrisch um \(\mu=\zeta(2)/2\). Deshalb gilt
$$p(\mu)=\Pr(Z\gt \mu)=\frac12,$$
und allgemeiner
$$p(\mu-u)+p(\mu+u)=1.$$
Insbesondere sind \(1\) und \(\zeta(2)-1\) Spiegelpunkte um \(\mu\), also haben beide die Wahrscheinlichkeit \(1/2\). Mit den feineren Rechenparametern ergibt die Implementierung für die eigentliche Aufgabe
$$p(0.5)\approx 0.565654540708545.$$
Die C++-, Python- und Java-Implementierungen berechnen alle dasselbe abgeschnittene Fourier-Integral. Die numerische Hauptarbeit geschieht in C++: Die Werte \(1/(2n^2)\) werden bis zum Cutoff vorab gespeichert, an jedem Simpson-Knoten wird das abgeschnittene Kosinusprodukt ausgewertet, und anschließend werden die Simpson-Gewichte aufsummiert.
Die inneren Knoten des Integrationsgitters werden auf mehrere Threads verteilt. Dadurch sinkt die Laufzeit auf Mehrkernsystemen, ohne dass sich das mathematische Ergebnis ändert. Der Punkt \(t=0\) wird separat über den Grenzwert behandelt, und die Anzahl der Teilintervalle wird bei Bedarf auf eine gerade Zahl angepasst.
Python und Java dienen als schlanke Hüllen um dieselbe kompilierte numerische Routine. Deshalb verwenden alle drei Sprachversionen dieselben Prüfungen, dieselbe Trunkierung und dieselbe Endausgabe, sofern die numerischen Parameter übereinstimmen.
Seien \(T\) die Integrationsobergrenze, \(h\) die Schrittweite und \(N\) der Produkt-Cutoff. Dann gibt es ungefähr \(T/h\) Simpson-Teilintervalle. Pro Integrationsknoten werden \(N\) Kosinusfaktoren multipliziert, also beträgt die Laufzeit
$$O\left(\frac{T}{h}\,N\right).$$
Der Speicherbedarf ist \(O(N)\) für die vorberechneten Koeffizienten sowie ein kleiner Zusatzspeicher für partielle Summen der Threads. Parallelisierung verbessert die Konstante, nicht die asymptotische Ordnung.
Rastgele seri
$$Z=\sum_{n\ge 1}\frac{B_n}{n^2},\qquad B_n\in\{0,1\},\qquad \Pr(B_n=0)=\Pr(B_n=1)=\tfrac12$$
için istenen değer
$$p(a)=\Pr(Z\gt a)$$
olup hedef nokta \(a=0.5\)'tir. \(Z\) her zaman \(0\) ile \(\zeta(2)=\pi^2/6\) arasında kalır, fakat seri sonsuz sayıda bağımsız ikili seçim içerdiği için doğrudan sayım yapılamaz. Çözüm bunun yerine seriyi merkezleyip karakteristik fonksiyonunu çıkarır ve olasılığı Fourier ters dönüşümüyle hesaplar.
Yöntem şu sırayı izler: önce seriyi ortalaması etrafında merkezle, sonra bağımsızlıktan karakteristik fonksiyonu kur, ardından olasılığı tek değişkenli bir integrale indir ve son olarak bu integrali sayısal olarak hesapla.
Her terimin beklenen değeri \(1/(2n^2)\) olduğundan
$$\mu=\mathbb E[Z]=\sum_{n\ge 1}\frac{1}{2n^2}=\frac{\zeta(2)}{2}=\frac{\pi^2}{12}$$
elde edilir. \(B_n=\frac{1+\varepsilon_n}{2}\) ve \(\varepsilon_n\in\{-1,+1\}\) yazarsak
$$Z=\mu+\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2}$$
olur. Böylece merkezlenmiş değişken
$$Y=Z-\mu=\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2}$$
sıfır etrafında simetriktir. Ayrıca tüm terimler negatif olmadığı için
$$a\lt 0 \implies p(a)=1,\qquad a\ge \zeta(2)\implies p(a)=0$$
sınırları hemen çıkar.
Sabit bir \(n\) için ilgili merkezlenmiş terimin karakteristik fonksiyonu
$$\mathbb E\left[e^{it\varepsilon_n/(2n^2)}\right]=\frac{e^{it/(2n^2)}+e^{-it/(2n^2)}}{2}=\cos\left(\frac{t}{2n^2}\right)$$
şeklindedir. Terimler bağımsız olduğundan tüm seri için
$$\varphi_Y(t)=\prod_{n\ge 1}\cos\left(\frac{t}{2n^2}\right)$$
elde edilir. Bu fonksiyonun reel ve çift olması, dağılımın simetrisini doğrudan yansıtır.
Şimdi
$$x=\mu-a$$
tanımını yapalım. O zaman \(p(a)=\Pr(Y\gt -x)\) olur. Gil-Pelaez tersleme formülü ile
$$p(a)=\frac12+\frac{1}{\pi}\int_0^{\infty}\frac{\sin(tx)}{t}\,\varphi_Y(t)\,dt$$
yazılır. Böylece sonsuz rastgele seri problemi, deterministik ve tek boyutlu bir integrale indirgenmiş olur.
Uygulamada sonsuz çarpım
$$\varphi_{Y,N}(t)=\prod_{n=1}^{N}\cos\left(\frac{t}{2n^2}\right)$$
ile yaklaşık alınır. Ayrıca integral \([0,\infty)\) yerine \([0,T]\) üzerinde hesaplanır. Bu yaklaşım mantıklıdır; çünkü büyük \(n\) için atılan kosinüs çarpanları \(1\)'e çok yakındır ve osilasyonlu integrand büyük \(t\) değerlerinde giderek daha az katkı yapar.
\(T=mh\) ve \(m\) çift ise
$$\int_0^T f(t)\,dt\approx \frac{h}{3}\left(f(0)+f(T)+4\sum_{j=1,3,\dots,m-1}f(jh)+2\sum_{j=2,4,\dots,m-2}f(jh)\right)$$
kullanılır; burada
$$f(t)=\frac{\sin(tx)}{t}\,\varphi_{Y,N}(t).$$
\(t=0\) noktasında görünen tekillik gerçek değildir; çünkü
$$\lim_{t\to 0}\frac{\sin(tx)}{t}=x.$$
Bu nedenle başlangıç noktası ayrı ve kararlı biçimde hesaplanabilir.
Her biti \(1-B_n\) ile değiştirirsek
$$\sum_{n\ge 1}\frac{1-B_n}{n^2}=\zeta(2)-Z$$
olur. Bu da dağılımın \(\mu=\zeta(2)/2\) etrafında simetrik olduğunu gösterir. Dolayısıyla
$$p(\mu)=\Pr(Z\gt \mu)=\frac12,$$
ve genel olarak
$$p(\mu-u)+p(\mu+u)=1.$$
Özellikle \(1\) ile \(\zeta(2)-1\) eşiklerinin ikisi de \(\mu\)'ya göre simetriktir; bu yüzden her ikisi için de olasılık \(1/2\)'dir. İnce ayarlı hesaplamada hedef değer
$$p(0.5)\approx 0.565654540708545$$
olarak elde edilir.
C++, Python ve Java uygulamaları aynı kesilmiş Fourier integralini hesaplar. Sayısal yükün tamamı C++ uygulamasındadır: \(1/(2n^2)\) katsayıları önceden hazırlanır, Simpson düğümlerinin her birinde kesilmiş kosinüs çarpımı hesaplanır ve sonra ağırlıklı toplam alınır.
İç Simpson düğümleri mevcut iş parçacıklarına bölünür. Böylece çok çekirdekli makinelerde süre kısalır, fakat elde edilen matematiksel değer değişmez. Başlangıç noktası yukarıdaki limit kullanılarak özel ele alınır ve Simpson yönteminin gereği olarak alt aralık sayısı çift tutulur.
Python ve Java uygulamaları ise aynı derlenmiş sayısal çekirdeğe yönlenen hafif arayüzlerdir. Bu yüzden üç dilde de kullanılan matematik, denetim noktaları ve sonuç aynıdır.
İntegral üst sınırı \(T\), adım genişliği \(h\), çarpım kesimi \(N\) olsun. Simpson alt aralık sayısı yaklaşık \(T/h\) olur. Her düğümde \(N\) adet kosinüs çarpıldığı için toplam zaman maliyeti
$$O\left(\frac{T}{h}\,N\right)$$
mertebesindedir. Bellek kullanımı, önceden hesaplanan katsayılar için \(O(N)\) ve iş parçacığı başına kısmi toplamlar için küçük bir ek maliyetten ibarettir. Paralellik sadece sabit katsayıyı iyileştirir.
Se considera la serie aleatoria
$$Z=\sum_{n\ge 1}\frac{B_n}{n^2},\qquad B_n\in\{0,1\},\qquad \Pr(B_n=0)=\Pr(B_n=1)=\tfrac12.$$
La cantidad buscada es
$$p(a)=\Pr(Z\gt a)$$
en el caso \(a=0.5\). Como \(0\le Z\le \zeta(2)=\pi^2/6\) y la serie depende de infinitas decisiones binarias independientes, no es viable enumerar todos los casos. La estrategia consiste en centrar la serie, obtener su función característica y recuperar la probabilidad mediante inversión de Fourier.
La idea central es transformar un problema de probabilidad infinita en una integral oscilatoria de una sola variable.
Cada término tiene esperanza \(1/(2n^2)\), luego
$$\mu=\mathbb E[Z]=\sum_{n\ge 1}\frac{1}{2n^2}=\frac{\zeta(2)}{2}=\frac{\pi^2}{12}.$$
Si escribimos \(B_n=\frac{1+\varepsilon_n}{2}\) con \(\varepsilon_n\in\{-1,+1\}\) equiprobable, entonces
$$Z=\mu+\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2}.$$
Por tanto, la variable centrada es
$$Y=Z-\mu=\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2},$$
que es simétrica alrededor de \(0\). Además, como todos los términos son no negativos y la suma máxima es \(\zeta(2)\), se obtiene inmediatamente
$$a\lt 0 \implies p(a)=1,\qquad a\ge \zeta(2)\implies p(a)=0.$$
Para un índice fijo \(n\),
$$\mathbb E\left[e^{it\varepsilon_n/(2n^2)}\right]=\frac{e^{it/(2n^2)}+e^{-it/(2n^2)}}{2}=\cos\left(\frac{t}{2n^2}\right).$$
La independencia de los sumandos implica entonces
$$\varphi_Y(t)=\prod_{n\ge 1}\cos\left(\frac{t}{2n^2}\right).$$
La función es real y par, exactamente lo que cabe esperar de una variable centrada y simétrica.
Definimos
$$x=\mu-a.$$
Entonces \(p(a)=\Pr(Y\gt -x)\). La inversión de Gil-Pelaez da
$$p(a)=\frac12+\frac{1}{\pi}\int_0^{\infty}\frac{\sin(tx)}{t}\,\varphi_Y(t)\,dt.$$
Así, el problema original queda reducido a calcular una integral real de una dimensión.
En la práctica se usa
$$\varphi_{Y,N}(t)=\prod_{n=1}^{N}\cos\left(\frac{t}{2n^2}\right)$$
en lugar del producto infinito. También se sustituye \([0,\infty)\) por \([0,T]\). Esto funciona bien porque los factores omitidos se acercan muy deprisa a \(1\) cuando \(n\) crece, y porque la parte oscilatoria del integrando tiene cada vez menos efecto más allá del corte elegido.
Si \(T=mh\) con \(m\) par, la aproximación utilizada es
$$\int_0^T f(t)\,dt\approx \frac{h}{3}\left(f(0)+f(T)+4\sum_{j=1,3,\dots,m-1}f(jh)+2\sum_{j=2,4,\dots,m-2}f(jh)\right),$$
donde
$$f(t)=\frac{\sin(tx)}{t}\,\varphi_{Y,N}(t).$$
En \(t=0\) se usa el límite
$$\lim_{t\to 0}\frac{\sin(tx)}{t}=x,$$
de modo que no queda ninguna singularidad real en el origen.
Si cada bit se reemplaza por \(1-B_n\), entonces
$$\sum_{n\ge 1}\frac{1-B_n}{n^2}=\zeta(2)-Z.$$
Eso muestra que la distribución de \(Z\) es simétrica alrededor de \(\mu=\zeta(2)/2\). Por tanto,
$$p(\mu)=\Pr(Z\gt \mu)=\frac12,$$
y en general
$$p(\mu-u)+p(\mu+u)=1.$$
En particular, \(1\) y \(\zeta(2)-1\) están a la misma distancia de \(\mu\), así que ambos umbrales tienen probabilidad \(1/2\). Con los parámetros finos de integración, el cálculo principal produce
$$p(0.5)\approx 0.565654540708545.$$
Las implementaciones en C++, Python y Java evalúan la misma integral de Fourier truncada. La versión en C++ realiza el trabajo numérico principal: precalcula los coeficientes \(1/(2n^2)\), evalúa el producto de cosenos truncado en cada nodo de Simpson y acumula la suma ponderada correspondiente.
Los nodos interiores de Simpson se reparten entre varios hilos cuando eso es posible. Así se reduce el tiempo real de ejecución sin modificar el valor matemático. El nodo \(t=0\) se trata mediante el límite anterior y el número de subintervalos se ajusta para que la regla de Simpson siempre use un conteo par.
Las versiones en Python y Java son envoltorios ligeros que delegan en el mismo núcleo numérico compilado. Por eso las tres variantes comparten el mismo método, las mismas comprobaciones de simetría y la misma salida final.
Si \(T\) es el límite superior de integración, \(h\) el paso y \(N\) el corte del producto, entonces el número de subintervalos de Simpson es del orden de \(T/h\). Cada evaluación del integrando requiere multiplicar \(N\) factores coseno, así que el coste temporal total es
$$O\left(\frac{T}{h}\,N\right).$$
El consumo de memoria es \(O(N)\) para los coeficientes precalculados, más un pequeño almacenamiento adicional para las sumas parciales por hilo. El paralelismo mejora el factor constante, no el orden asintótico.
设随机级数为
$$Z=\sum_{n\ge 1}\frac{B_n}{n^2},\qquad B_n\in\{0,1\},\qquad \Pr(B_n=0)=\Pr(B_n=1)=\tfrac12.$$
题目要求高精度计算尾概率
$$p(a)=\Pr(Z\gt a)$$
在 \(a=0.5\) 时的数值。由于 \(Z\) 的取值始终落在 \(0\) 与 \(\zeta(2)=\pi^2/6\) 之间，而且它依赖无穷多个独立的二元随机变量，直接枚举所有情况根本不可行。实现采用的路线是：先把级数中心化，再写出特征函数，最后用 Fourier 反演把目标概率还原成一个一维积分。
整个推导可以分成几个自然步骤：确定均值与对称性，利用独立性得到无穷乘积形式，再通过反演公式把概率问题转成数值积分问题。
每一项 \(B_n/n^2\) 的期望是 \(1/(2n^2)\)，因此
$$\mu=\mathbb E[Z]=\sum_{n\ge 1}\frac{1}{2n^2}=\frac{\zeta(2)}{2}=\frac{\pi^2}{12}.$$
把随机位写成
$$B_n=\frac{1+\varepsilon_n}{2},\qquad \varepsilon_n\in\{-1,+1\},$$
其中 \(\varepsilon_n\) 取 \(\pm1\) 的概率相同，则有
$$Z=\mu+\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2}.$$
于是中心化变量为
$$Y=Z-\mu=\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2},$$
它关于 \(0\) 对称。另一方面，由于每一项都非负，而最大值出现在所有 \(B_n=1\) 的情形，所以
$$0\le Z\le \zeta(2).$$
这立刻给出边界结论
$$a\lt 0 \implies p(a)=1,\qquad a\ge \zeta(2)\implies p(a)=0.$$
对固定的 \(n\)，中心化后单项的特征函数是
$$\mathbb E\left[e^{it\varepsilon_n/(2n^2)}\right]=\frac{e^{it/(2n^2)}+e^{-it/(2n^2)}}{2}=\cos\left(\frac{t}{2n^2}\right).$$
因为各项独立，所以整个级数的特征函数就是这些因子的乘积：
$$\varphi_Y(t)=\prod_{n\ge 1}\cos\left(\frac{t}{2n^2}\right).$$
这个函数是实值且偶对称的，正好对应 \(Y\) 的对称分布。也正因为如此，后面的反演公式可以化成纯实数的正弦积分，不需要另外处理虚部。
令
$$x=\mu-a.$$
那么
$$p(a)=\Pr(Z\gt a)=\Pr(Y\gt -x).$$
应用 Gil-Pelaez 反演公式可得
$$p(a)=\frac12+\frac{1}{\pi}\int_0^{\infty}\frac{\sin(tx)}{t}\,\varphi_Y(t)\,dt.$$
这一步非常关键，因为它把“无穷随机和超过阈值的概率”变成了“一个明确的一维积分”。从这里开始，问题的核心就从概率论转向数值积分。
实际计算时，不可能保留无限多个余弦因子，于是用
$$\varphi_{Y,N}(t)=\prod_{n=1}^{N}\cos\left(\frac{t}{2n^2}\right)$$
去近似 \(\varphi_Y(t)\)。同时，积分区间 \([0,\infty)\) 也被替换成 \([0,T]\)。这种处理之所以有效，是因为当 \(n\) 很大时，\(\frac{t}{2n^2}\) 极小，被省略的余弦因子会非常接近 \(1\)；另一方面，积分核本身具有振荡性，超过足够大的上限后，新增区间对结果的净贡献会迅速减弱。
若取 \(T=mh\)，其中 \(m\) 为偶数、\(h\) 为步长，则采用
$$\int_0^T f(t)\,dt\approx \frac{h}{3}\left(f(0)+f(T)+4\sum_{j=1,3,\dots,m-1}f(jh)+2\sum_{j=2,4,\dots,m-2}f(jh)\right),$$
其中
$$f(t)=\frac{\sin(tx)}{t}\,\varphi_{Y,N}(t).$$
在 \(t=0\) 处看起来有 \(\frac{0}{0}\) 的形式，但实际上
$$\lim_{t\to 0}\frac{\sin(tx)}{t}=x,$$
因此原点处是可去奇点，可以稳定地单独处理。这一点也是实现细节中必须特别注意的地方。
如果把每个随机位 \(B_n\) 都替换成 \(1-B_n\)，那么新的和变成
$$\sum_{n\ge 1}\frac{1-B_n}{n^2}=\zeta(2)-Z.$$
这说明 \(Z\) 的分布关于 \(\mu=\zeta(2)/2\) 对称，所以一定有
$$p(\mu)=\Pr(Z\gt \mu)=\frac12,$$
并且更一般地满足
$$p(\mu-u)+p(\mu+u)=1.$$
特别地，阈值 \(1\) 与 \(\zeta(2)-1\) 正好以 \(\mu\) 为中心对称，所以它们对应的尾概率都应等于 \(1/2\)。实现正是利用这些关系做数值校验。采用较细的积分参数时，目标值计算为
$$p(0.5)\approx 0.565654540708545,$$
这也是最终输出在四舍五入前所对应的高精度数值。
C++、Python 和 Java 三种实现最终都在计算同一个截断后的 Fourier 积分。真正的数值工作由 C++ 实现完成：先预计算直到截断点为止的 \(1/(2n^2)\) 系数，再在每个 Simpson 采样点上计算截断余弦乘积，最后把所有节点按 Simpson 权重累加起来。
为了加速，内部采样点会按区间划分给多个线程处理；端点则单独处理，其中 \(t=0\) 使用上面的极限值。实现还会自动保证 Simpson 公式使用偶数个子区间，这是该求积公式成立的基本条件。
Python 和 Java 版本只是轻量入口，它们把最终求值委托给同一个已编译的数值核心。因此三种语言在数学方法、截断策略、对称性校验以及最终概率上完全一致，只要使用相同的积分参数，就应得到相同结果。
设积分上限为 \(T\)，步长为 \(h\)，余弦乘积的截断长度为 \(N\)。则 Simpson 子区间数大约是 \(T/h\)。每个采样点都需要乘上 \(N\) 个余弦因子，所以总时间复杂度为
$$O\left(\frac{T}{h}\,N\right).$$
空间复杂度主要来自预计算的系数数组，因此是 \(O(N)\)。除此之外，只需要少量与线程数相关的部分和缓存。并行化只能改善运行常数，不改变渐近阶。
Рассматривается случайный ряд
$$Z=\sum_{n\ge 1}\frac{B_n}{n^2},\qquad B_n\in\{0,1\},\qquad \Pr(B_n=0)=\Pr(B_n=1)=\tfrac12.$$
Нужно найти хвостовую вероятность
$$p(a)=\Pr(Z\gt a)$$
для \(a=0.5\) с высокой точностью. Поскольку \(0\le Z\le \zeta(2)=\pi^2/6\), а ряд зависит от бесконечного числа независимых двоичных выборов, прямой перебор невозможен. Поэтому решение переходит к центрированной форме ряда, выписывает его характеристическую функцию и затем восстанавливает вероятность через обращение Фурье.
Основная схема такова: центрировать сумму, использовать независимость для получения бесконечного произведения, а затем свести задачу к вычислению одного осциллирующего интеграла.
Средний вклад одного слагаемого равен \(1/(2n^2)\), поэтому
$$\mu=\mathbb E[Z]=\sum_{n\ge 1}\frac{1}{2n^2}=\frac{\zeta(2)}{2}=\frac{\pi^2}{12}.$$
Запишем \(B_n=\frac{1+\varepsilon_n}{2}\), где \(\varepsilon_n\in\{-1,+1\}\) равновероятны. Тогда
$$Z=\mu+\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2}.$$
Следовательно, центрированная величина
$$Y=Z-\mu=\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2}$$
симметрична относительно нуля. Кроме того, все слагаемые неотрицательны, а максимальная сумма равна \(\zeta(2)\), так что
$$a\lt 0 \implies p(a)=1,\qquad a\ge \zeta(2)\implies p(a)=0.$$
Для фиксированного \(n\) имеем
$$\mathbb E\left[e^{it\varepsilon_n/(2n^2)}\right]=\frac{e^{it/(2n^2)}+e^{-it/(2n^2)}}{2}=\cos\left(\frac{t}{2n^2}\right).$$
Независимость слагаемых дает произведение:
$$\varphi_Y(t)=\prod_{n\ge 1}\cos\left(\frac{t}{2n^2}\right).$$
Эта функция вещественная и четная, что полностью согласуется с симметрией распределения \(Y\).
Положим
$$x=\mu-a.$$
Тогда \(p(a)=\Pr(Y\gt -x)\). По формуле Гила-Пелаэса получаем
$$p(a)=\frac12+\frac{1}{\pi}\int_0^{\infty}\frac{\sin(tx)}{t}\,\varphi_Y(t)\,dt.$$
Тем самым исходная вероятностная задача для бесконечного ряда превращается в задачу численного вычисления одномерного интеграла.
На практике используется приближение
$$\varphi_{Y,N}(t)=\prod_{n=1}^{N}\cos\left(\frac{t}{2n^2}\right).$$
Интервал интегрирования \([0,\infty)\) тоже заменяется на \([0,T]\). Это оправдано тем, что при больших \(n\) пропущенные косинусные множители почти равны \(1\), а колебательный интеграл за достаточно большим \(T\) дает все меньший чистый вклад.
Если \(T=mh\), где \(m\) четно, то применяется
$$\int_0^T f(t)\,dt\approx \frac{h}{3}\left(f(0)+f(T)+4\sum_{j=1,3,\dots,m-1}f(jh)+2\sum_{j=2,4,\dots,m-2}f(jh)\right),$$
где
$$f(t)=\frac{\sin(tx)}{t}\,\varphi_{Y,N}(t).$$
В точке \(t=0\) используется предел
$$\lim_{t\to 0}\frac{\sin(tx)}{t}=x,$$
поэтому видимая особенность в нуле является устранимой.
Если заменить каждый бит на \(1-B_n\), то
$$\sum_{n\ge 1}\frac{1-B_n}{n^2}=\zeta(2)-Z.$$
Значит, распределение \(Z\) симметрично относительно \(\mu=\zeta(2)/2\). Отсюда следует
$$p(\mu)=\Pr(Z\gt \mu)=\frac12,$$
а также более общее равенство
$$p(\mu-u)+p(\mu+u)=1.$$
В частности, точки \(1\) и \(\zeta(2)-1\) симметричны относительно \(\mu\), поэтому обе дают вероятность \(1/2\). При более точных параметрах интегрирования основное значение получается равным
$$p(0.5)\approx 0.565654540708545.$$
Реализации на C++, Python и Java вычисляют один и тот же усеченный Fourier-интеграл. Основная численная работа выполняется в C++: заранее вычисляются коэффициенты \(1/(2n^2)\) до выбранного усечения, затем в каждом узле схемы Симпсона считается усеченное произведение косинусов, после чего узлы суммируются с нужными весами.
Внутренние узлы сетки распределяются между доступными потоками. Это ускоряет выполнение на многоядерных системах, но не меняет математический результат. Точка \(t=0\) обрабатывается отдельно через предельное значение, а число подинтервалов при необходимости доводится до четного, как и требуется правилом Симпсона.
Версии на Python и Java являются легкими оболочками и передают вычисление тому же скомпилированному численному ядру. Поэтому все три реализации используют одинаковые проверки симметрии, одинаковые параметры усечения и выдают одно и то же значение при совпадающих настройках.
Пусть \(T\) - верхняя граница интегрирования, \(h\) - шаг сетки, а \(N\) - длина усеченного произведения. Тогда число подинтервалов порядка \(T/h\). В каждом узле нужно перемножить \(N\) косинусных множителей, поэтому суммарная временная сложность равна
$$O\left(\frac{T}{h}\,N\right).$$
По памяти требуется \(O(N)\) для предвычисленных коэффициентов и небольшой дополнительный буфер для частичных сумм по потокам. Распараллеливание улучшает только константу.
لدينا السلسلة العشوائية
$$Z=\sum_{n\ge 1}\frac{B_n}{n^2},\qquad B_n\in\{0,1\},\qquad \Pr(B_n=0)=\Pr(B_n=1)=\tfrac12.$$
والمطلوب هو حساب احتمال الذيل
$$p(a)=\Pr(Z\gt a)$$
عند \(a=0.5\) بدقة عالية. بما أن \(Z\) يبقى دائمًا بين \(0\) و\(\zeta(2)=\pi^2/6\)، وبما أن السلسلة تعتمد على عدد لا نهائي من الاختيارات الثنائية المستقلة، فإن التعداد المباشر غير عملي تمامًا. لذلك تعتمد الطريقة على تمركز السلسلة حول متوسطها، ثم اشتقاق الدالة المميزة، ثم استرجاع الاحتمال بواسطة عكس Fourier.
الفكرة الأساسية تمر عبر مراحل واضحة: نحسب المتوسط، نعيد كتابة السلسلة بصيغة متناظرة، نبني الدالة المميزة من الاستقلال، ثم نحوّل الاحتمال المطلوب إلى تكامل واحد يمكن تقريبه عدديًا.
القيمة المتوقعة لكل حد هي \(1/(2n^2)\)، ومن ثم
$$\mu=\mathbb E[Z]=\sum_{n\ge 1}\frac{1}{2n^2}=\frac{\zeta(2)}{2}=\frac{\pi^2}{12}.$$
إذا كتبنا
$$B_n=\frac{1+\varepsilon_n}{2},\qquad \varepsilon_n\in\{-1,+1\},$$
حيث تأخذ \(\varepsilon_n\) القيمتين \(\pm1\) باحتمالين متساويين، فإن
$$Z=\mu+\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2}.$$
إذن المتغير الممركز هو
$$Y=Z-\mu=\sum_{n\ge 1}\frac{\varepsilon_n}{2n^2},$$
وهو متناظر حول الصفر. كذلك بما أن كل الحدود غير سالبة وأكبر مجموع ممكن يساوي \(\zeta(2)\)، نحصل مباشرة على
$$a\lt 0 \implies p(a)=1,\qquad a\ge \zeta(2)\implies p(a)=0.$$
لكل \(n\) ثابت، تكون الدالة المميزة للحد الممركز
$$\mathbb E\left[e^{it\varepsilon_n/(2n^2)}\right]=\frac{e^{it/(2n^2)}+e^{-it/(2n^2)}}{2}=\cos\left(\frac{t}{2n^2}\right).$$
وبسبب استقلال الحدود يتفكك الناتج إلى حاصل ضرب لا نهائي:
$$\varphi_Y(t)=\prod_{n\ge 1}\cos\left(\frac{t}{2n^2}\right).$$
هذه الدالة حقيقية وزوجية، وهذا يعبّر مباشرة عن تناظر توزيع \(Y\).
نضع
$$x=\mu-a.$$
عندئذ يصبح
$$p(a)=\Pr(Y\gt -x).$$
وباستخدام صيغة Gil-Pelaez نحصل على
$$p(a)=\frac12+\frac{1}{\pi}\int_0^{\infty}\frac{\sin(tx)}{t}\,\varphi_Y(t)\,dt.$$
وهكذا تحوّل السؤال الاحتمالي الأصلي إلى تكامل حقيقي أحادي البعد يمكن التعامل معه عدديًا.
في التطبيق العملي لا يمكن الاحتفاظ بعدد لا نهائي من العوامل، لذلك يُستبدل الحاصل بـ
$$\varphi_{Y,N}(t)=\prod_{n=1}^{N}\cos\left(\frac{t}{2n^2}\right).$$
كما يُستبدل المجال \([0,\infty)\) بالمجال \([0,T]\). هذا التقريب معقول لأن عوامل جيب التمام المحذوفة تصبح قريبة جدًا من \(1\) عندما يكبر \(n\)، ولأن التكامل الاهتزازي يضيف مساهمة صافية أصغر فأصغر بعد حد علوي كبير بما يكفي.
إذا كان \(T=mh\) مع \(m\) زوجيًا، فالتقريب العددي هو
$$\int_0^T f(t)\,dt\approx \frac{h}{3}\left(f(0)+f(T)+4\sum_{j=1,3,\dots,m-1}f(jh)+2\sum_{j=2,4,\dots,m-2}f(jh)\right),$$
حيث
$$f(t)=\frac{\sin(tx)}{t}\,\varphi_{Y,N}(t).$$
وعند \(t=0\) لا توجد مشكلة حقيقية، لأن
$$\lim_{t\to 0}\frac{\sin(tx)}{t}=x.$$
لذلك يمكن التعامل مع نقطة البداية بطريقة مستقرة بدل اعتبارها مفردة حقيقية.
إذا استبدلنا كل بت \(B_n\) بـ \(1-B_n\)، فإن المجموع الجديد يساوي
$$\sum_{n\ge 1}\frac{1-B_n}{n^2}=\zeta(2)-Z.$$
وهذا يعني أن توزيع \(Z\) متناظر حول \(\mu=\zeta(2)/2\). لذلك
$$p(\mu)=\Pr(Z\gt \mu)=\frac12,$$
وبشكل أعم
$$p(\mu-u)+p(\mu+u)=1.$$
وبصورة خاصة، فإن العتبتين \(1\) و\(\zeta(2)-1\) متناظرتان حول \(\mu\)، ولذلك يعطي كل منهما احتمالًا يساوي \(1/2\). ومع الإعدادات العددية الأدق، تعطي الخوارزمية القيمة
$$p(0.5)\approx 0.565654540708545.$$
التنفيذات في C++ وPython وJava تحسب جميعها التكامل نفسه بعد قطع الحاصل اللانهائي. التنفيذ في C++ هو الذي ينفذ العمل العددي الفعلي: يهيئ مسبقًا معاملات \(1/(2n^2)\)، ثم يحسب حاصل ضرب جيوب التمام المقطوع عند كل عقدة من عقد Simpson، ثم يجمع القيم بالأوزان المناسبة.
تُوزَّع العقد الداخلية للتكامل على عدة خيوط عندما تكون المعالجة المتوازية متاحة. هذا يحسن زمن التنفيذ على الأجهزة متعددة الأنوية من دون أن يغيّر الناتج الرياضي. أما النقطة \(t=0\) فتعالج باستخدام النهاية السابقة، ويُعدَّل عدد الفواصل ليبقى زوجيًا كما تتطلب قاعدة Simpson.
أما نسختا Python وJava فهما مجرد واجهتين خفيفتين تفوضان الحساب إلى النواة العددية المترجمة نفسها. لذلك تشترك اللغات الثلاث في الطريقة الرياضية نفسها، وفي اختبارات التناظر نفسها، وفي الناتج النهائي نفسه عندما تتطابق الإعدادات العددية.
إذا كان حد التكامل الأعلى هو \(T\)، وكان طول الخطوة \(h\)، وكان طول قطع حاصل الضرب هو \(N\)، فإن عدد فواصل Simpson يكون تقريبًا \(T/h\). وكل عقدة تتطلب ضرب \(N\) من عوامل جيب التمام، لذا فإن الزمن الكلي هو
$$O\left(\frac{T}{h}\,N\right).$$
أما الذاكرة فتبلغ \(O(N)\) لتخزين المعاملات المحسوبة مسبقًا، مع زيادة صغيرة فقط لمجاميع جزئية مرتبطة بعدد الخيوط. التوازي يحسن العامل الثابت ولا يغير الرتبة التقاربية.