Problem 285: Pythagorean Odds

Project Euler Problem 285 Solution

EulerSolve provides an optimized solution for Project Euler Problem 285, Pythagorean Odds, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary On turn $k$, the random point is chosen uniformly as $$P=(1+ka,\ 1+kb),\qquad a,b\in[0,1],$$ so $P$ is uniform in the square $$Q_k=[1,k+1]\times[1,k+1].$$ The score is $k$ exactly when the distance from the origin rounds to $k$, i.e. $$k-\frac12 \le \|P\| < k+\frac12.$$ The program computes the expected total score up to $k=100000$. The final numeric answer is intentionally omitted here. Mathematical Approach 1) Expected score on one turn is an area. Since $P$ is uniform in a square of area $k^2$, the probability of scoring $k$ equals $$\frac{\operatorname{area}\bigl(Q_k\cap \mathcal A_k\bigr)}{k^2},$$ where $\mathcal A_k$ is the annulus $$\mathcal A_k=\left\{(x,y):k-\frac12\le \sqrt{x^2+y^2}<k+\frac12\right\}.$$ Therefore the expected contribution of turn $k$ is $$k\cdot \frac{\operatorname{area}(Q_k\cap \mathcal A_k)}{k^2} =\frac{\operatorname{area}(Q_k\cap \mathcal A_k)}{k}.$$ 2) Why the upper edges of the square do not matter. Every point in the outer disk $$x^2+y^2\le \left(k+\frac12\right)^2$$ has $x\le k+\frac12$ and $y\le k+\frac12$, which is strictly inside the square bounds $x,y\le k+1$. So the only relevant constraints from $Q_k$ are $$x\ge 1,\qquad y\ge 1.$$ This is the key simplification behind the whole solution. 3) Introduce a cumulative area primitive....

Detailed mathematical approach

Problem Summary

On a turn, Albert chooses a positive integer $k$, two real numbers $a$ and $b$ are drawn uniformly at random from the interval $[0,1]$, and the number

$$d=\sqrt{(ka+1)^2+(kb+1)^2}$$

is rounded to the nearest integer. If the result equals $k$, he scores $k$ points; otherwise he scores nothing. In the statement's example, $k=6$, $a=0.2$ and $b=0.85$ give $d^2=42.05$ and $d=6.484\ldots$, which rounds to $6$, so this turn scores $6$ points. If he plays one turn with each of $k=1,2,\dots,10$, his expected total score is $10.20914$ when rounded to five decimal places; we must find this rounded expected total for one turn with each of $k=1,2,\dots,10^5$. ^[1]

The whole solution in one sentence: the point $(ka+1,kb+1)$ is a uniform random point of a square, a turn scores exactly when this point lies in the ring between the circles of radii $k-\frac12$ and $k+\frac12$ around the origin, so the expected score of turn $k$ is the area of the part of that ring inside the square divided by $k$, and an elementary formula for a disc cut by two lines turns the expected total into $10^5$ explicit terms.

Mathematical Approach

We turn the game into geometry: Sections 1–4 reduce the expected score of a turn to an area, Sections 5 and 6 compute it, and Sections 7–10 check the result and study the size and the accuracy of the terms.

1. A random point in a square

Put $x=ka+1$ and $y=kb+1$. If $a$ is uniform on $[0,1]$, then $x$ is uniform on $[1,k+1]$: the map $a\mapsto 1+ka$ stretches the interval by the factor $k$ and shifts it by $1$. ^[2] The same holds for $y$, and the two numbers are chosen independently, so the point $(x,y)$ is uniform on the square

$$Q_k=[1,k+1]\times[1,k+1],$$

whose side is $k$, whose area is $k^2$ and whose lower left corner is $(1,1)$. By Pythagoras' theorem, $d$ is the distance from the origin $O=(0,0)$ to the point $(x,y)$, which explains the title. In the statement's example the point is $(2.2,\,6.1)$, and indeed $2.2^2+6.1^2=4.84+37.21=42.05$.

