Problem 285: Pythagorean Odds

Mathematical Foundation

Theorem (Winning Probability as Annular Area). The probability of winning round $k$ is $P_k = A(r_2, c) - A(r_1, c)$, where $r_1 = \frac{2k-1}{2k}$, $r_2 = \frac{2k+1}{2k}$, $c = \frac{1}{k}$, and

Proof. The winning condition $|k - \sqrt{(ka+1)^2+(kb+1)^2}| < \frac{1}{2}$ is equivalent to

Substituting $x = a + 1/k$, $y = b + 1/k$, the unit square $[0,1]^2$ in $(a,b)$ maps to $[c, 1+c]^2$ where $c = 1/k$, and the condition becomes $r_1^2 \le x^2 + y^2 \le r_2^2$, an annular region.

Since $r_2 = 1 + \frac{1}{2k} < 1 + c$ and $r_1 < r_2$, the outer disk lies within the square’s upper bounds. The binding constraints are only $x \ge c$ and $y \ge c$.

Define $A(r, c) = \text{Area}\{(x,y) : x \ge c,\, y \ge c,\, x^2 + y^2 \le r^2\}$. For $r \ge c\sqrt{2}$ (so the disk intersects the first-quadrant region beyond $(c,c)$):

Evaluating the integral using $\int \sqrt{r^2-x^2}\,dx = \frac{x\sqrt{r^2-x^2}}{2} + \frac{r^2}{2}\arcsin\frac{x}{r} + C$:

At $x = \sqrt{r^2 - c^2}$: $\sqrt{r^2-x^2} = c$, so the first term gives $\frac{c\sqrt{r^2-c^2}}{2}$, the second gives $\frac{r^2}{2}\arcsin\frac{\sqrt{r^2-c^2}}{r}$, and the third gives $-c\sqrt{r^2-c^2}$.

At $x = c$: $\sqrt{r^2-c^2}$ appears in the first term as $\frac{c\sqrt{r^2-c^2}}{2}$, the second gives $\frac{r^2}{2}\arcsin\frac{c}{r}$, and the third gives $-c^2$.

Combining:

Using $\arcsin\frac{\sqrt{r^2-c^2}}{r} = \frac{\pi}{2} - \arcsin\frac{c}{r}$:

The probability is $P_k = A(r_2, c) - A(r_1, c)$, and the expected score is $E = \sum_{k=1}^{N} k \cdot P_k$. $\quad\square$

Editorial

Albert picks a,b uniformly in [0,1]. Computes sqrt((ka+1)^2+(kb+1)^2), rounds to nearest integer. Scores k if result = k, else 0. Find expected total for k=1..100000, to 5 decimals. Winning condition: (k-1/2)^2 <= (ka+1)^2 + (kb+1)^2 < (k+1/2)^2 Substituting x=a+1/k, y=b+1/k: circle centered at origin, x in [1/k, 1+1/k], y same. Inner radius r1 = (k-1/2)/k, outer radius r2 = (k+1/2)/k. Area of {x>=c, y>=c, x^2+y^2 <= r^2} where c=1/k: For r^2 >= 2c^2: integral_c^sqrt(r^2-c^2) (sqrt(r^2-x^2) - c) dx = [xsqrt(r^2-x^2)/2 + r^2/2 * arcsin(x/r)] from c to sqrt(r^2-c^2) - c(sqrt(r^2-c^2)-c) Uses mpmath for sufficient precision.

Pseudocode

    total = 0.0
    For k from 1 to N:
        c = 1.0 / k
        r1 = (2*k - 1) / (2.0 * k)
        r2 = (2*k + 1) / (2.0 * k)
        P_k = area(r2, c) - area(r1, c)
        total += k * P_k
    Return round(total, 5)

    If r < c * sqrt(2) then
        Return 0.0
    Return (r^2 / 2) * (pi/2 - 2*arcsin(c/r)) - c*sqrt(r^2 - c^2) + c^2

Complexity Analysis

• Time: $O(N)$ where $N = 100000$. Each iteration involves $O(1)$ floating-point operations (including $\arcsin$ and $\sqrt{\cdot}$).
• Space: $O(1)$.

Answer

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

typedef long double ld;

// Area of {x >= c, y >= c, x^2+y^2 <= r^2}
// = integral_c^sqrt(r^2-c^2) (sqrt(r^2-x^2) - c) dx
// = [x*sqrt(r^2-x^2)/2 + r^2/2*arcsin(x/r)]_c^sqrt(r^2-c^2) - c*(sqrt(r^2-c^2)-c)
ld circle_quad_area(ld r, ld c) {
    if (r * r < 2 * c * c) return 0.0L;

    ld x_max = sqrtl(r * r - c * c);

    auto F = [&](ld x) -> ld {
        return x * sqrtl(r * r - x * x) / 2 + r * r / 2 * asinl(x / r);
    };

    ld area = F(x_max) - F(c) - c * (x_max - c);
    return max(area, (ld)0);
}

int main() {
    ld total = 0;
    for (int k = 1; k <= 100000; k++) {
        ld c = 1.0L / k;
        ld r1 = (k - 0.5L) / k;
        ld r2 = (k + 0.5L) / k;

        ld a2 = circle_quad_area(r2, c);
        ld a1 = circle_quad_area(r1, c);
        ld area = a2 - a1;

        total += k * area;
    }
    cout << fixed << setprecision(5) << (double)total << endl;
    return 0;
}

"""
Project Euler Problem 285: Pythagorean Odds

Albert picks a,b uniformly in [0,1]. Computes sqrt((ka+1)^2+(kb+1)^2), rounds to nearest
integer. Scores k if result = k, else 0. Find expected total for k=1..100000, to 5 decimals.

Winning condition: (k-1/2)^2 <= (ka+1)^2 + (kb+1)^2 < (k+1/2)^2

Substituting x=a+1/k, y=b+1/k: circle centered at origin, x in [1/k, 1+1/k], y same.
Inner radius r1 = (k-1/2)/k, outer radius r2 = (k+1/2)/k.

Area of {x>=c, y>=c, x^2+y^2 <= r^2} where c=1/k:
For r^2 >= 2c^2: integral_c^sqrt(r^2-c^2) (sqrt(r^2-x^2) - c) dx
= [x*sqrt(r^2-x^2)/2 + r^2/2 * arcsin(x/r)] from c to sqrt(r^2-c^2) - c*(sqrt(r^2-c^2)-c)

Uses mpmath for sufficient precision.
"""

from mpmath import mp, mpf, sqrt, asin

mp.dps = 25

def circle_quad_area(r, c):
    """Area of {x >= c, y >= c, x^2+y^2 <= r^2}."""
    if r * r < 2 * c * c:
        return mpf(0)
    x_max = sqrt(r * r - c * c)

    def F(x):
        return x * sqrt(r * r - x * x) / 2 + r * r / 2 * asin(x / r)

    return F(x_max) - F(c) - c * (x_max - c)

total = mpf(0)
for k in range(1, 100001):
    c = mpf(1) / k
    r1 = (mpf(2) * k - 1) / (2 * k)
    r2 = (mpf(2) * k + 1) / (2 * k)
    total += k * (circle_quad_area(r2, c) - circle_quad_area(r1, c))

print(f"{float(total):.5f}")