Problem 184: Triangles Containing the Origin

Project Euler Problem 184 Solution

EulerSolve provides an optimized solution for Project Euler Problem 184, Triangles Containing the Origin, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $$P_r=\{(x,y)\in\mathbb{Z}^2\setminus\{(0,0)\}:x^2+y^2\lt r^2\}.$$ For the actual Project Euler input, $r=105$. The task is to count unordered triangles whose three vertices are chosen from $P_{105}$ and whose interior contains the origin strictly, not merely on an edge or on a degenerate line. A direct approach would examine all $\binom{N}{3}$ triples of lattice points and run a point-in-triangle test on each one, where $N=|P_r|$. The implementations avoid that completely. They reorganize the point set by polar angle, group points that lie on the same ray from the origin, and then subtract the triples that are geometrically impossible instead of testing triangles one by one. Mathematical Approach From Triangles to Angular Intervals A nondegenerate triangle contains the origin if and only if the origin lies in the convex hull of its three vertices. Equivalently, a triangle misses the origin if there exists a line through the origin such that all three vertices lie in one closed half-plane determined by that line. In angular language, that means the three vertex directions fit inside some closed interval of length $\pi$....

Detailed mathematical approach

Problem Summary

Let $I_r$ be the set of points $(x,y)$ with integer coordinates and $x^2+y^2\lt r^2$, strictly inside the circle of radius $r$ centred at the origin. We must count the triangles with all three vertices in $I_{105}$ that contain the origin in their interior; a triangle is given by its three vertices. The statement gives three checkpoints: $I_2$ has nine points, the origin among them, and $8$ such triangles; $I_3$ gives $360$ and $I_5$ gives $10600$. ^[1]

The whole solution in one sentence: a triangle misses the origin exactly when some line through the origin has all three corners on one side of it or on it, so we count all triples of points and subtract the failing ones, which can be counted line by line from just the number of points on each line through the origin.

Mathematical Approach

We look at the smallest case, turn "the origin is inside" into a statement about directions, and count the failing triples one line through the origin at a time.

1. The smallest case by hand

For $r=2$ the nine points of $I_2$ are the origin, the axis points $(\pm1,0)$, $(0,\pm1)$ and the diagonal points $(\pm1,\pm1)$. A triangle with a corner at the origin has the origin on its boundary, so from now on the points are the $N$ points of $I_r$ other than the origin; $N=8$ for $r=2$. Three points on one line enclose nothing and never count either.

Of the $\binom83=56$ triples, where $\binom nk$ counts the ways to choose $k$ of $n$ objects, the eight good ones are the triangles $(1,0)$, $(-1,1)$, $(0,-1)$ and $(1,1)$, $(-1,0)$, $(1,-1)$ and their rotations by $90^\circ$, $180^\circ$ and $270^\circ$, the two shapes drawn in the statement. The other $48$ fail in one of two ways: the corners $(1,0)$, $(1,1)$, $(0,1)$ lie in one quarter of the plane, so the origin is outside, while the triangle $(1,0)$, $(0,1)$, $(-1,0)$ has the origin on its side from $(1,0)$ to $(-1,0)$. In both cases the $x$-axis has all three corners on one side of it or on it.

2. When the origin is inside a triangle

A line divides the plane into two open half-planes; a closed half-plane contains the line as well. ^[2] A triple of points is bad if some line through the origin has all three points in one of its closed half-planes, and good otherwise.

Claim. The origin lies in the interior of the triangle with corners $P,Q,R$ exactly when $\{P,Q,R\}$ is good.

