Project: A Small, Robust 2D Geometry Toolkit

Assemble the course's algorithms into one tested C++ program with a single coordinate contract, check every part against an independent reference, and learn when to reach for CGAL or Boost.Geometry instead.

Twelve lessons, one program

Each lesson in this course built one algorithm in its own self-contained program. That was the right way to learn them, and it hides the problem you meet when you put them together: the pieces were written with slightly different assumptions. The polygon lesson promised coordinates up to 10^6; the orientation test is exact up to 10^9; the in-circle test needs __int128 and stops at 5 * 10^8. The convex hull drops collinear points, while the triangulation must keep every point. Point-in-polygon reports "on the boundary" as a third answer; a library you call next to it may report false there. Each piece is correct on its own. Assembled carelessly, they disagree with each other on exactly the degenerate inputs that the course spent so much time on.

This lesson assembles a small 2D geometry toolkit in one file: exact predicates, convex hull, polygon area and point-in-polygon, segment intersection, k-d tree nearest neighbor, and Delaunay triangulation by edge flips. It comes with a test harness that checks every part against an independent reference on degenerate inputs, and with a picture that shows several parts working together. Then it looks at the libraries you would use instead of your own code (CGAL and Boost.Geometry), and ends with a checklist for writing geometry code that survives production data.

The full program is long, about 570 lines, and much of that is tests. On a first read, take in the contract in the next section, skim the algorithm parts (you have seen each of them), read the test harness carefully, and run the complete program.

The contract: one coordinate model, checked at the door

The toolkit makes one promise about its input and checks it in one place: integer coordinates with |x|, |y| <= kMaxCoord, with kMaxCoord = 10^6, and polygons of at most 10^6 vertices. That is the smallest of the bounds the earlier lessons needed, so every predicate is exact inside it:

cross and dist2 reach at most 8 * kMaxCoord^2, far below the long long limit (lesson 1).
twiceArea sums up to kMaxVertices terms of at most 2 * kMaxCoord^2 each, 2 * 10^18 in total, which still fits (lesson 5).
inCircle multiplies four coordinate differences; lesson 7 bounds its determinant by 192 * kMaxCoord^4, here 1.92 * 10^26, which __int128 holds with room to spare. Lesson 7's own bound, 5 * 10^8, would also do. __int128 is a GCC and Clang extension, not standard C++ (lesson 2).

Each of these is a static_assert, so raising kMaxCoord past what the arithmetic supports stops the build instead of producing wrong signs. If your data are real numbers, scale them to this grid once, as lesson 2 describes, and accept that the rounding happens there and only there.

The compile-time checks protect the code; they do not protect it from the data. That job belongs to the function that reads input: every point is checked against the bound before any predicate sees it, and a bad line is reported with its number. The program below is that function, fed from strings; on a real system you pass a file stream or std::cin.

main.cppc++23

#include <cstdlib>
#include <expected>
#include <format>
#include <limits>
#include <print>
#include <sstream>
#include <string>
#include <vector>

struct Point { long long x, y; };

constexpr long long kMaxCoord = 1'000'000;
constexpr long long kMaxVertices = 1'000'000;
constexpr long long kMaxLL = std::numeric_limits<long long>::max();
static_assert(8 * kMaxCoord * kMaxCoord < kMaxLL);                 // cross, dist2
static_assert(2 * kMaxCoord * kMaxCoord <= kMaxLL / kMaxVertices);  // twiceArea
using Wide = __int128;  // a GCC/Clang extension, not standard C++
constexpr Wide kWideMax = Wide((static_cast<unsigned __int128>(1) << 127) - 1);
static_assert(192 * Wide(kMaxCoord) * kMaxCoord * kMaxCoord * kMaxCoord <= kWideMax,
              "in-circle determinant must fit in __int128");

bool inRange(Point p) { return std::abs(p.x) <= kMaxCoord && std::abs(p.y) <= kMaxCoord; }

// One "x y" pair per line. Everything the predicates cannot handle exactly is
// rejected here, with a line number, instead of failing somewhere deep inside.
std::expected<std::vector<Point>, std::string> readPoints(std::istream& in) {
    std::vector<Point> pts;
    std::string line;
    for (int number = 1; std::getline(in, line); ++number) {
        std::istringstream fields(line);
        Point p{};
        std::string extra;
        if (!(fields >> p.x >> p.y) || fields >> extra)
            return std::unexpected(std::format("line {}: expected two integers", number));
        if (!inRange(p))
            return std::unexpected(std::format("line {}: ({}, {}) is out of range",
                                               number, p.x, p.y));
        pts.push_back(p);
    }
    if (pts.size() > static_cast<std::size_t>(kMaxVertices))
        return std::unexpected("too many points");
    return pts;
}

int main() {
    // On a real system this would be a file or std::cin.
    for (const char* text : {"0 0\n10 0\n5 7\n", "0 0\n2500000 1\n", "0 0\n1.5 2\n",
                             "0 0\n99999999999999999999 3\n"}) {
        std::istringstream in(text);
        auto pts = readPoints(in);
        if (pts) std::println("ok: {} points", pts->size());
        else std::println("rejected: {}", pts.error());
    }
}

Output

Press Run to compile and run this example.

ok: 3 points
rejected: line 2: (2500000, 1) is out of range
rejected: line 2: expected two integers
rejected: line 2: expected two integers

The last two inputs fail at parsing: 1.5 is not an integer, and 99999999999999999999 does not fit in a long long, so the stream sets its fail state. The second one parses but is out of range. In all three cases the problem is reported where it can be explained, not three function calls deep as an overflow or a wrong sign.

The algorithms, part by part

The toolkit has six parts. The first four, shown below, are the library; the fifth is the test harness, in the next section; the sixth, the drawing and main, appears only in the complete program. Each library part repeats, unchanged or nearly so, code from an earlier lesson, so the commentary is short.

Part 1: the contract and the predicates. Everything else calls these functions and nothing else does arithmetic on coordinates, except the shoelace sum and the k-d tree's distance bounds. operator==, lexLess and inBox come from lesson 3, onSegment from lesson 5.

// ---- Part 1: the coordinate model and the exact predicates ----------------
#include <algorithm>
#include <array>
#include <cstdlib>
#include <format>
#include <limits>
#include <map>
#include <numeric>
#include <print>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

struct Point { long long x, y; };

// The whole toolkit promises exact answers for |x|, |y| <= kMaxCoord and
// polygons of at most kMaxVertices vertices. Check input against it once.
constexpr long long kMaxCoord = 1'000'000;
constexpr long long kMaxVertices = 1'000'000;
constexpr long long kMaxLL = std::numeric_limits<long long>::max();
static_assert(8 * kMaxCoord * kMaxCoord < kMaxLL);                 // cross, dist2
static_assert(2 * kMaxCoord * kMaxCoord <= kMaxLL / kMaxVertices);  // twiceArea
using Wide = __int128;  // a GCC/Clang extension, not standard C++
constexpr Wide kWideMax = Wide((static_cast<unsigned __int128>(1) << 127) - 1);
static_assert(192 * Wide(kMaxCoord) * kMaxCoord * kMaxCoord * kMaxCoord <= kWideMax,
              "in-circle determinant must fit in __int128");

bool inRange(Point p) { return std::abs(p.x) <= kMaxCoord && std::abs(p.y) <= kMaxCoord; }
bool operator==(Point p, Point q) { return p.x == q.x && p.y == q.y; }
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }

long long cross(Point a, Point b, Point c) {  // twice the signed area of abc
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {  // -1 right, 0 collinear, +1 left
    long long v = cross(a, b, c);
    return (v > 0) - (v < 0);
}
// > 0 if d is strictly inside the circle through a, b, c (counter-clockwise),
// < 0 if strictly outside, 0 if on it. Exact for coordinates within kMaxCoord.
int inCircle(Point a, Point b, Point c, Point d) {
    long long adx = a.x - d.x, ady = a.y - d.y;
    long long bdx = b.x - d.x, bdy = b.y - d.y;
    long long cdx = c.x - d.x, cdy = c.y - d.y;
    long long aLift = adx * adx + ady * ady;
    long long bLift = bdx * bdx + bdy * bdy;
    long long cLift = cdx * cdx + cdy * cdy;
    Wide det = Wide(aLift) * (bdx * cdy - cdx * bdy)
             + Wide(bLift) * (cdx * ady - adx * cdy)
             + Wide(cLift) * (adx * bdy - bdx * ady);
    return (det > 0) - (det < 0);
}
long long dist2(Point a, Point b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
// For p collinear with a and b: is p within the segment's bounding box?
bool inBox(Point a, Point b, Point p) {
    return std::min(a.x, b.x) <= p.x && p.x <= std::max(a.x, b.x) &&
           std::min(a.y, b.y) <= p.y && p.y <= std::max(a.y, b.y);
}
bool onSegment(Point a, Point b, Point p) {
    return cross(a, b, p) == 0 && inBox(a, b, p);
}

Part 2: hull, polygons and segments. convexHull is Andrew's monotone chain from lesson 4, unchanged: counter-clockwise from the lexicographically smallest point, duplicates removed, and a keepCollinear parameter for the collinear policy (dropped by default). twiceArea and crossingTest come from lesson 5, with the three-valued Where result; intersects is lesson 3's closed-segment test. Skip this part on a first read.

// ---- Part 2: convex hull, polygons, segments --------------------------------
// Andrew's monotone chain. Returns the hull counter-clockwise, starting at the
// lexicographically smallest point, without repeating it at the end.
std::vector<Point> convexHull(std::vector<Point> pts, bool keepCollinear = false) {
    std::ranges::sort(pts, lexLess);
    auto dups = std::ranges::unique(pts);
    pts.erase(dups.begin(), dups.end());
    if (pts.size() < 3) return pts;
    auto onLine = [&](Point p) { return orientation(pts.front(), pts.back(), p) == 0; };
    if (std::ranges::all_of(pts, onLine))  // no area: the hull is a segment
        return keepCollinear ? pts : std::vector{pts.front(), pts.back()};
    std::vector<Point> hull;
    for (int pass = 0; pass < 2; ++pass) {  // lower hull, then upper hull
        const size_t start = hull.size();
        for (Point p : pts) {
            while (hull.size() >= start + 2) {
                int turn = orientation(hull[hull.size() - 2], hull.back(), p);
                if (turn > 0 || (turn == 0 && keepCollinear)) break;
                hull.pop_back();  // not a left turn: the middle point is not a vertex
            }
            hull.push_back(p);
        }
        hull.pop_back();  // the last point of this chain starts the next one
        std::ranges::reverse(pts);
    }
    return hull;
}

// Twice the signed area: positive for counter-clockwise polygons.
long long twiceArea(const std::vector<Point>& poly) {
    long long sum = 0;
    for (std::size_t i = 0, n = poly.size(); i < n; ++i) {
        Point p = poly[i], q = poly[(i + 1) % n];
        sum += p.x * q.y - q.x * p.y;
    }
    return sum;
}

enum class Where { Outside, Inside, Boundary };

// Crossing number with the half-open rule; the boundary is reported as such.
Where crossingTest(const std::vector<Point>& poly, Point p) {
    bool inside = false;
    for (std::size_t i = 0, n = poly.size(); i < n; ++i) {
        Point a = poly[i], b = poly[(i + 1) % n];
        if (onSegment(a, b, p)) return Where::Boundary;
        if ((a.y > p.y) != (b.y > p.y)) {
            int side = orientation(a, b, p);
            if (b.y > a.y ? side > 0 : side < 0) inside = !inside;
        }
    }
    return inside ? Where::Inside : Where::Outside;
}

struct Segment { Point a, b; };

bool intersects(Segment s, Segment t) {  // closed segments: touching counts
    int d1 = orientation(t.a, t.b, s.a), d2 = orientation(t.a, t.b, s.b);
    int d3 = orientation(s.a, s.b, t.a), d4 = orientation(s.a, s.b, t.b);
    if (d1 * d2 < 0 && d3 * d4 < 0) return true;
    return (d1 == 0 && inBox(t.a, t.b, s.a)) || (d2 == 0 && inBox(t.a, t.b, s.b)) ||
           (d3 == 0 && inBox(s.a, s.b, t.a)) || (d4 == 0 && inBox(s.a, s.b, t.b));
}

Part 3: the k-d tree. Lesson 8's tree in one array, with the median of each range as the node and the cell test for pruning. The one change: nearest returns the point, not only its squared distance, because the drawing needs it. Skip this part on a first read too.

// ---- Part 3: k-d tree for nearest-neighbor queries ---------------------------
bool lessOn(int axis, Point a, Point b) {
    return axis == 0 ? (a.x != b.x ? a.x < b.x : a.y < b.y)
                     : (a.y != b.y ? a.y < b.y : a.x < b.x);
}

struct Box { long long xlo, ylo, xhi, yhi; };

long long boxDist2(Point q, Box b) {  // from q to the closest point of b
    long long dx = std::max({b.xlo - q.x, 0LL, q.x - b.xhi});
    long long dy = std::max({b.ylo - q.y, 0LL, q.y - b.yhi});
    return dx * dx + dy * dy;
}

// Implicit tree in one array: the median of each range is its root.
class KdTree {
public:
    explicit KdTree(std::vector<Point> p) : pts_(std::move(p)) {
        build(0, pts_.size(), 0);
    }

    // A nearest point to q (the tree must not be empty).
    Point nearest(Point q) const {
        long long best = std::numeric_limits<long long>::max();
        Point found = pts_.front();
        Box world{-kMaxCoord, -kMaxCoord, kMaxCoord, kMaxCoord};
        search(0, pts_.size(), 0, world, q, best, found);
        return found;
    }

private:
    void build(std::size_t lo, std::size_t hi, int axis) {
        if (hi - lo <= 1) return;
        std::size_t mid = lo + (hi - lo) / 2;
        std::nth_element(pts_.begin() + lo, pts_.begin() + mid, pts_.begin() + hi,
                         [axis](Point a, Point b) { return lessOn(axis, a, b); });
        build(lo, mid, 1 - axis);
        build(mid + 1, hi, 1 - axis);
    }
    void search(std::size_t lo, std::size_t hi, int axis, Box cell, Point q,
                long long& best, Point& found) const {
        if (lo >= hi || boxDist2(q, cell) >= best) return;  // nothing here can win
        std::size_t mid = lo + (hi - lo) / 2;
        Point p = pts_[mid];
        if (dist2(p, q) < best) {
            best = dist2(p, q);
            found = p;
        }
        Box left = cell, right = cell;
        if (axis == 0) left.xhi = right.xlo = p.x;
        else left.yhi = right.ylo = p.y;
        // Search the side that contains q first: it usually shrinks best quickly.
        if (lessOn(axis, q, p)) {
            search(lo, mid, 1 - axis, left, q, best, found);
            search(mid + 1, hi, 1 - axis, right, q, best, found);
        } else {
            search(mid + 1, hi, 1 - axis, right, q, best, found);
            search(lo, mid, 1 - axis, left, q, best, found);
        }
    }
    std::vector<Point> pts_;
};

Part 4: Delaunay by flips. The monotone-chain sweep produces some triangulation, the half-edge mesh from lesson 11 stores it, and Lawson's flips turn it into the Delaunay triangulation. triangles lists each face once. The input must not contain duplicate points.

One addition guards against a failure the earlier lessons did not show. With an exact inCircle, Lawson's algorithm performs at most n(n-1)/2 flips. When we planted an off-by-one bug in inCircle while testing this toolkit, the flips never stopped: the program ran until we killed it. So delaunay counts its flips and throws std::logic_error when it exceeds the bound. A loop whose termination depends on a predicate being right should fail loudly, not hang.

// ---- Part 4: Delaunay triangulation by edge flips (half-edge mesh) ----------
using Triangle = std::array<int, 3>;  // vertex indices, counter-clockwise

// Monotone-chain sweep: every point popped off a hull chain becomes a triangle.
std::vector<Triangle> sweepTriangulate(const std::vector<Point>& p) {
    std::vector<int> order(p.size());
    std::iota(order.begin(), order.end(), 0);
    std::ranges::sort(order, [&](int i, int j) { return lexLess(p[i], p[j]); });
    std::vector<Triangle> tris;
    std::vector<int> lo, up;
    for (int q : order) {
        while (lo.size() > 1 && orientation(p[lo.end()[-2]], p[lo.back()], p[q]) < 0) {
            tris.push_back({lo.end()[-2], q, lo.back()});
            lo.pop_back();
        }
        while (up.size() > 1 && orientation(p[up.end()[-2]], p[up.back()], p[q]) > 0) {
            tris.push_back({up.end()[-2], up.back(), q});
            up.pop_back();
        }
        lo.push_back(q);
        up.push_back(q);
    }
    return tris;
}

struct HalfEdge { int origin, twin, next, face; };  // face -1: outside

struct Mesh {
    std::vector<HalfEdge> he;
    std::vector<int> out;  // per vertex: one outgoing half-edge
    int dest(int h) const { return he[he[h].twin].origin; }
    int prev(int h) const { return he[he[h].next].next; }  // triangles only
};

Mesh build(int vertexCount, const std::vector<Triangle>& tris) {
    Mesh m;
    m.out.assign(vertexCount, -1);
    std::map<std::pair<int, int>, int> directed;
    for (int f = 0; f < static_cast<int>(tris.size()); ++f)
        for (int k = 0; k < 3; ++k) {
            int a = tris[f][k], b = tris[f][(k + 1) % 3];
            directed[{a, b}] = 3 * f + k;
            m.he.push_back({a, -1, 3 * f + (k + 1) % 3, f});
            m.out[a] = 3 * f + k;
        }
    std::map<int, int> boundaryOut;
    const int interior = static_cast<int>(m.he.size());
    for (int h = 0; h < interior; ++h) {
        int a = m.he[h].origin, b = m.he[m.he[h].next].origin;
        if (auto it = directed.find({b, a}); it != directed.end()) {
            m.he[h].twin = it->second;
            continue;
        }
        int g = static_cast<int>(m.he.size());
        m.he.push_back({b, h, -1, -1});
        m.he[h].twin = g;
        m.out[b] = boundaryOut[b] = g;
    }
    for (int g = interior; g < static_cast<int>(m.he.size()); ++g)
        m.he[g].next = boundaryOut.at(m.dest(g));
    return m;
}

// Replace edge a-b (triangles a b c and b a d) by edge c-d.
void flip(Mesh& m, int h) {
    int t = m.he[h].twin;
    int h1 = m.he[h].next, h2 = m.he[h1].next;  // b->c, c->a
    int t1 = m.he[t].next, t2 = m.he[t1].next;  // a->d, d->b
    int a = m.he[h].origin, b = m.he[t].origin;
    int c = m.he[h2].origin, d = m.he[t2].origin;
    int f0 = m.he[h].face, f1 = m.he[t].face;
    m.he[h] = {d, t, h2, f0};
    m.he[h2].next = t1;
    m.he[t1].next = h;
    m.he[t1].face = f0;
    m.he[t] = {c, h, t2, f1};
    m.he[t2].next = h1;
    m.he[h1].next = t;
    m.he[h1].face = f1;
    if (m.out[a] == h) m.out[a] = t1;
    if (m.out[b] == t) m.out[b] = h1;
}

bool illegal(const Mesh& m, const std::vector<Point>& p, int h) {
    int t = m.he[h].twin;
    if (m.he[h].face < 0 || m.he[t].face < 0) return false;
    int a = m.he[h].origin, b = m.he[t].origin;
    int c = m.he[m.prev(h)].origin, d = m.he[m.prev(t)].origin;
    return inCircle(p[a], p[b], p[c], p[d]) > 0;
}

// The Delaunay triangulation of distinct points: sweep, then Lawson's flips.
Mesh delaunay(const std::vector<Point>& p) {
    Mesh m = build(static_cast<int>(p.size()), sweepTriangulate(p));
    std::vector<int> stack(m.he.size());
    std::iota(stack.begin(), stack.end(), 0);
    // With an exact inCircle, Lawson needs at most n(n-1)/2 flips (lesson 11).
    // More means a broken predicate, which would otherwise loop forever.
    long long flipsLeft = static_cast<long long>(p.size()) * (p.size() - 1) / 2;
    while (!stack.empty()) {
        int h = stack.back();
        stack.pop_back();
        if (!illegal(m, p, h)) continue;
        if (flipsLeft-- == 0) throw std::logic_error("delaunay: too many flips");
        flip(m, h);
        int t = m.he[h].twin;
        for (int g : {m.he[h].next, m.prev(h), m.he[t].next, m.prev(t)})
            stack.push_back(g);
    }
    return m;
}

// The triangles of a mesh, each once, as vertex index triples.
std::vector<Triangle> triangles(const Mesh& m) {
    std::vector<Triangle> tris;
    for (int h = 0; h < static_cast<int>(m.he.size()); ++h)
        if (m.he[h].face >= 0 && h < m.he[h].next && h < m.prev(h))
            tris.push_back({m.he[h].origin, m.dest(h), m.he[m.prev(h)].origin});
    return tris;
}

A test harness that does not trust the code it tests

A test that recomputes the answer with the same idea shares the same bugs. Each test in the harness uses a different formulation of the question, and each runs on inputs built to be degenerate: coordinates between 0 and 2 or 0 and 4 in a quarter of the cases each, so that duplicates, collinear triples and cocircular quadruples are common rather than rare.

Part	Independent reference
`inCircle`	distance to the circumcenter, built exactly as a fraction
`convexHull`	the definition: strictly convex, contains every input point, vertices are input points
`crossingTest`	star-shaped polygons, tiled by a fan of triangles from their center
`intersects`	Cramer's rule on the parametric form, on every pair over a 4 by 4 grid
`KdTree`	brute force over all points, comparing distances, not points
`delaunay`	empty circles by brute force; boundary and triangle count against `convexHull`; area equal to the hull's; nearest neighbors among Delaunay neighbors

The hull test deserves a sentence. Instead of comparing with a second hull algorithm, it checks the three properties that define the hull: a polygon that turns strictly left at every vertex, contains every input point, and uses only input points as vertices is the convex hull, and nothing else is. The Delaunay test also crosses module boundaries: the number of outside half-edges must equal convexHull(p, true).size(), the triangle count must be 2n - 2 minus that number, the triangles' areas must add up to the area of convexHull(p), and every point's nearest neighbor, found by brute force, must be at the same distance as its nearest Delaunay neighbor, found by walking the one-ring in the half-edge mesh. That second check exercises the per-vertex handles that flip must keep up to date. The run helper counts an exception as a failed case, so the flip guard in delaunay shows up as a failure count instead of ending the run.

// ---- Part 5: property tests against independent references ------------------
struct Lcg {
    unsigned long long s;
    unsigned next() {
        s = s * 6364136223846793005ULL + 1442695040888963407ULL;
        return unsigned(s >> 33);
    }
    int range(int lo, int hi) {  // inclusive
        return lo + int(next() % unsigned(hi - lo + 1));
    }
};

std::vector<Point> randomPoints(Lcg& rng, int n, int range, bool distinct) {
    std::vector<Point> p;
    for (int i = 0; i < n; ++i) p.push_back({rng.range(0, range), rng.range(0, range)});
    if (distinct) {
        std::ranges::sort(p, lexLess);
        auto dups = std::ranges::unique(p);
        p.erase(dups.begin(), dups.end());
    }
    return p;
}

// Small ranges make ties, duplicates and collinear points common.
int pickRange(int i) { return std::array{2, 4, 10, 1000}[i % 4]; }

int sign(Wide v) { return (v > 0) - (v < 0); }

// In-circle from the circumcenter, built exactly as a fraction U / D.
int inCircleByCenter(Point a, Point b, Point c, Point d) {
    Wide bx = b.x - a.x, by = b.y - a.y, cx = c.x - a.x, cy = c.y - a.y;
    Wide D = 2 * (bx * cy - by * cx);  // nonzero: abc is not collinear
    Wide ux = cy * (bx * bx + by * by) - by * (cx * cx + cy * cy);
    Wide uy = bx * (cx * cx + cy * cy) - cx * (bx * bx + by * by);
    Wide dx = D * (d.x - a.x) - ux, dy = D * (d.y - a.y) - uy;
    return sign(ux * ux + uy * uy - dx * dx - dy * dy) * sign(D);  // a is on the circle
}

bool testPredicates(Lcg& rng, int i) {
    auto p = randomPoints(rng, 4, pickRange(i), false);
    Point a = p[0], b = p[1], c = p[2], d = p[3];
    if (orientation(a, b, c) != -orientation(b, a, c)) return false;
    if (orientation(a, b, c) != orientation(b, c, a)) return false;
    if (orientation(a, b, c) == 0) return true;  // no circle through a, b, c
    return inCircle(a, b, c, d) == inCircleByCenter(a, b, c, d);
}

// A convex polygon of input points that turns strictly left everywhere and has
// every input point inside or on it is, by definition, the convex hull.
bool testHull(Lcg& rng, int i) {
    auto p = randomPoints(rng, rng.range(0, 12), pickRange(i), false);
    auto h = convexHull(p);
    const std::size_t k = h.size();
    for (Point v : h)
        if (std::ranges::find(p, v) == p.end()) return false;
    if (k == 0) return p.empty();
    if (k == 1) return std::ranges::all_of(p, [&](Point q) { return q == h[0]; });
    if (k == 2)
        return h[0] != h[1] &&
               std::ranges::all_of(p, [&](Point q) { return onSegment(h[0], h[1], q); });
    for (std::size_t j = 0; j < k; ++j) {
        if (orientation(h[j], h[(j + 1) % k], h[(j + 2) % k]) <= 0) return false;
        for (Point q : p)
            if (orientation(h[j], h[(j + 1) % k], q) < 0) return false;
    }
    return true;
}

// A star-shaped polygon around (0, 0); the reference answer comes from its fan
// of triangles, which tile it exactly.
bool testPointInPolygon(Lcg& rng, int i) {
    const Point dirs[] = {{1, 0}, {2, 1}, {1, 1}, {1, 2}, {0, 1}, {-1, 2},
        {-1, 1}, {-2, 1}, {-1, 0}, {-2, -1}, {-1, -1}, {-1, -2}, {0, -1}, {1, -2},
        {1, -1}, {2, -1}};  // 16 directions in counter-clockwise order
    std::vector<Point> poly;
    for (int j = 0; j < 16; ++j) {
        int r = rng.range(1, 1 + i % 3);
        if (j % 4 == 0 || rng.range(0, 1))  // keep the axes: gaps stay below 180 degrees
            poly.push_back({r * dirs[j].x, r * dirs[j].y});
    }
    const std::size_t n = poly.size();
    for (long long x = -7; x <= 7; ++x)
        for (long long y = -7; y <= 7; ++y) {
            Point q{x, y};
            Where want = Where::Outside;
            for (std::size_t j = 0; j < n; ++j) {
                Point a = poly[j], b = poly[(j + 1) % n];
                if (onSegment(a, b, q)) { want = Where::Boundary; break; }
                if (orientation({0, 0}, a, q) >= 0 && orientation(a, b, q) >= 0 &&
                    orientation(b, {0, 0}, q) >= 0)
                    want = Where::Inside;
            }
            if (crossingTest(poly, q) != want) return false;
        }
    return true;
}

// Segments by Cramer's rule: a + t (b - a) = c + u (d - c) with t, u in [0, 1].
bool intersectByParameters(Point a, Point b, Point c, Point d) {
    long long rx = b.x - a.x, ry = b.y - a.y, sx = d.x - c.x, sy = d.y - c.y;
    long long qx = c.x - a.x, qy = c.y - a.y;
    long long den = rx * sy - ry * sx;
    if (den != 0) {
        long long t = qx * sy - qy * sx, u = qx * ry - qy * rx;
        if (den < 0) den = -den, t = -t, u = -u;
        return 0 <= t && t <= den && 0 <= u && u <= den;
    }
    if (rx == 0 && ry == 0)  // ab is a single point: swap the roles
        return sx == 0 && sy == 0 ? a == c : intersectByParameters(c, d, a, b);
    if (qx * ry - qy * rx != 0) return false;  // parallel, on different lines
    long long t0 = qx * rx + qy * ry, t1 = t0 + sx * rx + sy * ry;  // project onto r
    long long len2 = rx * rx + ry * ry;  // ab covers [0, len2] on that axis
    return std::max(std::min(t0, t1), 0LL) <= std::min(std::max(t0, t1), len2);
}

bool testNearest(Lcg& rng, int i) {
    auto p = randomPoints(rng, rng.range(1, 40), pickRange(i), false);
    KdTree tree(p);
    for (int k = 0; k < 10; ++k) {
        Point q{rng.range(-2, pickRange(i) + 2), rng.range(-2, pickRange(i) + 2)};
        long long brute = std::numeric_limits<long long>::max();
        for (Point v : p) brute = std::min(brute, dist2(v, q));
        if (dist2(tree.nearest(q), q) != brute) return false;  // ties: compare distances
    }
    return true;
}

bool testDelaunay(Lcg& rng, int i) {
    auto p = randomPoints(rng, rng.range(1, 30), pickRange(i), true);
    const int n = static_cast<int>(p.size());
    Mesh m = delaunay(p);
    auto tris = triangles(m);
    long long area = 0;
    for (auto [a, b, c] : tris) {
        if (orientation(p[a], p[b], p[c]) <= 0) return false;
        for (Point q : p)
            if (inCircle(p[a], p[b], p[c], q) > 0) return false;  // empty circles
        area += cross(p[a], p[b], p[c]);
    }
    auto corners = convexHull(p), boundary = convexHull(p, true);
    if (corners.size() < 3) return tris.empty();  // all points collinear
    const int t = static_cast<int>(tris.size()), h = static_cast<int>(boundary.size());
    const int outside = static_cast<int>(m.he.size()) - 3 * t;  // outside half-edges
    if (outside != h || t != 2 * n - 2 - h || area != twiceArea(corners)) return false;
    // A nearest neighbor of every point is one of its Delaunay neighbors.
    for (int v = 0; v < n; ++v) {
        long long ring = std::numeric_limits<long long>::max();
        long long all = ring;
        int h = m.out[v];
        do {
            ring = std::min(ring, dist2(p[v], p[m.dest(h)]));
            h = m.he[m.he[h].twin].next;
        } while (h != m.out[v]);
        for (int w = 0; w < n; ++w)
            if (w != v) all = std::min(all, dist2(p[v], p[w]));
        if (ring != all) return false;
    }
    return true;
}

bool runAllTests() {
    Lcg rng{42};
    bool allPassed = true;
    auto run = [&](const char* name, int total, auto test) {
        int passed = 0;
        for (int i = 0; i < total; ++i) {
            try {
                passed += test(i);
            } catch (const std::exception&) {  // a failed check, not a crash
            }
        }
        std::println("{:<30} {:>6}/{}", name, passed, total);
        allPassed = allPassed && passed == total;
    };
    run("orientation and in-circle", 20000,
        [&](int i) { return testPredicates(rng, i); });
    run("convex hull", 4000, [&](int i) { return testHull(rng, i); });
    run("point in polygon (225 each)", 300,
        [&](int i) { return testPointInPolygon(rng, i); });
    run("all segment pairs, 4x4 grid", 65536, [](int i) {  // every pair, exhaustively
        auto g = [](int k) { return Point{k % 4, k / 4}; };
        Point a = g(i & 15), b = g(i >> 4 & 15), c = g(i >> 8 & 15), d = g(i >> 12);
        return intersects({a, b}, {c, d}) == intersectByParameters(a, b, c, d);
    });
    run("k-d tree nearest (10 each)", 2000, [&](int i) { return testNearest(rng, i); });
    run("Delaunay by flips", 2000, [&](int i) { return testDelaunay(rng, i); });
    return allPassed;
}

The complete toolkit

Here is the whole program: the five parts above, then the drawing helper from earlier lessons, an overview drawing, and main. Run it.

main.cppc++23

// ---- Part 1: the coordinate model and the exact predicates ----------------
#include <algorithm>
#include <array>
#include <cstdlib>
#include <format>
#include <limits>
#include <map>
#include <numeric>
#include <print>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

struct Point { long long x, y; };

// The whole toolkit promises exact answers for |x|, |y| <= kMaxCoord and
// polygons of at most kMaxVertices vertices. Check input against it once.
constexpr long long kMaxCoord = 1'000'000;
constexpr long long kMaxVertices = 1'000'000;
constexpr long long kMaxLL = std::numeric_limits<long long>::max();
static_assert(8 * kMaxCoord * kMaxCoord < kMaxLL);                 // cross, dist2
static_assert(2 * kMaxCoord * kMaxCoord <= kMaxLL / kMaxVertices);  // twiceArea
using Wide = __int128;  // a GCC/Clang extension, not standard C++
constexpr Wide kWideMax = Wide((static_cast<unsigned __int128>(1) << 127) - 1);
static_assert(192 * Wide(kMaxCoord) * kMaxCoord * kMaxCoord * kMaxCoord <= kWideMax,
              "in-circle determinant must fit in __int128");

bool inRange(Point p) { return std::abs(p.x) <= kMaxCoord && std::abs(p.y) <= kMaxCoord; }
bool operator==(Point p, Point q) { return p.x == q.x && p.y == q.y; }
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }

long long cross(Point a, Point b, Point c) {  // twice the signed area of abc
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {  // -1 right, 0 collinear, +1 left
    long long v = cross(a, b, c);
    return (v > 0) - (v < 0);
}
// > 0 if d is strictly inside the circle through a, b, c (counter-clockwise),
// < 0 if strictly outside, 0 if on it. Exact for coordinates within kMaxCoord.
int inCircle(Point a, Point b, Point c, Point d) {
    long long adx = a.x - d.x, ady = a.y - d.y;
    long long bdx = b.x - d.x, bdy = b.y - d.y;
    long long cdx = c.x - d.x, cdy = c.y - d.y;
    long long aLift = adx * adx + ady * ady;
    long long bLift = bdx * bdx + bdy * bdy;
    long long cLift = cdx * cdx + cdy * cdy;
    Wide det = Wide(aLift) * (bdx * cdy - cdx * bdy)
             + Wide(bLift) * (cdx * ady - adx * cdy)
             + Wide(cLift) * (adx * bdy - bdx * ady);
    return (det > 0) - (det < 0);
}
long long dist2(Point a, Point b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
// For p collinear with a and b: is p within the segment's bounding box?
bool inBox(Point a, Point b, Point p) {
    return std::min(a.x, b.x) <= p.x && p.x <= std::max(a.x, b.x) &&
           std::min(a.y, b.y) <= p.y && p.y <= std::max(a.y, b.y);
}
bool onSegment(Point a, Point b, Point p) {
    return cross(a, b, p) == 0 && inBox(a, b, p);
}

// ---- Part 2: convex hull, polygons, segments --------------------------------
// Andrew's monotone chain. Returns the hull counter-clockwise, starting at the
// lexicographically smallest point, without repeating it at the end.
std::vector<Point> convexHull(std::vector<Point> pts, bool keepCollinear = false) {
    std::ranges::sort(pts, lexLess);
    auto dups = std::ranges::unique(pts);
    pts.erase(dups.begin(), dups.end());
    if (pts.size() < 3) return pts;
    auto onLine = [&](Point p) { return orientation(pts.front(), pts.back(), p) == 0; };
    if (std::ranges::all_of(pts, onLine))  // no area: the hull is a segment
        return keepCollinear ? pts : std::vector{pts.front(), pts.back()};
    std::vector<Point> hull;
    for (int pass = 0; pass < 2; ++pass) {  // lower hull, then upper hull
        const size_t start = hull.size();
        for (Point p : pts) {
            while (hull.size() >= start + 2) {
                int turn = orientation(hull[hull.size() - 2], hull.back(), p);
                if (turn > 0 || (turn == 0 && keepCollinear)) break;
                hull.pop_back();  // not a left turn: the middle point is not a vertex
            }
            hull.push_back(p);
        }
        hull.pop_back();  // the last point of this chain starts the next one
        std::ranges::reverse(pts);
    }
    return hull;
}

// Twice the signed area: positive for counter-clockwise polygons.
long long twiceArea(const std::vector<Point>& poly) {
    long long sum = 0;
    for (std::size_t i = 0, n = poly.size(); i < n; ++i) {
        Point p = poly[i], q = poly[(i + 1) % n];
        sum += p.x * q.y - q.x * p.y;
    }
    return sum;
}

enum class Where { Outside, Inside, Boundary };

// Crossing number with the half-open rule; the boundary is reported as such.
Where crossingTest(const std::vector<Point>& poly, Point p) {
    bool inside = false;
    for (std::size_t i = 0, n = poly.size(); i < n; ++i) {
        Point a = poly[i], b = poly[(i + 1) % n];
        if (onSegment(a, b, p)) return Where::Boundary;
        if ((a.y > p.y) != (b.y > p.y)) {
            int side = orientation(a, b, p);
            if (b.y > a.y ? side > 0 : side < 0) inside = !inside;
        }
    }
    return inside ? Where::Inside : Where::Outside;
}

struct Segment { Point a, b; };

bool intersects(Segment s, Segment t) {  // closed segments: touching counts
    int d1 = orientation(t.a, t.b, s.a), d2 = orientation(t.a, t.b, s.b);
    int d3 = orientation(s.a, s.b, t.a), d4 = orientation(s.a, s.b, t.b);
    if (d1 * d2 < 0 && d3 * d4 < 0) return true;
    return (d1 == 0 && inBox(t.a, t.b, s.a)) || (d2 == 0 && inBox(t.a, t.b, s.b)) ||
           (d3 == 0 && inBox(s.a, s.b, t.a)) || (d4 == 0 && inBox(s.a, s.b, t.b));
}

// ---- Part 3: k-d tree for nearest-neighbor queries ---------------------------
bool lessOn(int axis, Point a, Point b) {
    return axis == 0 ? (a.x != b.x ? a.x < b.x : a.y < b.y)
                     : (a.y != b.y ? a.y < b.y : a.x < b.x);
}

struct Box { long long xlo, ylo, xhi, yhi; };

long long boxDist2(Point q, Box b) {  // from q to the closest point of b
    long long dx = std::max({b.xlo - q.x, 0LL, q.x - b.xhi});
    long long dy = std::max({b.ylo - q.y, 0LL, q.y - b.yhi});
    return dx * dx + dy * dy;
}

// Implicit tree in one array: the median of each range is its root.
class KdTree {
public:
    explicit KdTree(std::vector<Point> p) : pts_(std::move(p)) {
        build(0, pts_.size(), 0);
    }

    // A nearest point to q (the tree must not be empty).
    Point nearest(Point q) const {
        long long best = std::numeric_limits<long long>::max();
        Point found = pts_.front();
        Box world{-kMaxCoord, -kMaxCoord, kMaxCoord, kMaxCoord};
        search(0, pts_.size(), 0, world, q, best, found);
        return found;
    }

private:
    void build(std::size_t lo, std::size_t hi, int axis) {
        if (hi - lo <= 1) return;
        std::size_t mid = lo + (hi - lo) / 2;
        std::nth_element(pts_.begin() + lo, pts_.begin() + mid, pts_.begin() + hi,
                         [axis](Point a, Point b) { return lessOn(axis, a, b); });
        build(lo, mid, 1 - axis);
        build(mid + 1, hi, 1 - axis);
    }
    void search(std::size_t lo, std::size_t hi, int axis, Box cell, Point q,
                long long& best, Point& found) const {
        if (lo >= hi || boxDist2(q, cell) >= best) return;  // nothing here can win
        std::size_t mid = lo + (hi - lo) / 2;
        Point p = pts_[mid];
        if (dist2(p, q) < best) {
            best = dist2(p, q);
            found = p;
        }
        Box left = cell, right = cell;
        if (axis == 0) left.xhi = right.xlo = p.x;
        else left.yhi = right.ylo = p.y;
        // Search the side that contains q first: it usually shrinks best quickly.
        if (lessOn(axis, q, p)) {
            search(lo, mid, 1 - axis, left, q, best, found);
            search(mid + 1, hi, 1 - axis, right, q, best, found);
        } else {
            search(mid + 1, hi, 1 - axis, right, q, best, found);
            search(lo, mid, 1 - axis, left, q, best, found);
        }
    }
    std::vector<Point> pts_;
};

// ---- Part 4: Delaunay triangulation by edge flips (half-edge mesh) ----------
using Triangle = std::array<int, 3>;  // vertex indices, counter-clockwise

// Monotone-chain sweep: every point popped off a hull chain becomes a triangle.
std::vector<Triangle> sweepTriangulate(const std::vector<Point>& p) {
    std::vector<int> order(p.size());
    std::iota(order.begin(), order.end(), 0);
    std::ranges::sort(order, [&](int i, int j) { return lexLess(p[i], p[j]); });
    std::vector<Triangle> tris;
    std::vector<int> lo, up;
    for (int q : order) {
        while (lo.size() > 1 && orientation(p[lo.end()[-2]], p[lo.back()], p[q]) < 0) {
            tris.push_back({lo.end()[-2], q, lo.back()});
            lo.pop_back();
        }
        while (up.size() > 1 && orientation(p[up.end()[-2]], p[up.back()], p[q]) > 0) {
            tris.push_back({up.end()[-2], up.back(), q});
            up.pop_back();
        }
        lo.push_back(q);
        up.push_back(q);
    }
    return tris;
}

struct HalfEdge { int origin, twin, next, face; };  // face -1: outside

struct Mesh {
    std::vector<HalfEdge> he;
    std::vector<int> out;  // per vertex: one outgoing half-edge
    int dest(int h) const { return he[he[h].twin].origin; }
    int prev(int h) const { return he[he[h].next].next; }  // triangles only
};

Mesh build(int vertexCount, const std::vector<Triangle>& tris) {
    Mesh m;
    m.out.assign(vertexCount, -1);
    std::map<std::pair<int, int>, int> directed;
    for (int f = 0; f < static_cast<int>(tris.size()); ++f)
        for (int k = 0; k < 3; ++k) {
            int a = tris[f][k], b = tris[f][(k + 1) % 3];
            directed[{a, b}] = 3 * f + k;
            m.he.push_back({a, -1, 3 * f + (k + 1) % 3, f});
            m.out[a] = 3 * f + k;
        }
    std::map<int, int> boundaryOut;
    const int interior = static_cast<int>(m.he.size());
    for (int h = 0; h < interior; ++h) {
        int a = m.he[h].origin, b = m.he[m.he[h].next].origin;
        if (auto it = directed.find({b, a}); it != directed.end()) {
            m.he[h].twin = it->second;
            continue;
        }
        int g = static_cast<int>(m.he.size());
        m.he.push_back({b, h, -1, -1});
        m.he[h].twin = g;
        m.out[b] = boundaryOut[b] = g;
    }
    for (int g = interior; g < static_cast<int>(m.he.size()); ++g)
        m.he[g].next = boundaryOut.at(m.dest(g));
    return m;
}

// Replace edge a-b (triangles a b c and b a d) by edge c-d.
void flip(Mesh& m, int h) {
    int t = m.he[h].twin;
    int h1 = m.he[h].next, h2 = m.he[h1].next;  // b->c, c->a
    int t1 = m.he[t].next, t2 = m.he[t1].next;  // a->d, d->b
    int a = m.he[h].origin, b = m.he[t].origin;
    int c = m.he[h2].origin, d = m.he[t2].origin;
    int f0 = m.he[h].face, f1 = m.he[t].face;
    m.he[h] = {d, t, h2, f0};
    m.he[h2].next = t1;
    m.he[t1].next = h;
    m.he[t1].face = f0;
    m.he[t] = {c, h, t2, f1};
    m.he[t2].next = h1;
    m.he[h1].next = t;
    m.he[h1].face = f1;
    if (m.out[a] == h) m.out[a] = t1;
    if (m.out[b] == t) m.out[b] = h1;
}

bool illegal(const Mesh& m, const std::vector<Point>& p, int h) {
    int t = m.he[h].twin;
    if (m.he[h].face < 0 || m.he[t].face < 0) return false;
    int a = m.he[h].origin, b = m.he[t].origin;
    int c = m.he[m.prev(h)].origin, d = m.he[m.prev(t)].origin;
    return inCircle(p[a], p[b], p[c], p[d]) > 0;
}

// The Delaunay triangulation of distinct points: sweep, then Lawson's flips.
Mesh delaunay(const std::vector<Point>& p) {
    Mesh m = build(static_cast<int>(p.size()), sweepTriangulate(p));
    std::vector<int> stack(m.he.size());
    std::iota(stack.begin(), stack.end(), 0);
    // With an exact inCircle, Lawson needs at most n(n-1)/2 flips (lesson 11).
    // More means a broken predicate, which would otherwise loop forever.
    long long flipsLeft = static_cast<long long>(p.size()) * (p.size() - 1) / 2;
    while (!stack.empty()) {
        int h = stack.back();
        stack.pop_back();
        if (!illegal(m, p, h)) continue;
        if (flipsLeft-- == 0) throw std::logic_error("delaunay: too many flips");
        flip(m, h);
        int t = m.he[h].twin;
        for (int g : {m.he[h].next, m.prev(h), m.he[t].next, m.prev(t)})
            stack.push_back(g);
    }
    return m;
}

// The triangles of a mesh, each once, as vertex index triples.
std::vector<Triangle> triangles(const Mesh& m) {
    std::vector<Triangle> tris;
    for (int h = 0; h < static_cast<int>(m.he.size()); ++h)
        if (m.he[h].face >= 0 && h < m.he[h].next && h < m.prev(h))
            tris.push_back({m.he[h].origin, m.dest(h), m.he[m.prev(h)].origin});
    return tris;
}

// ---- Part 5: property tests against independent references ------------------
struct Lcg {
    unsigned long long s;
    unsigned next() {
        s = s * 6364136223846793005ULL + 1442695040888963407ULL;
        return unsigned(s >> 33);
    }
    int range(int lo, int hi) {  // inclusive
        return lo + int(next() % unsigned(hi - lo + 1));
    }
};

std::vector<Point> randomPoints(Lcg& rng, int n, int range, bool distinct) {
    std::vector<Point> p;
    for (int i = 0; i < n; ++i) p.push_back({rng.range(0, range), rng.range(0, range)});
    if (distinct) {
        std::ranges::sort(p, lexLess);
        auto dups = std::ranges::unique(p);
        p.erase(dups.begin(), dups.end());
    }
    return p;
}

// Small ranges make ties, duplicates and collinear points common.
int pickRange(int i) { return std::array{2, 4, 10, 1000}[i % 4]; }

int sign(Wide v) { return (v > 0) - (v < 0); }

// In-circle from the circumcenter, built exactly as a fraction U / D.
int inCircleByCenter(Point a, Point b, Point c, Point d) {
    Wide bx = b.x - a.x, by = b.y - a.y, cx = c.x - a.x, cy = c.y - a.y;
    Wide D = 2 * (bx * cy - by * cx);  // nonzero: abc is not collinear
    Wide ux = cy * (bx * bx + by * by) - by * (cx * cx + cy * cy);
    Wide uy = bx * (cx * cx + cy * cy) - cx * (bx * bx + by * by);
    Wide dx = D * (d.x - a.x) - ux, dy = D * (d.y - a.y) - uy;
    return sign(ux * ux + uy * uy - dx * dx - dy * dy) * sign(D);  // a is on the circle
}

bool testPredicates(Lcg& rng, int i) {
    auto p = randomPoints(rng, 4, pickRange(i), false);
    Point a = p[0], b = p[1], c = p[2], d = p[3];
    if (orientation(a, b, c) != -orientation(b, a, c)) return false;
    if (orientation(a, b, c) != orientation(b, c, a)) return false;
    if (orientation(a, b, c) == 0) return true;  // no circle through a, b, c
    return inCircle(a, b, c, d) == inCircleByCenter(a, b, c, d);
}

// A convex polygon of input points that turns strictly left everywhere and has
// every input point inside or on it is, by definition, the convex hull.
bool testHull(Lcg& rng, int i) {
    auto p = randomPoints(rng, rng.range(0, 12), pickRange(i), false);
    auto h = convexHull(p);
    const std::size_t k = h.size();
    for (Point v : h)
        if (std::ranges::find(p, v) == p.end()) return false;
    if (k == 0) return p.empty();
    if (k == 1) return std::ranges::all_of(p, [&](Point q) { return q == h[0]; });
    if (k == 2)
        return h[0] != h[1] &&
               std::ranges::all_of(p, [&](Point q) { return onSegment(h[0], h[1], q); });
    for (std::size_t j = 0; j < k; ++j) {
        if (orientation(h[j], h[(j + 1) % k], h[(j + 2) % k]) <= 0) return false;
        for (Point q : p)
            if (orientation(h[j], h[(j + 1) % k], q) < 0) return false;
    }
    return true;
}

// A star-shaped polygon around (0, 0); the reference answer comes from its fan
// of triangles, which tile it exactly.
bool testPointInPolygon(Lcg& rng, int i) {
    const Point dirs[] = {{1, 0}, {2, 1}, {1, 1}, {1, 2}, {0, 1}, {-1, 2},
        {-1, 1}, {-2, 1}, {-1, 0}, {-2, -1}, {-1, -1}, {-1, -2}, {0, -1}, {1, -2},
        {1, -1}, {2, -1}};  // 16 directions in counter-clockwise order
    std::vector<Point> poly;
    for (int j = 0; j < 16; ++j) {
        int r = rng.range(1, 1 + i % 3);
        if (j % 4 == 0 || rng.range(0, 1))  // keep the axes: gaps stay below 180 degrees
            poly.push_back({r * dirs[j].x, r * dirs[j].y});
    }
    const std::size_t n = poly.size();
    for (long long x = -7; x <= 7; ++x)
        for (long long y = -7; y <= 7; ++y) {
            Point q{x, y};
            Where want = Where::Outside;
            for (std::size_t j = 0; j < n; ++j) {
                Point a = poly[j], b = poly[(j + 1) % n];
                if (onSegment(a, b, q)) { want = Where::Boundary; break; }
                if (orientation({0, 0}, a, q) >= 0 && orientation(a, b, q) >= 0 &&
                    orientation(b, {0, 0}, q) >= 0)
                    want = Where::Inside;
            }
            if (crossingTest(poly, q) != want) return false;
        }
    return true;
}

// Segments by Cramer's rule: a + t (b - a) = c + u (d - c) with t, u in [0, 1].
bool intersectByParameters(Point a, Point b, Point c, Point d) {
    long long rx = b.x - a.x, ry = b.y - a.y, sx = d.x - c.x, sy = d.y - c.y;
    long long qx = c.x - a.x, qy = c.y - a.y;
    long long den = rx * sy - ry * sx;
    if (den != 0) {
        long long t = qx * sy - qy * sx, u = qx * ry - qy * rx;
        if (den < 0) den = -den, t = -t, u = -u;
        return 0 <= t && t <= den && 0 <= u && u <= den;
    }
    if (rx == 0 && ry == 0)  // ab is a single point: swap the roles
        return sx == 0 && sy == 0 ? a == c : intersectByParameters(c, d, a, b);
    if (qx * ry - qy * rx != 0) return false;  // parallel, on different lines
    long long t0 = qx * rx + qy * ry, t1 = t0 + sx * rx + sy * ry;  // project onto r
    long long len2 = rx * rx + ry * ry;  // ab covers [0, len2] on that axis
    return std::max(std::min(t0, t1), 0LL) <= std::min(std::max(t0, t1), len2);
}

bool testNearest(Lcg& rng, int i) {
    auto p = randomPoints(rng, rng.range(1, 40), pickRange(i), false);
    KdTree tree(p);
    for (int k = 0; k < 10; ++k) {
        Point q{rng.range(-2, pickRange(i) + 2), rng.range(-2, pickRange(i) + 2)};
        long long brute = std::numeric_limits<long long>::max();
        for (Point v : p) brute = std::min(brute, dist2(v, q));
        if (dist2(tree.nearest(q), q) != brute) return false;  // ties: compare distances
    }
    return true;
}

bool testDelaunay(Lcg& rng, int i) {
    auto p = randomPoints(rng, rng.range(1, 30), pickRange(i), true);
    const int n = static_cast<int>(p.size());
    Mesh m = delaunay(p);
    auto tris = triangles(m);
    long long area = 0;
    for (auto [a, b, c] : tris) {
        if (orientation(p[a], p[b], p[c]) <= 0) return false;
        for (Point q : p)
            if (inCircle(p[a], p[b], p[c], q) > 0) return false;  // empty circles
        area += cross(p[a], p[b], p[c]);
    }
    auto corners = convexHull(p), boundary = convexHull(p, true);
    if (corners.size() < 3) return tris.empty();  // all points collinear
    const int t = static_cast<int>(tris.size()), h = static_cast<int>(boundary.size());
    const int outside = static_cast<int>(m.he.size()) - 3 * t;  // outside half-edges
    if (outside != h || t != 2 * n - 2 - h || area != twiceArea(corners)) return false;
    // A nearest neighbor of every point is one of its Delaunay neighbors.
    for (int v = 0; v < n; ++v) {
        long long ring = std::numeric_limits<long long>::max();
        long long all = ring;
        int h = m.out[v];
        do {
            ring = std::min(ring, dist2(p[v], p[m.dest(h)]));
            h = m.he[m.he[h].twin].next;
        } while (h != m.out[v]);
        for (int w = 0; w < n; ++w)
            if (w != v) all = std::min(all, dist2(p[v], p[w]));
        if (ring != all) return false;
    }
    return true;
}

bool runAllTests() {
    Lcg rng{42};
    bool allPassed = true;
    auto run = [&](const char* name, int total, auto test) {
        int passed = 0;
        for (int i = 0; i < total; ++i) {
            try {
                passed += test(i);
            } catch (const std::exception&) {  // a failed check, not a crash
            }
        }
        std::println("{:<30} {:>6}/{}", name, passed, total);
        allPassed = allPassed && passed == total;
    };
    run("orientation and in-circle", 20000,
        [&](int i) { return testPredicates(rng, i); });
    run("convex hull", 4000, [&](int i) { return testHull(rng, i); });
    run("point in polygon (225 each)", 300,
        [&](int i) { return testPointInPolygon(rng, i); });
    run("all segment pairs, 4x4 grid", 65536, [](int i) {  // every pair, exhaustively
        auto g = [](int k) { return Point{k % 4, k / 4}; };
        Point a = g(i & 15), b = g(i >> 4 & 15), c = g(i >> 8 & 15), d = g(i >> 12);
        return intersects({a, b}, {c, d}) == intersectByParameters(a, b, c, d);
    });
    run("k-d tree nearest (10 each)", 2000, [&](int i) { return testNearest(rng, i); });
    run("Delaunay by flips", 2000, [&](int i) { return testDelaunay(rng, i); });
    return allPassed;
}

// ---- Part 6: the overview drawing and main ----------------------------------
class Svg {
public:
    Svg(double width, double height) : w_(width), h_(height) {}
    void line(Point a, Point b, const char* color = "#1f2937", double width = 1.5) {
        add("<line x1=\"{}\" y1=\"{}\" x2=\"{}\" y2=\"{}\" stroke=\"{}\" stroke-width=\"{}\"/>",
            a.x, fy(a.y), b.x, fy(b.y), color, width);
    }
    void dot(Point p, const char* color = "#1f2937", double r = 3) {
        add("<circle cx=\"{}\" cy=\"{}\" r=\"{}\" fill=\"{}\"/>", p.x, fy(p.y), r, color);
    }
    void polygon(const std::vector<Point>& pts, const char* stroke = "#1f2937",
                 const char* fill = "none") {
        std::string list;
        for (Point p : pts) list += std::format("{},{} ", p.x, fy(p.y));
        add("<polygon points=\"{}\" stroke=\"{}\" fill=\"{}\" stroke-width=\"1.5\"/>", list, stroke, fill);
    }
    void text(Point p, const std::string& s, double size = 8) {
        add("<text x=\"{}\" y=\"{}\" font-size=\"{}\">{}</text>", p.x, fy(p.y), size, s);
    }
    void print() const {
        std::println("<svg viewBox=\"0 0 {} {}\">", w_, h_);
        std::print("{}", body_);
        std::println("</svg>");
    }

private:
    double fy(double y) const { return h_ - y; }  // SVG's y axis points down
    template <class... Args>
    void add(std::format_string<Args...> fmt, Args&&... args) {
        body_ += std::format(fmt, std::forward<Args>(args)...);
        body_ += '\n';
    }
    double w_, h_;
    std::string body_;
};

void drawOverview() {
    Lcg rng{11};
    auto p = randomPoints(rng, 60, 270, true);
    for (Point& v : p) v = {v.x + 15, v.y / 2 + 20};  // fit a 300 x 180 picture
    Svg svg(300, 180);
    auto hull = convexHull(p);
    svg.polygon(hull, "#2563eb", "rgba(37,99,235,0.15)");
    for (auto [a, b, c] : triangles(delaunay(p))) {
        svg.line(p[a], p[b], "#1f2937", 0.6);
        svg.line(p[b], p[c], "#1f2937", 0.6);
        svg.line(p[c], p[a], "#1f2937", 0.6);
    }
    for (Point v : p) svg.dot(v, "#1f2937", 2);
    Point q{150, 100};
    Point nn = KdTree(p).nearest(q);
    svg.line(q, nn, "#dc2626", 1.5);
    svg.dot(nn, "#16a34a", 3.5);
    svg.dot(q, "#dc2626", 3.5);
    const char* where = crossingTest(hull, q) == Where::Inside ? "inside" : "not inside";
    svg.text({15, 8}, std::format("query ({}, {}) is {} the hull; nearest point ({}, {})",
                                  q.x, q.y, where, nn.x, nn.y));
    svg.print();
}

constexpr bool kDrawOverview = false;  // set to true to draw the picture too

int main() {
    bool ok = runAllTests();
    std::println("all tests passed: {}", ok);
    if (ok && kDrawOverview) drawOverview();
}

Output

Press Run to compile and run this example.

orientation and in-circle       20000/20000
convex hull                      4000/4000
point in polygon (225 each)       300/300
all segment pairs, 4x4 grid     65536/65536
k-d tree nearest (10 each)       2000/2000
Delaunay by flips                2000/2000
all tests passed: true

Every count is a pass count, and every test passes: 20,000 predicate cases, 4,000 hulls, 300 polygons with 225 query points each, all 65,536 ordered pairs of segments with endpoints on a 4 by 4 grid (degenerate segments included), 2,000 point sets with 10 nearest-neighbor queries each, and 2,000 Delaunay triangulations. The whole run takes well under a second.

Now set kDrawOverview to true and run it again. The program then also draws 60 seeded points with their convex hull (shaded), their Delaunay triangulation, and one query point in red joined to its nearest neighbor in green, with a caption that reports the point-in-polygon answer for the query against the hull. Every element of the picture comes from a different part of the toolkit; if any two parts disagreed about a convention, this is where you would see it.

When a property test fails, shrink the input

Property tests find bugs; they do not explain them. A failing case with 18 random points is hard to debug, but most of those points are irrelevant. Shrinking removes the noise: repeatedly delete a point or move a coordinate towards zero, and keep each change as long as the test still fails. What remains is usually small enough to understand by hand.

The program below plants a realistic bug in lesson 4's monotone chain: the upper pass keeps collinear points while the lower pass drops them, so the two halves follow different policies. It runs the strict-convexity property on 1,000 seeded inputs, three in four with coordinates up to 1,000 and one in four on a 5 by 5 grid, and then shrinks the first failure.

main.cppc++23

#include <algorithm>
#include <format>
#include <print>
#include <string>
#include <vector>

struct Point { long long x, y; };

struct Lcg {
    unsigned long long s;
    unsigned next() {
        s = s * 6364136223846793005ULL + 1442695040888963407ULL;
        return unsigned(s >> 33);
    }
    int range(int lo, int hi) {  // inclusive
        return lo + int(next() % unsigned(hi - lo + 1));
    }
};

long long cross(Point a, Point b, Point c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
int orientation(Point a, Point b, Point c) {
    long long v = cross(a, b, c);
    return (v > 0) - (v < 0);
}

bool operator==(Point p, Point q) { return p.x == q.x && p.y == q.y; }
bool lexLess(Point p, Point q) { return p.x < q.x || (p.x == q.x && p.y < q.y); }

// Lesson 4's monotone chain with a planted bug: the upper pass keeps collinear
// points although the hull should drop them, so the two chains disagree.
std::vector<Point> buggyHull(std::vector<Point> pts) {
    std::ranges::sort(pts, lexLess);
    auto dups = std::ranges::unique(pts);
    pts.erase(dups.begin(), dups.end());
    if (pts.size() < 3) return pts;
    auto onLine = [&](Point p) { return orientation(pts.front(), pts.back(), p) == 0; };
    if (std::ranges::all_of(pts, onLine))  // no area: the hull is a segment
        return std::vector{pts.front(), pts.back()};
    std::vector<Point> hull;
    for (int pass = 0; pass < 2; ++pass) {  // lower hull, then upper hull
        const size_t start = hull.size();
        for (Point p : pts) {
            while (hull.size() >= start + 2) {
                int turn = orientation(hull[hull.size() - 2], hull.back(), p);
                if (turn > 0 || (turn == 0 && pass == 1)) break;  // the bug
                hull.pop_back();
            }
            hull.push_back(p);
        }
        hull.pop_back();
        std::ranges::reverse(pts);
    }
    return hull;
}

// The property: every turn of the hull is a strict left turn.
bool strictlyConvex(const std::vector<Point>& p) {
    auto h = buggyHull(p);
    for (std::size_t j = 0, k = h.size(); k >= 3 && j < k; ++j)
        if (orientation(h[j], h[(j + 1) % k], h[(j + 2) % k]) <= 0) return false;
    return true;
}

// Make a failing input smaller while it keeps failing: drop points, then pull
// coordinates towards zero. What is left is a counterexample a human can read.
std::vector<Point> shrink(std::vector<Point> p) {
    for (bool progress = true; progress;) {
        progress = false;
        for (std::size_t i = 0; i < p.size(); ++i) {
            auto smaller = p;
            smaller.erase(smaller.begin() + static_cast<std::ptrdiff_t>(i));
            if (!strictlyConvex(smaller)) { p = smaller; progress = true; --i; }
        }
        for (Point& v : p)
            for (long long* c : {&v.x, &v.y})
                while (*c != 0) {
                    long long old = *c;
                    *c /= 2;
                    if (strictlyConvex(p)) { *c = old; break; }
                    progress = true;
                }
    }
    return p;
}

std::string show(const std::vector<Point>& p) {
    std::string s;
    for (Point v : p) s += std::format("{}({}, {})", s.empty() ? "" : " ", v.x, v.y);
    return s;
}

int main() {
    Lcg rng{1};
    int wideFails = 0, gridFails = 0;
    std::vector<Point> firstFailure;
    for (int i = 0; i < 1000; ++i) {
        std::vector<Point> p;
        for (int n = rng.range(3, 25); n > 0; --n)
            p.push_back({rng.range(0, 1000), rng.range(0, 1000)});
        bool grid = i % 4 == 3;  // every fourth case on a 5 x 5 grid: ties are common
        if (grid)
            for (Point& v : p) v = {v.x % 5, v.y % 5};
        if (strictlyConvex(p)) continue;
        ++(grid ? gridFails : wideFails);
        if (firstFailure.empty()) firstFailure = p;
    }
    std::println("failures: {} of 750 wide cases, {} of 250 grid cases", wideFails,
                 gridFails);
    std::println("first failure, {} points:", firstFailure.size());
    std::println("{}", show(firstFailure));
    std::println("shrunk: {}", show(shrink(firstFailure)));
}

Output

Press Run to compile and run this example.

failures: 0 of 750 wide cases, 146 of 250 grid cases
first failure, 18 points:
(0, 0) (4, 1) (1, 4) (0, 4) (1, 1) (4, 1) (1, 1) (2, 3) (3, 2) (4, 2) (4, 4) (3, 4) (2, 4) (3, 3) (4, 2) (3, 3) (3, 2) (1, 1)
shrunk: (2, 4) (0, 4) (1, 4) (0, 0)

Two results matter. First, none of the 750 wide-range inputs exposes the bug, and more than half of the grid inputs do: random coordinates over a large range almost never put three hull points on one line, so a test suite without small-range generators would pass forever. Second, the shrunk counterexample has four points: three collinear on the top edge and one below them to give the hull an area (an all-collinear input takes a separate branch and would not fail). That points straight at the collinear policy. Libraries such as RapidCheck for C++ and Hypothesis for Python build shrinking in; a hand-written loop like this one is enough to start.

CGAL in practice: kernels and exact predicates

CGAL, the Computational Geometry Algorithms Library, is the reference implementation of most of this course and a great deal more. Its central design idea is the one from lesson 2: separate predicates (decisions, which must be exact) from constructions (new coordinates, which may be rounded). A CGAL kernel bundles the geometric types (Point_2, Segment_2, Triangle_2 and so on), the predicates on them and the constructions, for one choice of number type. Algorithms are templates over the kernel, so the same convex hull code runs with any of them.

The CGAL 6.2.1 manual (checked in October 2026) describes the predefined kernels this way:

CGAL::Simple_cartesian<double>: plain double coordinates. Fast and easy to debug, with exactly the failures lesson 2 measured.
CGAL::Exact_predicates_inexact_constructions_kernel: "exact geometric predicates, but geometric constructions may be inexact due to round-off errors". Coordinates are double; predicates are filtered (CGAL's filtered kernels use interval arithmetic) and fall back to exact arithmetic when the filter cannot decide.
CGAL::Exact_predicates_exact_constructions_kernel: exact constructions as well, so constructed points such as segment intersections are exact. Slower, and needed when constructed points are fed back into predicates.

Exact predicates with inexact constructions is the usual default because most algorithms only decide: a convex hull or a Delaunay triangulation outputs input points, so it never constructs anything, and it gets exact answers at close to double speed. You pay for exact constructions only when your algorithm builds new points and then reasons about them, as boolean operations on polygons do.

The following function uses the default kernel for a hull, a Delaunay triangulation and a nearest-site query. It needs CGAL, which the playground does not have, so it is not runnable here. Every name in it was checked against the CGAL 6.2.1 reference manual.

// Requires CGAL 6 (header-only) and its dependencies (Boost headers, GMP, MPFR).
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/convex_hull_2.h>
#include <iterator>
#include <print>
#include <vector>

using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Point_2 = K::Point_2;
using Delaunay = CGAL::Delaunay_triangulation_2<K>;

void summarize(const std::vector<Point_2>& points, const Point_2& query) {
    // A predicate: exact, even for nearly collinear doubles.
    if (points.size() >= 3 &&
        CGAL::orientation(points[0], points[1], points[2]) == CGAL::COLLINEAR)
        std::println("the first three points are collinear");

    std::vector<Point_2> hull;  // the extreme points, counter-clockwise
    CGAL::convex_hull_2(points.begin(), points.end(), std::back_inserter(hull));

    Delaunay dt(points.begin(), points.end());
    double area = 0;  // a construction: rounded, and fine for reporting
    for (Delaunay::Face_handle f : dt.finite_face_handles())
        area += dt.triangle(f).area();
    std::println("{} hull vertices; {} vertices and {} triangles, total area {}",
                 hull.size(), dt.number_of_vertices(), dt.number_of_faces(), area);

    Delaunay::Vertex_handle site = dt.nearest_vertex(query);
    std::println("nearest site: ({}, {})", dt.point(site).x(), dt.point(site).y());
}

Note what CGAL decides for you that this course made explicit. number_of_faces counts finite triangles only, because CGAL closes the triangulation with an "infinite vertex" instead of outside half-edges. convex_hull_2 returns extreme points. Inserting a point that coincides with an existing vertex returns that vertex and leaves the triangulation unchanged. Read the manual page of each function for those policies, exactly as you would read the comments in this toolkit.

Meshes and the Half-Edge Data Structure