2. When a turn scores

A number rounds to $k$ exactly when it lies between $k-\frac12$ and $k+\frac12$. So turn $k$ scores when

$$k-\tfrac12\le d\lt k+\tfrac12,$$

that is, when $(x,y)$ lies in the ring $R_k$ between the circles around $O$ with radii $k-\frac12$ and $k+\frac12$, an annulus. ^[3] A point exactly on one of the circles has probability $0$, so the rounding of halves does not matter. In the example, $d=6.484\ldots$ lies between $5.5$ and $6.5$, so $(2.2,\,6.1)$ lies in $R_6$.

3. Probability is area

For a uniform random point of a region, the chance of landing in a part of it is the fraction of the area that this part covers, the basic principle of geometric probability. ^[4] Hence turn $k$ scores with probability

$$p_k=\frac{\operatorname{area}(R_k\cap Q_k)}{k^2}.$$

The score of turn $k$ is $k$ with probability $p_k$ and $0$ otherwise, so its expected value is $e_k=k\,p_k$. The expected value of a sum is the sum of the expected values, whether or not the terms are independent (linearity of expectation). ^[5] For the turns $k=1,\dots,K$ the expected total score is therefore

$$\begin{aligned}E_K&=\sum_{k=1}^{K}e_k\\&=\sum_{k=1}^{K}\frac{\operatorname{area}(R_k\cap Q_k)}{k}.\end{aligned}$$

4. The far edges of the square never matter

Every point of the ring $R_k$ lies at distance less than $k+\frac12$ from $O$, so both of its coordinates are less than $k+\frac12$ and hence less than $k+1$. On $R_k$ the conditions $x\le k+1$ and $y\le k+1$ of the square therefore hold automatically, and only $x\ge1$ and $y\ge1$ remain: for $k\ge2$ the ring crosses the square from its bottom edge to its left edge and never reaches the top or the right edge. For $r\ge0$ let $A(r)$ be the area of the set of points with $x\ge1$, $y\ge1$ and $x^2+y^2\le r^2$, the part of the disc of radius $r$ around $O$ in this quarter plane. Removing the part of the inner disc from the part of the outer disc gives

$$E_K=\sum_{k=1}^{K}\frac{A(k+\frac12)-A(k-\frac12)}{k}.$$

5. A disc cut by two lines

The point of the quarter plane closest to $O$ is its corner $(1,1)$, at distance $\sqrt2\approx1.414$. So $A(r)=0$ for $r\le\sqrt2$. For $r\ge\sqrt2$ we start from the quarter disc

$$D=\{(x,y):x\ge0,\ y\ge0,\ x^2+y^2\le r^2\}$$

of area $\pi r^2/4$ and remove its points with $x\lt1$ or $y\lt1$. These points form two strips, the part of $D$ with $x\lt1$ and the part with $y\lt1$, mirror images of each other in the line $y=x$ with the same area $s(r)$. They overlap in the unit square $[0,1)\times[0,1)$, which lies inside $D$ because all its points are closer to $O$ than $\sqrt2\le r$. By inclusion–exclusion, the area of the union of two sets is the sum of their areas minus the area of their overlap, ^[6] so the two strips together cover the area $2s(r)-1$, and

$$A(r)=\frac{\pi r^2}{4}-2s(r)+1\qquad(r\ge\sqrt2).$$

6. The strip: a triangle and a sector

The circle $x^2+y^2=r^2$ meets the line $x=1$ at $P=(1,\sqrt{r^2-1})$. The segment from $O$ to $P$ cuts the strip $0\le x\le1$ of $D$ into two pieces:

