Problem 252: Convex Holes

Project Euler Problem 252 Solution

EulerSolve provides an optimized solution for Project Euler Problem 252, Convex Holes, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The problem generates a deterministic set of 500 integer points and asks for the largest-area convex polygon whose interior contains no other generated point. In computational geometry this is a largest convex hole . The point set is not arbitrary: it comes from the recurrence $$s_0=290797,\qquad s_n=s_{n-1}^2\bmod 50515093,$$ $$t_n=(s_n\bmod 2000)-1000,$$ and the points are $$P_i=(t_{2i-1},t_{2i})\qquad(1\le i\le 500).$$ The code does not enumerate polygons directly. Instead, it reduces the problem to fast emptiness tests for triangles and a dynamic program over angularly ordered vertices. Mathematical Approach Point Generation The recurrence above produces the first three points $$P_1=(527,144),\qquad P_2=(-488,732),\qquad P_3=(-454,-947).$$ The oriented doubled area of this triangle is $$\operatorname{cross}(P_1,P_2,P_3)=1684193,$$ so its ordinary area is $842096.5$. This already explains why the implementation stores doubled areas as integers: all coordinates are integral, so every polygon area is an integer multiple of $1/2$. Cross Product and Doubled Area For three points $A,B,C$, define $$\operatorname{cross}(A,B,C)=(B_x-A_x)(C_y-A_y)-(B_y-A_y)(C_x-A_x).$$ This quantity is positive for a left turn, negative for a right turn, and zero for collinearity....

Detailed mathematical approach

Problem Summary

A pseudo-random generator produces points with integer coordinates: $S_0=290797$, $S_{n+1}=S_n^2 \bmod 50515093$ and $T_n=(S_n \bmod 2000)-1000$, and the $k$-th point is $P_k=(T_{2k-1},T_{2k})$; the first points are $(527,144)$, $(-488,732)$ and $(-454,-947)$. A convex hole of a set of points is a convex polygon whose vertices are points of the set and whose interior contains no point of the set; other points of the set may lie on its perimeter. For the first $20$ points the largest convex hole is a heptagon of area $1049694.5$. We must find the largest area of a convex hole of the first $500$ points, written with one digit after the decimal point. ^[1]

The whole solution in one sentence: cut every convex hole from its lowest vertex into a fan of triangles, which must have no point inside and no point on the diagonals between them, and for each possible lowest vertex build the best fan triangle by triangle over the points sorted by angle, keeping for every possible last edge only the largest area.

Mathematical Approach

We make areas exact, fix which polygons count, reduce the hole condition to tests on triangles and diagonals, and assemble the triangles by dynamic programming.

1. The points

Since $S_n\bmod 2000$ lies between $0$ and $1999$, every coordinate is an integer from $-1000$ to $999$. The $500$ points are all different but not in general position: $66$ triples of them lie on a common line ($62$ lines through exactly three of the points and one line through four), so $65$ segments between two of the points contain a third point. The first $20$ points have no three on a line.

2. Turns and doubled areas

For three points $a,b,c$ put

$$\begin{aligned}\operatorname{cross}(a,b,c)={}&(b_x-a_x)(c_y-a_y)\\&-(b_y-a_y)(c_x-a_x).\end{aligned}$$

This number is positive when $c$ lies to the left of the directed line from $a$ to $b$, so that the path $a\to b\to c$ turns counterclockwise, negative when $c$ lies to the right, and zero when the points are collinear; ^[2] its absolute value is twice the area of the triangle $abc$. ^[3] It is unchanged by a cyclic rotation of $a,b,c$ and changes sign when two of them are swapped. For example $\operatorname{cross}(P_1,P_2,P_3)=1684193$: the first three points form a counterclockwise triangle of area $842096.5$.

By the shoelace formula, twice the area of a polygon with vertices $p_1,\dots,p_r$ in counterclockwise order is $\sum_i(x_iy_{i+1}-x_{i+1}y_i)$ with $p_{r+1}=p_1$. ^[3] With integer coordinates this is an integer, so every area ends in $.0$ or $.5$, which is why one digit after the decimal point suffices. The programs therefore work only with doubled areas, exact integers of at most $2\cdot1999^2\lt8\cdot10^6$, since all points lie in a square of side $1999$.

3. What a hole may contain

Only the interior must be free of points. Take $A=(0,0)$, $B=(6,0)$, $C=(6,4)$, $D=(0,4)$ and $E=(3,2)$. The rectangle $ABCD$ is not a hole, because its centre $E$ lies in its interior, but the triangle $ABC$ is a hole of area $12$, although $E$ lies on its side $CA$, and no hole of these points is larger.

A vertex with an angle of $180^\circ$, in the middle of a straight piece of the boundary, can be removed without changing the region; it then becomes a point on the perimeter, which is allowed. So we may assume that every hole is strictly convex: listed counterclockwise, any three consecutive vertices $p,q,s$ satisfy $\operatorname{cross}(p,q,s)\gt0$. Then no three vertices are collinear, and every vertex lies strictly to the left of each directed edge that does not contain it.

4. The anchor and the angular order

