Detailed mathematical approach

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