Proof. If all three corners lie in a closed half-plane $K$ of a line $m$ through the origin, then $K$, which contains the segment between any two of its points, ^[2] contains the whole triangle; but points arbitrarily close to the origin lie on the other side of $m$, outside $K$, so the origin is not an interior point. Conversely, suppose the origin is not interior. If $P,Q,R$ lie on one line $\ell$, the parallel to $\ell$ through the origin, or $\ell$ itself, has all of $\ell$ on one closed side. Otherwise the interior consists of the points that lie, for each side of the triangle, strictly on the same side of its line as the opposite corner, ^[3] so for some side, say $QR$, the origin is not strictly on the side of $P$. The parallel $m$ to $QR$ through the origin then lies beyond the line $QR$ or on it, and the closed half-plane of $m$ that contains $P$ also contains the line $QR$. In both cases the triple is bad. $\square$

3. Directions, rays and lines

Seen from the origin, every point $P$ has a direction, its polar angle $\varphi_P$, measured anticlockwise from the positive $x$-axis and taken modulo a full turn $2\pi$. ^[4] A closed half-plane of a line through the origin holds, apart from the origin, exactly the points whose directions lie in a closed half-turn from some $\psi$ to $\psi+\pi$; so a triple is bad exactly when its three directions fit into one closed half-turn. The points with directions strictly between $\varphi_P$ and $\varphi_P+\pi$ are ahead of $P$: an observer at the origin facing $P$ sees them strictly to the left. In coordinates, $Q$ is ahead of $P$ exactly when

$$P\times Q=x_Py_Q-y_Px_Q\gt0,$$

the determinant of the two vectors, which is positive for an anticlockwise turn of less than a half-turn and zero for parallel vectors. ^[5] As $Q\times P=-P\times Q$, two points are never ahead of each other.

Several points can share a ray. Dividing an integer point of a ray by the greatest common divisor (gcd) of its coordinates gives a point $(u,v)$ of the ray with $\gcd(u,v)=1$, and the integer points of the ray are its multiples $k(u,v)$, $k\ge1$: if $(x,y)=t(u,v)$ with a real $t\gt0$, then $t=sx+wy$ for the integers $s,w$ with $su+wv=1$ from Bézout's identity, ^[6] so $t$ is a whole number. The ray thus holds $c$ points of the disk, the number of $k\ge1$ with $k^2(u^2+v^2)\lt r^2$; for $r=3$ this is $2$ for $(1,0)$ and $(1,1)$, and $1$ for $(2,1)$, as $(4,2)$ is too far out. A line through the origin consists of two opposite rays.

4. What one line sees

The half-turn $P\mapsto-P$ about the origin maps the disk and the integer points onto themselves, with three consequences:

$N$ is even, as the map pairs each point with a different one, and we write $H=N/2$.
The two rays of a line $L$ through the origin hold the same number $c$ of points, so $L$ holds $2c$.
The map swaps the two open sides of $L$, so the $N-2c$ points off $L$ split evenly: each open side holds $H-c$ points. In particular, exactly $H-c$ points are ahead of each point of $L$.

For $r=2$ every line has $c=1$, and the $H-c=3$ points ahead of $(1,0)$ are $(1,1)$, $(0,1)$ and $(-1,1)$.

5. Four kinds of bad triples on one line

Fix a line $L$ through the origin with $c$ points on each ray. We attach four kinds of triples to it:

Kind A: a point $P$ of $L$ and two points ahead of $P$, possibly on one ray. With $2c$ choices for $P$ and $\binom{H-c}{2}$ for the pair, there are $2c\binom{H-c}{2}=c(H-c)(H-c-1)$.
Kind B: one point on each ray of $L$ and a third point off $L$, so $c^2(N-2c)$ triples. The side between the first two passes through the origin.
Kind C: two points on one ray of $L$ and a third point ahead of them, so $2\binom c2(H-c)=c(c-1)(H-c)$ triples.
Kind D: three points of $L$, so $\binom{2c}{3}$ triples.

All of them are bad: in kinds A and C all points lie in the closed half-turn that starts at the direction of $P$ or of the doubled ray, in kind B all lie in the closed half-plane of $L$ that holds the third point, and in kind D all lie on $L$.

6. Every bad triple is counted exactly once