Call a point $q$ above a point $A$ if $q_y\gt A_y$, or if $q_y=A_y$ and $q_x\gt A_x$. Every polygon has exactly one vertex with all other vertices above it: its lowest vertex, the leftmost one on ties. We call it the anchor of the polygon; the Graham scan for convex hulls starts from the same point. ^[2] Seen from $A$, a point above $A$ has a direction angle $\theta$ in $[0^\circ,180^\circ)$, and for two such points $p$ and $q$

$$\operatorname{cross}(A,p,q)=|Ap|\,|Aq|\,\sin(\theta_q-\theta_p).$$

As $\theta_q-\theta_p$ lies strictly between $-180^\circ$ and $180^\circ$, the sign of $\operatorname{cross}(A,p,q)$ is the sign of $\theta_q-\theta_p$, so the points above $A$ can be sorted by angle in exact integer arithmetic; points on a common ray from $A$ are ordered by distance.

Let $A=w_0,w_1,\dots,w_{r-1}$ be the vertices of a strictly convex polygon with anchor $A$, in counterclockwise order. By Section 3, $A$ lies strictly to the left of each edge $w_iw_{i+1}$ with $1\le i\le r-2$, so $\operatorname{cross}(A,w_i,w_{i+1})\gt0$. The angles of $w_1,\dots,w_{r-1}$ therefore increase strictly, and in the sorted list of the points above $A$ the vertices appear in their order along the boundary.

5. Cutting a hole into a fan

The diagonals from $A$ to $w_2,\dots,w_{r-2}$ cut the polygon into the $r-2$ triangles $T_i=Aw_iw_{i+1}$, a fan triangulation. ^[4] A point of the polygon other than $A$ lies on a ray from $A$ between the rays to $w_1$ and $w_{r-1}$, and so it lies

strictly inside some $T_i$: in the interior;
at a vertex, on an edge $w_iw_{i+1}$, or on $Aw_1$ or $Aw_{r-1}$: on the perimeter;
strictly between $A$ and $w_i$ for some $2\le i\le r-2$: in the interior, since the line through two non-adjacent vertices of a strictly convex polygon meets its boundary only at these two vertices.

Hole test. A strictly convex polygon with anchor $A$ is a hole if and only if (a) no point of the set lies strictly inside any triangle $T_i$, and (b) no point of the set lies strictly between $A$ and $w_i$ for any $2\le i\le r-2$.

Condition (b) cannot be dropped. The fan of the rectangle of Section 3 consists of $ABC$ and $ACD$; neither has a point strictly inside, since $E$ lies on their common side, yet $E$ lies on the diagonal $AC$ and therefore inside the rectangle.

6. Testing one triangle

The interior of a counterclockwise triangle $abc$ is the intersection of the three open half-planes to the left of $ab$, $bc$ and $ca$. Let

$$L(a,b)=\{\,p:\ \operatorname{cross}(a,b,p)\gt0\,\}$$

be the set of given points strictly to the left of the directed line from $a$ to $b$. The triangle has no point strictly inside if and only if

$$L(a,b)\cap L(b,c)\cap L(c,a)=\varnothing.$$

A point on a side has cross product $0$ with that side and is not counted, exactly as the statement requires: in the rectangle, $\operatorname{cross}(C,A,E)=0$, so the triangle $ABC$ passes the test.

7. Sets as rows of bits

A set of points is stored as a row of $n$ bits, bit $p$ being $1$ when point $p$ belongs to the set, and intersections become bitwise ANDs, one machine word of $64$ bits at a time. ^[5] The program stores $L(a,b)$ for every ordered pair as $\lceil n/64\rceil=8$ words. Filling all rows tests each of the $n(n-1)(n-2)=124\,251\,000$ ordered triples once; afterwards a triangle test looks at no more than $8$ word positions, with two ANDs per position, and stops at the first position where the three rows share a bit.

8. Building the fan edge by edge

Fix an anchor $A$ and let $v_0,v_1,\dots,v_{m-1}$ be the points above $A$ in angular order. For $j\lt k$ let $F(j,k)$ be the largest doubled area of a strictly convex hole with anchor $A$ whose last two vertices, in counterclockwise order, are $v_j$ and $v_k$, if such a hole exists. Write $t=\operatorname{cross}(A,v_j,v_k)$ for twice the area of the last fan triangle. Such a hole exists only if $t\gt0$ and the triangle $Av_jv_k$ has no point strictly inside, and then

$$F(j,k)=t+\max\bigl(0,\ \max_h F(h,j)\bigr).$$

Here $h$ runs over the indices for which $F(h,j)$ exists and $\operatorname{cross}(v_h,v_j,v_k)\gt0$, and these predecessors are used only if no point lies strictly between $A$ and $v_j$; otherwise $F(j,k)=t$. The term $0$ stands for the triangle $Av_jv_k$ on its own, and a maximum over no predecessors counts as $0$. Since $h\lt j$, all values $F(h,j)$ are known when the pairs are processed with $j$ increasing. The answer is half the largest value of $F$ over all anchors and pairs. This is dynamic programming: a best hole is a best smaller hole plus one triangle, and the same smaller holes serve many larger ones. ^[6]

9. Why the recurrence is right

