Problem Summary

The task asks for a rope-loop probability \(P(n)\) after \(n-1\) random joining operations. The implementations do not estimate this with Monte Carlo simulation. Instead, they keep the exact distribution of the remaining continuous parameter and update it symbolically from one step to the next.

For the actual Project Euler instance the quantity of interest is \(P(80)\). The central idea is that every intermediate configuration can be reduced to a discrete state index and a continuous parameter \(r \in [0,1]\), so the whole process becomes an exact dynamic program on polynomial densities.

Mathematical Approach

Write \(m=n-1\). After \(k\) joining steps, the implementations represent the process by densities

$$f_k(s,r), \qquad -k \le s \le k,\quad 0 \le r \le 1,$$

where \(s\) is the discrete state and \(r\) is the remaining continuous parameter. The desired probability is the total mass of the terminal state \(s=0\) after \(m\) steps.

Step 1: Initialize the State Distribution

At the start there is only one possible discrete state, so the initial density is uniform:

$$f_0(0,r)=1,\qquad f_0(s,r)=0 \text{ for } s\ne 0.$$

Because \(\int_0^1 1\,dr=1\), this is already a properly normalized probability distribution.

Step 2: Derive the Three Transition Kernels

Suppose a current state contributes density \(P(t)\) in the old variable \(t\). The next step can send mass to the same discrete state, to the state \(s+1\), or to the state \(s-1\). The three kernels encoded by the implementations are

$$K_0(r,t)=1-|r-t|,$$

$$K_+(r,t)=\max(r-t,0),$$

$$K_-(r,t)=\max(t-r,0).$$

Therefore one source density \(P\) contributes

$$f_{k+1}(s,r)\mathrel{+}= \int_0^1 K_0(r,t)P(t)\,dt,$$

$$f_{k+1}(s+1,r)\mathrel{+}= \int_0^1 K_+(r,t)P(t)\,dt=\int_0^r (r-t)P(t)\,dt,$$

$$f_{k+1}(s-1,r)\mathrel{+}= \int_0^1 K_-(r,t)P(t)\,dt=\int_r^1 (t-r)P(t)\,dt.$$

The identity

$$K_0(r,t)+K_+(r,t)+K_-(r,t)=1$$

shows immediately that total probability mass is preserved.

Step 3: Rewrite the Kernels with Antiderivatives

If

$$P(t)=\sum_{i=0}^{d} c_i t^i,$$

define the two prefix integrals and the two full moments

$$A(r)=\int_0^r P(t)\,dt,\qquad B(r)=\int_0^r tP(t)\,dt,$$

$$T=\int_0^1 P(t)\,dt,\qquad T_1=\int_0^1 tP(t)\,dt.$$

Then the three transition contributions become

$$\int_0^r (r-t)P(t)\,dt = rA(r)-B(r),$$

$$\int_r^1 (t-r)P(t)\,dt = \bigl(T_1-B(r)\bigr)-r\bigl(T-A(r)\bigr),$$

$$\int_0^1 (1-|r-t|)P(t)\,dt = (1-r)A(r)+B(r)+(1+r)\bigl(T-A(r)\bigr)-\bigl(T_1-B(r)\bigr).$$

These are exactly the polynomial combinations constructed by the C++, Python, and Java implementations.

Step 4: Why Polynomial Densities Are Enough

The initial density has degree \(0\). If \(P\) has degree \(d\), then \(A\) has degree at most \(d+1\), \(B\) has degree at most \(d+2\), and each transition formula uses only linear combinations of \(A\), \(B\), \(rA\), and \(r(T-A)\). So one step increases the degree by at most \(2\).

After \(k\) steps every density therefore has degree at most \(2k\), and after all \(m=n-1\) steps the bound is

$$\deg f_m(s,\cdot)\le 2m.$$

That is why the implementations can store each density as a fixed coefficient array of length \(2m+1\).

Step 5: Worked Example for \(n=3\)

Here \(m=2\). Starting from \(f_0(0,r)=1\), one transition gives

$$f_1(0,r)=\frac{1}{2}+r-r^2,$$

$$f_1(1,r)=\frac{r^2}{2},$$

$$f_1(-1,r)=\frac{(1-r)^2}{2}=\frac{1}{2}-r+\frac{r^2}{2}.$$

Integrating over \(r\in[0,1]\) yields masses

$$\int_0^1 f_1(0,r)\,dr=\frac{2}{3},\qquad \int_0^1 f_1(1,r)\,dr=\frac{1}{6},\qquad \int_0^1 f_1(-1,r)\,dr=\frac{1}{6},$$

which already sum to \(1\).

Applying the transition a second time, the central density becomes

$$f_2(0,r)=\frac{11}{24}+\frac{1}{2}r-\frac{1}{4}r^2-\frac{1}{2}r^3+\frac{1}{4}r^4.$$

Hence

$$P(3)=\int_0^1 f_2(0,r)\,dr=\frac{11}{20}.$$

The implementations also agree on the further checkpoint

$$P(5)=\frac{15619}{36288}.$$

Step 6: Extract the Final Probability

After all \(m=n-1\) steps, the answer is simply

$$P(n)=\int_0^1 f_m(0,r)\,dr.$$

The normalization identity

$$\sum_{s=-m}^{m}\int_0^1 f_m(s,r)\,dr=1$$

is a strong internal consistency check and is explicitly verified by the high-precision implementation.

How the Code Works

The C++, Python, and Java implementations keep one table of polynomial coefficients for the current step and one for the next step. For each occupied discrete state they compute the two prefix antiderivatives \(A\) and \(B\), together with the full integrals \(T\) and \(T_1\). Because integrating a monomial only divides its coefficient by \(i+1\) or \(i+2\), every update is exact at the coefficient level.

Those polynomial pieces are then added to the three destination states according to the formulas above. After the final step, the implementation integrates the polynomial attached to state \(0\) over \([0,1]\) and prints the result. The C++ version additionally checks the exact small cases \(P(3)=11/20\), \(P(5)=15619/36288\), and normalization, while the Python and Java versions apply the same recurrence to produce the final decimal value.

Complexity Analysis

There are \(m=n-1\) stages. At stage \(k\), at most \(2k+1=O(n)\) discrete states are active, and each density uses \(O(n)\) coefficients. Processing one state requires only a constant number of passes over those coefficients, so one stage costs \(O(n^2)\) arithmetic and the full run costs \(O(n^3)\).

The coefficient tables store \(O(n)\) states with \(O(n)\) coefficients each, so the memory usage is \(O(n^2)\). High-precision decimal arithmetic is used because the exact rational values quickly develop large denominators.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=807
  2. Probability density function: Wikipedia - Probability density function
  3. Dynamic programming: Wikipedia - Dynamic programming
  4. Continuous uniform distribution: Wikipedia - Continuous uniform distribution
  5. Integral: Wikipedia - Integral

Problemzusammenfassung

Gesucht ist eine Seil-Schleifen-Wahrscheinlichkeit \(P(n)\) nach \(n-1\) zufälligen Verknüpfungen. Die Implementierungen schätzen diese Größe nicht per Monte-Carlo-Simulation, sondern verfolgen die exakte Verteilung eines verbleibenden kontinuierlichen Parameters symbolisch von Schritt zu Schritt.