Claim. Every bad triple belongs to exactly one line and one kind, and no good triple belongs to any.

Proof. All four kinds are bad, so good triples belong to none. For a bad triple, take the first of these cases that applies.

All three points lie on one line through the origin: kind D for that line, the only line through the origin containing two of them.
Two points lie on opposite rays of a line $L$. The third is off $L$, as case 1 failed, so the triple is of kind B for $L$. A second opposite pair would share a point with the first and lie on $L$ too, so $L$ is unique.
Two points share a ray, and the third point $R$ is off its line. If $R$ is ahead of the two, the triple is of kind C; otherwise the two are ahead of $R$, and it is of kind A with $P=R$.
Three different directions, no two opposite. They fit into a closed half-turn from $\psi$ to $\psi+\pi$; we measure all three directions in this interval and let $P$ be the point that comes first. The other two directions come after $\varphi_P$, not after $\psi+\pi\le\varphi_P+\pi$, and not at $\varphi_P+\pi$, so both points are ahead of $P$: kind A.

No triple has two kinds or two lines. Kinds A, B and C have a point on $L$ and one off $L$, so, unlike kind D, their points never lie on one line through the origin. In kinds A and C all directions lie in a half-turn without its end, so, unlike kind B, no two points are opposite. In kind A, $P$ is the only point with both others ahead of it, as two points are never ahead of each other, which fixes $P$ and its line; in kind C no point has both others ahead, as the two on the shared ray are not ahead of each other and the third is ahead of them. So kinds A and C never overlap, and the shared ray fixes the line of kind C. $\square$

7. The formula

Subtracting the bad triples from all $\binom N3$ triples gives the number of triangles that contain the origin:

$$G(r)=\binom N3-\sum_L f(c_L),$$

$$\begin{aligned}f(c)&=c(H-c)(H-c-1)\\&\quad+c^2(N-2c)\\&\quad+c(c-1)(H-c)+\binom{2c}{3},\end{aligned}$$

where $L$ runs over the lines through the origin that contain points of the disk and $c_L$ is the number of points on each ray of $L$. For $r=2$ there are four lines, the axes and the diagonals, each with $c=1$: kind A gives $1\cdot3\cdot2=6$ triples per line, kind B gives $1\cdot6=6$, kinds C and D give none, so $f(1)=12$ and $G(2)=56-4\cdot12=8$, as in the statement. The two failing examples of Section 1 are of kinds A and B for the $x$-axis.

8. Worked example: $r=3$

The points with $x^2+y^2\lt9$ are exactly those of the square $-2\le x,y\le2$, since there $x^2+y^2\le8$. So $N=24$, $H=12$ and there are $\binom{24}3=2024$ triples. The points lie on eight lines: the two axes and the two diagonals with $c=2$, and the four lines through $(2,1)$, $(1,2)$, $(-1,2)$, $(-2,1)$ with $c=1$. On a line with $c=2$, such as the $x$-axis:

kind A: $2\cdot10\cdot9=180$, as each open side holds $H-c=10$ points;
kind B: $2^2\cdot20=80$, one point on each half of the axis and one of the $20$ points off it;
kind C: $2\cdot1\cdot10=20$, the pair $(1,0),(2,0)$ with a point above the axis, or $(-1,0),(-2,0)$ with a point below it;
kind D: $\binom43=4$, three of the four axis points.

This makes $f(2)=284$. A line with $c=1$ has $1\cdot11\cdot10=110$ triples of kind A, $1\cdot22=22$ of kind B and none of kinds C and D, so $f(1)=132$. Hence

$$\begin{aligned}G(3)&=2024-4\cdot284-4\cdot132\\&=2024-1664=360,\end{aligned}$$

the value in the statement. A program that sorts each of the $2024$ triples into good ones and the four kinds confirms these numbers on every line; in total it finds $1160$ triples of kind A, $408$ of kind B, $80$ of kind C, $16$ of kind D and $360$ good ones. The same formula gives $G(4)=2768$, from $N=44$ points on $16$ lines.

