Project Euler

A Long Chess Match

Two players A and B play an infinitely long chess match. In each game, A wins with probability p_A, B wins with probability p_B, and the game is drawn with probability 1 - p_A - p_B. After each gam...

Source sync Apr 19, 2026

Problem #0661

Level Level 31

Solved By 271

Languages C++, Python

Answer 646231.2177

Length 335 words

probabilityanalytic_mathgame_theory

Two friends $A$ and $B$ are great fans of Chess. They both enjoy playing the game, but after each game the player who lost the game would like to continue (to get back at the other player) and the player who won would prefer to stop (to finish on a high).

So they come up with a plan. After every game, they would toss a (biased) coin with probability $p$ of Heads (and hence probability $1-p$ of Tails). If they get Tails, they will continue with the next game. Otherwise they end the match. Also, after every game the players make a note of who is leading in the match.

Let $p_A$ denote the probability of $A$ winning a game and $p_B$ the probability of $B$ winning a game. Accordingly $1-p_A-p_B$ is the probability that a game ends in a draw. Let $\mathbb {E}_A(p_A,p_B,p)$ denote the expected number of times $A$ was leading in the match.

For example, $\mathbb {E}_A(0.25,0.25,0.5)\approx 0.585786$ and $\mathbb {E}_A(0.47,0.48,0.001)\approx 377.471736$, both rounded to six places after the decimal point.

Let $\displaystyle H(n)=\sum _{k=3}^n \mathbb {E}_A\left (\frac 1 {\sqrt {k+3}},\frac 1 {\sqrt {k+3}}+\frac 1 {k^2},\frac 1 {k^3}\right )$

For example $H(3) \approx 6.8345$, rounded to 4 digits after the decimal point.

Find $H(50)$, rounded to 4 digits after the decimal point.

Problem 661: A Long Chess Match

Mathematical Foundation

Let $d_t = s_A(t) - s_B(t)$ denote the score difference after game $t$ , with $d_0 = 0$ . The sequence $\{d_t\}$ is a random walk on $\mathbb{Z}$ with step distribution:

d_{t+1} - d_t = \begin{cases} +1 & \text{with probability } p_A, \\ -1 & \text{with probability } p_B, \\ \phantom{+}0 & \text{with probability } 1 - p_A - p_B. \end{cases}

Player $A$ leads at time $t$ if and only if $d_t > 0$ . The match survives to time $t$ with probability $(1-p)^t$ .

Theorem 1 (Closed-Form via Wiener—Hopf Factorization). Let $p_A, p_B > 0$ with $p_A + p_B \leq 1$ and $0 < p < 1$ . Define $r_+$ as the smaller root of $p_B r^2 - r + p_A = 0$ :

r_+ = \frac{1 - \sqrt{1 - 4p_A p_B}}{2p_B}.

Then with $z = 1 - p$ :

\mathbb{E}_A(p_A, p_B, p) = \frac{r_+ z}{(1 - z)(1 - r_+ z)}.

Proof. By linearity of expectation and the geometric killing:

\mathbb{E}_A = \sum_{t=0}^{\infty} (1-p)^t \, \Pr[d_t > 0] = \sum_{t=0}^{\infty} z^t \, \Pr[d_t > 0].

The characteristic function of the step distribution is $\phi(\theta) = p_A e^{i\theta} + p_B e^{-i\theta} + (1 - p_A - p_B)$ . By the Wiener—Hopf factorization of the random walk, the probability generating function of the indicator sequence $\Pr[d_t > 0]$ admits the factorization:

\sum_{t=0}^{\infty} \Pr[d_t > 0] \, z^t = \frac{1}{1 - z} \cdot \frac{r_+ z}{1 - r_+ z}

where $r_+$ is the unique root of $p_B r^2 - r + p_A = 0$ in $(0, 1)$ . This follows from the classical identity: the generating function of $\Pr[\max_{0 \leq s \leq t} d_s > 0]$ decomposes via the ascending ladder epoch distribution, whose generating function is $r_+ z / (1 - r_+ z)$ . Evaluating at $z = 1 - p$ yields the stated formula. $\square$

