Problem 906: A Collective Decision

Project Euler Problem 906 Solution

EulerSolve provides an optimized solution for Project Euler Problem 906, A Collective Decision, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Problem 906 studies two independent random preference orders on the same set of $n$ alternatives and asks for a probability $P(n)$; the final evaluation uses $n=20000$. The implementations show that the whole computation can be organized around one distinguished alternative $x$. If $x$ is in position $a+1$ in the first order and in position $b+1$ in the second, then among the other $N=n-1$ alternatives there are $a$ above $x$ in the first order and $b$ above $x$ in the second. The remaining alternatives split into three blocks: those above $x$ only in the first order, those above $x$ only in the second order, and those below $x$ in both orders. The event counted by the implementations requires the two upper blocks to be disjoint, and then it adds an auxiliary contribution coming from the common lower block. Averaging that contribution over the $n^2$ equally likely rank pairs of $x$ yields the desired probability. Mathematical Approach Let $N=n-1$. For a fixed alternative $x$, the two permutations induce a pair of rank coordinates. The solution turns the original probability question into a clean double sum over those coordinates. Fix one alternative and its two rank coordinates Write $y\succ_1 x$ when the first preference order places $y$ above $x$, and $y\succ_2 x$ for the second order....

Detailed mathematical approach

Problem Summary

Problem 906 studies two independent random preference orders on the same set of $n$ alternatives and asks for a probability $P(n)$; the final evaluation uses $n=20000$. The implementations show that the whole computation can be organized around one distinguished alternative $x$. If $x$ is in position $a+1$ in the first order and in position $b+1$ in the second, then among the other $N=n-1$ alternatives there are $a$ above $x$ in the first order and $b$ above $x$ in the second.

The remaining alternatives split into three blocks: those above $x$ only in the first order, those above $x$ only in the second order, and those below $x$ in both orders. The event counted by the implementations requires the two upper blocks to be disjoint, and then it adds an auxiliary contribution coming from the common lower block. Averaging that contribution over the $n^2$ equally likely rank pairs of $x$ yields the desired probability.

Mathematical Approach

Let $N=n-1$. For a fixed alternative $x$, the two permutations induce a pair of rank coordinates. The solution turns the original probability question into a clean double sum over those coordinates.

Fix one alternative and its two rank coordinates

Write $y\succ_1 x$ when the first preference order places $y$ above $x$, and $y\succ_2 x$ for the second order. Define

$$A=\{y:y\succ_1 x\},\qquad B=\{y:y\succ_2 x\}.$$

If $|A|=a$ and $|B|=b$, then $x$ has rank $a+1$ in the first order and rank $b+1$ in the second. Because the two orders are independent and each rank of $x$ is uniform on $\{1,\dots,n\}$, every pair $(a,b)\in\{0,\dots,N\}^2$ occurs with probability $1/n^2$.

Why the upper sets must be disjoint

The implemented formula contributes only when no alternative lies above $x$ in both orders. In set language that condition is

$$A\cap B=\varnothing.$$

For fixed sizes $a$ and $b$, the set $B$ is a uniformly random $b$-subset of the $N$ alternatives other than $x$. Exactly $\binom{N-a}{b}$ such subsets avoid the $a$ elements already placed in $A$. Therefore

$$q(a,b)=\Pr(A\cap B=\varnothing\mid |A|=a,\ |B|=b)=\frac{\binom{N-a}{b}}{\binom{N}{b}}.$$

This immediately explains the triangular summation region in the code: if $a+b>N$, then disjoint sets of those sizes cannot exist, so the contribution is zero. It also explains the multiplicative update used inside the inner loop:

$$q(a,0)=1,\qquad q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

The unanimous lower block

Once $A$ and $B$ are disjoint, every remaining alternative is below $x$ in both orders. Denote that common lower block by

$$C=\{y:x\succ_1 y\ \text{and}\ x\succ_2 y\},\qquad |C|=m=N-a-b.$$

The implementations attach to this block the auxiliary quantity

$$R_m=\sum_{c=0}^{m}\prod_{t=0}^{c-1}\frac{m-t}{N-t}=\sum_{c=0}^{m}\frac{(m)_c}{(N)_c}=\sum_{c=0}^{m}\frac{\binom{m}{c}}{\binom{N}{c}}.$$