9. Telling rays apart in floating point

The programs find the rays without gcds: they sort the angles in floating point and group angles that differ by at most $\varepsilon=10^{-8}$. Points on one ray have exactly the same angle, so their computed angles differ only by the rounding error of atan2. For Java's Math.atan2 this is at most 2 units in the last place, less than $10^{-15}$, since doubles below $4$ lie at most $2^{-51}$ apart; ^[7] ^[8] C++ and Python call the platform's C library, whose accuracy the language standards leave open. Two points $v,w$ on different, not opposite rays make an angle $\delta$ with $0\lt\delta\lt\pi$. The parallelogram they span has area $|v\times w|=|v|\,|w|\sin\delta$, and $v\times w$ is a nonzero integer, so $|v|\,|w|\sin\delta\ge1$. With $|v|,|w|\lt r$ and $\delta\ge\sin\delta$ this gives

$$\delta\gt\frac1{r^2}.$$

For $r=105$ the bound is about $9.07\cdot10^{-5}$, some $9000$ times larger than $\varepsilon$, and the smallest gap between neighbouring rays there is about $9.33\cdot10^{-5}$. So with any atan2 accurate to far better than $10^{-8}$, two angles at $r=105$ fall into one group exactly when their points share a ray.

How the Code Works

The C++ program ^[9] starts in main: parse_arguments accepts --radius=R, whose digits parse_int_after_prefix reads, with $R\ge2$ required, and --skip-checkpoints. Unless the latter is given, run_checkpoints compares solve(2), solve(3), solve(5) with $8$, $360$, $10600$ and stops the program with exit code $2$ on a mismatch. Then it prints solve(R), with $R=105$ by default.

solve scans the square $-r\le x,y\le r$, skips the origin and the points with $x^2+y^2\ge r^2$, stores std::atan2(y, x), the polar angle in $(-\pi,\pi]$ with $+\pi$ on the negative $x$-axis, ^[4] and sorts, so the points of each ray become neighbours. A guard returns $0$ for fewer than three points. Next it appends a copy of the angles shifted by $2\pi$, as a sweep past the angle $\pi$ would need; nothing reads this copy.

The first half $=H$ sorted angles hold exactly one ray of every line: of two opposite directions exactly one lies in $(-\pi,0]$, so the rays in this interval hold $\sum_L c_L=H$ points (Section 4), all with smaller angles than the rest. The loop over $i\lt H$ thus meets every line once, through its ray in the lower half-plane or on the positive $x$-axis. Its inner while loop counts the angles within kEps $=10^{-8}$ of angles[i], which by Section 9 is the $c$ of that ray, kept in cnt. The total ans starts at $\binom n3=n(n-1)(n-2)/6$ and loses $f(c)$ per group in three statements: kind A first, then kinds B and C together, then kind D as $d(d-1)(d-2)/6$ with doubled $=d=2c$.

In C++ the triangle counts are 64-bit long long values, while radius, the coordinates and total are int; at $r=105$ the largest product, $n(n-1)(n-2)$, is about $4.1\cdot10^{13}$. Beyond $R\ge2$, nothing checks --radius: from $r=818$ on this product exceeds $2^{63}-1$, and from $r=32768$ on the int value $x^2+y^2$ overflows. Python does the same in one function with $r=105$ fixed and no checkpoints, with the unused list angles2 as its copy and unbounded integers. Java sorts a List<Double> with Collections.sort, without the copy or checkpoints; as $4.1\cdot10^{13}$ exceeds the int limit $2^{31}-1$, it casts $n$ to long before multiplying and declares cnt as long.

Complexity and Verification