The triangle alone has doubled area $t$. Now take a strictly convex hole $A,w_1,\dots,w_{r-1}$ with $r\ge4$, $w_{r-2}=v_j$ and $w_{r-1}=v_k$, and remove its last vertex. The rest is again strictly convex: the other turns are unchanged, and by cyclic rotation the new turns at $v_j$ and at $A$ are $\operatorname{cross}(A,w_{r-3},v_j)$ and $\operatorname{cross}(A,w_1,v_j)$, positive because the angles increase. Its fan triangles and diagonals are among those of the larger hole, so it is a hole ending with $v_h=w_{r-3}$ and $v_j$. As the turn at $v_j$ is a left turn and $Av_j$ is a diagonal of the larger hole, $h$ is admissible, and the doubled area is that of the smaller hole plus $t$. Hence $F(j,k)$ is at most the right-hand side.

Conversely, take an admissible $h$ and a hole of doubled area $F(h,j)$ ending with $v_h,v_j$, with first vertex $w_1$ after $A$, and append $v_k$. The new turns are $\operatorname{cross}(v_h,v_j,v_k)\gt0$ at $v_j$, $\operatorname{cross}(v_j,v_k,A)=t\gt0$ at $v_k$ and $\operatorname{cross}(A,w_1,v_k)\gt0$ at $A$, because the angle grows from $w_1$ to $v_k$ by less than $180^\circ$. The edges between consecutive vertices other than $A$ lie in non-overlapping angular sectors, and the two edges at $A$ lie on the outermost rays, so non-adjacent edges never meet: the polygon is simple, and a simple polygon with only strict left turns is strictly convex. Its fan triangles and diagonals, the old ones together with $Av_jv_k$ and $Av_j$, all pass the hole test. Hence equality holds. Every strictly convex hole is found from its own anchor, so the maximum over all anchors is the answer.

10. A worked example

Take $A=(0,0)$, $B=(4,0)$, $C=(4,3)$, $D=(2,3)$ and $E=(2,4)$. Their convex hull $ABCE$ is not a hole, since $D$ lies inside it. The other points are above $A$ in the angular order $B,C,D,E$, because $\operatorname{cross}(A,B,C)=12$, $\operatorname{cross}(A,C,D)=6$ and $\operatorname{cross}(A,D,E)=2$ are positive. The six pairs, in the order of the program:

$(B,C)$: $t=12$, empty, $F=12$; $(B,D)$: $t=12$, empty, $F=12$.
$(B,E)$: $t=16$, but $D$ lies strictly inside $ABE$, so the pair is rejected.
$(C,D)$: $t=6$, empty; the predecessor $(B,C)$ makes the left turn $\operatorname{cross}(B,C,D)=6$ and no point lies strictly between $A$ and $C$, so $F=12+6=18$.
$(C,E)$: $t=10$, but $D$ lies strictly inside $ACE$: rejected.
$(D,E)$: $t=2$, empty; both predecessors give right turns, $\operatorname{cross}(B,D,E)=-2$ and $\operatorname{cross}(C,D,E)=-2$, so $F=2$.

The best value $18$ belongs to the trapezoid $ABCD$, whose horizontal sides of lengths $4$ and $2$ lie $3$ apart, so its area is $(4+2)\cdot3/2=9$. The other anchors give values of at most $8$ (area $4$, from $B$ for the quadrilateral $BCED$), so the largest hole has area $9$, as a check of all subsets confirms.

11. The example of the statement

For the first $20$ points the largest value of $F$ belongs to the heptagon of the statement, the only hole of that area. Its anchor is $P_3=(-454,-947)$, the lowest of the twenty points, and its other vertices in counterclockwise order are $P_9=(913,-860)$, $P_{20}=(273,-175)$, $P_{12}=(-60,138)$, $P_{19}=(-477,-23)$, $P_{14}=(-646,-422)$ and $P_{15}=(-665,-726)$. The doubled areas of its five fan triangles add up to

$$\begin{aligned}&992075+484627+389011\\&{}+165333+68343=2099389,\end{aligned}$$

which gives the area $1049694.5$ of the statement. These points have $1393$ holes, from $456$ triangles to $2$ octagons, so the largest area does not come with the most vertices.

How the Code Works

The C++ program ^[9] runs run_validation_checkpoints (next section), then generate_points(500) and solve_max_area2, and prints the result with print_area_with_single_decimal: the integer part of half the doubled area, followed by $.5$ if the doubled area is odd and by $.0$ otherwise. generate_points runs the recurrence in 64-bit integers, since the square of a value below $50515093$ stays below $2^{52}$.

solve_max_area2 first calls precompute_geometry. For every ordered pair $(a,b)$ it fills a row of left_bits whose bit $p$ is set when $\operatorname{cross}(a,b,p)\gt0$, the set $L(a,b)$ of Section 6, and the flag has_inner_collinear, which is set when another point lies on the segment from $a$ to $b$; point_on_segment_inclusive decides this from a zero cross product and the bounding box of the segment. triangle_is_empty_strict ANDs the rows of $(a,b)$, $(b,c)$ and $(c,a)$ word by word and stops at the first nonzero word. Then solve_for_anchor runs for every point as anchor, and the largest result is returned.

solve_for_anchor implements Sections 4 and 8. It collects the points above the anchor in verts and sorts them with angle_cmp: by the sign of the cross product, then by distance, then by index. The table dp has $m^2$ entries set to $-1$, and incoming[k] lists the indices $j$ for which the pair $(j,k)$ has an entry. For $j\lt k$, a pair is skipped when tri_area2, the value $t$, is not positive or when the triangle is not empty. Otherwise current starts at $t$, and if can_extend_over_j holds, meaning that no point lies between the anchor and $v_j$, every predecessor $h$ in incoming[j] with $\operatorname{cross}(v_h,v_j,v_k)\gt0$ is tried. The value is stored, $j$ is appended to incoming[k], and the best value is updated.

