Problem 807: Loops of Ropes

Project Euler Problem 807 Solution

EulerSolve provides an optimized solution for Project Euler Problem 807, Loops of Ropes, with C++, Python, Java, and a step-by-step mathematical explanation.

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$....

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

$$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

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

Problem 807 source code

C++

#include <cassert>
#include <iomanip>
#include <iostream>
#include <vector>
#include <boost/multiprecision/cpp_dec_float.hpp>

using namespace std;

namespace {

using BigFloat = boost::multiprecision::cpp_dec_float_50;
using Poly = vector<BigFloat>;

struct Result {
    BigFloat prob_zero;
    BigFloat total_mass;
};

static Result compute_probability(int n) {
    assert(n > 1);
    const int m = n - 1;
    const int max_deg = 2 * m;
    const int offset = m;

    auto make_zero_poly = [&]() { return Poly(max_deg + 1, BigFloat(0)); };

    vector<Poly> curr(2 * m + 1, make_zero_poly());
    curr[offset][0] = BigFloat(1);

    for (int step = 1; step <= m; ++step) {
        vector<Poly> next(2 * m + 1, make_zero_poly());

        for (int s = -(step - 1); s <= (step - 1); ++s) {
            const Poly& P = curr[offset + s];

            BigFloat total = 0;
            BigFloat total_r = 0;
            for (int i = 0; i <= max_deg; ++i) {
                if (P[i] == 0) continue;
                total += P[i] / BigFloat(i + 1);
                total_r += P[i] / BigFloat(i + 2);
            }

            Poly A(max_deg + 1, BigFloat(0));
            Poly B(max_deg + 1, BigFloat(0));
            for (int i = 0; i <= max_deg; ++i) {
                if (P[i] == 0) continue;
                if (i + 1 <= max_deg) A[i + 1] = P[i] / BigFloat(i + 1);
                if (i + 2 <= max_deg) B[i + 2] = P[i] / BigFloat(i + 2);
            }

            Poly TminusA(max_deg + 1, BigFloat(0));
            TminusA[0] = total;
            for (int i = 1; i <= max_deg; ++i) TminusA[i] = -A[i];

            Poly& same = next[offset + s];
            Poly& plus = next[offset + s + 1];
            Poly& minus = next[offset + s - 1];

            for (int i = 0; i <= max_deg; ++i) {
                BigFloat i1 = A[i] + B[i];
                if (i > 0) i1 -= A[i - 1];

                BigFloat i2 = TminusA[i] + B[i];
                if (i > 0) i2 += TminusA[i - 1];
                if (i == 0) i2 -= total_r;

                same[i] += i1 + i2;

                BigFloat jplus = -B[i];
                if (i > 0) jplus += A[i - 1];
                plus[i] += jplus;

                BigFloat jminus = -B[i];
                if (i == 0) jminus += total_r;
                if (i > 0) jminus -= TminusA[i - 1];
                minus[i] += jminus;
            }
        }

        curr.swap(next);
    }

    auto integrate_0_1 = [&](const Poly& P) {
        BigFloat res = 0;
        for (int i = 0; i <= max_deg; ++i) {
            if (P[i] == 0) continue;
            res += P[i] / BigFloat(i + 1);
        }
        return res;
    };

    Result out;
    out.prob_zero = integrate_0_1(curr[offset]);

    BigFloat total_mass = 0;
    for (int s = -m; s <= m; ++s) {
        total_mass += integrate_0_1(curr[offset + s]);
    }
    out.total_mass = total_mass;
    return out;
}

static void validate() {
    using boost::multiprecision::abs;
    const BigFloat eps = BigFloat("1e-22");

    const Result r3 = compute_probability(3);
    const Result r5 = compute_probability(5);

    const BigFloat expected3 = BigFloat(11) / BigFloat(20);
    const BigFloat expected5 = BigFloat(15619) / BigFloat(36288);

    if (abs(r3.prob_zero - expected3) > eps) {
        cerr << "Validation failed for P(3)." << '\n';
        std::exit(1);
    }
    if (abs(r5.prob_zero - expected5) > eps) {
        cerr << "Validation failed for P(5)." << '\n';
        std::exit(1);
    }
    if (abs(r5.total_mass - BigFloat(1)) > eps) {
        cerr << "Validation failed for normalization." << '\n';
        std::exit(1);
    }
}

}

int main() {
    using boost::multiprecision::abs;
    validate();

    const Result r80 = compute_probability(80);
    const BigFloat eps = BigFloat("1e-22");
    if (abs(r80.total_mass - BigFloat(1)) > eps) {
        cerr << "Normalization check failed for n=80." << '\n';
        return 1;
    }

    cout << fixed << setprecision(10) << r80.prob_zero << '\n';
    return 0;
}

Python

import decimal