Für den eigentlichen Project-Euler-Fall wird \(P(80)\) berechnet. Der entscheidende Punkt ist, dass sich jede Zwischenkonfiguration auf einen diskreten Zustandsindex und einen kontinuierlichen Parameter \(r \in [0,1]\) reduzieren lässt. Dadurch wird das Problem zu einer exakten dynamischen Programmierung auf Polynomdichten.

Mathematischer Ansatz

Setze \(m=n-1\). Nach \(k\) Verknüpfungen beschreiben die Implementierungen den Prozess durch Dichten

$$f_k(s,r), \qquad -k \le s \le k,\quad 0 \le r \le 1,$$

wobei \(s\) der diskrete Zustand und \(r\) der verbleibende kontinuierliche Parameter ist. Die gesuchte Wahrscheinlichkeit ist die Gesamtmasse des Endzustands \(s=0\) nach \(m\) Schritten.

Schritt 1: Anfangsverteilung festlegen

Zu Beginn gibt es nur einen möglichen diskreten Zustand, also ist die Anfangsdichte gleichmäßig:

$$f_0(0,r)=1,\qquad f_0(s,r)=0 \text{ for } s\ne 0.$$

Da \(\int_0^1 1\,dr=1\) gilt, ist diese Verteilung bereits korrekt normiert.

Schritt 2: Die drei Übergangskerne bestimmen

Nehme an, ein aktueller Zustand trägt die Dichte \(P(t)\) in der alten Variablen \(t\). Im nächsten Schritt kann Masse im gleichen diskreten Zustand bleiben, nach \(s+1\) gehen oder nach \(s-1\) gehen. Die in den Implementierungen kodierten Kerne sind

$$K_0(r,t)=1-|r-t|,$$

$$K_+(r,t)=\max(r-t,0),$$

$$K_-(r,t)=\max(t-r,0).$$

Daraus folgen die Beiträge

$$f_{k+1}(s,r)\mathrel{+}= \int_0^1 K_0(r,t)P(t)\,dt,$$

$$f_{k+1}(s+1,r)\mathrel{+}= \int_0^1 K_+(r,t)P(t)\,dt=\int_0^r (r-t)P(t)\,dt,$$

$$f_{k+1}(s-1,r)\mathrel{+}= \int_0^1 K_-(r,t)P(t)\,dt=\int_r^1 (t-r)P(t)\,dt.$$

Die Identität

$$K_0(r,t)+K_+(r,t)+K_-(r,t)=1$$

zeigt sofort, dass die gesamte Wahrscheinlichkeitsmasse erhalten bleibt.

Schritt 3: Die Kerne mit Stammfunktionen umschreiben

Sei

$$P(t)=\sum_{i=0}^{d} c_i t^i.$$

Dann definieren wir zwei Präfixintegrale und zwei Gesamtmomente:

$$A(r)=\int_0^r P(t)\,dt,\qquad B(r)=\int_0^r tP(t)\,dt,$$

$$T=\int_0^1 P(t)\,dt,\qquad T_1=\int_0^1 tP(t)\,dt.$$

Die drei Übergänge lassen sich dann schreiben als

$$\int_0^r (r-t)P(t)\,dt = rA(r)-B(r),$$

$$\int_r^1 (t-r)P(t)\,dt = \bigl(T_1-B(r)\bigr)-r\bigl(T-A(r)\bigr),$$

$$\int_0^1 (1-|r-t|)P(t)\,dt = (1-r)A(r)+B(r)+(1+r)\bigl(T-A(r)\bigr)-\bigl(T_1-B(r)\bigr).$$

Genau diese Polynomkombinationen bauen die C++-, Python- und Java-Implementierungen auf.

Schritt 4: Warum Polynome als Zustandsraum genügen

Die Anfangsdichte hat Grad \(0\). Hat \(P\) Grad \(d\), so hat \(A\) höchstens Grad \(d+1\), \(B\) höchstens Grad \(d+2\), und jede Übergangsformel benutzt nur lineare Kombinationen von \(A\), \(B\), \(rA\) und \(r(T-A)\). Ein Schritt erhöht den Grad also um höchstens \(2\).

Nach \(k\) Schritten hat jede Dichte daher Grad höchstens \(2k\), und nach allen \(m=n-1\) Schritten gilt

$$\deg f_m(s,\cdot)\le 2m.$$

Darum kann jede Dichte als festes Koeffizientenarray der Länge \(2m+1\) gespeichert werden.

Schritt 5: Durchgerechnetes Beispiel für \(n=3\)

Hier ist \(m=2\). Ausgehend von \(f_0(0,r)=1\) liefert ein Übergang

$$f_1(0,r)=\frac{1}{2}+r-r^2,$$

$$f_1(1,r)=\frac{r^2}{2},$$

$$f_1(-1,r)=\frac{(1-r)^2}{2}=\frac{1}{2}-r+\frac{r^2}{2}.$$

Nach Integration über \(r\in[0,1]\) erhält man

$$\int_0^1 f_1(0,r)\,dr=\frac{2}{3},\qquad \int_0^1 f_1(1,r)\,dr=\frac{1}{6},\qquad \int_0^1 f_1(-1,r)\,dr=\frac{1}{6},$$

also bereits eine normierte Verteilung.

Nach dem zweiten Übergang lautet die zentrale Dichte

$$f_2(0,r)=\frac{11}{24}+\frac{1}{2}r-\frac{1}{4}r^2-\frac{1}{2}r^3+\frac{1}{4}r^4.$$

Damit folgt

$$P(3)=\int_0^1 f_2(0,r)\,dr=\frac{11}{20}.$$

Ein weiterer Kontrollwert der Implementierung ist

$$P(5)=\frac{15619}{36288}.$$

Schritt 6: Die gesuchte Wahrscheinlichkeit auslesen

Nach allen \(m=n-1\) Schritten ist die Antwort schlicht

$$P(n)=\int_0^1 f_m(0,r)\,dr.$$

Die Normalisierungsidentität

$$\sum_{s=-m}^{m}\int_0^1 f_m(s,r)\,dr=1$$

ist eine starke Plausibilitätsprüfung und wird in der hochpräzisen Implementierung explizit kontrolliert.

Wie der Code arbeitet

Die C++-, Python- und Java-Implementierungen halten eine Tabelle von Polynomkoeffizienten für den aktuellen Schritt und eine zweite für den nächsten Schritt. Für jeden belegten diskreten Zustand werden die beiden Präfix-Stammfunktionen \(A\) und \(B\) sowie die Gesamtintegrale \(T\) und \(T_1\) berechnet. Da das Integrieren eines Monoms nur eine Division des Koeffizienten durch \(i+1\) oder \(i+2\) erfordert, erfolgen die Updates auf Koeffizientenebene exakt.

Diese Polynomstücke werden dann gemäß den obigen Formeln auf die drei Zielzustände verteilt. Nach dem letzten Schritt integriert die Implementierung das Polynom im Zustand \(0\) über \([0,1]\) und gibt das Ergebnis aus. Die C++-Version prüft zusätzlich die exakten Kleinwerte \(P(3)=11/20\), \(P(5)=15619/36288\) sowie die Normierung, während die Python- und Java-Versionen dieselbe Rekurrenz zur Erzeugung des Endwerts verwenden.