The Java program mirrors these functions with long arithmetic, a long[][] array of rows and a boolean[][] array of flags. It runs no checkpoints, and it prints half the doubled area with String.format and one decimal, which is exact because the doubled area is an integer far below $2^{53}$. The Python file has no geometry of its own: it compiles the C++ source with clang++ or g++ when no up-to-date compiled copy exists, runs it and prints its output, so it shows exactly the result of the C++ program.

Complexity and Verification

Let $n=500$. The cache tests $n(n-1)(n-2)$ ordered triples, $O(n^3)$ time, and its rows take $n^2\lceil n/64\rceil=2\cdot10^6$ words, $16$ MB in C++ (more in Java, where each row is a separate array), plus $n^2$ flags. In the anchor searches every triple of points is examined exactly once, at its lowest point, so the pair loops make $\binom n3=20\,708\,500$ visits. Of these, $66$ are collinear, and the other $20\,708\,434$ tests look at $29\,213\,511$ word positions, about $1.4$ per test, since a large triangle usually contains a point found in the first word. Only $487\,135$ triangles are empty; they are exactly the entries of the tables. Sorting adds $O(n^2\log n)$ time, and the table of one anchor has at most $499^2$ entries of 8 bytes, about 2 MB, plus one index in incoming per entry.

For $n$ points chosen independently and uniformly at random in a convex region, the expected number of empty triangles is at most $4\binom n2=2n(n-1)$, and this is tight up to a lower-order term. ^[8] Our points are only pseudo-random, yet their $487\,135$ empty triangles are close to the bound $499\,000$ for $n=500$. The predecessor loop makes $3\,638\,725$ steps, and $817\,758$ of them pass the turn test, exactly the number of holes with four vertices, since each such hole arises from one successful step. Each predecessor step uses four points, the anchor and $v_h,v_j,v_k$, so there are at most $\binom n4$ steps, and points in convex position, where every triangle is empty, reach this bound; together with the triangle tests the worst case is $O(n^4)$. For these $500$ points the C++ and Java programs finish in well under a second on a desktop computer. Dobkin, Edelsbrunner and Overmars list all empty triangles with an optimal algorithm and extend it to empty convex polygons. ^[7]

All three programs print the same value, and the checkpoints pass: the generator gives the three points of the statement, the fast method agrees on the first $10$ points with brute_force_max_area2, which tries every subset in convex position, and the first $20$ points give the doubled area $2099389$. Further checks, made for this article and not part of the published solutions:

A separate brute force over the convex hulls of all subsets agrees with the program on the first $n$ points for every $n$ from $3$ to $20$.
On all $466$ subsets of at least three points of the grid $\{0,1,2\}^2$, and on $30\,000$ random sets of $3$ to $11$ distinct points with coordinates from $0$ to $6$, $21\,876$ of them with three points on a line, solve_max_area2 agrees with that brute force every time.
An enumeration of all $2\,903\,576$ holes of the $500$ points, with its own angular sort, direct point-in-triangle tests and areas from the shoelace formula, gives the same maximum; the hole with the most vertices has $13$.
Without condition (b) of Section 5 the program would report $24$ instead of $12$ for the rectangle of Section 3. For the $500$ points of the problem the condition happens not to change the maximum.

Footnotes and References

The references below supply the background facts used above; the hole test, the recurrence and all counts are derived or computed in this article.

Project Euler 252 — Convex Holes. The official statement defines convex holes (no given point in the interior, further points allowed on the perimeter), gives the generator and the largest area 1049694.5 for the first 20 points, a heptagon, and asks for the largest area for the first 500 points with one decimal digit.
Wikipedia — Graham scan. The algorithm starts from the lowest point, the leftmost one on ties, sorts the others by angle around it, and classifies turns by the sign of (x₂ − x₁)(y₃ − y₁) − (y₂ − y₁)(x₃ − x₁).
Wikipedia — Shoelace formula. Twice the area of a simple polygon is the cyclic sum of xᵢyᵢ₊₁ − xᵢ₊₁yᵢ, positive for counterclockwise order, and half of each term is the oriented area of a triangle with one vertex at the origin.
Wikipedia — Fan triangulation. One vertex is joined to all other vertices; every convex polygon can be triangulated in this way, and n vertices give n − 2 triangles.
Wikipedia — Bit array. A bit array represents a set, a 1-bit marking membership, and the intersection of two such sets is a bitwise AND taken one machine word at a time.
Wikipedia — Dynamic programming. A problem is broken into simpler subproblems; the method applies when optimal solutions combine optimal solutions of subproblems and the same subproblems recur.
D. Dobkin, H. Edelsbrunner and M. Overmars — Searching for empty convex polygons, Algorithmica 5 (1990), 561–571. Finding empty triangles is related to visibility in a star-shaped polygon; this yields an optimal algorithm for all empty triangles, extended to empty convex polygons and a largest empty convex subset.
M. Balko, M. Scheucher and P. Valtr — Holes and islands in random point sets (arXiv:2003.00909). A k-hole of a point set in general position is an empty convex polygon with k vertices from the set. For n uniform random points in a convex body of volume 1 in dimension d, the expected number of empty simplices is at most 2ᵈ⁻¹ · d! · C(n, d), tight in the plane up to a lower-order term.
C++ — Euler252.cpp. The linked immutable C++ revision contains precompute_geometry, solve_for_anchor and the checks of run_validation_checkpoints. It is the implementation record for this article.

