The local solutions model the pie-eating process through an expectation function \(E(x)\). Here \(x \ge 1\) is the current scaled state, one step is spent immediately, and then the state is multiplied by a random factor \(W \in [0,1]\) whose density is \(2(1-w)\). The process stops as soon as the state falls below 1, so the boundary condition is \(E(u)=0\) for \(u \lt 1\).
Therefore the programming task is not to simulate random trials, but to derive and evaluate the exact formula for \(E(x)\), especially at the default input \(x=40\).
All three implementations use the same closed form:
$$E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3},\qquad x\ge1.$$
Condition on the first random multiplier \(W\). We always spend one step now, and after that the expected remaining work is \(E(xW)\). Averaging over the density \(2(1-w)\) gives the renewal equation used in the C++ comments:
$$E(x)=1+2\int_0^1 (1-w)\,E(xw)\,dw.$$
Because \(E(xw)=0\) whenever \(xw \lt 1\), only the interval \(w \ge 1/x\) contributes when \(x \ge 1\). So the same equation can be rewritten as
$$E(x)=1+2\int_{1/x}^1 (1-w)\,E(xw)\,dw.$$
Set \(t=xw\). Then \(dw=dt/x\) and \(w=t/x\), which converts the integral into a Volterra-type expression with fixed lower limit:
$$E(x)=1+\frac{2}{x^2}\int_1^x (x-t)\,E(t)\,dt.$$
This form is much easier to differentiate because the unknown function now appears as \(E(t)\) instead of \(E(xw)\).
Define
$$I(x)=\int_1^x (x-t)\,E(t)\,dt.$$
The integral equation becomes
$$x^2\bigl(E(x)-1\bigr)=2I(x).$$
Differentiate once:
$$2x\bigl(E(x)-1\bigr)+x^2E'(x)=2\int_1^x E(t)\,dt.$$
Differentiate a second time:
$$2\bigl(E(x)-1\bigr)+4xE'(x)+x^2E''(x)=2E(x).$$
The \(E(x)\) terms cancel, leaving the differential equation
$$x^2E''(x)+4xE'(x)=2.$$
Let \(Y(x)=E'(x)\). Then
$$x^2Y'(x)+4xY(x)=2.$$
After division by \(x^2\), this is a first-order linear equation:
$$Y'(x)+\frac{4}{x}Y(x)=\frac{2}{x^2}.$$
The integrating factor is \(x^4\), so
$$\frac{d}{dx}\left(x^4Y(x)\right)=2x^2.$$
Integrating gives
$$x^4Y(x)=\frac{2}{3}x^3+C_1,$$
hence
$$E'(x)=\frac{2}{3x}+\frac{C_1}{x^4}.$$
Integrating once more,
$$E(x)=\frac{2}{3}\ln(x)-\frac{C_1}{3x^3}+C_2.$$
The logarithm is the dominant term, while the homogeneous part contributes the \(x^{-3}\) correction seen in the code.
At \(x=1\), the integral part vanishes, so
$$E(1)=1.$$
Now evaluate the once-differentiated identity at \(x=1\):
$$2\bigl(E(1)-1\bigr)+E'(1)=0.$$
Since \(E(1)=1\), this yields
$$E'(1)=0.$$
Substitute into the general solution:
$$0=E'(1)=\frac{2}{3}+C_1,$$
so \(C_1=-\frac{2}{3}\). Then
$$1=E(1)=0+\frac{2}{9}+C_2,$$
which gives \(C_2=\frac{7}{9}\). Therefore
$$\boxed{E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3}.}$$
The C++ program verifies the formula at three checkpoints, and they match exactly:
$$E(1)=1.0000000000,$$
$$E(2)=1.2676536759,$$
$$E(7.5)=2.1215732071.$$
For the default Project Euler input, the same formula gives
$$E(40)=3.2370342194.$$
As \(x\) grows, the correction term decays rapidly, so the expectation behaves like
$$E(x)=\frac{2}{3}\ln(x)+\frac{7}{9}+O(x^{-3}).$$
The C++ solution stores the mathematics in expected_steps(x) and evaluates the closed form with long double. It also supports an optional --x=... argument and a --skip-checkpoints switch, while the helper nearly_equal validates the known values \(E(1)\), \(E(2)\), and \(E(7.5)\).
The Python and Java versions intentionally strip the implementation down to the essential formula. Each computes the same expression and formats the value at \(x=40\) with 10 digits after the decimal point. None of the three languages performs simulation or numerical integration at runtime.
Once the closed form is known, evaluation needs one logarithm and a constant number of arithmetic operations. Therefore the running time is \(O(1)\) and the memory usage is \(O(1)\). This is far better than Monte Carlo sampling or iterative quadrature, both of which would need many repeated evaluations to approach the same precision.
Die lokalen Lösungen modellieren den Pie-Prozess durch eine Erwartungswertfunktion \(E(x)\). Dabei beschreibt \(x \ge 1\) den aktuellen skalierten Zustand, ein Schritt wird sofort verbraucht, und danach wird der Zustand mit einem Zufallsfaktor \(W \in [0,1]\) multipliziert, dessen Dichte \(2(1-w)\) ist. Sobald der Zustand unter 1 fällt, endet der Prozess, also gilt als Randbedingung \(E(u)=0\) für \(u \lt 1\).
Die eigentliche Programmieraufgabe besteht daher nicht in einer Simulation, sondern in der exakten Herleitung und Auswertung von \(E(x)\), insbesondere für den Standardwert \(x=40\).
Alle drei Implementierungen verwenden dieselbe geschlossene Formel:
$$E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3},\qquad x\ge1.$$
Konditioniert man auf den ersten Zufallsfaktor \(W\), dann wird sofort ein Schritt gezahlt und der erwartete Restaufwand ist \(E(xW)\). Das Mitteln gegen die Dichte \(2(1-w)\) liefert die in den C++-Kommentaren genannte Erneuerungsgleichung:
$$E(x)=1+2\int_0^1 (1-w)\,E(xw)\,dw.$$
Da \(E(xw)=0\) gilt, sobald \(xw \lt 1\), trägt für \(x \ge 1\) nur der Bereich \(w \ge 1/x\) bei. Also auch
$$E(x)=1+2\int_{1/x}^1 (1-w)\,E(xw)\,dw.$$
Mit der Substitution \(t=xw\) gilt \(dw=dt/x\) und \(w=t/x\). Dadurch wird die Gleichung zu einer Volterra-artigen Integralform mit fester Untergrenze:
$$E(x)=1+\frac{2}{x^2}\int_1^x (x-t)\,E(t)\,dt.$$
Diese Form ist deutlich günstiger, weil die unbekannte Funktion nun als \(E(t)\) statt als \(E(xw)\) vorkommt.
Definiere
$$I(x)=\int_1^x (x-t)\,E(t)\,dt.$$
Dann lautet die Integralgleichung
$$x^2\bigl(E(x)-1\bigr)=2I(x).$$
Einmal ableiten ergibt
$$2x\bigl(E(x)-1\bigr)+x^2E'(x)=2\int_1^x E(t)\,dt.$$
Noch einmal ableiten liefert
$$2\bigl(E(x)-1\bigr)+4xE'(x)+x^2E''(x)=2E(x).$$
Die \(E(x)\)-Terme heben sich weg, und es bleibt
$$x^2E''(x)+4xE'(x)=2.$$
Setze \(Y(x)=E'(x)\). Dann gilt
$$x^2Y'(x)+4xY(x)=2.$$
Nach Division durch \(x^2\) erhält man die lineare Gleichung erster Ordnung
$$Y'(x)+\frac{4}{x}Y(x)=\frac{2}{x^2}.$$
Der Integrationsfaktor ist \(x^4\), also
$$\frac{d}{dx}\left(x^4Y(x)\right)=2x^2.$$
Integration ergibt
$$x^4Y(x)=\frac{2}{3}x^3+C_1,$$
also
$$E'(x)=\frac{2}{3x}+\frac{C_1}{x^4}.$$
Noch einmal integrieren:
$$E(x)=\frac{2}{3}\ln(x)-\frac{C_1}{3x^3}+C_2.$$
Der logarithmische Anteil dominiert, während die homogene Lösung genau den im Code sichtbaren \(x^{-3}\)-Korrekturterm liefert.
Für \(x=1\) verschwindet der Integralteil, also
$$E(1)=1.$$
Setzt man \(x=1\) in die einmal abgeleitete Gleichung ein, erhält man
$$2\bigl(E(1)-1\bigr)+E'(1)=0.$$
Wegen \(E(1)=1\) folgt daraus
$$E'(1)=0.$$
Nun in die allgemeine Lösung einsetzen:
$$0=E'(1)=\frac{2}{3}+C_1,$$
also \(C_1=-\frac{2}{3}\). Außerdem
$$1=E(1)=0+\frac{2}{9}+C_2,$$
daher \(C_2=\frac{7}{9}\). Somit ergibt sich
$$\boxed{E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3}.}$$
Die im C++-Programm verwendeten Kontrollwerte stimmen exakt:
$$E(1)=1.0000000000,$$
$$E(2)=1.2676536759,$$
$$E(7.5)=2.1215732071.$$
Für den Standardwert des Project-Euler-Problems ergibt sich
$$E(40)=3.2370342194.$$
Für große \(x\) bleibt im Wesentlichen
$$E(x)=\frac{2}{3}\ln(x)+\frac{7}{9}+O(x^{-3}).$$
Die C++-Datei kapselt die Mathematik in expected_steps(x) und wertet die geschlossene Formel mit long double aus. Zusätzlich gibt es ein optionales Argument --x=..., einen Schalter --skip-checkpoints und die Hilfsfunktion nearly_equal, um die bekannten Werte \(E(1)\), \(E(2)\) und \(E(7.5)\) zu prüfen.
Die Python- und Java-Versionen sind bewusst kompakter. Beide berechnen exakt denselben Ausdruck und formatieren den Wert bei \(x=40\) mit 10 Nachkommastellen. Keine der Implementierungen führt zur Laufzeit eine Simulation oder numerische Integration aus.
Ist die geschlossene Formel einmal hergeleitet, braucht die Auswertung nur einen Logarithmus und konstant viele arithmetische Operationen. Damit beträgt die Laufzeit \(O(1)\) und der Speicherverbrauch \(O(1)\). Gegenüber Monte-Carlo-Simulation oder iterativer Quadratur ist das sowohl schneller als auch genauer.
Yerel çözümler pie sürecini \(E(x)\) adlı bir beklenen değer fonksiyonu ile modelliyor. Burada \(x \ge 1\) mevcut ölçeklenmiş durumu gösterir; bir adım hemen harcanır ve ardından durum, yoğunluğu \(2(1-w)\) olan \(W \in [0,1]\) rastgele çarpanı ile çarpılır. Durum 1'in altına düştüğünde süreç biter; dolayısıyla sınır koşulu \(u \lt 1\) için \(E(u)=0\) olur.
Bu nedenle programın işi rastgele denemeler simüle etmek değil, \(E(x)\) için kapalı formu türetip özellikle varsayılan \(x=40\) değerinde hesaplamaktır.
Üç dildeki çözüm de aynı formülü kullanır:
$$E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3},\qquad x\ge1.$$
İlk rastgele çarpan \(W\)'ye koşullayalım. Şu andaki bir adımı kesin harcarız; kalan beklenen iş yükü ise \(E(xW)\) olur. Bunu yoğunluk \(2(1-w)\) ile ortaladığımızda C++ yorumunda verilen yenilenme denklemi elde edilir:
$$E(x)=1+2\int_0^1 (1-w)\,E(xw)\,dw.$$
\(xw \lt 1\) olduğunda \(E(xw)=0\) olduğu için, \(x \ge 1\) durumunda yalnızca \(w \ge 1/x\) aralığı katkı verir. Bu yüzden eşdeğer biçim
$$E(x)=1+2\int_{1/x}^1 (1-w)\,E(xw)\,dw.$$
\(t=xw\) değişken dönüşümünü yapalım. O zaman \(dw=dt/x\) ve \(w=t/x\) olur. Denklem sabit alt sınırlı, Volterra tipinde bir integrale dönüşür:
$$E(x)=1+\frac{2}{x^2}\int_1^x (x-t)\,E(t)\,dt.$$
Bu biçim daha kullanışlıdır; çünkü bilinmeyen fonksiyon artık \(E(xw)\) yerine \(E(t)\) olarak görünür.
Şunu tanımlayalım:
$$I(x)=\int_1^x (x-t)\,E(t)\,dt.$$
Bu durumda integral denklem
$$x^2\bigl(E(x)-1\bigr)=2I(x)$$
haline gelir. Bir kez türev alırsak
$$2x\bigl(E(x)-1\bigr)+x^2E'(x)=2\int_1^x E(t)\,dt$$
elde edilir. Bir kez daha türev alınca
$$2\bigl(E(x)-1\bigr)+4xE'(x)+x^2E''(x)=2E(x)$$
olur. \(E(x)\) terimleri sadeleşir ve geriye
$$x^2E''(x)+4xE'(x)=2$$
kalır.
\(Y(x)=E'(x)\) yazalım. Böylece
$$x^2Y'(x)+4xY(x)=2$$
olur. \(x^2\)'ye böldüğümüzde birinci mertebeden lineer denklem elde ederiz:
$$Y'(x)+\frac{4}{x}Y(x)=\frac{2}{x^2}.$$
İntegrasyon çarpanı \(x^4\) olduğundan
$$\frac{d}{dx}\left(x^4Y(x)\right)=2x^2$$
yazılır. İntegral alınca
$$x^4Y(x)=\frac{2}{3}x^3+C_1,$$
dolayısıyla
$$E'(x)=\frac{2}{3x}+\frac{C_1}{x^4}.$$
Bir kez daha integre edersek
$$E(x)=\frac{2}{3}\ln(x)-\frac{C_1}{3x^3}+C_2$$
bulunur. Kodda görülen \(x^{-3}\) düzeltme terimi tam olarak homojen çözümden gelir.
\(x=1\) için integral kısmı yok olur, yani
$$E(1)=1.$$
Tek türev alınmış bağıntıda \(x=1\) yazarsak
$$2\bigl(E(1)-1\bigr)+E'(1)=0$$
elde edilir. \(E(1)=1\) olduğundan
$$E'(1)=0$$
çıkar. Şimdi genel çözüme yerleştirelim:
$$0=E'(1)=\frac{2}{3}+C_1,$$
buradan \(C_1=-\frac{2}{3}\). Ayrıca
$$1=E(1)=0+\frac{2}{9}+C_2,$$
olduğu için \(C_2=\frac{7}{9}\). Sonuç olarak
$$\boxed{E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3}.}$$
C++ dosyasındaki kontrol değerleri bu formülle birebir uyuşur:
$$E(1)=1.0000000000,$$
$$E(2)=1.2676536759,$$
$$E(7.5)=2.1215732071.$$
Varsayılan Project Euler girdisi için de
$$E(40)=3.2370342194$$
elde edilir. Büyük \(x\) için davranış
$$E(x)=\frac{2}{3}\ln(x)+\frac{7}{9}+O(x^{-3})$$
şeklindedir.
C++ çözümü tüm matematiği expected_steps(x) fonksiyonunda toplar ve kapalı formu long double ile hesaplar. Ayrıca isteğe bağlı --x=... argümanı, --skip-checkpoints seçeneği ve \(E(1)\), \(E(2)\), \(E(7.5)\) değerlerini doğrulayan nearly_equal yardımcı fonksiyonu vardır.
Python ve Java sürümleri bilinçli olarak daha kısadır. İkisi de aynı ifadeyi hesaplar ve \(x=40\) sonucunu 10 ondalık basamakla biçimlendirir. Çalışma zamanında hiçbir dilde simülasyon ya da sayısal integral yapılmaz.
Kapalı form elde edildikten sonra hesaplama yalnızca bir logaritma ve sabit sayıda aritmetik işlem gerektirir. Bu yüzden zaman karmaşıklığı \(O(1)\), bellek karmaşıklığı \(O(1)\) olur. Monte Carlo ya da yinelemeli nümerik integrasyonla karşılaştırıldığında hem daha hızlı hem de daha doğrudur.
Las soluciones locales describen el proceso de comer la tarta mediante una función de esperanza \(E(x)\). Aquí \(x \ge 1\) representa el estado escalado actual, se consume un paso de inmediato y luego el estado se multiplica por un factor aleatorio \(W \in [0,1]\) con densidad \(2(1-w)\). El proceso termina en cuanto el estado cae por debajo de 1, así que la condición de frontera es \(E(u)=0\) para \(u \lt 1\).
Por tanto, la tarea no consiste en simular muchas trayectorias aleatorias, sino en derivar y evaluar la expresión exacta de \(E(x)\), sobre todo para el valor por defecto \(x=40\).
Las tres implementaciones usan la misma forma cerrada:
$$E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3},\qquad x\ge1.$$
Condicionamos respecto al primer multiplicador aleatorio \(W\). Siempre gastamos un paso ahora y después el trabajo esperado restante es \(E(xW)\). Promediando con la densidad \(2(1-w)\) se obtiene la ecuación de renovación citada en el comentario de C++:
$$E(x)=1+2\int_0^1 (1-w)\,E(xw)\,dw.$$
Como \(E(xw)=0\) cuando \(xw \lt 1\), para \(x \ge 1\) solo contribuye el intervalo \(w \ge 1/x\). Así, también podemos escribir
$$E(x)=1+2\int_{1/x}^1 (1-w)\,E(xw)\,dw.$$
Tomamos el cambio de variable \(t=xw\). Entonces \(dw=dt/x\) y \(w=t/x\). La integral pasa a una forma de tipo Volterra con límite inferior fijo:
$$E(x)=1+\frac{2}{x^2}\int_1^x (x-t)\,E(t)\,dt.$$
Esta versión es mucho más cómoda porque la incógnita aparece como \(E(t)\) y no como \(E(xw)\).
Definimos
$$I(x)=\int_1^x (x-t)\,E(t)\,dt.$$
Entonces la ecuación integral se reescribe como
$$x^2\bigl(E(x)-1\bigr)=2I(x).$$
Derivando una vez:
$$2x\bigl(E(x)-1\bigr)+x^2E'(x)=2\int_1^x E(t)\,dt.$$
Derivando una segunda vez:
$$2\bigl(E(x)-1\bigr)+4xE'(x)+x^2E''(x)=2E(x).$$
Los términos con \(E(x)\) se cancelan y queda
$$x^2E''(x)+4xE'(x)=2.$$
Sea \(Y(x)=E'(x)\). Entonces
$$x^2Y'(x)+4xY(x)=2.$$
Dividiendo por \(x^2\) obtenemos una ecuación lineal de primer orden:
$$Y'(x)+\frac{4}{x}Y(x)=\frac{2}{x^2}.$$
El factor integrante es \(x^4\), de modo que
$$\frac{d}{dx}\left(x^4Y(x)\right)=2x^2.$$
Integrando:
$$x^4Y(x)=\frac{2}{3}x^3+C_1,$$
y por tanto
$$E'(x)=\frac{2}{3x}+\frac{C_1}{x^4}.$$
Una integración adicional produce
$$E(x)=\frac{2}{3}\ln(x)-\frac{C_1}{3x^3}+C_2.$$
La parte principal es logarítmica, y la solución homogénea explica el término correctivo \(x^{-3}\).
En \(x=1\) la parte integral desaparece, así que
$$E(1)=1.$$
Si evaluamos en \(x=1\) la identidad derivada una vez, obtenemos
$$2\bigl(E(1)-1\bigr)+E'(1)=0.$$
Como \(E(1)=1\), resulta
$$E'(1)=0.$$
Sustituyendo en la solución general:
$$0=E'(1)=\frac{2}{3}+C_1,$$
luego \(C_1=-\frac{2}{3}\). Además,
$$1=E(1)=0+\frac{2}{9}+C_2,$$
de donde \(C_2=\frac{7}{9}\). En consecuencia,
$$\boxed{E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3}.}$$
Los puntos de control del archivo C++ se recuperan exactamente:
$$E(1)=1.0000000000,$$
$$E(2)=1.2676536759,$$
$$E(7.5)=2.1215732071.$$
Para la entrada por defecto del problema,
$$E(40)=3.2370342194.$$
Además, para \(x\) grande se tiene la expansión
$$E(x)=\frac{2}{3}\ln(x)+\frac{7}{9}+O(x^{-3}).$$
La versión en C++ concentra toda la matemática en expected_steps(x) y evalúa la fórmula con long double. También admite el argumento opcional --x=..., el interruptor --skip-checkpoints y una función auxiliar nearly_equal para comprobar los valores conocidos \(E(1)\), \(E(2)\) y \(E(7.5)\).
Las versiones en Python y Java son deliberadamente mínimas. Ambas calculan exactamente la misma expresión y formatean el resultado para \(x=40\) con 10 decimales. Ninguna implementación hace integración numérica ni simulación estocástica durante la ejecución.
Una vez conocida la forma cerrada, la evaluación requiere un logaritmo y un número constante de operaciones aritméticas. Por eso el coste temporal es \(O(1)\) y el coste espacial también es \(O(1)\). Esto supera con claridad a un enfoque de Monte Carlo o a una cuadratura iterativa, que necesitarían muchas más operaciones para alcanzar precisión comparable.
本仓库中的解法把“Eating Pie”过程写成一个期望函数 \(E(x)\)。其中 \(x \ge 1\) 表示当前的缩放状态；先必定消耗一步，然后状态再乘上一个随机因子 \(W \in [0,1]\)，其密度为 \(2(1-w)\)。一旦状态降到 1 以下，过程立即结束，因此边界条件是 \(E(u)=0\)（当 \(u \lt 1\) 时）。
所以程序并不是做随机模拟，而是先把 \(E(x)\) 的精确公式推出来，再直接代入，默认求的是 \(x=40\)。
C++、Python 和 Java 三个版本使用完全相同的闭式：
$$E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3},\qquad x\ge1.$$
对第一次随机乘子 \(W\) 分类讨论。当前这一步一定会发生，之后剩余的期望步数是 \(E(xW)\)。按密度 \(2(1-w)\) 取平均，就得到代码注释中的更新方程：
$$E(x)=1+2\int_0^1 (1-w)\,E(xw)\,dw.$$
由于当 \(xw \lt 1\) 时有 \(E(xw)=0\)，所以对 \(x \ge 1\) 来说，真正有贡献的只是 \(w \ge 1/x\) 这一段。于是可改写为
$$E(x)=1+2\int_{1/x}^1 (1-w)\,E(xw)\,dw.$$
令 \(t=xw\)，则 \(dw=dt/x\)，而 \(w=t/x\)。这样积分变成具有固定下限的 Volterra 型形式：
$$E(x)=1+\frac{2}{x^2}\int_1^x (x-t)\,E(t)\,dt.$$
这个形式更容易处理，因为未知函数已经从 \(E(xw)\) 变成了 \(E(t)\)。
定义
$$I(x)=\int_1^x (x-t)\,E(t)\,dt.$$
则原式变为
$$x^2\bigl(E(x)-1\bigr)=2I(x).$$
先对 \(x\) 求一次导数：
$$2x\bigl(E(x)-1\bigr)+x^2E'(x)=2\int_1^x E(t)\,dt.$$
再求一次导数：
$$2\bigl(E(x)-1\bigr)+4xE'(x)+x^2E''(x)=2E(x).$$
两边的 \(E(x)\) 项消掉后，得到
$$x^2E''(x)+4xE'(x)=2.$$
记 \(Y(x)=E'(x)\)。那么
$$x^2Y'(x)+4xY(x)=2.$$
两边除以 \(x^2\)，可得一阶线性方程
$$Y'(x)+\frac{4}{x}Y(x)=\frac{2}{x^2}.$$
积分因子为 \(x^4\)，因此
$$\frac{d}{dx}\left(x^4Y(x)\right)=2x^2.$$
积分后得到
$$x^4Y(x)=\frac{2}{3}x^3+C_1,$$
也就是
$$E'(x)=\frac{2}{3x}+\frac{C_1}{x^4}.$$
再积分一次：
$$E(x)=\frac{2}{3}\ln(x)-\frac{C_1}{3x^3}+C_2.$$
因此主导项是对数项，而代码中的 \(x^{-3}\) 小修正项来自齐次解。
当 \(x=1\) 时，积分项消失，所以
$$E(1)=1.$$
把 \(x=1\) 代入一次求导后的等式：
$$2\bigl(E(1)-1\bigr)+E'(1)=0.$$
由于 \(E(1)=1\)，可知
$$E'(1)=0.$$
代回通解：
$$0=E'(1)=\frac{2}{3}+C_1,$$
所以 \(C_1=-\frac{2}{3}\)。再由
$$1=E(1)=0+\frac{2}{9}+C_2$$
得到 \(C_2=\frac{7}{9}\)。最终有
$$\boxed{E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3}.}$$
C++ 文件中的三个校验点与公式完全一致：
$$E(1)=1.0000000000,$$
$$E(2)=1.2676536759,$$
$$E(7.5)=2.1215732071.$$
默认输入下的结果是
$$E(40)=3.2370342194.$$
当 \(x\) 很大时，期望值的渐近式为
$$E(x)=\frac{2}{3}\ln(x)+\frac{7}{9}+O(x^{-3}).$$
C++ 实现把数学核心放在 expected_steps(x) 中，并使用 long double 计算闭式。它还支持可选参数 --x=...、跳过校验的 --skip-checkpoints，以及用 nearly_equal 检查 \(E(1)\)、\(E(2)\)、\(E(7.5)\) 三个已知值。
Python 和 Java 版本则刻意保持极简，只保留同一个公式，并把 \(x=40\) 的结果格式化为 10 位小数。三个版本在运行时都没有做随机模拟、数值积分或动态规划。
闭式一旦得到，求值只需要一次对数运算和常数次加减乘除。因此时间复杂度是 \(O(1)\)，空间复杂度也是 \(O(1)\)。与蒙特卡洛采样或迭代积分相比，这种做法既更快，也能直接给出高精度结果。
Локальные решения описывают процесс Eating Pie через функцию математического ожидания \(E(x)\). Здесь \(x \ge 1\) обозначает текущее нормированное состояние; один шаг тратится сразу, после чего состояние умножается на случайный множитель \(W \in [0,1]\) с плотностью \(2(1-w)\). Как только состояние становится меньше 1, процесс прекращается, поэтому граничное условие имеет вид \(E(u)=0\) при \(u \lt 1\).
Следовательно, задача программы состоит не в случайном моделировании, а в точном выводе формулы для \(E(x)\) и ее вычислении, прежде всего при значении по умолчанию \(x=40\).
Все три реализации используют одну и ту же замкнутую формулу:
$$E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3},\qquad x\ge1.$$
Зафиксируем первый случайный множитель \(W\). Один шаг совершается немедленно, а оставшееся ожидание равно \(E(xW)\). Усреднение по плотности \(2(1-w)\) дает уравнение восстановления, прямо упомянутое в комментарии C++:
$$E(x)=1+2\int_0^1 (1-w)\,E(xw)\,dw.$$
Так как \(E(xw)=0\) при \(xw \lt 1\), для \(x \ge 1\) ненулевая часть интеграла начинается только с \(w=1/x\). Поэтому можно переписать уравнение так:
$$E(x)=1+2\int_{1/x}^1 (1-w)\,E(xw)\,dw.$$
Сделаем замену \(t=xw\). Тогда \(dw=dt/x\) и \(w=t/x\), поэтому получаем интеграл типа Вольтерры с фиксированным нижним пределом:
$$E(x)=1+\frac{2}{x^2}\int_1^x (x-t)\,E(t)\,dt.$$
Эта форма удобнее, потому что неизвестная функция входит как \(E(t)\), а не как \(E(xw)\).
Обозначим
$$I(x)=\int_1^x (x-t)\,E(t)\,dt.$$
Тогда интегральное уравнение принимает вид
$$x^2\bigl(E(x)-1\bigr)=2I(x).$$
После первого дифференцирования получаем
$$2x\bigl(E(x)-1\bigr)+x^2E'(x)=2\int_1^x E(t)\,dt.$$
После второго:
$$2\bigl(E(x)-1\bigr)+4xE'(x)+x^2E''(x)=2E(x).$$
Члены с \(E(x)\) сокращаются, и остается
$$x^2E''(x)+4xE'(x)=2.$$
Положим \(Y(x)=E'(x)\). Тогда
$$x^2Y'(x)+4xY(x)=2.$$
После деления на \(x^2\) имеем линейное уравнение первого порядка:
$$Y'(x)+\frac{4}{x}Y(x)=\frac{2}{x^2}.$$
Интегрирующий множитель равен \(x^4\), следовательно
$$\frac{d}{dx}\left(x^4Y(x)\right)=2x^2.$$
Интегрируя, получаем
$$x^4Y(x)=\frac{2}{3}x^3+C_1,$$
то есть
$$E'(x)=\frac{2}{3x}+\frac{C_1}{x^4}.$$
Еще одно интегрирование дает
$$E(x)=\frac{2}{3}\ln(x)-\frac{C_1}{3x^3}+C_2.$$
Главный вклад логарифмический, а однородная часть порождает поправку порядка \(x^{-3}\), которая и видна в коде.
При \(x=1\) интегральная часть исчезает, значит
$$E(1)=1.$$
Подставим \(x=1\) в уже продифференцированное равенство:
$$2\bigl(E(1)-1\bigr)+E'(1)=0.$$
Так как \(E(1)=1\), получаем
$$E'(1)=0.$$
Теперь подставим это в общее решение:
$$0=E'(1)=\frac{2}{3}+C_1,$$
откуда \(C_1=-\frac{2}{3}\). Далее
$$1=E(1)=0+\frac{2}{9}+C_2,$$
поэтому \(C_2=\frac{7}{9}\). В итоге
$$\boxed{E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3}.}$$
Контрольные точки из C++ совпадают с формулой:
$$E(1)=1.0000000000,$$
$$E(2)=1.2676536759,$$
$$E(7.5)=2.1215732071.$$
Для стандартного входа Project Euler получаем
$$E(40)=3.2370342194.$$
При больших \(x\) действует асимптотика
$$E(x)=\frac{2}{3}\ln(x)+\frac{7}{9}+O(x^{-3}).$$
В версии на C++ вся математика собрана в функции expected_steps(x), которая вычисляет замкнутую формулу в типе long double. Также поддерживаются параметр --x=..., флаг --skip-checkpoints и функция nearly_equal для проверки известных значений \(E(1)\), \(E(2)\) и \(E(7.5)\).
Версии на Python и Java намеренно минималистичны: они вычисляют ту же формулу и форматируют результат для \(x=40\) с 10 знаками после запятой. Ни одна из реализаций не выполняет численное интегрирование или моделирование во время работы.
После вывода закрытой формулы вычисление требует одного логарифма и константного числа арифметических операций. Следовательно, время работы равно \(O(1)\), а память тоже \(O(1)\). Это существенно лучше, чем статистическое моделирование или повторная численная квадратура.
الحلول المحلية تصف عملية Eating Pie بواسطة دالة توقع \(E(x)\). تمثل \(x \ge 1\) الحالة المعيارية الحالية؛ نستهلك خطوة واحدة فورًا، ثم تُضرب الحالة في عامل عشوائي \(W \in [0,1]\) كثافته \(2(1-w)\). وعندما تصبح الحالة أقل من 1 تنتهي العملية، لذلك يكون شرط الحد هو \(E(u)=0\) عندما \(u \lt 1\).
إذًا المطلوب برمجيًا ليس تنفيذ محاكاة عشوائية كثيرة، بل اشتقاق الصيغة الدقيقة لـ \(E(x)\) ثم تقييمها مباشرة، وخصوصًا عند القيمة الافتراضية \(x=40\).
جميع تطبيقات C++ و Python و Java تستخدم الصيغة المغلقة نفسها:
$$E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3},\qquad x\ge1.$$
نشرط على أول عامل عشوائي \(W\). نحن ندفع خطوة واحدة الآن دائمًا، وبعدها يصبح التوقع المتبقي هو \(E(xW)\). بأخذ المتوسط بالنسبة إلى الكثافة \(2(1-w)\) نحصل على معادلة التجديد المذكورة في تعليق C++:
$$E(x)=1+2\int_0^1 (1-w)\,E(xw)\,dw.$$
وبما أن \(E(xw)=0\) عندما \(xw \lt 1\)، فإن الجزء الفعلي من التكامل عند \(x \ge 1\) يبدأ من \(w=1/x\). لذا نستطيع أيضًا كتابة
$$E(x)=1+2\int_{1/x}^1 (1-w)\,E(xw)\,dw.$$
نأخذ التبديل \(t=xw\)، وعندها \(dw=dt/x\) و \(w=t/x\). فيتحول التكامل إلى صيغة من نوع Volterra ذات حد سفلي ثابت:
$$E(x)=1+\frac{2}{x^2}\int_1^x (x-t)\,E(t)\,dt.$$
هذه الصيغة أسهل بكثير لأن الدالة المجهولة تظهر الآن على صورة \(E(t)\) بدلًا من \(E(xw)\).
لنعرّف
$$I(x)=\int_1^x (x-t)\,E(t)\,dt.$$
فتصبح المعادلة التكاملية
$$x^2\bigl(E(x)-1\bigr)=2I(x).$$
بالتفاضل مرة واحدة نحصل على
$$2x\bigl(E(x)-1\bigr)+x^2E'(x)=2\int_1^x E(t)\,dt.$$
وبالتفاضل مرة ثانية نحصل على
$$2\bigl(E(x)-1\bigr)+4xE'(x)+x^2E''(x)=2E(x).$$
تُختزل حدود \(E(x)\)، ويبقى
$$x^2E''(x)+4xE'(x)=2.$$
لنكتب \(Y(x)=E'(x)\). عندئذ
$$x^2Y'(x)+4xY(x)=2.$$
وبعد القسمة على \(x^2\) نحصل على معادلة خطية من الرتبة الأولى:
$$Y'(x)+\frac{4}{x}Y(x)=\frac{2}{x^2}.$$
عامل التكامل هو \(x^4\)، وبالتالي
$$\frac{d}{dx}\left(x^4Y(x)\right)=2x^2.$$
بالتكامل:
$$x^4Y(x)=\frac{2}{3}x^3+C_1,$$
ومن ثم
$$E'(x)=\frac{2}{3x}+\frac{C_1}{x^4}.$$
وبتكامل جديد نحصل على
$$E(x)=\frac{2}{3}\ln(x)-\frac{C_1}{3x^3}+C_2.$$
إذن الحد الرئيسي لوغاريتمي، أما الحد \(x^{-3}\) الظاهر في الكود فيأتي من الحل المتجانس.
عندما \(x=1\) يختفي الجزء التكاملي، لذلك
$$E(1)=1.$$
وبالتعويض عن \(x=1\) في المعادلة بعد التفاضل مرة واحدة نجد
$$2\bigl(E(1)-1\bigr)+E'(1)=0.$$
ومادام \(E(1)=1\)، ينتج
$$E'(1)=0.$$
نعوض الآن في الحل العام:
$$0=E'(1)=\frac{2}{3}+C_1,$$
ومنها \(C_1=-\frac{2}{3}\). كذلك
$$1=E(1)=0+\frac{2}{9}+C_2,$$
فتكون \(C_2=\frac{7}{9}\). إذن
$$\boxed{E(x)=\frac{7}{9}+\frac{2}{3}\ln(x)+\frac{2}{9x^3}.}$$
نقاط التحقق الموجودة في ملف C++ تتطابق تمامًا مع الصيغة:
$$E(1)=1.0000000000,$$
$$E(2)=1.2676536759,$$
$$E(7.5)=2.1215732071.$$
وللقيمة الافتراضية في المسألة نحصل على
$$E(40)=3.2370342194.$$
أما عندما يكبر \(x\)، فالسلوك التقاربي هو
$$E(x)=\frac{2}{3}\ln(x)+\frac{7}{9}+O(x^{-3}).$$
نسخة C++ تجمع الرياضيات كلها داخل الدالة expected_steps(x) وتقيّم الصيغة باستخدام long double. كما تدعم الوسيط الاختياري --x=... والمفتاح --skip-checkpoints، وتستخدم الدالة nearly_equal للتحقق من القيم \(E(1)\) و \(E(2)\) و \(E(7.5)\).
أما نسختا Python و Java فهما أبسط عمدًا: كل منهما يحسب التعبير نفسه ويعرض نتيجة \(x=40\) بدقة 10 منازل عشرية. لا توجد في أي نسخة محاكاة عشوائية أو تكامل عددي أثناء التنفيذ.
بعد معرفة الصيغة المغلقة، لا يلزم سوى عملية لوغاريتم واحدة وعدد ثابت من عمليات الجمع والطرح والضرب والقسمة. لذلك يكون الزمن \(O(1)\) والذاكرة \(O(1)\). وهذا أفضل بكثير من المحاكاة الاحتمالية أو التكامل العددي التكراري من حيث السرعة والدقة.