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

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. Geometrically,

$$|\operatorname{cross}(A,B,C)|=2\,[ABC],$$

where $[ABC]$ denotes the Euclidean area of triangle $ABC$. Using doubled area avoids all floating-point roundoff during the DP.

Left-of-Edge Sets and Empty Triangles

For every directed edge $(a,b)$, the code precomputes the set

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

that is, all generated points lying strictly to the left of the oriented line from $a$ to $b$.

If triangle $(a,b,c)$ is oriented counterclockwise, then a point $p$ lies strictly inside it exactly when

$$p\in L(a,b)\cap L(b,c)\cap L(c,a).$$

Therefore the triangle is strictly empty if and only if

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

This criterion is ideal for bitsets: emptiness reduces to a few word-wise AND operations.

The word “strictly” matters. The problem forbids points in the interior of the polygon, but points on the boundary are allowed. That is why the test uses $\gt 0$ instead of $\ge 0$.

Why Anchors Are Enough

Every convex polygon has a unique lowest vertex; if several vertices have the same $y$-coordinate, the leftmost among them is unique. Call this vertex the anchor $A$.

Then every other vertex $v$ of the polygon must satisfy

$$v_y>A_y\quad\text{or}\quad(v_y=A_y\text{ and }v_x>A_x).$$

So once $A$ is fixed, every admissible polygon using $A$ lies entirely in the upper half-plane relative to the horizontal through $A$. The code uses exactly this fact to build the candidate list for each anchor.

It then sorts those candidates by polar angle around $A$. For a convex polygon containing $A$ as lowest-leftmost vertex, this angular order is precisely the boundary order of the remaining vertices.

Dynamic Programming on Convex Chains

Fix an anchor $A$, and let the angle-sorted candidates be

$$v_0,v_1,\dots,v_{m-1}.$$

The DP state $dp[j,k]$ stores the largest doubled area of an empty convex chain

$$A\to \cdots \to v_j \to v_k$$

whose final edge is $(v_j,v_k)$. The base case is the triangle $A,v_j,v_k$, whose doubled area contribution is

$$\Delta(A,v_j,v_k)=\operatorname{cross}(A,v_j,v_k).$$

If $\Delta\le 0$, the triple is not counterclockwise and cannot be part of the chain.

The transition is

$$dp[j,k]=\max\Big(\Delta(A,v_j,v_k),\max_h\{dp[h,j]+\Delta(A,v_j,v_k)\}\Big),$$

subject to three constraints:

1. The triangle $(A,v_j,v_k)$ must be strictly empty.

2. The turn at $v_j$ must stay convex:

$$\operatorname{cross}(v_h,v_j,v_k)\gt 0.$$

3. The segment $(A,v_j)$ may not contain another generated point in its open interior when it is used as an internal diagonal.

Why the Collinearity Check Matters

The code precomputes a Boolean flag has_inner_collinear[a][b] that records whether some generated point lies on the open segment from $a$ to $b$.

This is not needed for boundary edges of the final polygon, because boundary points are allowed. But when the DP extends a chain through $v_j$, the segment $(A,v_j)$ becomes an internal diagonal of the anchor fan triangulation. If another point lay on that open segment, it would lie in the polygon interior, which is forbidden. Hence such transitions are blocked.

Why the DP Computes the Correct Maximum

Take any empty convex polygon and choose its lowest-leftmost vertex as anchor $A$. Its remaining vertices appear in increasing angular order around $A$, and the polygon can be triangulated into the fan

$$ (A,v_{i_1},v_{i_2}),\ (A,v_{i_2},v_{i_3}),\ \dots .$$

Because the polygon interior is empty, each of these triangles is strictly empty. Because the polygon is convex, consecutive triples make only left turns. Therefore every valid polygon corresponds to a valid DP chain, and the chain area is exactly the sum of its triangle contributions.

Conversely, any DP chain satisfying the emptiness, left-turn, and diagonal constraints reconstructs an empty convex polygon with anchor $A$. So the DP searches exactly the feasible family and maximizes its area.

Validation Checkpoint

The C++ code includes a checkpoint for the first 20 generated points. It verifies that the maximum doubled area is

$$2099389,$$

which corresponds to the area

$$1049694.5.$$

This is useful both as a correctness test and as a reminder that half-integer areas naturally arise.

How the Code Works

The implementation has four main stages. First, it generates the 500 points from the pseudo-random recurrence. Second, it precomputes for every directed pair $(a,b)$ both the left-of-edge bitset $L(a,b)$ and the flag telling whether the open segment $ab$ contains a generated point. Third, for each anchor $A$, it filters the candidate vertices above/right of $A$, sorts them by angle, and runs the chain DP described above. Finally, it takes the best doubled area among all anchors and prints it as either .0 or .5.

Complexity Analysis

Let $n=500$. The bitset table contains one bitset for each ordered pair of vertices, so its size is

$$O\!\left(n^2\cdot \left\lceil \frac{n}{64}\right\rceil\right)=O\!\left(\frac{n^3}{64}\right)$$

machine words, plus an $O(n^2)$ table for the collinearity flags.

The geometric precomputation scans every triple $(a,b,p)$, so its arithmetic cost is $O(n^3)$, but the later emptiness tests become very fast because they are bit-parallel.

For a fixed anchor with $m$ visible candidates, the DP is quadratic in the state space and cubic in the worst case when all transitions are possible:

$$O(m^3).$$

In practice, the emptiness tests, orientation pruning, and incoming-edge filtering discard many transitions. The outer loop over anchors is parallelized across threads.

Footnotes and References

Problem page: https://projecteuler.net/problem=252
Cross products and orientation tests: cp-algorithms — Basic Geometry
Shoelace formula / polygon area: Wikipedia — Shoelace formula
Convex hull and related structures: Wikipedia — Convex hull
Bit array / bitset background: Wikipedia — Bit array

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());
    }
}