Problem 252 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;

constexpr int kTargetPointCount = 500;
constexpr i64 kMod = 50'515'093;

struct Point {
    int x;
    int y;

    bool operator==(const Point& other) const {
        return x == other.x && y == other.y;
    }
};

inline i64 cross(const Point& a, const Point& b, const Point& c) {
    return static_cast<i64>(b.x - a.x) * (c.y - a.y) -
           static_cast<i64>(b.y - a.y) * (c.x - a.x);
}

inline i64 dist2(const Point& a, const Point& b) {
    const i64 dx = static_cast<i64>(a.x) - b.x;
    const i64 dy = static_cast<i64>(a.y) - b.y;
    return dx * dx + dy * dy;
}

std::vector<Point> generate_points(int count) {
    std::vector<int> t_values;
    t_values.reserve(static_cast<std::size_t>(count) * 2);

    i64 s = 290'797;
    for (int i = 0; i < count * 2; ++i) {
        s = (s * s) % kMod;
        t_values.push_back(static_cast<int>(s % 2000) - 1000);
    }

    std::vector<Point> points;
    points.reserve(count);
    for (int i = 0; i < count; ++i) {
        points.push_back(Point{t_values[2 * i], t_values[2 * i + 1]});
    }

    return points;
}

bool point_on_segment_inclusive(const Point& a, const Point& b, const Point& p) {
    if (cross(a, b, p) != 0) {
        return false;
    }
    const int min_x = std::min(a.x, b.x);
    const int max_x = std::max(a.x, b.x);
    const int min_y = std::min(a.y, b.y);
    const int max_y = std::max(a.y, b.y);
    return p.x >= min_x && p.x <= max_x && p.y >= min_y && p.y <= max_y;
}

struct GeometryCache {
    int n = 0;
    int word_count = 0;
    std::vector<u64> left_bits;   // left_bits[(a*n + b)*word_count + w]
    std::vector<std::uint8_t> has_inner_collinear; // has inner point on open segment a-b
};

GeometryCache precompute_geometry(const std::vector<Point>& points, unsigned threads) {
    GeometryCache cache;
    cache.n = static_cast<int>(points.size());
    cache.word_count = (cache.n + 63) / 64;
    cache.left_bits.assign(static_cast<std::size_t>(cache.n) * cache.n * cache.word_count, 0);
    cache.has_inner_collinear.assign(static_cast<std::size_t>(cache.n) * cache.n, 0);

    if (threads == 0) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0) {
            threads = 1;
        }
    }
    threads = std::max<unsigned>(1, std::min<unsigned>(threads, static_cast<unsigned>(cache.n)));

    auto worker = [&](int a_begin, int a_end) {
        for (int a = a_begin; a < a_end; ++a) {
            for (int b = 0; b < cache.n; ++b) {
                if (a == b) {
                    continue;
                }

                u64* bits = &cache.left_bits[(static_cast<std::size_t>(a) * cache.n + b) *
                                             cache.word_count];
                bool has_inner = false;

                for (int p = 0; p < cache.n; ++p) {
                    if (p == a || p == b) {
                        continue;
                    }

                    const i64 cr = cross(points[a], points[b], points[p]);
                    if (cr > 0) {
                        bits[p >> 6] |= (1ULL << (p & 63));
                    } else if (!has_inner && cr == 0 &&
                               point_on_segment_inclusive(points[a], points[b], points[p])) {
                        has_inner = true;
                    }
                }

                cache.has_inner_collinear[static_cast<std::size_t>(a) * cache.n + b] =
                    static_cast<std::uint8_t>(has_inner);
            }
        }
    };

    if (threads == 1 || cache.n < 64) {
        worker(0, cache.n);
        return cache;
    }

    std::vector<std::thread> pool;
    pool.reserve(threads);

    const int chunk = (cache.n + static_cast<int>(threads) - 1) / static_cast<int>(threads);
    for (unsigned tid = 0; tid < threads; ++tid) {
        const int begin = static_cast<int>(tid) * chunk;
        const int end = std::min(cache.n, begin + chunk);
        if (begin >= end) {
            continue;
        }
        pool.emplace_back(worker, begin, end);
    }

    for (std::thread& th : pool) {
        th.join();
    }

    return cache;
}

inline bool triangle_is_empty_strict(int a,
                                     int b,
                                     int c,
                                     const GeometryCache& cache) {
    const u64* ab = &cache.left_bits[(static_cast<std::size_t>(a) * cache.n + b) *
                                     cache.word_count];
    const u64* bc = &cache.left_bits[(static_cast<std::size_t>(b) * cache.n + c) *
                                     cache.word_count];
    const u64* ca = &cache.left_bits[(static_cast<std::size_t>(c) * cache.n + a) *
                                     cache.word_count];

    for (int w = 0; w < cache.word_count; ++w) {
        if ((ab[w] & bc[w] & ca[w]) != 0) {
            return false;
        }
    }
    return true;
}