the right triangle with vertices $O$, $(1,0)$ and $P$, with legs $1$ and $\sqrt{r^2-1}$ and area $\frac12\sqrt{r^2-1}$;
the circular sector between the segment $OP$ and the positive $y$-axis. Its angle $\varphi$ at $O$ satisfies $\sin\varphi=1/r$, because $P$ lies at distance $1$ from the $y$-axis and at distance $r$ from $O$. So, in radians, $\varphi=\arcsin(1/r)$. The sector is the fraction $\varphi/(2\pi)$ of the disc of area $\pi r^2$, so its area is $r^2\varphi/2$. ^[7]

Adding the two pieces and inserting $s(r)$ into Section 5 gives the closed formula

$$s(r)=\frac{\sqrt{r^2-1}+r^2\arcsin(1/r)}{2},$$

$$A(r)=r^2\Bigl(\frac{\pi}{4}-\arcsin\frac1r\Bigr)-\sqrt{r^2-1}+1.$$

At $r=\sqrt2$ we have $\arcsin(1/\sqrt2)=\pi/4$ and $\sqrt{r^2-1}=1$, so the formula gives $0$ and joins the zero branch without a jump. The case distinction is necessary: for $r\lt1$ the arcsine of $1/r$ is undefined, and for $1\le r\lt\sqrt2$ the formula is wrong (about $0.2146$ at $r=1$), because the unit square no longer lies inside $D$. Among the radii $k\pm\frac12$, only $\frac12$ is below $\sqrt2$.

7. Worked example: the first two turns

Turn $1$ uses the radii $\frac12$ and $\frac32$, and turn $2$ the radii $\frac32$ and $\frac52$. Here $A(\frac12)=0$, and Sections 5 and 6 give:

$r=\frac32$: the quarter disc $9\pi/16\approx1.7671459$, the triangle $\frac12\sqrt{5/4}\approx0.5590170$ and the sector $\frac98\arcsin\frac23$, about $0.8209436$, so $s(\frac32)\approx1.3799606$ and $A(\frac32)\approx1.7671459-2.7599212$${}+1\approx0.0072247$;
$r=\frac52$: the quarter disc $25\pi/16\approx4.9087385$, the triangle $\frac12\sqrt{21/4}\approx1.1456439$ and the sector $\frac{25}8\arcsin\frac25$, about $1.2859901$, so $s(\frac52)\approx2.4316341$ and $A(\frac52)\approx1.0454704$.

The square $Q_1=[1,2]\times[1,2]$ has area $1$, and only a tiny piece at its corner $(1,1)$ lies within distance $1.5$ of $O$, so $p_1=e_1\approx0.0072247$. On turn $2$ the part of the ring has the area $1.0454704-0.0072247=1.0382457$, so $p_2$ is this area divided by $4$, about $0.2595614$, and $e_2$ is the area divided by $2$, about $0.5191229$. Together the first two turns are worth $0.5263475$ points on average.

8. All ten turns of the statement

The same computation for $k=1,\dots,10$ gives these expected scores $e_k$, rounded to seven decimals:

$k=1,\dots,5$: $0.0072247$, $0.5191229$, $0.8906989$, $1.0653446$, $1.1680518$;
$k=6,\dots,10$: $1.2358889$, $1.2840960$, $1.3201381$, $1.3481129$, $1.3704606$.

Their sum is $E_{10}\approx10.2091394$, which rounds to the checkpoint $10.20914$ of the statement. The probability $p_k$ is largest at $k=3$, about $0.2969$, and decreases from there on, while $e_k=k\,p_k$ keeps growing; both were checked numerically for all $k\le10^5$, and the next section finds the limit of $e_k$.

9. How large the terms become

For large $r$, the arcsine series $\arcsin u=u+u^3/6+\cdots$ ^[8] and the expansion of the square root give the following, where $O(r^{-3})$ stands for a remainder of at most a constant times $r^{-3}$:

$$\begin{aligned}r^2\arcsin(1/r)&=r+\tfrac{1}{6r}+O(r^{-3}),\\ \sqrt{r^2-1}&=r-\tfrac{1}{2r}+O(r^{-3}),\end{aligned}$$