The term for a fixed $c$ has a direct counting interpretation: $\binom{m}{c}/\binom{N}{c}$ is the probability that a uniformly chosen $c$-subset of the $N$ alternatives other than $x$ lies entirely inside the unanimous lower block $C$. The total $R_m$ adds those probabilities for every feasible subset size $c$.

Collapsing the auxiliary sum

The direct sum for $R_m$ is already good enough for computation, but the mathematics becomes clearer after simplifying it. Start from

$$\frac{\binom{m}{c}}{\binom{N}{c}}=\frac{1}{\binom{N}{m}}\binom{N-c}{m-c}.$$

Substituting this into the previous formula gives

$$R_m=\frac{1}{\binom{N}{m}}\sum_{c=0}^{m}\binom{N-c}{m-c}.$$

Now set $j=m-c$. The sum becomes

$$\sum_{j=0}^{m}\binom{N-m+j}{j},$$

which is a standard hockey-stick sum, so

$$R_m=\frac{\binom{N+1}{m}}{\binom{N}{m}}=\frac{N+1}{N+1-m}.$$

Since $m=N-a-b$, this can be rewritten in the especially convenient form

$$R_{N-a-b}=\frac{n}{a+b+1}.$$

The final probability formula

Putting the two pieces together, the probability computed by the C++, Python, and Java implementations is

$$\boxed{P(n)=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a} q(a,b)\,R_{N-a-b}=\frac{1}{n^2}\sum_{a=0}^{N}\sum_{b=0}^{N-a}\frac{\binom{N-a}{b}}{\binom{N}{b}}\cdot\frac{n}{a+b+1}.}$$

Everything in the program is a fast way to evaluate this triangular double sum without ever forming large factorials or binomial coefficients explicitly.

Worked Example: $n=3$

Take $n=3$, so $N=2$. Then

$$R_0=1,\qquad R_1=1+\frac12=\frac32,\qquad R_2=1+1+1=3.$$

The admissible pairs $(a,b)$ satisfy $a+b\le 2$. Their contributions are

$$\begin{aligned} (0,0)&:\ q=1,\ m=2,\ qR_m=3,\\ (0,1)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (0,2)&:\ q=1,\ m=0,\ qR_m=1,\\ (1,0)&:\ q=1,\ m=1,\ qR_m=\frac32,\\ (1,1)&:\ q=\frac12,\ m=0,\ qR_m=\frac12,\\ (2,0)&:\ q=1,\ m=0,\ qR_m=1. \end{aligned}$$

The total is $17/2$, and dividing by $n^2=9$ gives

$$P(3)=\frac{17/2}{9}=\frac{17}{18}.$$

This small case shows all the moving parts of the full solution: the disjointness factor, the lower-block sum, and the final average over rank pairs.

How the Code Works

Precomputation of reciprocal factors and lower-block values

The implementations first build a table of reciprocals $1,2,\dots,N$. That allows every ratio in the double sum to be updated by multiplication instead of by repeated division. They also precompute all values of $R_m$ for $m=0,\dots,N$. For a fixed $m$, the term with $c=0$ is $1$, and the next term is obtained from the previous one by multiplying by $(m-c+1)/(N-c+1)$. Summing those terms produces the required lower-block table.

The triangular sweep over $(a,b)$

After the $R_m$ table is available, the main work is the triangular loop over all $a$ and $b$ with $0\le a\le N$ and $0\le b\le N-a$. For each fixed $a$, the implementation starts from $q(a,0)=1$ and updates the disjointness factor with

$$q(a,b+1)=q(a,b)\frac{N-a-b}{N-b}.$$

At each lattice point $(a,b)$, it looks up $R_{N-a-b}$, adds $q(a,b)R_{N-a-b}$ to the running total, and finally divides the accumulated sum by $n^2$.

Numerical accuracy and compiled-language optimizations

The calculation is purely floating-point. The compiled implementations therefore use compensated summation to reduce accumulation error, and the C++ version also keeps the main total in extended precision. The compiled versions split the outer $a$-range into blocks and process those blocks in parallel, while the Python implementation follows the same mathematics serially.

Complexity Analysis

Let $N=n-1$. Building the reciprocal table costs $O(N)$ time and $O(N)$ memory. Precomputing all $R_m$ values with the direct inner products used in the implementations costs $O(N^2)$ arithmetic operations. The main triangular sum over $(a,b)$ is also $O(N^2)$.