Komplexitätsanalyse

Es gibt \(m=n-1\) Stufen. In Stufe \(k\) sind höchstens \(2k+1=O(n)\) diskrete Zustände aktiv, und jede Dichte besitzt \(O(n)\) Koeffizienten. Ein Zustand wird mit nur konstant vielen Durchläufen über diese Koeffizienten verarbeitet. Eine Stufe kostet daher \(O(n^2)\) Operationen, und der Gesamtlauf benötigt \(O(n^3)\).

Die Koeffiziententabellen speichern \(O(n)\) Zustände mit jeweils \(O(n)\) Koeffizienten, also \(O(n^2)\) Speicher. Hochpräzise Dezimalarithmetik ist sinnvoll, weil die exakten rationalen Werte schnell große Nenner entwickeln.

Fußnoten und Quellen

  1. Problemseite: https://projecteuler.net/problem=807
  2. Wahrscheinlichkeitsdichtefunktion: Wikipedia - Probability density function
  3. Dynamische Programmierung: Wikipedia - Dynamic programming
  4. Stetige Gleichverteilung: Wikipedia - Continuous uniform distribution
  5. Integral: Wikipedia - Integral

Problem Özeti

Sorulan büyüklük, \(n-1\) rastgele bağlama adımından sonra ortaya çıkan ip-ilmek olasılığı \(P(n)\)'dir. İmplementasyonlar bu değeri Monte Carlo simülasyonuyla tahmin etmez; bunun yerine elde kalan sürekli parametrenin tam dağılımını adım adım sembolik olarak taşır.

Project Euler örneğinde hedeflenen değer \(P(80)\)'dir. Temel fikir, her ara konfigürasyonun bir ayrık durum indisi ile \(r \in [0,1]\) aralığında bir sürekli parametreye indirgenebilmesidir. Böylece süreç, polinom yoğunluklar üzerinde çalışan tam bir dinamik program haline gelir.

Matematiksel Yaklaşım

\(m=n-1\) yazalım. \(k\) bağlama adımından sonra implementasyonlar süreci

$$f_k(s,r), \qquad -k \le s \le k,\quad 0 \le r \le 1,$$

yoğunluklarıyla temsil eder. Burada \(s\) ayrık durum, \(r\) ise elde kalan sürekli parametredir. İstenen olasılık, \(m\) adım sonunda \(s=0\) durumunda kalan toplam kütledir.

Adım 1: Başlangıç dağılımını kur

Başlangıçta yalnızca tek bir ayrık durum mümkündür; dolayısıyla başlangıç yoğunluğu düzgündür:

$$f_0(0,r)=1,\qquad f_0(s,r)=0 \text{ for } s\ne 0.$$

\(\int_0^1 1\,dr=1\) olduğu için bu zaten normalize edilmiş bir olasılık dağılımıdır.

Adım 2: Üç geçiş çekirdeğini çıkar

Eski değişken \(t\) için bir kaynak yoğunluk \(P(t)\) alalım. Bir sonraki adımda kütle aynı ayrık durumda kalabilir, \(s+1\)'e gidebilir ya da \(s-1\)'e gidebilir. Implementasyonların kodladığı üç çekirdek şunlardır:

$$K_0(r,t)=1-|r-t|,$$

$$K_+(r,t)=\max(r-t,0),$$

$$K_-(r,t)=\max(t-r,0).$$

Buna göre bir kaynak yoğunluğun katkıları

$$f_{k+1}(s,r)\mathrel{+}= \int_0^1 K_0(r,t)P(t)\,dt,$$

$$f_{k+1}(s+1,r)\mathrel{+}= \int_0^1 K_+(r,t)P(t)\,dt=\int_0^r (r-t)P(t)\,dt,$$

$$f_{k+1}(s-1,r)\mathrel{+}= \int_0^1 K_-(r,t)P(t)\,dt=\int_r^1 (t-r)P(t)\,dt$$

olur. Ayrıca

$$K_0(r,t)+K_+(r,t)+K_-(r,t)=1$$

eşitliği toplam olasılık kütlesinin korunduğunu doğrudan gösterir.

Adım 3: Çekirdekleri antitürevlerle yaz

Eğer

$$P(t)=\sum_{i=0}^{d} c_i t^i$$

ise iki önek integral ve iki tam moment tanımlayalım:

$$A(r)=\int_0^r P(t)\,dt,\qquad B(r)=\int_0^r tP(t)\,dt,$$

$$T=\int_0^1 P(t)\,dt,\qquad T_1=\int_0^1 tP(t)\,dt.$$

Böylece üç katkı şu hale gelir:

$$\int_0^r (r-t)P(t)\,dt = rA(r)-B(r),$$

$$\int_r^1 (t-r)P(t)\,dt = \bigl(T_1-B(r)\bigr)-r\bigl(T-A(r)\bigr),$$

$$\int_0^1 (1-|r-t|)P(t)\,dt = (1-r)A(r)+B(r)+(1+r)\bigl(T-A(r)\bigr)-\bigl(T_1-B(r)\bigr).$$

C++, Python ve Java implementasyonlarının kurduğu polinom kombinasyonları tam olarak bunlardır.

Adım 4: Neden polinomlar kapalı bir sınıf oluşturur

Başlangıç yoğunluğunun derecesi \(0\)'dır. \(P\)'nin derecesi \(d\) ise \(A\)'nın derecesi en fazla \(d+1\), \(B\)'nin derecesi en fazla \(d+2\) olur; geçiş formülleri de yalnızca \(A\), \(B\), \(rA\) ve \(r(T-A)\) gibi doğrusal birleşimler kullanır. Bu yüzden tek bir adım dereceyi en fazla \(2\) artırır.

Dolayısıyla \(k\) adımdan sonra her yoğunluk en fazla \(2k\) derecelidir; bütün süreç sonunda

$$\deg f_m(s,\cdot)\le 2m$$

elde edilir. Bu nedenle her yoğunluğu uzunluğu \(2m+1\) olan sabit bir katsayı dizisiyle saklamak yeterlidir.

Adım 5: \(n=3\) için çalışılmış örnek

Burada \(m=2\)'dir. \(f_0(0,r)=1\) başlangıcından bir geçiş sonra

$$f_1(0,r)=\frac{1}{2}+r-r^2,$$

$$f_1(1,r)=\frac{r^2}{2},$$

$$f_1(-1,r)=\frac{(1-r)^2}{2}=\frac{1}{2}-r+\frac{r^2}{2}$$

elde edilir. Bunların \([0,1]\) üzerindeki integralleri

$$\int_0^1 f_1(0,r)\,dr=\frac{2}{3},\qquad \int_0^1 f_1(1,r)\,dr=\frac{1}{6},\qquad \int_0^1 f_1(-1,r)\,dr=\frac{1}{6}$$

olup toplamları yine \(1\)'dir.

Geçiş bir kez daha uygulandığında merkez yoğunluk

$$f_2(0,r)=\frac{11}{24}+\frac{1}{2}r-\frac{1}{4}r^2-\frac{1}{2}r^3+\frac{1}{4}r^4$$

olur. Buradan