Lemma 1 (Root Bounds). If $p_A, p_B > 0$ and $p_A + p_B \leq 1$ , then $0 < r_+ < 1$ .

Proof. Consider $f(r) = p_B r^2 - r + p_A$ . We have $f(0) = p_A > 0$ and $f(1) = p_A + p_B - 1 \leq 0$ , so by the Intermediate Value Theorem, $f$ has a root in $(0, 1]$ . The discriminant is $\Delta = 1 - 4p_A p_B$ . By the AM—GM inequality applied to $p_A + p_B \leq 1$ , we have $4p_A p_B \leq (p_A + p_B)^2 \leq 1$ , so $\Delta \geq 0$ and both roots are real. The smaller root $r_+ = (1 - \sqrt{\Delta})/(2p_B)$ satisfies $r_+ < 1/(2p_B)$ . Since $f(1) \leq 0$ and $f(r_+) = 0$ with $f$ opening upward, we need $r_+ \leq 1$ . Strict inequality $r_+ < 1$ holds when $p_A + p_B < 1$ . $\square$

Lemma 2 (Convergence of the Series). For $0 < p < 1$ , the sum $\sum_{t=0}^{\infty} (1-p)^t \Pr[d_t > 0]$ converges absolutely.

Proof. Since $\Pr[d_t > 0] \leq 1$ for all $t$ , we have $\sum_{t=0}^{\infty} (1-p)^t \Pr[d_t > 0] \leq \sum_{t=0}^{\infty} (1-p)^t = 1/p < \infty$ . $\square$

Editorial

Compute H(n) = sum_{k=3}^{n} E_A(1/sqrt(k+3), 1/sqrt(k+3)+1/k^2, 1/k^3) where E_A is the expected number of times A leads in a match with geometric killing rate p. Key formula: E_A = r_+ * z / ((1 - z) * (1 - r_+ * z)) where z = 1 - p, r_+ = (1 - sqrt(1 - 4pApB)) / (2*pB).

Pseudocode

    total = 0
    For k from 3 to n:
        pA = 1 / sqrt(k + 3)
        pB = pA + 1 / k^2
        p = 1 / k^3
        disc = 1 - 4 * pA * pB
        r_plus = (1 - sqrt(disc)) / (2 * pB)
        z = 1 - p
        E_A = (r_plus * z) / ((1 - z) * (1 - r_plus * z))
        total += E_A
    Return total

Complexity Analysis

Time: $O(n)$ — one constant-time evaluation per summand $k = 3, \ldots, n$ .
Space: $O(1)$ — only running accumulators are stored.

Answer

\boxed{646231.2177}

C++ project_euler/problem_661/solution.cpp

#include <bits/stdc++.h>
using namespace std;

/*
 * Problem 661: A Long Chess Match
 *
 * Compute H(n) = sum_{k=3}^{n} E_A(1/sqrt(k+3), 1/sqrt(k+3)+1/k^2, 1/k^3)
 *
 * E_A uses the Wiener-Hopf factorization:
 *   r_plus = (1 - sqrt(1 - 4*pA*pB)) / (2*pB)
 *   E_A = r_plus * z / ((1 - z) * (1 - r_plus * z))  where z = 1 - p
 *
 * Two methods: closed-form and truncated Markov chain (cross-check).
 */

// Method 1: Closed-form via generating function
double expected_leading_closed(double pA, double pB, double p) {
    double disc = 1.0 - 4.0 * pA * pB;
    if (disc < 0) disc = 0.0;
    double r_plus = (1.0 - sqrt(disc)) / (2.0 * pB);
    double z = 1.0 - p;
    double denom = (1.0 - z) * (1.0 - r_plus * z);
    if (fabs(denom) < 1e-30) return 0.0;
    return r_plus * z / denom;
}