and hence, by Section 6,

$$A(r)=\frac{\pi r^2}{4}-2r+1+\frac{1}{3r}+O(r^{-3}).$$

So the part of the ring on turn $k$ has the area $\pi k/2-2+O(k^{-2})$. By the area formula of an annulus, ^[3] the whole ring has the area

$$\pi\bigl((k+\tfrac12)^2-(k-\tfrac12)^2\bigr)=2\pi k,$$

its quarter in the first quadrant has the area $\pi k/2$, and each of the two strips of width $1$ along the axes cuts away almost a unit square. In the difference, the term $1/(3r)$ contributes

$$\frac13\Bigl(\frac{1}{k+\frac12}-\frac{1}{k-\frac12}\Bigr)=-\frac{1}{3(k^2-\frac14)},$$

and the two remainders contribute $O(k^{-3})$, so after the division by $k$

$$e_k=\frac{\pi}{2}-\frac{2}{k}-\frac{1}{3k^3}+O(k^{-4}).$$

For $k=10$ the right-hand side without the error term is $1.3704630$, against $e_{10}\approx1.3704606$. The probability of scoring falls roughly like $\pi/(2k)$, while the expected score of a turn approaches $\pi/2\approx1.5708$; the term $-2/k$ is the price of the shift by $1$ in $ka+1$ and $kb+1$.

10. Five decimals in double precision

The code subtracts $A(k-\frac12)$ from $A(k+\frac12)$ in floating point. At $k=10^5$ both values are about $7.85\cdot10^9$, their difference only about $\pi k/2$, so about five of the roughly sixteen significant digits of a double cancel, an effect called catastrophic cancellation. ^[9] Compared with a 35-digit recomputation, the values of $A$ from the three programs on the Mac used here are off by at most $1.8\cdot10^{-6}$, each term by at most $2\cdot10^{-11}$ and all terms together by at most $2.8\cdot10^{-7}$; the $10^5$ additions, each rounded by at most half a unit in the last place of the running total, add at most $6.5\cdot10^{-7}$. So the double sum is within $9.3\cdot10^{-7}$ of the exact value even if all errors point the same way; the actual difference is at most $6.3\cdot10^{-10}$. The exact value lies about $1.7\cdot10^{-6}$ from the nearest point where its fifth decimal would change, so the rounding is correct, although the worst-case margin is less than a factor of two. Subtracting the two terms $\pi r^2/4$ by hand, which gives exactly $\pi k/2$, would remove the largest cancellation.

How the Code Works

All three programs follow the same plan; the C++ program ^[10] adds options and a self-check. area_inside_square_and_circle(r), in Java areaInsideSquareAndCircle, is the function $A$ of Sections 5 and 6: it returns $0$ when $r\le\sqrt2$, compared with the constant kSqrt2, and otherwise computes t, the square root of $r^2-1$, takes $\pi$ as acos(-1) and returns $r^2(\pi/4-\arcsin(1/r))-t+1$. solve adds (outer - inner) / k to expected for $k=1,\dots,K$, where outer is $A(k+\frac12)$ and inner is $A(k-\frac12)$; this is the sum of Section 4. Since inner repeats the previous turn's outer, almost every value is computed twice; reusing it would halve work that takes milliseconds anyway.

The guard $r\le\sqrt2$ is more than a shortcut: for the inner radius $\frac12$ of turn $1$ the formula would take the square root of $-3/4$ and the arcsine of $2$, which give NaN in C++ and Java and would spoil the sum, and raise a ValueError in Python.

The C++ main reads --max-k=N (default $10^5$) and --skip-checkpoints and rejects anything else and $K\lt1$. parse_int_after_prefix accumulates the digits in an int without an overflow check, so an out-of-range value is undefined behaviour (the tested builds turned --max-k=4294967298 into two turns), and the largest int as limit would overflow the final ++k in solve; neither affects $K=10^5$. Unless skipped, run_checkpoints rounds solve(10) to five decimals with std::round and stops the program with exit code 2 if this differs from $10.20914$ by more than $10^{-12}$. The total is converted to double and printed with five decimals.