i64 solve_for_anchor(int anchor,
                     const std::vector<Point>& points,
                     const GeometryCache& cache) {
    const int n = static_cast<int>(points.size());

    std::vector<int> verts;
    verts.reserve(n - 1);
    for (int idx = 0; idx < n; ++idx) {
        if (idx == anchor) {
            continue;
        }
        if (points[idx].y > points[anchor].y ||
            (points[idx].y == points[anchor].y && points[idx].x > points[anchor].x)) {
            verts.push_back(idx);
        }
    }

    const int m = static_cast<int>(verts.size());
    if (m < 2) {
        return 0;
    }

    auto angle_cmp = [&](int lhs, int rhs) {
        const i64 cr = cross(points[anchor], points[lhs], points[rhs]);
        if (cr != 0) {
            return cr > 0;
        }
        const i64 d_lhs = dist2(points[anchor], points[lhs]);
        const i64 d_rhs = dist2(points[anchor], points[rhs]);
        if (d_lhs != d_rhs) {
            return d_lhs < d_rhs;
        }
        return lhs < rhs;
    };
    std::sort(verts.begin(), verts.end(), angle_cmp);

    std::vector<i64> dp(static_cast<std::size_t>(m) * m, -1);
    std::vector<std::vector<int>> incoming(m);

    i64 best = 0;

    for (int j = 0; j < m - 1; ++j) {
        const int vj = verts[j];
        const bool can_extend_over_j =
            !cache.has_inner_collinear[static_cast<std::size_t>(anchor) * n + vj];

        for (int k = j + 1; k < m; ++k) {
            const int vk = verts[k];
            const i64 tri_area2 = cross(points[anchor], points[vj], points[vk]);
            if (tri_area2 <= 0) {
                continue;
            }

            if (!triangle_is_empty_strict(anchor, vj, vk, cache)) {
                continue;
            }

            i64 current = tri_area2; // base triangle (anchor, vj, vk)

            if (can_extend_over_j) {
                for (const int h : incoming[j]) {
                    const int vh = verts[h];
                    if (cross(points[vh], points[vj], points[vk]) <= 0) {
                        continue;
                    }
                    const i64 candidate = dp[static_cast<std::size_t>(h) * m + j] + tri_area2;
                    if (candidate > current) {
                        current = candidate;
                    }
                }
            }

            dp[static_cast<std::size_t>(j) * m + k] = current;
            incoming[k].push_back(j);
            if (current > best) {
                best = current;
            }
        }
    }

    return best;
}

i64 solve_max_area2(const std::vector<Point>& points,
                    bool allow_multithreading,
                    unsigned requested_threads = 0) {
    if (points.size() < 3) {
        return 0;
    }

    unsigned threads = requested_threads;
    if (threads == 0) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0) {
            threads = 1;
        }
    }

    if (!allow_multithreading || points.size() < 120) {
        threads = 1;
    }

    const GeometryCache cache = precompute_geometry(points, threads);

    if (threads == 1) {
        i64 best = 0;
        for (int anchor = 0; anchor < static_cast<int>(points.size()); ++anchor) {
            const i64 local = solve_for_anchor(anchor, points, cache);
            if (local > best) {
                best = local;
            }
        }
        return best;
    }

    std::atomic<int> next_anchor{0};
    std::vector<i64> partial_best(threads, 0);
    std::vector<std::thread> workers;
    workers.reserve(threads);

    for (unsigned tid = 0; tid < threads; ++tid) {
        workers.emplace_back([&, tid]() {
            i64 local_best = 0;
            while (true) {
                const int anchor = next_anchor.fetch_add(1, std::memory_order_relaxed);
                if (anchor >= static_cast<int>(points.size())) {
                    break;
                }
                const i64 value = solve_for_anchor(anchor, points, cache);
                if (value > local_best) {
                    local_best = value;
                }
            }
            partial_best[tid] = local_best;
        });
    }

    for (std::thread& worker : workers) {
        worker.join();
    }

    i64 best = 0;
    for (const i64 value : partial_best) {
        if (value > best) {
            best = value;
        }
    }
    return best;
}

std::vector<int> convex_hull_indices(std::vector<int> ids, const std::vector<Point>& points) {
    std::sort(ids.begin(), ids.end(), [&](int lhs, int rhs) {
        if (points[lhs].x != points[rhs].x) {
            return points[lhs].x < points[rhs].x;
        }
        if (points[lhs].y != points[rhs].y) {
            return points[lhs].y < points[rhs].y;
        }
        return lhs < rhs;
    });

    std::vector<int> lower;
    for (const int idx : ids) {
        while (lower.size() >= 2 &&
               cross(points[lower[lower.size() - 2]], points[lower.back()], points[idx]) <= 0) {
            lower.pop_back();
        }
        lower.push_back(idx);
    }

    std::vector<int> upper;
    for (int i = static_cast<int>(ids.size()) - 1; i >= 0; --i) {
        const int idx = ids[i];
        while (upper.size() >= 2 &&
               cross(points[upper[upper.size() - 2]], points[upper.back()], points[idx]) <= 0) {
            upper.pop_back();
        }
        upper.push_back(idx);
    }

    if (lower.size() <= 1 || upper.size() <= 1) {
        return {};
    }

    lower.pop_back();
    upper.pop_back();
    lower.insert(lower.end(), upper.begin(), upper.end());

    if (lower.size() < 3) {
        return {};
    }
    return lower;
}