$$P(3)=\int_0^1 f_2(0,r)\,dr=\frac{11}{20}$$

sonucu çıkar. Implementasyonların uyuştuğu bir başka kontrol değeri de

$$P(5)=\frac{15619}{36288}$$

değeridir.

Adım 6: Son olasılığı çıkar

Bütün \(m=n-1\) adımlar tamamlandığında cevap doğrudan

$$P(n)=\int_0^1 f_m(0,r)\,dr$$

şeklindedir. Ayrıca

$$\sum_{s=-m}^{m}\int_0^1 f_m(s,r)\,dr=1$$

özdeşliği güçlü bir iç tutarlılık kontrolüdür ve yüksek hassasiyetli implementasyon tarafından açıkça doğrulanır.

Kod Nasıl Çalışır

C++, Python ve Java implementasyonları mevcut adım için bir polinom katsayı tablosu ve sonraki adım için ikinci bir tablo tutar. Her dolu ayrık durum için iki önek antitürev \(A\) ve \(B\) ile tam integraller \(T\) ve \(T_1\) hesaplanır. Bir monomu entegre etmek yalnızca katsayısını \(i+1\) ya da \(i+2\) ile bölmek anlamına geldiğinden, tüm güncellemeler katsayı seviyesinde tam biçimde yapılır.

Bu polinom parçaları daha sonra yukarıdaki formüllere göre üç hedef duruma eklenir. Son adımın ardından implementasyon, durum \(0\)'a ait polinomu \([0,1]\) üzerinde integre eder ve sonucu yazdırır. C++ sürümü ayrıca \(P(3)=11/20\), \(P(5)=15619/36288\) ve normalizasyon kontrollerini yapar; Python ve Java sürümleri aynı bağıntıyı kullanarak nihai ondalık değeri üretir.

Karmaşıklık Analizi

\(m=n-1\) adet aşama vardır. \(k\). aşamada en fazla \(2k+1=O(n)\) ayrık durum aktiftir ve her yoğunluk \(O(n)\) katsayı taşır. Bir durumun işlenmesi bu katsayılar üzerinde yalnızca sabit sayıda geçiş gerektirir. Bu yüzden tek bir aşama \(O(n^2)\), tüm süreç ise \(O(n^3)\) aritmetik işlem maliyetine sahiptir.

Katsayı tabloları \(O(n)\) durum ve her biri için \(O(n)\) katsayı sakladığından bellek kullanımı \(O(n^2)\)'dir. Kesin rasyonel değerlerin paydaları hızla büyüdüğü için yüksek hassasiyetli ondalık aritmetik tercih edilir.

Dipnotlar ve Kaynaklar

  1. Problem sayfası: https://projecteuler.net/problem=807
  2. Olasılık yoğunluk fonksiyonu: Wikipedia - Probability density function
  3. Dinamik programlama: Wikipedia - Dynamic programming
  4. Sürekli düzgün dağılım: Wikipedia - Continuous uniform distribution
  5. İntegral: Wikipedia - Integral

Resumen del Problema

La cantidad pedida es una probabilidad \(P(n)\) asociada a \(n-1\) uniones aleatorias de cuerdas. Las implementaciones no la aproximan con simulación Monte Carlo. En cambio, conservan de forma exacta la distribución del parámetro continuo restante y la actualizan simbólicamente en cada paso.

En el caso real de Project Euler se necesita \(P(80)\). La observación clave es que cada configuración intermedia puede comprimirse en un índice de estado discreto y en un parámetro continuo \(r \in [0,1]\). Así, todo el proceso se convierte en una programación dinámica exacta sobre densidades polinómicas.

Enfoque Matemático

Escribamos \(m=n-1\). Después de \(k\) uniones, las implementaciones describen el proceso mediante densidades

$$f_k(s,r), \qquad -k \le s \le k,\quad 0 \le r \le 1,$$

donde \(s\) es el estado discreto y \(r\) es el parámetro continuo restante. La probabilidad buscada es la masa total del estado terminal \(s=0\) tras \(m\) pasos.

Paso 1: Fijar la distribución inicial

Al principio solo existe un estado discreto posible, por lo que la densidad inicial es uniforme:

$$f_0(0,r)=1,\qquad f_0(s,r)=0 \text{ for } s\ne 0.$$

Como \(\int_0^1 1\,dr=1\), la distribución ya está normalizada.

Paso 2: Obtener los tres núcleos de transición

Supongamos que un estado actual aporta una densidad \(P(t)\) en la variable antigua \(t\). En el siguiente paso, la masa puede quedarse en el mismo estado discreto, pasar a \(s+1\) o pasar a \(s-1\). Los tres núcleos codificados por las implementaciones son

$$K_0(r,t)=1-|r-t|,$$

$$K_+(r,t)=\max(r-t,0),$$

$$K_-(r,t)=\max(t-r,0).$$

Por tanto, las contribuciones son

$$f_{k+1}(s,r)\mathrel{+}= \int_0^1 K_0(r,t)P(t)\,dt,$$

$$f_{k+1}(s+1,r)\mathrel{+}= \int_0^1 K_+(r,t)P(t)\,dt=\int_0^r (r-t)P(t)\,dt,$$

$$f_{k+1}(s-1,r)\mathrel{+}= \int_0^1 K_-(r,t)P(t)\,dt=\int_r^1 (t-r)P(t)\,dt.$$

La identidad

$$K_0(r,t)+K_+(r,t)+K_-(r,t)=1$$

muestra de inmediato que la masa total de probabilidad se conserva.

Paso 3: Reescribir los núcleos con primitivas

Si

$$P(t)=\sum_{i=0}^{d} c_i t^i,$$

definimos dos integrales prefijo y dos momentos globales:

$$A(r)=\int_0^r P(t)\,dt,\qquad B(r)=\int_0^r tP(t)\,dt,$$

$$T=\int_0^1 P(t)\,dt,\qquad T_1=\int_0^1 tP(t)\,dt.$$

Entonces las tres contribuciones quedan como

$$\int_0^r (r-t)P(t)\,dt = rA(r)-B(r),$$

$$\int_r^1 (t-r)P(t)\,dt = \bigl(T_1-B(r)\bigr)-r\bigl(T-A(r)\bigr),$$

$$\int_0^1 (1-|r-t|)P(t)\,dt = (1-r)A(r)+B(r)+(1+r)\bigl(T-A(r)\bigr)-\bigl(T_1-B(r)\bigr).$$

Estas son exactamente las combinaciones polinómicas que construyen las implementaciones en C++, Python y Java.

Paso 4: Por qué basta trabajar con polinomios

La densidad inicial tiene grado \(0\). Si \(P\) tiene grado \(d\), entonces \(A\) tiene grado como mucho \(d+1\), \(B\) tiene grado como mucho \(d+2\), y cada fórmula de transición usa solo combinaciones lineales de \(A\), \(B\), \(rA\) y \(r(T-A)\). Por ello, un paso aumenta el grado a lo sumo en \(2\).

Después de \(k\) pasos, cada densidad tiene entonces grado a lo sumo \(2k\), y tras todos los \(m=n-1\) pasos se obtiene

$$\deg f_m(s,\cdot)\le 2m.$$