In C++, solve(max_k) returns a long double. Python's solve(max_k=100000) has the limit as a default parameter and Java's solve() fixes maxK at $10^5$; both return the total as a string with five decimals and have no options or checkpoint. Java uses Math.PI, the same double as acos(-1), and Locale.US, which fixes the decimal point.

The C++ code computes in long double: 80-bit extended precision with GCC on x86, but the same 64-bit format as double on 64-bit ARM under Windows and macOS, ^[11] including the Mac used here, where Java's double and Python's float have this format as well. Equal formats do not guarantee equal results: Java reproduces Python's values of $A$ bit for bit, while the tested C++ builds, which fuse multiplications with the following subtractions into single instructions, differ in the last bits for about a quarter of the radii. Section 10 covers all three on this Mac; an 80-bit build has smaller error bounds but was not tested.

Complexity and Verification

The program evaluates $A$ twice per turn at a constant cost, among others one square root of $r^2-1$ and one arcsine, so the time is $O(K)$ and the memory $O(1)$, as only the running total is kept. The C++ build runs in a few milliseconds, Python and Java in well under a second. Unless --skip-checkpoints is given, the C++ program checks before the main sum that $E_{10}$, rounded to five decimals, is $10.20914$. For this article the following was also verified; these checks are not part of the published solutions:

the formula of Section 6 agrees with numerical integration of $\sqrt{r^2-x^2}-1$ for $r=\frac32,\frac52,\frac72,\frac{11}2,\frac{13}2,\frac{21}2$ in all $15$ digits compared;
the area of $R_k\cap Q_k$, integrated numerically over the whole square including its far edges, equals $A(k+\frac12)-A(k-\frac12)$ to within $10^{-29}$ for every $k\le10$;
a seeded simulation of the game, $10^7$ random turns for each $k\le10$, gives $10.21035\pm0.00233$ (one standard error), consistent with $E_{10}$;
the error figures of Section 10 compare the values of $A$ from each of the three programs with a 35-digit recomputation;
C++ (built with Clang and with GCC), Python and Java print the same five decimals for $K=10^5$, and the C++ program prints $10.20914$ for --max-k=10 and $0.52635$ for --max-k=2, the two turns of Section 7.

Footnotes and References

The references below supply the background facts used above. The reduction to the function $A$, its closed formula and the expansion of the terms are derived in this article; the error figures of Section 10 are measurements of the tested programs plus a bound for the additions.

Project Euler 285 — Pythagorean Odds. The official statement: for a chosen k and two numbers a, b uniform in [0, 1], k points are scored if the square root of (ka + 1)² + (kb + 1)², rounded to the nearest integer, equals k. The example k = 6, a = 0.2, b = 0.85 gives the sum of squares 42.05; the expected total for k = 1, …, 10 rounds to 10.20914, and the task asks for k = 1, …, 10⁵.
Wikipedia — Continuous uniform distribution. If u is sampled from the standard uniform distribution on [0, 1], then c + (d − c)u is uniform on [c, d]; with c = 1 and d = k + 1 this is the map of Section 1.
Wikipedia — Annulus (mathematics). An annulus is the region between two concentric circles; with outer radius R and inner radius r its area is π(R² − r²), as used in Section 9.
Wikipedia — Geometric probability. The article credits Newton with the fundamental principle that chance is proportional to area fraction, the rule used in Section 3.
Wikipedia — Expected value. Expectation is linear: E[X + Y] = E[X] + E[Y] for any random variables X and Y whenever the right-hand side is defined, so no independence is needed.
Wikipedia — Inclusion–exclusion principle. For two sets, |A ∪ B| = |A| + |B| − |A ∩ B|; the article notes that the identities hold for any measure and sets of finite measure, hence for areas (Section 5).
Wikipedia — Circular sector. A sector of radius r with central angle θ in radians has area r²θ/2; Section 6 uses it with θ = arcsin(1/r).
Wikipedia — Inverse trigonometric functions. The series arcsin z = z + (1/2)z³/3 + (1·3/(2·4))z⁵/5 + ⋯ holds for |z| ≤ 1; Section 9 uses its first two terms.
Wikipedia — Catastrophic cancellation. Subtracting good approximations of two nearby numbers may give a very bad approximation of their difference; the relative error is inversely proportional to the true difference (Section 10).
C++ — Euler285.cpp. The linked immutable C++ revision contains area_inside_square_and_circle, solve and the checkpoint in run_checkpoints. It is the implementation record for this article.
Wikipedia — Long double. With GCC, long double is the 80-bit extended precision format on x86 processors; with Microsoft Visual C++ on x86 and on 64-bit ARM under Windows and macOS it is simply a synonym for double.