i64 polygon_area2(const std::vector<int>& hull, const std::vector<Point>& points) {
    i64 twice_area = 0;
    for (std::size_t i = 0; i < hull.size(); ++i) {
        const Point& a = points[hull[i]];
        const Point& b = points[hull[(i + 1) % hull.size()]];
        twice_area += static_cast<i64>(a.x) * b.y - static_cast<i64>(a.y) * b.x;
    }
    if (twice_area < 0) {
        twice_area = -twice_area;
    }
    return twice_area;
}

bool point_strictly_inside_convex_polygon(const Point& p,
                                          const std::vector<int>& hull,
                                          const std::vector<Point>& points) {
    bool on_boundary = false;

    for (std::size_t i = 0; i < hull.size(); ++i) {
        const Point& a = points[hull[i]];
        const Point& b = points[hull[(i + 1) % hull.size()]];
        const i64 cr = cross(a, b, p);
        if (cr < 0) {
            return false;
        }
        if (cr == 0 && point_on_segment_inclusive(a, b, p)) {
            on_boundary = true;
        }
    }

    return !on_boundary;
}

i64 brute_force_max_area2(const std::vector<Point>& points) {
    const int n = static_cast<int>(points.size());
    if (n >= 20) {
        return -1; // safety guard: not intended for larger N.
    }

    i64 best = 0;
    const int total_masks = 1 << n;
    for (int mask = 0; mask < total_masks; ++mask) {
        if (__builtin_popcount(static_cast<unsigned>(mask)) < 3) {
            continue;
        }

        std::vector<int> subset;
        subset.reserve(n);
        for (int i = 0; i < n; ++i) {
            if (mask & (1 << i)) {
                subset.push_back(i);
            }
        }

        const std::vector<int> hull = convex_hull_indices(subset, points);
        if (hull.size() != subset.size()) {
            continue;
        }

        const i64 area2 = polygon_area2(hull, points);
        if (area2 <= best) {
            continue;
        }

        bool ok = true;
        for (int p = 0; p < n; ++p) {
            if (mask & (1 << p)) {
                continue;
            }
            if (point_strictly_inside_convex_polygon(points[p], hull, points)) {
                ok = false;
                break;
            }
        }

        if (ok) {
            best = area2;
        }
    }

    return best;
}

bool run_validation_checkpoints() {
    {
        const std::vector<Point> p = generate_points(3);
        const std::vector<Point> expected = {
            {527, 144},
            {-488, 732},
            {-454, -947},
        };
        if (p != expected) {
            std::cerr << "Generator checkpoint failed.\n";
            return false;
        }
    }

    {
        const std::vector<Point> p = generate_points(10);
        const i64 fast = solve_max_area2(p, false, 1);
        const i64 brute = brute_force_max_area2(p);
        if (fast != brute) {
            std::cerr << "Brute-force checkpoint failed for N=10: got " << fast
                      << ", expected " << brute << "\n";
            return false;
        }
    }

    {
        const std::vector<Point> p = generate_points(20);
        constexpr i64 kExpectedArea2 = 2'099'389; // 1049694.5
        const i64 got = solve_max_area2(p, false, 1);
        if (got != kExpectedArea2) {
            std::cerr << "Sample checkpoint failed for N=20: got " << got
                      << ", expected " << kExpectedArea2 << "\n";
            return false;
        }
    }

    return true;
}

void print_area_with_single_decimal(i64 area2) {
    std::cout << (area2 / 2) << ((area2 & 1) ? ".5" : ".0") << '\n';
}

} // namespace

int main() {
    if (!run_validation_checkpoints()) {
        return 1;
    }

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) {
        threads = 4;
    }

    const std::vector<Point> points = generate_points(kTargetPointCount);
    const i64 answer_area2 = solve_max_area2(points, true, threads);
    print_area_with_single_decimal(answer_area2);
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""

    answer_candidates = []
    equal_candidates = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answer_candidates.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equal_candidates.append(m2.group(1).strip())

    if answer_candidates:
        return answer_candidates[-1]
    if equal_candidates:
        return equal_candidates[-1]
    return lines[-1]


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = subprocess.check_output([str(binary)], text=True)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.*;
import java.util.concurrent.*;

public class Euler252 {
    static final int TARGET_POINT_COUNT = 500;
    static final long MOD = 50515093L;

    static class Point {
        int x, y;

        Point(int x, int y) {
            this.x = x;
            this.y = y;
        }
    }

    static long cross(Point a, Point b, Point c) {
        return (long) (b.x - a.x) * (c.y - a.y) - (long) (b.y - a.y) * (c.x - a.x);
    }

    static long dist2(Point a, Point b) {
        long dx = (long) a.x - b.x;
        long dy = (long) a.y - b.y;
        return dx * dx + dy * dy;
    }

    static List<Point> generatePoints(int count) {
        long s = 290797;
        List<Integer> tValues = new ArrayList<>(count * 2);
        for (int i = 0; i < count * 2; ++i) {
            s = (s * s) % MOD;
            tValues.add((int) (s % 2000) - 1000);
        }
        List<Point> points = new ArrayList<>(count);
        for (int i = 0; i < count; ++i) {
            points.add(new Point(tValues.get(2 * i), tValues.get(2 * i + 1)));
        }
        return points;
    }