Therefore the overall running time is $O(N^2)$ and the memory usage is $O(N)$. Parallel execution reduces wall-clock time for the compiled versions, but it does not change the asymptotic order.

Footnotes and References

Problem page: Project Euler 906
Total order: Wikipedia - Total order
Binomial coefficient: Wikipedia - Binomial coefficient
Falling and rising factorials: Wikipedia - Falling and rising factorials
Hypergeometric distribution: Wikipedia - Hypergeometric distribution
Hockey-stick identity: Wikipedia - Hockey-stick identity
Kahan summation algorithm: Wikipedia - Kahan summation algorithm

Problem 906 source code

C++

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <thread>
#include <vector>

struct KahanSum {
    long double sum = 0.0L;
    long double c = 0.0L;

    void add(long double value) {
        long double y = value - c;
        long double t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }
};

std::vector<long double> build_reciprocals(int n) {
    std::vector<long double> inv(n + 1, 0.0L);
    for (int i = 1; i <= n; ++i) {
        inv[i] = 1.0L / static_cast<long double>(i);
    }
    return inv;
}

std::vector<long double> build_g_values(int n, const std::vector<long double>& inv) {
    const int N = n - 1;
    std::vector<long double> G(N + 1, 0.0L);
    G[0] = 1.0L;
    for (int m = 1; m <= N; ++m) {
        long double term = 1.0L;
        long double sum = 1.0L;
        for (int c = 1; c <= m; ++c) {
            term *= static_cast<long double>(m - c + 1) * inv[N - c + 1];
            sum += term;
        }
        G[m] = sum;
    }
    return G;
}

long double compute_sum_sequential(int n,
                                   const std::vector<long double>& G,
                                   const std::vector<long double>& inv) {
    const int N = n - 1;
    KahanSum total;
    for (int a = 0; a <= N; ++a) {
        long double f = 1.0L;
        const int max_b = N - a;
        for (int b = 0; b <= max_b; ++b) {
            const int m = N - a - b;
            total.add(f * G[m]);
            if (b < max_b) {
                f *= static_cast<long double>(N - a - b) * inv[N - b];
            }
        }
    }
    return total.sum;
}

long double compute_sum_parallel(int n,
                                 const std::vector<long double>& G,
                                 const std::vector<long double>& inv,
                                 int threads) {
    const int N = n - 1;
    if (threads <= 1 || N < 200) {
        return compute_sum_sequential(n, G, inv);
    }
    const int max_threads = std::min(threads, N + 1);
    const int block = (N + 1 + max_threads - 1) / max_threads;
    std::vector<long double> partial(max_threads, 0.0L);
    std::vector<std::thread> pool;
    pool.reserve(max_threads);
    for (int t = 0; t < max_threads; ++t) {
        const int start = t * block;
        const int end = std::min(N + 1, start + block);
        if (start >= end) {
            continue;
        }
        pool.emplace_back([&, start, end, t]() {
            KahanSum local;
            for (int a = start; a < end; ++a) {
                long double f = 1.0L;
                const int max_b = N - a;
                for (int b = 0; b <= max_b; ++b) {
                    const int m = N - a - b;
                    local.add(f * G[m]);
                    if (b < max_b) {
                        f *= static_cast<long double>(N - a - b) * inv[N - b];
                    }
                }
            }
            partial[t] = local.sum;
        });
    }
    for (auto& th : pool) {
        th.join();
    }
    KahanSum total;
    for (long double value : partial) {
        total.add(value);
    }
    return total.sum;
}

long double compute_probability(int n, int threads) {
    if (n <= 1) {
        return 1.0L;
    }
    const int N = n - 1;
    const std::vector<long double> inv = build_reciprocals(N);
    const std::vector<long double> G = build_g_values(n, inv);
    const long double sum = compute_sum_parallel(n, G, inv, threads);
    return sum / (static_cast<long double>(n) * static_cast<long double>(n));
}

bool approx_equal(long double a, long double b, long double tol) {
    return std::fabs(a - b) <= tol;
}