Por eso las implementaciones pueden guardar cada densidad en un arreglo fijo de coeficientes de longitud \(2m+1\).

Paso 5: Ejemplo trabajado para \(n=3\)

Aquí \(m=2\). Partiendo de \(f_0(0,r)=1\), una transición produce

$$f_1(0,r)=\frac{1}{2}+r-r^2,$$

$$f_1(1,r)=\frac{r^2}{2},$$

$$f_1(-1,r)=\frac{(1-r)^2}{2}=\frac{1}{2}-r+\frac{r^2}{2}.$$

Al integrar en \(r\in[0,1]\) aparecen las masas

$$\int_0^1 f_1(0,r)\,dr=\frac{2}{3},\qquad \int_0^1 f_1(1,r)\,dr=\frac{1}{6},\qquad \int_0^1 f_1(-1,r)\,dr=\frac{1}{6},$$

que ya suman \(1\).

Tras aplicar la transición una segunda vez, la densidad central pasa a ser

$$f_2(0,r)=\frac{11}{24}+\frac{1}{2}r-\frac{1}{4}r^2-\frac{1}{2}r^3+\frac{1}{4}r^4.$$

Así se obtiene

$$P(3)=\int_0^1 f_2(0,r)\,dr=\frac{11}{20}.$$

Las implementaciones también coinciden en el punto de control adicional

$$P(5)=\frac{15619}{36288}.$$

Paso 6: Extraer la probabilidad final

Después de los \(m=n-1\) pasos completos, la respuesta es

$$P(n)=\int_0^1 f_m(0,r)\,dr.$$

La identidad de normalización

$$\sum_{s=-m}^{m}\int_0^1 f_m(s,r)\,dr=1$$

es una comprobación interna muy fuerte y la implementación de alta precisión la verifica explícitamente.

Cómo Funciona el Código

Las implementaciones en C++, Python y Java mantienen una tabla de coeficientes polinómicos para el paso actual y otra para el siguiente. Para cada estado discreto activo calculan las dos primitivas prefijo \(A\) y \(B\), junto con las integrales completas \(T\) y \(T_1\). Como integrar un monomio solo exige dividir su coeficiente entre \(i+1\) o \(i+2\), cada actualización es exacta al nivel de coeficientes.

Luego esas piezas polinómicas se añaden a los tres estados de destino siguiendo las fórmulas anteriores. Al final, la implementación integra el polinomio asociado al estado \(0\) en \([0,1]\) y muestra el resultado. La versión en C++ además comprueba los valores exactos \(P(3)=11/20\), \(P(5)=15619/36288\) y la normalización, mientras que las versiones en Python y Java aplican la misma recurrencia para producir el valor decimal final.

Análisis de Complejidad

Hay \(m=n-1\) etapas. En la etapa \(k\) pueden estar activos como máximo \(2k+1=O(n)\) estados discretos, y cada densidad usa \(O(n)\) coeficientes. Procesar un estado requiere solo un número constante de recorridos sobre esos coeficientes, de modo que una etapa cuesta \(O(n^2)\) y toda la ejecución cuesta \(O(n^3)\).

Las tablas de coeficientes almacenan \(O(n)\) estados con \(O(n)\) coeficientes cada uno, así que la memoria es \(O(n^2)\). Se emplea aritmética decimal de alta precisión porque los valores racionales exactos desarrollan denominadores grandes con rapidez.

Notas y Referencias

  1. Página del problema: https://projecteuler.net/problem=807
  2. Función de densidad de probabilidad: Wikipedia - Probability density function
  3. Programación dinámica: Wikipedia - Dynamic programming
  4. Distribución uniforme continua: Wikipedia - Continuous uniform distribution
  5. Integral: Wikipedia - Integral

问题概述

题目要求计算一个与绳索随机连接过程有关的概率 \(P(n)\)。整个过程一共进行 \(n-1\) 次随机连接。实现并没有采用 Monte Carlo 模拟,而是精确追踪每一步剩余连续参数的概率密度,因此得到的是确定性的高精度结果,而不是带随机误差的估计值。

在实际的 Project Euler 问题中,需要求的是 \(P(80)\)。关键观察是:虽然原始过程看起来是组合对象和连续量混合在一起的随机演化,但每个中间状态都可以压缩为一个离散状态指标和一个位于 \(r \in [0,1]\) 的连续参数。这样一来,整个问题就可以转写成在多项式密度上的精确动态规划。

数学方法

记 \(m=n-1\)。在完成 \(k\) 次连接之后,实现把系统表示为一族密度函数

$$f_k(s,r), \qquad -k \le s \le k,\quad 0 \le r \le 1,$$

其中 \(s\) 表示离散状态,\(r\) 表示仍然保留下来的连续参数。题目所需的最终概率,就是做完全部 \(m\) 步之后终态 \(s=0\) 的总概率质量。

步骤 1:建立初始分布

开始时只有一个可能的离散状态,因此初始密度是均匀的:

$$f_0(0,r)=1,\qquad f_0(s,r)=0 \text{ for } s\ne 0.$$

因为 \(\int_0^1 1\,dr=1\),所以这已经是一个正确归一化的概率分布。

步骤 2:写出三个转移核

设旧变量 \(t\) 上有一个来源密度 \(P(t)\)。下一步中,这部分概率质量可以留在原来的离散状态,也可以流向 \(s+1\) 或 \(s-1\)。实现所对应的三个核函数是

$$K_0(r,t)=1-|r-t|,$$

$$K_+(r,t)=\max(r-t,0),$$

$$K_-(r,t)=\max(t-r,0).$$

于是,一个来源密度对下一层的贡献分别为

$$f_{k+1}(s,r)\mathrel{+}= \int_0^1 K_0(r,t)P(t)\,dt,$$

$$f_{k+1}(s+1,r)\mathrel{+}= \int_0^1 K_+(r,t)P(t)\,dt=\int_0^r (r-t)P(t)\,dt,$$

$$f_{k+1}(s-1,r)\mathrel{+}= \int_0^1 K_-(r,t)P(t)\,dt=\int_r^1 (t-r)P(t)\,dt.$$

而且

$$K_0(r,t)+K_+(r,t)+K_-(r,t)=1$$

这个恒等式立刻说明总概率质量在转移过程中保持不变。这也是代码中归一化检查成立的根本原因。

步骤 3:把转移核改写成前缀积分

$$P(t)=\sum_{i=0}^{d} c_i t^i,$$

定义两个前缀积分和两个整体矩:

$$A(r)=\int_0^r P(t)\,dt,\qquad B(r)=\int_0^r tP(t)\,dt,$$

$$T=\int_0^1 P(t)\,dt,\qquad T_1=\int_0^1 tP(t)\,dt.$$

那么三个转移贡献分别可以写成

$$\int_0^r (r-t)P(t)\,dt = rA(r)-B(r),$$

$$\int_r^1 (t-r)P(t)\,dt = \bigl(T_1-B(r)\bigr)-r\bigl(T-A(r)\bigr),$$

$$\int_0^1 (1-|r-t|)P(t)\,dt = (1-r)A(r)+B(r)+(1+r)\bigl(T-A(r)\bigr)-\bigl(T_1-B(r)\bigr).$$