    static boolean pointOnSegmentInclusive(Point a, Point b, Point p) {
        if (cross(a, b, p) != 0)
            return false;
        int minX = Math.min(a.x, b.x);
        int maxX = Math.max(a.x, b.x);
        int minY = Math.min(a.y, b.y);
        int maxY = Math.max(a.y, b.y);
        return p.x >= minX && p.x <= maxX && p.y >= minY && p.y <= maxY;
    }

    static class GeometryCache {
        int n;
        long[][] leftBits;
        boolean[][] hasInner;

        GeometryCache(int n) {
            this.n = n;
            int words = (n + 63) / 64;
            leftBits = new long[n * n][words];
            hasInner = new boolean[n][n];
        }
    }

    static GeometryCache precomputeGeometry(List<Point> points) {
        int n = points.size();
        GeometryCache cache = new GeometryCache(n);

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<?>> futures = new ArrayList<>();

        for (int a = 0; a < n; ++a) {
            final int fa = a;
            futures.add(executor.submit(() -> {
                for (int b = 0; b < n; ++b) {
                    if (fa == b)
                        continue;

                    long[] bits = cache.leftBits[fa * n + b];
                    boolean inner = false;

                    for (int p = 0; p < n; ++p) {
                        if (p == fa || p == b)
                            continue;
                        long cr = cross(points.get(fa), points.get(b), points.get(p));
                        if (cr > 0) {
                            bits[p >> 6] |= (1L << (p & 63));
                        } else if (!inner && cr == 0
                                && pointOnSegmentInclusive(points.get(fa), points.get(b), points.get(p))) {
                            inner = true;
                        }
                    }
                    cache.hasInner[fa][b] = inner;
                }
            }));
        }

        for (Future<?> f : futures) {
            try {
                f.get();
            } catch (Exception e) {
            }
        }
        executor.shutdown();
        return cache;
    }

    static boolean triangleIsEmptyStrict(int a, int b, int c, GeometryCache cache) {
        long[] ab = cache.leftBits[a * cache.n + b];
        long[] bc = cache.leftBits[b * cache.n + c];
        long[] ca = cache.leftBits[c * cache.n + a];

        for (int w = 0; w < ab.length; ++w) {
            if ((ab[w] & bc[w] & ca[w]) != 0)
                return false;
        }
        return true;
    }

    static long solveForAnchor(int anchor, List<Point> points, GeometryCache cache) {
        int n = points.size();
        List<Integer> verts = new ArrayList<>();

        for (int idx = 0; idx < n; ++idx) {
            if (idx == anchor)
                continue;
            if (points.get(idx).y > points.get(anchor).y ||
                    (points.get(idx).y == points.get(anchor).y && points.get(idx).x > points.get(anchor).x)) {
                verts.add(idx);
            }
        }

        int m = verts.size();
        if (m < 2)
            return 0;

        verts.sort((lhs, rhs) -> {
            long cr = cross(points.get(anchor), points.get(lhs), points.get(rhs));
            if (cr != 0)
                return cr > 0 ? -1 : 1;
            long dLhs = dist2(points.get(anchor), points.get(lhs));
            long dRhs = dist2(points.get(anchor), points.get(rhs));
            if (dLhs != dRhs)
                return dLhs < dRhs ? -1 : 1;
            return Integer.compare(lhs, rhs);
        });

        long[] dp = new long[m * m];
        Arrays.fill(dp, -1);
        List<List<Integer>> incoming = new ArrayList<>(m);
        for (int i = 0; i < m; i++)
            incoming.add(new ArrayList<>());

        long best = 0;

        for (int j = 0; j < m - 1; ++j) {
            int vj = verts.get(j);
            boolean canExtendOverJ = !cache.hasInner[anchor][vj];

            for (int k = j + 1; k < m; ++k) {
                int vk = verts.get(k);
                long triArea2 = cross(points.get(anchor), points.get(vj), points.get(vk));
                if (triArea2 <= 0)
                    continue;

                if (!triangleIsEmptyStrict(anchor, vj, vk, cache))
                    continue;

                long current = triArea2;

                if (canExtendOverJ) {
                    for (int h : incoming.get(j)) {
                        int vh = verts.get(h);
                        if (cross(points.get(vh), points.get(vj), points.get(vk)) <= 0)
                            continue;
                        long candidate = dp[h * m + j] + triArea2;
                        if (candidate > current) {
                            current = candidate;
                        }
                    }
                }

                dp[j * m + k] = current;
                incoming.get(k).add(j);
                if (current > best)
                    best = current;
            }
        }
        return best;
    }

    public static String solve() {
        List<Point> points = generatePoints(TARGET_POINT_COUNT);
        GeometryCache cache = precomputeGeometry(points);

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<Long>> futures = new ArrayList<>();

        for (int anchor = 0; anchor < points.size(); ++anchor) {
            final int fAnchor = anchor;
            futures.add(executor.submit(() -> solveForAnchor(fAnchor, points, cache)));
        }

        long best = 0;
        for (Future<Long> f : futures) {
            try {
                long val = f.get();
                if (val > best)
                    best = val;
            } catch (Exception e) {
            }
        }
        executor.shutdown();

        double ans = best / 2.0;
        return String.format(Locale.US, "%.1f", ans);
    }

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