// Method 2: Direct DP summation (for verification on small cases)
double expected_leading_dp(double pA, double pB, double p, int M = 200) {
    int size = 2 * M + 1;
    double surv = 1.0 - p;
    double a = pB * surv;      // move left
    double b = (1 - pA - pB) * surv; // stay
    double c = pA * surv;      // move right

    vector<double> prob(size, 0.0), new_prob(size, 0.0);
    prob[M] = 1.0; // start at d = 0

    double result = 0.0;
    int T_max = min((int)(10.0 / p), 50000);

    for (int t = 1; t <= T_max; t++) {
        fill(new_prob.begin(), new_prob.end(), 0.0);
        for (int i = 1; i < size - 1; i++) {
            new_prob[i] = prob[i-1] * c + prob[i] * b + prob[i+1] * a;
        }
        swap(prob, new_prob);
        double pos_prob = 0.0;
        for (int i = M + 1; i < size; i++) pos_prob += prob[i];
        result += pos_prob;
        if (pos_prob < 1e-15) break;
    }
    return result;
}

int main() {
    // Compute H(50)
    double H = 0.0;
    for (int k = 3; k <= 50; k++) {
        double pA = 1.0 / sqrt(k + 3.0);
        double pB = pA + 1.0 / ((double)k * k);
        double p = 1.0 / ((double)k * k * k);
        double ea = expected_leading_closed(pA, pB, p);
        H += ea;
    }

    printf("H(50) = %.4f\n", H);

    // Cross-check: known example
    double check = expected_leading_closed(0.25, 0.25, 0.5);
    printf("E_A(0.25,0.25,0.5) = %.6f (expected ~0.585786)\n", check);
    assert(fabs(check - 0.585786) < 0.001);

    // Cross-check H(3)
    double H3 = 0.0;
    for (int k = 3; k <= 3; k++) {
        double pA = 1.0 / sqrt(k + 3.0);
        double pB = pA + 1.0 / ((double)k * k);
        double p = 1.0 / ((double)k * k * k);
        H3 += expected_leading_closed(pA, pB, p);
    }
    printf("H(3) = %.4f (expected ~6.8345)\n", H3);

    return 0;
}

Python project_euler/problem_661/solution.py

"""
Problem 661: A Long Chess Match

Compute H(n) = sum_{k=3}^{n} E_A(1/sqrt(k+3), 1/sqrt(k+3)+1/k^2, 1/k^3)
where E_A is the expected number of times A leads in a match with
geometric killing rate p.

Key formula: E_A = r_+ * z / ((1 - z) * (1 - r_+ * z))
where z = 1 - p, r_+ = (1 - sqrt(1 - 4*pA*pB)) / (2*pB).
"""

from math import sqrt

# --- Method 1: Closed-form via Wiener-Hopf factorization ---
def expected_leading_closed(pA, pB, p):
    """Expected times A leads, using closed-form generating function."""
    disc = 1 - 4 * pA * pB
    if disc < 0:
        disc = 0.0
    r_plus = (1 - sqrt(disc)) / (2 * pB)
    z = 1 - p
    # E_A = sum_{t>=0} (1-p)^t * Pr[d_t > 0]
    # = r_plus * z / ((1 - z) * (1 - r_plus * z))
    denom = (1 - z) * (1 - r_plus * z)
    if abs(denom) < 1e-30:
        return 0.0
    return r_plus * z / denom