At $r=105$ the scan visits $(2r+1)^2=44521$ grid points and keeps $N=34608$, ^[10] close to the area $\pi\cdot105^2\approx34636$. Sorting costs $O(N\log N)$ and the counting loop $O(N)$; it reads $H=17304$ angles in $10500$ groups, one per line. As $N\le(2r+1)^2$, the method takes $O(r^2\log r)$ time and $O(r^2)$ memory: $2N$ doubles in C++ because of the unused copy, two lists in Python, $N$ boxed Double objects in Java. On the test machine each program, start-up included, finished in under $0.1$ seconds, whereas testing all triples would take about $7\cdot10^{12}$ point-in-triangle tests.

Unless --skip-checkpoints is given, the C++ program checks the three values of the statement before solving. For this article we also checked, outside the published programs:

a brute force over all triples of $I_r$, the origin included, agrees with Section 7 for every $r$ from $2$ to $30$; its exact integer test asks whether $P\times Q$, $Q\times R$, $R\times P$ are all positive or all negative, as the barycentric coordinates of the origin are proportional to $Q\times R$, $R\times P$, $P\times Q$; ^[3]
for $2\le r\le7$, sorting every triple into the kinds of Section 5 confirms the claim of Section 6 and, on every line, the four terms of $f(c)$;
an implementation of the formula that finds the rays with gcds agrees with the C++ program for every $r$ from $2$ to $104$ and with all three published programs at $r=105$;
at $r=105$ the angles on each ray are bit-for-bit identical within each language (Java's values differ from those of C++ and Python by one unit in the last place for $6020$ of the $34608$ points), and in all three languages every group of the inner loop is exactly one gcd ray.

Footnotes and References

The references below supply the background facts used above. The four kinds of bad triples, the proof that they count every bad triple once and the floating-point bound are derived in this article.

Project Euler 184 — Triangles Containing the Origin. The official statement defines I_r by x² + y² < r², lists the nine points of I_2, gives 8, 360 and 10600 triangles with the origin in their interior for I_2, I_3 and I_5, and asks for I_105.
Wikipedia — Half-space (geometry). In the plane a half-space is a half-plane, open or closed; a closed one is an open one together with its boundary line, and every half-space is convex.
Wikipedia — Barycentric coordinate system. Barycentric coordinates with respect to a triangle are ratios of signed areas; a point is interior when all three lie strictly between 0 and 1, and a coordinate is negative beyond the edge line opposite its vertex.
Wikipedia — Atan2. atan2(y, x) is the angle between the positive x-axis and the ray from the origin to (x, y), in the interval (−π, π]; for y = 0 and x < 0 it is π.
Wikipedia — Determinant. For a 2 × 2 matrix, |ad − bc| is the area of the parallelogram spanned by the columns; the sign is positive for an anticlockwise rotation from the first column to the second and zero for parallel columns.
Wikipedia — Bézout's identity. For integers a and b with greatest common divisor d there are integers x and y with ax + by = d.
Wikipedia — Double-precision floating-point format. A binary64 number has a 53-bit significand, and representable numbers between 2ⁿ and 2ⁿ⁺¹ are 2ⁿ⁻⁵² apart.
Java — Math.atan2 requires a result within 2 ulps (units in the last place) of the exact one; Python — math consists mostly of thin wrappers around the platform C math library.
C++ — Euler184.cpp. The linked immutable C++ revision contains parse_arguments, solve and run_checkpoints. It is the implementation record for this article.
OEIS A051132. Number of integer pairs (x, y) with x² + y² < n², the origin included: 0, 1, 9, 25, 45, 69, …, and a(105) = 34609 in its b-file.

Problem 184 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = long long;

struct Options {
    int radius = 105;
    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, "--radius=", options.radius)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.radius >= 2;
}