bool run_validation() {
    bool ok = true;
    const long double p3 = compute_probability(3, 1);
    const long double p3_expected = 17.0L / 18.0L;
    if (!approx_equal(p3, p3_expected, 1e-12L)) {
        std::cerr << "Validation failed for n=3: got " << std::setprecision(15)
                  << p3 << " expected " << p3_expected
                  << "\n";
        ok = false;
    }
    const long double p10 = compute_probability(10, 1);
    const long double p10_expected = 0.6760292265L;
    if (!approx_equal(p10, p10_expected, 1e-9L)) {
        std::cerr << "Validation failed for n=10: got " << std::setprecision(15)
                  << p10 << " expected " << p10_expected
                  << "\n";
        ok = false;
    }
    if (ok) {
        std::cerr << "Validation checkpoints passed (n=3, n=10).\n";
    }
    return ok;
}

int main(int argc, char** argv) {
    int n = 20000;
    int threads = static_cast<int>(std::thread::hardware_concurrency());
    if (threads <= 0) {
        threads = 1;
    }
    if (argc >= 2) {
        n = std::stoi(argv[1]);
    }
    if (argc >= 3) {
        threads = std::max(1, std::stoi(argv[2]));
    }

    const bool ok = run_validation();
    if (!ok) {
        return 1;
    }

    const long double result = compute_probability(n, threads);
    if (result < -1e-12L || result > 1.0L + 1e-12L) {
        std::cerr << "Probability out of range: " << std::setprecision(15)
                  << result << "\n";
    }
    std::cout << std::fixed << std::setprecision(10) << result << "\n";
    return 0;
}

Python

def solve():
    n = 20000
    N = n - 1
    inv = [0.0]*(N+1)
    for i in range(1, N+1): inv[i] = 1.0/i

    G = [0.0]*(N+1); G[0] = 1.0
    for m in range(1, N+1):
        term = 1.0; s = 1.0
        for c in range(1, m+1):
            term *= (m-c+1)*inv[N-c+1]; s += term
        G[m] = s

    total = 0.0
    for a in range(N+1):
        f = 1.0; mb = N - a
        for b in range(mb+1):
            m = N - a - b; total += f * G[m]
            if b < mb: f *= (N-a-b)*inv[N-b]

    return f'{total / (n*n):.10f}'

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

Java

public class Euler906 {

    static class KahanSum {
        double sum = 0.0;
        double c = 0.0;

        void add(double value) {
            double y = value - c;
            double t = sum + y;
            c = (t - sum) - y;
            sum = t;
        }
    }

    static double[] buildReciprocals(int n) {
        double[] inv = new double[n + 1];
        for (int i = 1; i <= n; ++i) {
            inv[i] = 1.0 / i;
        }
        return inv;
    }

    static double[] buildGValues(int n, double[] inv) {
        int N = n - 1;
        double[] G = new double[N + 1];
        G[0] = 1.0;
        for (int m = 1; m <= N; ++m) {
            double term = 1.0;
            double sum = 1.0;
            for (int c = 1; c <= m; ++c) {
                term *= (m - c + 1) * inv[N - c + 1];
                sum += term;
            }
            G[m] = sum;
        }
        return G;
    }

    public static String solve() {
        int n = 20000;
        if (n <= 1)
            return "1.0000000000";

        int N = n - 1;
        double[] inv = buildReciprocals(N);
        double[] G = buildGValues(n, inv);

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
        double[] partial = new double[threads];
        Thread[] pool = new Thread[threads];

        int maxThreads = Math.min(threads, N + 1);
        int block = (N + 1 + maxThreads - 1) / maxThreads;

        for (int t = 0; t < maxThreads; ++t) {
            final int tIdx = t;
            final int start = t * block;
            final int end = Math.min(N + 1, start + block);

            if (start >= end)
                continue;

            pool[t] = new Thread(() -> {
                KahanSum local = new KahanSum();
                for (int a = start; a < end; ++a) {
                    double f = 1.0;
                    int maxB = N - a;
                    for (int b = 0; b <= maxB; ++b) {
                        int m = N - a - b;
                        local.add(f * G[m]);
                        if (b < maxB) {
                            f *= (N - a - b) * inv[N - b];
                        }
                    }
                }
                partial[tIdx] = local.sum;
            });
            pool[t].start();
        }

        try {
            for (int t = 0; t < maxThreads; t++) {
                if (pool[t] != null)
                    pool[t].join();
            }
        } catch (InterruptedException e) {
            e.printStackTrace();
        }

        KahanSum total = new KahanSum();
        for (double value : partial) {
            total.add(value);
        }

        double result = total.sum / ((double) n * n);
        return String.format(java.util.Locale.US, "%.10f", result);
    }

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