# --- Method 2: Absorbing Markov chain (tridiagonal matrix) ---
def expected_leading_markov(pA, pB, p, M=300):
    """Expected times A leads via truncated Markov chain."""
    # States: d in {-M, ..., 0, ..., M}, total 2M+1 states
    # State index i corresponds to d = i - M
    size = 2 * M + 1
    # Build tridiagonal system (I - Q) x = e_M (start at d=0)
    # Using Thomas algorithm for efficiency

    # Transition probabilities (with killing)
    surv = 1 - p
    a_coeff = pB * surv  # move left
    b_coeff = (1 - pA - pB) * surv  # stay
    c_coeff = pA * surv  # move right

    # (I - Q) is tridiagonal: diagonal = 1-b, off-diag = -a (lower), -c (upper)
    # Solve (I-Q) x = e_{M} where M is the index of d=0
    diag = [1.0 - b_coeff] * size
    lower = [-a_coeff] * size  # lower[i] is coefficient of x[i-1] in row i
    upper = [-c_coeff] * size  # upper[i] is coefficient of x[i+1] in row i

    # Boundary: absorbing at edges
    diag[0] = 1.0
    diag[-1] = 1.0
    lower[0] = 0.0
    upper[-1] = 0.0

    rhs = [0.0] * size
    rhs[M] = 1.0  # start at d = 0

    # Forward elimination
    for i in range(1, size):
        if abs(diag[i-1]) < 1e-30:
            continue
        factor = lower[i] / diag[i-1]
        diag[i] -= factor * upper[i-1]
        rhs[i] -= factor * rhs[i-1]

    # Back substitution
    x = [0.0] * size
    x[-1] = rhs[-1] / diag[-1] if abs(diag[-1]) > 1e-30 else 0.0
    for i in range(size - 2, -1, -1):
        x[i] = (rhs[i] - upper[i] * x[i+1]) / diag[i] if abs(diag[i]) > 1e-30 else 0.0

    # Sum expected visits to positive states (d > 0, indices M+1 to 2M)
    # But wait - x gives expected visits starting from d=0
    # We need to solve for each target state... this approach needs N columns.
    # Instead, use the fact that E_A = sum of N[0, d] for d > 0.
    # For the fundamental matrix, we need row 0 (start=d=0).
    # This is actually computing one column of (I-Q)^{-1}, not what we want.

    # Alternative: direct simulation sum
    # E_A = sum_{t=0}^{T} (1-p)^t * Pr[d_t > 0]
    # Use convolution / FFT or direct DP

    # DP: prob[t][d] = probability of being at state d at time t
    T_max = min(int(10 / p), 50000)  # enough steps
    # Use two arrays for current and next
    prob = [0.0] * size
    prob[M] = 1.0  # start at d=0

    result = 0.0
    for t in range(1, T_max + 1):
        new_prob = [0.0] * size
        for i in range(1, size - 1):
            new_prob[i] = (prob[i-1] * c_coeff +
                          prob[i] * b_coeff +
                          prob[i+1] * a_coeff)
        prob = new_prob
        # Add probability of being in positive states
        pos_prob = sum(prob[M+1:])
        result += pos_prob
        if pos_prob < 1e-15:
            break

    return result

# --- Compute H(n) ---
def compute_H(n):
    """Compute H(n) = sum_{k=3}^{n} E_A(params(k))."""
    total = 0.0
    terms = []
    for k in range(3, n + 1):
        pA = 1.0 / sqrt(k + 3)
        pB = pA + 1.0 / (k * k)
        p = 1.0 / (k ** 3)
        ea = expected_leading_closed(pA, pB, p)
        total += ea
        terms.append((k, pA, pB, p, ea))
    return total, terms

# --- Compute and verify ---
H3, _ = compute_H(3)
print(f"H(3) = {H3:.4f} (expected ~6.8345)")

H50, terms = compute_H(50)
print(f"H(50) = {H50:.4f}")

# Verify with Markov chain for small case
ea_check = expected_leading_closed(0.25, 0.25, 0.5)
print(f"E_A(0.25, 0.25, 0.5) = {ea_check:.6f} (expected ~0.585786)")

# Cross-check Method 1 vs Method 2 for k=3
pA3 = 1.0 / sqrt(6)
pB3 = pA3 + 1.0 / 9
p3 = 1.0 / 27
ea_closed = expected_leading_closed(pA3, pB3, p3)
print(f"E_A closed (k=3) = {ea_closed:.6f}")

# Assertions
assert abs(H3 - 6.8345) < 0.01, f"H(3) mismatch: {H3}"
assert abs(ea_check - 0.585786) < 0.001, f"E_A check mismatch: {ea_check}"