这正是 C++、Python 和 Java 实现中实际构造出来的多项式组合。由于每个单项式积分时只需要把系数除以 \(i+1\) 或 \(i+2\),所以整套更新在系数层面是完全精确的。

步骤 4:为什么始终停留在多项式空间内

初始密度的次数是 \(0\)。如果 \(P\) 的次数为 \(d\),那么 \(A\) 的次数至多为 \(d+1\),\(B\) 的次数至多为 \(d+2\),而转移公式只用到 \(A\)、\(B\)、\(rA\) 与 \(r(T-A)\) 这类线性组合。因此,一次转移最多只会让次数增加 \(2\)。

于是做完 \(k\) 步以后,每个密度的次数都不会超过 \(2k\)。特别地,在全部 \(m=n-1\) 步之后有

$$\deg f_m(s,\cdot)\le 2m.$$

这就解释了为什么实现可以为每个状态预留长度固定为 \(2m+1\) 的系数数组,而不需要动态改变表示结构。

步骤 5:\(n=3\) 的完整小例子

当 \(n=3\) 时,\(m=2\)。从 \(f_0(0,r)=1\) 出发,做一次转移后得到

$$f_1(0,r)=\frac{1}{2}+r-r^2,$$

$$f_1(1,r)=\frac{r^2}{2},$$

$$f_1(-1,r)=\frac{(1-r)^2}{2}=\frac{1}{2}-r+\frac{r^2}{2}.$$

把它们在 \([0,1]\) 上积分,可得对应概率质量

$$\int_0^1 f_1(0,r)\,dr=\frac{2}{3},\qquad \int_0^1 f_1(1,r)\,dr=\frac{1}{6},\qquad \int_0^1 f_1(-1,r)\,dr=\frac{1}{6}.$$

三项之和仍然等于 \(1\),说明一步之后分布依然完全归一化。

再做第二次转移,中心状态的密度变成

$$f_2(0,r)=\frac{11}{24}+\frac{1}{2}r-\frac{1}{4}r^2-\frac{1}{2}r^3+\frac{1}{4}r^4.$$

因此

$$P(3)=\int_0^1 f_2(0,r)\,dr=\frac{11}{20}.$$

实现还与另一个精确检查值一致:

$$P(5)=\frac{15619}{36288}.$$

步骤 6:提取最终答案

完成全部 \(m=n-1\) 步之后,答案就是

$$P(n)=\int_0^1 f_m(0,r)\,dr.$$

同时还应满足归一化恒等式

$$\sum_{s=-m}^{m}\int_0^1 f_m(s,r)\,dr=1.$$

高精度实现会显式检查这一点,因此不仅输出最终概率,也顺带验证整个递推过程没有丢失或重复计算概率质量。

代码如何工作

C++、Python 和 Java 实现都维护两张多项式系数表,一张对应当前步,一张对应下一步。对每一个当前非零的离散状态,它们先计算前缀积分 \(A\)、\(B\),再计算整体积分 \(T\)、\(T_1\)。这些量随后按照上面的三个公式分别加到相同状态、上一个状态和下一个状态的多项式中。

因为每一步都只是在有限长度的系数数组上做加法、减法和按小整数除法,所以实现非常稳定,而且不会引入数值积分误差。全部步骤结束后,程序只需要把状态 \(0\) 对应的最终多项式在 \([0,1]\) 上积分,就得到 \(P(n)\)。其中 C++ 版本还会额外核对 \(P(3)=11/20\)、\(P(5)=15619/36288\) 以及总质量是否为 \(1\);Python 和 Java 版本则使用同样的递推直接输出最终十进制结果。

复杂度分析

总共有 \(m=n-1\) 个阶段。第 \(k\) 个阶段最多有 \(2k+1=O(n)\) 个活跃的离散状态,而每个状态的密度用 \(O(n)\) 个系数表示。处理一个状态只需要对这些系数做常数次线性扫描,因此单个阶段的代价是 \(O(n^2)\),整体代价是 \(O(n^3)\)。

存储方面,需要为 \(O(n)\) 个状态各保存 \(O(n)\) 个系数,所以空间复杂度是 \(O(n^2)\)。由于精确有理数的分母会迅速变大,程序采用高精度十进制类型来保证最终结果的小数位稳定可靠。

注释与参考资料

  1. 题目页面: https://projecteuler.net/problem=807
  2. 概率密度函数: Wikipedia - Probability density function
  3. 动态规划: Wikipedia - Dynamic programming
  4. 连续均匀分布: Wikipedia - Continuous uniform distribution
  5. 积分: Wikipedia - Integral

Краткое описание задачи

Нужно вычислить вероятность \(P(n)\), связанную с процессом случайного соединения веревок за \(n-1\) шагов. Реализации не используют моделирование Монте-Карло. Вместо этого они точно переносят распределение оставшегося непрерывного параметра от шага к шагу, поэтому итоговое значение получается детерминированным и высокоточным.

Для самой задачи Project Euler требуется \(P(80)\). Ключевое наблюдение состоит в том, что каждую промежуточную конфигурацию можно сжать до дискретного индекса состояния и непрерывного параметра \(r \in [0,1]\). После этого задача превращается в точное динамическое программирование на полиномиальных плотностях.

Математический подход

Положим \(m=n-1\). После \(k\) соединений процесс описывается плотностями

$$f_k(s,r), \qquad -k \le s \le k,\quad 0 \le r \le 1,$$

где \(s\) обозначает дискретное состояние, а \(r\) - оставшийся непрерывный параметр. Искомая вероятность равна полной массе конечного состояния \(s=0\) после \(m\) шагов.

Шаг 1: Задать начальное распределение

В начале возможен только один дискретный статус, поэтому начальная плотность равномерна:

$$f_0(0,r)=1,\qquad f_0(s,r)=0 \text{ for } s\ne 0.$$

Так как \(\int_0^1 1\,dr=1\), это уже корректно нормированное распределение вероятности.

Шаг 2: Выписать три ядра перехода

Пусть некоторое текущее состояние имеет плотность \(P(t)\) по старой переменной \(t\). На следующем шаге масса может остаться в том же дискретном состоянии, перейти в \(s+1\) или в \(s-1\). Три ядра, зашитые в реализации, имеют вид

$$K_0(r,t)=1-|r-t|,$$

$$K_+(r,t)=\max(r-t,0),$$

$$K_-(r,t)=\max(t-r,0).$$

Отсюда получаются вклады

$$f_{k+1}(s,r)\mathrel{+}= \int_0^1 K_0(r,t)P(t)\,dt,$$

$$f_{k+1}(s+1,r)\mathrel{+}= \int_0^1 K_+(r,t)P(t)\,dt=\int_0^r (r-t)P(t)\,dt,$$

$$f_{k+1}(s-1,r)\mathrel{+}= \int_0^1 K_-(r,t)P(t)\,dt=\int_r^1 (t-r)P(t)\,dt.$$

Тождество

$$K_0(r,t)+K_+(r,t)+K_-(r,t)=1$$

немедленно показывает, что полная вероятность сохраняется.

Шаг 3: Переписать переходы через первообразные

Если

$$P(t)=\sum_{i=0}^{d} c_i t^i,$$

то введем две префиксные интегральные функции и два полных момента:

$$A(r)=\int_0^r P(t)\,dt,\qquad B(r)=\int_0^r tP(t)\,dt,$$

$$T=\int_0^1 P(t)\,dt,\qquad T_1=\int_0^1 tP(t)\,dt.$$

Тогда три вклада записываются так:

$$\int_0^r (r-t)P(t)\,dt = rA(r)-B(r),$$

$$\int_r^1 (t-r)P(t)\,dt = \bigl(T_1-B(r)\bigr)-r\bigl(T-A(r)\bigr),$$

$$\int_0^1 (1-|r-t|)P(t)\,dt = (1-r)A(r)+B(r)+(1+r)\bigl(T-A(r)\bigr)-\bigl(T_1-B(r)\bigr).$$

Именно эти полиномиальные комбинации строят реализации на C++, Python и Java.

Шаг 4: Почему достаточно полиномов

Начальная плотность имеет степень \(0\). Если степень \(P\) равна \(d\), то степень \(A\) не превышает \(d+1\), степень \(B\) не превышает \(d+2\), а формулы перехода используют только линейные комбинации \(A\), \(B\), \(rA\) и \(r(T-A)\). Значит, один шаг увеличивает степень не более чем на \(2\).

После \(k\) шагов каждая плотность имеет степень не более \(2k\), а после всех \(m=n-1\) шагов выполняется

$$\deg f_m(s,\cdot)\le 2m.$$

Именно поэтому можно хранить каждую плотность в фиксированном массиве коэффициентов длины \(2m+1\).

Шаг 5: Разобранный пример для \(n=3\)

Здесь \(m=2\). Если начать с \(f_0(0,r)=1\), то после одного перехода получаем

$$f_1(0,r)=\frac{1}{2}+r-r^2,$$

$$f_1(1,r)=\frac{r^2}{2},$$

$$f_1(-1,r)=\frac{(1-r)^2}{2}=\frac{1}{2}-r+\frac{r^2}{2}.$$

Интегрирование по \(r\in[0,1]\) дает массы

$$\int_0^1 f_1(0,r)\,dr=\frac{2}{3},\qquad \int_0^1 f_1(1,r)\,dr=\frac{1}{6},\qquad \int_0^1 f_1(-1,r)\,dr=\frac{1}{6},$$

и их сумма снова равна \(1\).

После второго перехода центральная плотность становится

$$f_2(0,r)=\frac{11}{24}+\frac{1}{2}r-\frac{1}{4}r^2-\frac{1}{2}r^3+\frac{1}{4}r^4.$$

Следовательно,

$$P(3)=\int_0^1 f_2(0,r)\,dr=\frac{11}{20}.$$

Реализации также совпадают с дополнительной проверкой

$$P(5)=\frac{15619}{36288}.$$

Шаг 6: Извлечь окончательную вероятность

После всех \(m=n-1\) шагов ответ равен

$$P(n)=\int_0^1 f_m(0,r)\,dr.$$

Дополнительно должно выполняться нормировочное тождество

$$\sum_{s=-m}^{m}\int_0^1 f_m(s,r)\,dr=1.$$

Это сильная внутренняя проверка правильности, и высокоточная реализация контролирует ее явно.

Как Работает Код

Реализации на C++, Python и Java хранят одну таблицу полиномиальных коэффициентов для текущего шага и вторую для следующего. Для каждого активного дискретного состояния они вычисляют две префиксные первообразные \(A\) и \(B\), а также полные интегралы \(T\) и \(T_1\). Поскольку интегрирование монома сводится к делению коэффициента на \(i+1\) или \(i+2\), обновления выполняются точно на уровне коэффициентов.

Затем эти полиномиальные части добавляются в три целевых состояния по формулам выше. После последнего шага реализация интегрирует полином, относящийся к состоянию \(0\), по отрезку \([0,1]\) и печатает результат. Версия на C++ дополнительно проверяет точные значения \(P(3)=11/20\), \(P(5)=15619/36288\) и нормировку; версии на Python и Java используют ту же рекуррентную схему для получения итогового десятичного значения.

Анализ Сложности

Имеется \(m=n-1\) этапов. На этапе \(k\) активно не более \(2k+1=O(n)\) дискретных состояний, и каждая плотность хранится с \(O(n)\) коэффициентами. Обработка одного состояния требует лишь постоянного числа проходов по этим коэффициентам. Поэтому один этап стоит \(O(n^2)\), а весь алгоритм - \(O(n^3)\).

Таблицы коэффициентов содержат \(O(n)\) состояний по \(O(n)\) коэффициентов в каждом, то есть требуют \(O(n^2)\) памяти. Используется десятичная арифметика повышенной точности, потому что точные рациональные значения быстро приобретают большие знаменатели.

Примечания и Ссылки

  1. Страница задачи: https://projecteuler.net/problem=807
  2. Плотность вероятности: Wikipedia - Probability density function
  3. Динамическое программирование: Wikipedia - Dynamic programming
  4. Непрерывное равномерное распределение: Wikipedia - Continuous uniform distribution
  5. Интеграл: Wikipedia - Integral

ملخص المسألة

المطلوب هو حساب احتمال \(P(n)\) المرتبط بعملية عشوائية لربط الحبال عبر \(n-1\) خطوة. لا تعتمد التطبيقات على المحاكاة من نوع Monte Carlo، بل تحتفظ بالتوزيع الدقيق للمعامل المستمر المتبقي وتحدّثه رمزيا من خطوة إلى أخرى، ولذلك تكون النتيجة نهائية وحتمية وعالية الدقة.

في حالة Project Euler الفعلية نحتاج إلى \(P(80)\). الفكرة الأساسية هي أن كل حالة وسيطة يمكن اختزالها إلى فهرس حالة متقطع ومعامل مستمر \(r \in [0,1]\). عندها تصبح المسألة كلها برنامجا ديناميكيا دقيقا يعمل على كثافات متعددة الحدود.

المنهج الرياضي

لنكتب \(m=n-1\). بعد \(k\) خطوات ربط، تمثل التطبيقات العملية بواسطة كثافات

$$f_k(s,r), \qquad -k \le s \le k,\quad 0 \le r \le 1,$$

حيث إن \(s\) هو الحالة المتقطعة و\(r\) هو المعامل المستمر المتبقي. والاحتمال المطلوب هو الكتلة الكلية للحالة النهائية \(s=0\) بعد إتمام \(m\) خطوات.

الخطوة 1: بناء التوزيع الابتدائي

في البداية توجد حالة متقطعة واحدة فقط، ولذلك تكون الكثافة الابتدائية منتظمة:

$$f_0(0,r)=1,\qquad f_0(s,r)=0 \text{ for } s\ne 0.$$

وبما أن \(\int_0^1 1\,dr=1\)، فهذا بالفعل توزيع احتمالي مُطبّع.

الخطوة 2: اشتقاق نوى الانتقال الثلاث

لنفترض أن لدينا كثافة مصدر \(P(t)\) في المتغير القديم \(t\). في الخطوة التالية يمكن للكتلة أن تبقى في الحالة نفسها، أو تنتقل إلى \(s+1\)، أو تنتقل إلى \(s-1\). النوى الثلاث المشفرة في التطبيقات هي