Problem 285 source code

C++

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <string>

namespace {

struct Options {
    int max_k = 100000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--max-k=", options.max_k)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.max_k >= 1;
}

long double area_inside_square_and_circle(const long double r) {
    static const long double kSqrt2 = std::sqrt(2.0L);
    if (r <= kSqrt2) {
        return 0.0L;
    }
    const long double t = std::sqrt(r * r - 1.0L);
    const long double pi = std::acos(-1.0L);
    return r * r * (pi / 4.0L - std::asin(1.0L / r)) - t + 1.0L;
}

long double solve(const int max_k) {
    long double expected = 0.0L;
    for (int k = 1; k <= max_k; ++k) {
        const long double outer = area_inside_square_and_circle(static_cast<long double>(k) + 0.5L);
        const long double inner = area_inside_square_and_circle(static_cast<long double>(k) - 0.5L);
        expected += (outer - inner) / static_cast<long double>(k);
    }
    return expected;
}

bool run_checkpoints() {
    const long double sample = solve(10);
    const long double rounded_sample = std::round(sample * 100000.0L) / 100000.0L;
    if (std::fabsl(rounded_sample - 10.20914L) > 1e-12L) {
        std::cerr << "Checkpoint failed for 10-turn sample" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    const long double answer = solve(options.max_k);
    std::cout << std::fixed << std::setprecision(5) << static_cast<double>(answer) << '\n';
    return 0;
}

Python

import math

def area_inside_square_and_circle(r):
    kSqrt2 = math.sqrt(2.0)
    if r <= kSqrt2:
        return 0.0
    t = math.sqrt(r * r - 1.0)
    pi = math.acos(-1.0)
    return r * r * (pi / 4.0 - math.asin(1.0 / r)) - t + 1.0

def solve(max_k=100000):
    expected = 0.0
    for k in range(1, max_k + 1):
        outer = area_inside_square_and_circle(k + 0.5)
        inner = area_inside_square_and_circle(k - 0.5)
        expected += (outer - inner) / float(k)
    return f"{expected:.5f}"

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

Java

import java.util.Locale;

public class Euler285 {
    static double areaInsideSquareAndCircle(double r) {
        double kSqrt2 = Math.sqrt(2.0);
        if (r <= kSqrt2) {
            return 0.0;
        }
        double t = Math.sqrt(r * r - 1.0);
        double pi = Math.PI;
        return r * r * (pi / 4.0 - Math.asin(1.0 / r)) - t + 1.0;
    }

    public static String solve() {
        int maxK = 100000;
        double expected = 0.0;
        for (int k = 1; k <= maxK; ++k) {
            double outer = areaInsideSquareAndCircle(k + 0.5);
            double inner = areaInsideSquareAndCircle(k - 0.5);
            expected += (outer - inner) / (double) k;
        }
        return String.format(Locale.US, "%.5f", expected);
    }

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