i64 solve(const int radius) {
    constexpr double kPi = 3.141592653589793238462643383279502884;
    constexpr double kEps = 1e-8;

    std::vector<double> angles;
    for (int x = -radius; x <= radius; ++x) {
        for (int y = -radius; y <= radius; ++y) {
            if ((x == 0 && y == 0) || x * x + y * y >= radius * radius) {
                continue;
            }
            angles.push_back(std::atan2(static_cast<double>(y), static_cast<double>(x)));
        }
    }

    std::sort(angles.begin(), angles.end());
    const int total = static_cast<int>(angles.size());
    if (total < 3) {
        return 0;
    }

    angles.reserve(static_cast<std::size_t>(2 * total));
    for (int i = 0; i < total; ++i) {
        angles.push_back(angles[static_cast<std::size_t>(i)] + 2.0 * kPi);
    }

    const i64 n = total;
    const i64 half = n / 2;
    i64 ans = n * (n - 1) * (n - 2) / 6;

    for (i64 i = 0; i < half;) {
        i64 cnt = 1;
        while (i + cnt < half &&
               angles[static_cast<std::size_t>(i + cnt)] - angles[static_cast<std::size_t>(i)] <= kEps) {
            ++cnt;
        }

        ans -= (half - cnt) * (half - cnt - 1) * cnt;
        ans -= cnt * cnt * (n - 2 * cnt) + cnt * (cnt - 1) * (half - cnt);

        const i64 doubled = 2 * cnt;
        ans -= doubled * (doubled - 1) * (doubled - 2) / 6;
        i += cnt;
    }

    return ans;
}

bool run_checkpoints() {
    if (solve(2) != 8LL) {
        std::cerr << "Checkpoint failed for r=2" << '\n';
        return false;
    }
    if (solve(3) != 360LL) {
        std::cerr << "Checkpoint failed for r=3" << '\n';
        return false;
    }
    if (solve(5) != 10600LL) {
        std::cerr << "Checkpoint failed for r=5" << '\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;
    }

    std::cout << solve(options.radius) << '\n';
    return 0;
}

Python

# Problem 184: Triangles containing the origin in a circle of radius 105.
import math

def solve():
    R = 105
    EPS = 1e-8
    PI = math.pi
    angles = []
    for x in range(-R, R+1):
        for y in range(-R, R+1):
            if (x==0 and y==0) or x*x+y*y >= R*R: continue
            angles.append(math.atan2(y, x))
    angles.sort()
    n = len(angles)
    # Duplicate for circular handling
    angles2 = angles + [a + 2*PI for a in angles]
    # Count triangles containing origin
    ans = n*(n-1)*(n-2)//6
    half = n // 2
    # Subtract collinear/half-plane cases
    i = 0
    while i < half:
        cnt = 1
        while i+cnt < half and angles[i+cnt] - angles[i] <= EPS:
            cnt += 1
        ans -= (half-cnt)*(half-cnt-1)*cnt
        ans -= cnt*cnt*(n-2*cnt) + cnt*(cnt-1)*(half-cnt)
        doubled = 2*cnt
        ans -= doubled*(doubled-1)*(doubled-2)//6
        i += cnt
    print(ans)

solve()

Java

import java.util.*;

public class Euler184 {
    public static void main(String[] args) {
        int R = 105;
        double PI = Math.PI, EPS = 1e-8;
        List<Double> angles = new ArrayList<>();
        for (int x = -R; x <= R; x++)
            for (int y = -R; y <= R; y++) {
                if ((x == 0 && y == 0) || x * x + y * y >= R * R)
                    continue;
                angles.add(Math.atan2(y, x));
            }
        Collections.sort(angles);
        int n = angles.size();
        int half = n / 2;
        long ans = (long) n * (n - 1) * (n - 2) / 6;
        int i = 0;
        while (i < half) {
            long cnt = 1;
            while (i + cnt < half && angles.get((int) (i + cnt)) - angles.get(i) <= EPS)
                cnt++;
            ans -= (half - cnt) * (half - cnt - 1) * cnt;
            ans -= cnt * cnt * (n - 2 * cnt) + cnt * (cnt - 1) * (half - cnt);
            long d = 2 * cnt;
            ans -= d * (d - 1) * (d - 2) / 6;
            i += (int) cnt;
        }
        System.out.println(ans);
    }
}