def compute_probability(n):
    decimal.getcontext().prec = 60
    BigFloat = decimal.Decimal
    
    m = n - 1
    max_deg = 2 * m
    offset = m
    
    def make_zero_poly():
        return [BigFloat(0) for _ in range(max_deg + 1)]
        
    curr = [make_zero_poly() for _ in range(2 * m + 1)]
    curr[offset][0] = BigFloat(1)
    
    for step in range(1, m + 1):
        next_poly = [make_zero_poly() for _ in range(2 * m + 1)]
        
        for s in range(-(step - 1), step):
            P = curr[offset + s]
            
            total = BigFloat(0)
            total_r = BigFloat(0)
            
            for i in range(max_deg + 1):
                if P[i] == 0: continue
                total += P[i] / BigFloat(i + 1)
                total_r += P[i] / BigFloat(i + 2)
                
            A = make_zero_poly()
            B = make_zero_poly()
            
            for i in range(max_deg + 1):
                if P[i] == 0: continue
                if i + 1 <= max_deg: A[i + 1] = P[i] / BigFloat(i + 1)
                if i + 2 <= max_deg: B[i + 2] = P[i] / BigFloat(i + 2)
                
            TminusA = make_zero_poly()
            TminusA[0] = total
            for i in range(1, max_deg + 1):
                TminusA[i] = -A[i]
                
            same = next_poly[offset + s]
            plus = next_poly[offset + s + 1]
            minus = next_poly[offset + s - 1]
            
            for i in range(max_deg + 1):
                i1 = A[i] + B[i]
                if i > 0: i1 -= A[i - 1]
                
                i2 = TminusA[i] + B[i]
                if i > 0: i2 += TminusA[i - 1]
                if i == 0: i2 -= total_r
                
                same[i] += i1 + i2
                
                jplus = -B[i]
                if i > 0: jplus += A[i - 1]
                plus[i] += jplus
                
                jminus = -B[i]
                if i == 0: jminus += total_r
                if i > 0: jminus -= TminusA[i - 1]
                minus[i] += jminus
                
        curr = next_poly
        
    def integrate_0_1(P):
        res = BigFloat(0)
        for i in range(max_deg + 1):
            if P[i] == 0: continue
            res += P[i] / BigFloat(i + 1)
        return res
        
    prob_zero = integrate_0_1(curr[offset])
    return prob_zero

def solve():
    ans = compute_probability(80)
    # The precision required is 10 decimal places
    return "{:.10f}".format(ans)

if __name__ == "__main__":
    print(solve())

Java

import java.math.BigDecimal;
import java.math.MathContext;

public class Euler807 {

    static class Result {
        BigDecimal probZero;
        BigDecimal totalMass;
    }

    static Result computeProbability(int n) {
        MathContext mc = new MathContext(60);

        int m = n - 1;
        int maxDeg = 2 * m;
        int offset = m;

        BigDecimal[][] curr = new BigDecimal[2 * m + 1][maxDeg + 1];
        for (int i = 0; i < 2 * m + 1; i++) {
            for (int j = 0; j <= maxDeg; j++) {
                curr[i][j] = BigDecimal.ZERO;
            }
        }
        curr[offset][0] = BigDecimal.ONE;

        for (int step = 1; step <= m; ++step) {
            BigDecimal[][] next = new BigDecimal[2 * m + 1][maxDeg + 1];
            for (int i = 0; i < 2 * m + 1; i++) {
                for (int j = 0; j <= maxDeg; j++) {
                    next[i][j] = BigDecimal.ZERO;
                }
            }

            for (int s = -(step - 1); s <= (step - 1); ++s) {
                BigDecimal[] P = curr[offset + s];

                BigDecimal total = BigDecimal.ZERO;
                BigDecimal totalR = BigDecimal.ZERO;

                for (int i = 0; i <= maxDeg; ++i) {
                    if (P[i].signum() == 0)
                        continue;
                    total = total.add(P[i].divide(new BigDecimal(i + 1), mc), mc);
                    totalR = totalR.add(P[i].divide(new BigDecimal(i + 2), mc), mc);
                }

                BigDecimal[] A = new BigDecimal[maxDeg + 1];
                BigDecimal[] B = new BigDecimal[maxDeg + 1];
                for (int i = 0; i <= maxDeg; i++) {
                    A[i] = BigDecimal.ZERO;
                    B[i] = BigDecimal.ZERO;
                }

                for (int i = 0; i <= maxDeg; ++i) {
                    if (P[i].signum() == 0)
                        continue;
                    if (i + 1 <= maxDeg)
                        A[i + 1] = P[i].divide(new BigDecimal(i + 1), mc);
                    if (i + 2 <= maxDeg)
                        B[i + 2] = P[i].divide(new BigDecimal(i + 2), mc);
                }

                BigDecimal[] TminusA = new BigDecimal[maxDeg + 1];
                TminusA[0] = total;
                for (int i = 1; i <= maxDeg; ++i)
                    TminusA[i] = A[i].negate();

                BigDecimal[] same = next[offset + s];
                BigDecimal[] plus = next[offset + s + 1];
                BigDecimal[] minus = next[offset + s - 1];

                for (int i = 0; i <= maxDeg; ++i) {
                    BigDecimal i1 = A[i].add(B[i]);
                    if (i > 0)
                        i1 = i1.subtract(A[i - 1]);

                    BigDecimal i2 = TminusA[i].add(B[i]);
                    if (i > 0)
                        i2 = i2.add(TminusA[i - 1]);
                    if (i == 0)
                        i2 = i2.subtract(totalR);

                    same[i] = same[i].add(i1).add(i2);

                    BigDecimal jplus = B[i].negate();
                    if (i > 0)
                        jplus = jplus.add(A[i - 1]);
                    plus[i] = plus[i].add(jplus);

                    BigDecimal jminus = B[i].negate();
                    if (i == 0)
                        jminus = jminus.add(totalR);
                    if (i > 0)
                        jminus = jminus.subtract(TminusA[i - 1]);
                    minus[i] = minus[i].add(jminus);
                }
            }

            curr = next;
        }

        Result out = new Result();
        BigDecimal res01 = BigDecimal.ZERO;
        for (int i = 0; i <= maxDeg; ++i) {
            BigDecimal val = curr[offset][i];
            if (val.signum() == 0)
                continue;
            res01 = res01.add(val.divide(new BigDecimal(i + 1), mc), mc);
        }
        out.probZero = res01;

        return out;
    }

    public static String solve() {
        Result r80 = computeProbability(80);
        return String.format(java.util.Locale.US, "%.10f", r80.probZero);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}