$$K_0(r,t)=1-|r-t|,$$

$$K_+(r,t)=\max(r-t,0),$$

$$K_-(r,t)=\max(t-r,0).$$

ومن ثم تكون المساهمات

$$f_{k+1}(s,r)\mathrel{+}= \int_0^1 K_0(r,t)P(t)\,dt,$$

$$f_{k+1}(s+1,r)\mathrel{+}= \int_0^1 K_+(r,t)P(t)\,dt=\int_0^r (r-t)P(t)\,dt,$$

$$f_{k+1}(s-1,r)\mathrel{+}= \int_0^1 K_-(r,t)P(t)\,dt=\int_r^1 (t-r)P(t)\,dt.$$

والهوية

$$K_0(r,t)+K_+(r,t)+K_-(r,t)=1$$

توضح مباشرة أن الكتلة الاحتمالية الكلية محفوظة.

الخطوة 3: إعادة كتابة الانتقالات باستعمال التكاملات التراكمية

إذا كانت

$$P(t)=\sum_{i=0}^{d} c_i t^i,$$

فنعرّف تكاملين تراكميين وعزمين كليين:

$$A(r)=\int_0^r P(t)\,dt,\qquad B(r)=\int_0^r tP(t)\,dt,$$

$$T=\int_0^1 P(t)\,dt,\qquad T_1=\int_0^1 tP(t)\,dt.$$

عندئذ تصبح المساهمات الثلاث

$$\int_0^r (r-t)P(t)\,dt = rA(r)-B(r),$$

$$\int_r^1 (t-r)P(t)\,dt = \bigl(T_1-B(r)\bigr)-r\bigl(T-A(r)\bigr),$$

$$\int_0^1 (1-|r-t|)P(t)\,dt = (1-r)A(r)+B(r)+(1+r)\bigl(T-A(r)\bigr)-\bigl(T_1-B(r)\bigr).$$

هذه هي بالضبط التركيبات متعددة الحدود التي تبنيها تطبيقات ++C وPython وJava.

الخطوة 4: لماذا تبقى الكثافات متعددة حدود

درجة الكثافة الابتدائية هي \(0\). وإذا كانت درجة \(P\) تساوي \(d\)، فإن درجة \(A\) لا تتجاوز \(d+1\)، ودرجة \(B\) لا تتجاوز \(d+2\)، كما أن صيغ الانتقال لا تستخدم إلا تراكيب خطية من \(A\) و\(B\) و\(rA\) و\(r(T-A)\). لذلك فإن كل خطوة تزيد الدرجة بمقدار لا يتجاوز \(2\).

وبعد \(k\) خطوات تصبح درجة كل كثافة على الأكثر \(2k\)، وبعد جميع الخطوات \(m=n-1\) نحصل على

$$\deg f_m(s,\cdot)\le 2m.$$

ولهذا السبب تستطيع التطبيقات تخزين كل كثافة في مصفوفة معاملات ثابتة الطول مقدارها \(2m+1\).

الخطوة 5: مثال محلول عندما \(n=3\)

هنا \(m=2\). إذا بدأنا من \(f_0(0,r)=1\)، فإن خطوة انتقال واحدة تعطي

$$f_1(0,r)=\frac{1}{2}+r-r^2,$$

$$f_1(1,r)=\frac{r^2}{2},$$

$$f_1(-1,r)=\frac{(1-r)^2}{2}=\frac{1}{2}-r+\frac{r^2}{2}.$$

وبتكامل هذه الكثافات على الفترة \([0,1]\) نحصل على الكتل

$$\int_0^1 f_1(0,r)\,dr=\frac{2}{3},\qquad \int_0^1 f_1(1,r)\,dr=\frac{1}{6},\qquad \int_0^1 f_1(-1,r)\,dr=\frac{1}{6},$$

ومجموعها يساوي \(1\) كما يجب.

وبعد تطبيق الانتقال مرة ثانية تصبح كثافة الحالة المركزية

$$f_2(0,r)=\frac{11}{24}+\frac{1}{2}r-\frac{1}{4}r^2-\frac{1}{2}r^3+\frac{1}{4}r^4.$$

إذن

$$P(3)=\int_0^1 f_2(0,r)\,dr=\frac{11}{20}.$$

كما تتفق التطبيقات أيضا على قيمة فحص إضافية هي

$$P(5)=\frac{15619}{36288}.$$

الخطوة 6: استخراج الاحتمال النهائي

بعد إكمال جميع الخطوات \(m=n-1\)، يكون الجواب ببساطة

$$P(n)=\int_0^1 f_m(0,r)\,dr.$$

ويجب كذلك أن تتحقق هوية التطبيع

$$\sum_{s=-m}^{m}\int_0^1 f_m(s,r)\,dr=1.$$

وهذا فحص داخلي قوي، وتقوم النسخة ذات الدقة العالية بالتحقق منه صراحة.

كيف يعمل الكود

تحتفظ تطبيقات ++C وPython وJava بجدول معاملات متعدد الحدود للخطوة الحالية وجدول آخر للخطوة التالية. ولكل حالة متقطعة غير صفرية تحسب التكاملين التراكميين \(A\) و\(B\)، ثم التكاملين الكليين \(T\) و\(T_1\). وبما أن تكامل حد أحادي لا يحتاج إلا إلى قسمة معامله على \(i+1\) أو \(i+2\)، فإن التحديثات كلها تتم بدقة كاملة على مستوى المعاملات.

ثم تُضاف هذه الأجزاء متعددة الحدود إلى الحالات الهدف الثلاث حسب الصيغ السابقة. وبعد الخطوة الأخيرة يقوم التنفيذ بتكامل متعدد الحدود المرتبط بالحالة \(0\) على \([0,1]\) ويطبع الناتج. كما أن نسخة ++C تتحقق أيضا من القيم الدقيقة \(P(3)=11/20\) و\(P(5)=15619/36288\) ومن شرط التطبيع، بينما تستخدم نسختا Python وJava العلاقة نفسها لإنتاج القيمة العشرية النهائية.

تحليل التعقيد

لدينا \(m=n-1\) مراحل. في المرحلة \(k\) يوجد على الأكثر \(2k+1=O(n)\) حالات متقطعة فعالة، وكل كثافة تخزن بواسطة \(O(n)\) معاملات. ومعالجة حالة واحدة تحتاج فقط إلى عدد ثابت من المرور على هذه المعاملات، ولذلك تكون كلفة المرحلة الواحدة \(O(n^2)\) وتكون الكلفة الكلية \(O(n^3)\).

أما الذاكرة فتبلغ \(O(n^2)\) لأن جداول المعاملات تخزن \(O(n)\) حالات، ولكل حالة \(O(n)\) معاملات. وتُستخدم الحسابات العشرية عالية الدقة لأن القيم الكسرية الدقيقة تكتسب مقامات كبيرة بسرعة.

حواشٍ ومراجع

  1. صفحة المسألة: https://projecteuler.net/problem=807
  2. دالة كثافة الاحتمال: Wikipedia - Probability density function
  3. البرمجة الديناميكية: Wikipedia - Dynamic programming
  4. التوزيع المنتظم المستمر: Wikipedia - Continuous uniform distribution
  5. التكامل: Wikipedia - Integral