Problem 514: Geoboard Shapes

Project Euler Problem 514 Solution

EulerSolve provides an optimized solution for Project Euler Problem 514, Geoboard Shapes, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Consider the lattice square $\{0,\dots,N\}^2$, which contains $(N+1)^2$ grid points. Each point is selected independently with probability $$p=\frac{1}{N+1},\qquad q=1-p.$$ Let $H$ be the convex hull of the selected points. The goal is to compute $\mathbb{E}[\operatorname{Area}(H)]$. A direct enumeration of all $2^{(N+1)^2}$ subsets is impossible, so the solution rewrites the expected area as a sum over possible hull edges and then evaluates that sum direction by direction. Mathematical Approach The central idea is that area can be decomposed into oriented edge contributions. Once that is done, linearity of expectation allows us to analyze one candidate hull edge at a time and sum the results. Step 1: Rewrite area as a sum over oriented hull edges If the hull vertices in counterclockwise order are $v_0,\dots,v_{m-1}$, then the shoelace formula gives $$2\,\operatorname{Area}(H)=\sum_{t=0}^{m-1} v_t\times v_{t+1},\qquad v_t\times v_{t+1}=x_t y_{t+1}-y_t x_{t+1}.$$ Taking expectation and using linearity, we may sum over all ordered lattice-point pairs $(u,v)$: $$2\,\mathbb{E}[\operatorname{Area}(H)]=\sum_{u,v}(u\times v)\,\mathbb{P}(u\to v\text{ appears as a counterclockwise hull edge}).$$ So the problem becomes a probability question: for a fixed oriented segment $u\to v$, when does it contribute as an actual edge of the convex hull?...

Detailed mathematical approach

Problem Summary

Consider the lattice square $\{0,\dots,N\}^2$, which contains $(N+1)^2$ grid points. Each point is selected independently with probability

$$p=\frac{1}{N+1},\qquad q=1-p.$$

Let $H$ be the convex hull of the selected points. The goal is to compute $\mathbb{E}[\operatorname{Area}(H)]$. A direct enumeration of all $2^{(N+1)^2}$ subsets is impossible, so the solution rewrites the expected area as a sum over possible hull edges and then evaluates that sum direction by direction.

Mathematical Approach

The central idea is that area can be decomposed into oriented edge contributions. Once that is done, linearity of expectation allows us to analyze one candidate hull edge at a time and sum the results.

Step 1: Rewrite area as a sum over oriented hull edges

If the hull vertices in counterclockwise order are $v_0,\dots,v_{m-1}$, then the shoelace formula gives

$$2\,\operatorname{Area}(H)=\sum_{t=0}^{m-1} v_t\times v_{t+1},\qquad v_t\times v_{t+1}=x_t y_{t+1}-y_t x_{t+1}.$$

Taking expectation and using linearity, we may sum over all ordered lattice-point pairs $(u,v)$:

$$2\,\mathbb{E}[\operatorname{Area}(H)]=\sum_{u,v}(u\times v)\,\mathbb{P}(u\to v\text{ appears as a counterclockwise hull edge}).$$

So the problem becomes a probability question: for a fixed oriented segment $u\to v$, when does it contribute as an actual edge of the convex hull?

Step 2: Group candidate edges by primitive direction and supporting line

Any lattice segment has direction proportional to a primitive vector $d=(dx,dy)$ with $\gcd(|dx|,|dy|)=1$. Fix such a direction. All lattice points on a line parallel to $d$ satisfy

$$s(x,y)=dy\,x-dx\,y=\text{constant}.$$

Therefore the implementation groups candidate edges by primitive direction and by the value of $s$, which identifies one supporting line. On a fixed line, ordering the lattice points by repeated steps of $d$ gives

$$P_0,P_1,\dots,P_{k-1}.$$

Every possible hull edge with that slope is some pair $P_i\to P_j$ with $0\le i\lt j\lt k$.

Step 3: Compute the probability that one pair is the hull edge

Fix one supporting line and one pair $P_i\to P_j$. For this segment to be the hull edge whose interior lies on the left, four independent conditions must hold:

$$\text{(a) }P_i\text{ and }P_j\text{ are selected},$$

$$\text{(b) every point strictly to the right of the line is unselected},$$

$$\text{(c) every collinear point before }P_i\text{ or after }P_j\text{ is unselected},$$

$$\text{(d) at least one point strictly to the left of the line is selected}.$$

If $L$ is the number of lattice points strictly on the left, $R$ the number strictly on the right, and

$$O=i+(k-1-j)$$

the number of collinear points outside the segment, then independence gives

$$\mathbb{P}(P_i\to P_j\text{ is that hull edge})=p^2 q^{R} q^{O}(1-q^{L}).$$

Points lying between $P_i$ and $P_j$ on the same line are allowed to be selected; they remain on the boundary segment and do not change the area contribution. The factor $1-q^L$ also removes the degenerate case in which all selected points are collinear and the area is zero.

Step 4: Sum all pairs on one supporting line

For a fixed line, $L$ and $R$ depend only on the line, while $O=i+(k-1-j)$ depends on the chosen pair. Hence the total expected oriented contribution from one line is

$$p^2 q^{R}(1-q^{L})\sum_{0\le i\lt j\lt k}(P_i\times P_j)\,q^{i}q^{k-1-j}.$$

This is exactly why the implementation loops over all ordered pairs on a line: the cross product $P_i\times P_j$ supplies the geometric term, and the powers of $q$ encode the probability that these two points are the extreme selected points on that line.

Step 5: Sum over all primitive directions

Now sum the previous expression over every supporting line of every primitive direction. Each genuine convex-hull edge is counted exactly once: its slope determines the primitive direction, and the counterclockwise orientation is the unique orientation for which the hull lies on the left and the exterior lies on the right.

Therefore

$$\boxed{\mathbb{E}[\operatorname{Area}(H)]=\frac12\sum_{\text{primitive }d}\ \sum_{\ell\parallel d} p^2 q^{R(\ell)}(1-q^{L(\ell)})\sum_{0\le i\lt j\lt k_\ell}(P_i\times P_j)\,q^{i}q^{k_\ell-1-j}.}$$

This is the formula evaluated by the program.

Worked Example: $N=1$

For $N=1$, the grid is the unit square with four corner points and

$$p=\frac12.$$

The hull has area $1$ only when all four corners are selected, which happens with probability $1/16$. It has area $1/2$ when exactly three of the four corners are selected, and there are four such configurations, each with probability $1/16$. All remaining configurations have area $0$.

Hence

$$\mathbb{E}[\operatorname{Area}(H)]=1\cdot\frac{1}{16}+\frac12\cdot 4\cdot\frac{1}{16}=\frac{3}{16}=0.18750,$$

which matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. They first compute $p$, $q$, and a table of powers $q^0,q^1,\dots,q^{(N+1)^2}$, because every probability factor is assembled from these values. Next they enumerate every primitive direction $(dx,dy)$ inside the square.

For one direction, the implementation scans all lattice points and evaluates $s(x,y)=dy\,x-dx\,y$. This simultaneously gives the population of each parallel line and identifies the boundary point from which that line should be traversed so that every line is visited exactly once. Prefix sums over the line populations then give, for each line, how many points lie strictly to its left and strictly to its right.

After that, each line is walked in order, its points are stored as $P_0,\dots,P_{k-1}$, and every pair $P_i,P_j$ contributes

$$(P_i\times P_j)\,q^i q^{k-1-j}$$

to the line sum. Multiplying by the line factor

$$p^2 q^R(1-q^L)$$

turns that geometric sum into an expected oriented-area contribution. Adding all directions and dividing by $2$ yields the final expected area. Implementations with native parallel support split the direction set across workers, but the mathematical result is identical.

Complexity Analysis

The number of primitive directions with $|dx|,|dy|\le N$ is $O(N^2)$. For a fixed direction, building line statistics over the $(N+1)^2$ lattice points costs $O(N^2)$, while the pair accumulation on one line of length $k$ costs $O(k^2)$. Summed over all lines of that direction, $\sum k=(N+1)^2$ and the worst case is $\sum k^2=O(N^3)$, so one direction costs $O(N^3)$ time in the worst case.

Consequently the total work is $O(N^5)$. The memory usage is $O(N^2)$: the power table, the per-direction line counts, and the temporary storage for a single line are all quadratic or smaller. Parallel execution reduces wall-clock time but does not change the asymptotic bounds.

Footnotes and References

Problem page: https://projecteuler.net/problem=514
The oriented-area identity for polygons, often called the shoelace formula.
Convex hulls, supporting lines, and the fact that a hull edge is determined by one empty half-plane.
Primitive lattice directions, i.e. coprime integer step vectors, which enumerate all lattice slopes without duplication.

Problem 514 source code

C++

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <thread>
#include <vector>

using namespace std;

struct Direction {
    int dx;
    int dy;
};

struct Point {
    int x;
    int y;
};

static vector<Direction> build_directions(int N) {
    vector<Direction> dirs;
    dirs.reserve(4 * N * N + 4);
    for (int dx = -N; dx <= N; ++dx) {
        for (int dy = -N; dy <= N; ++dy) {
            if (dx == 0 && dy == 0) continue;
            int g = std::gcd(std::abs(dx), std::abs(dy));
            if (g != 1) continue;
            dirs.push_back({dx, dy});
        }
    }
    return dirs;
}

static long double process_direction(int N,
                                     const Direction& dir,
                                     const vector<long double>& pow_q,
                                     long double p,
                                     int total_points) {
    // For a line with direction (dx, dy), each oriented edge contributes via
    // P(edge) = p^2 * q^R * q^O * (1 - q^L), with R/L counts of points on the
    // right/left half-planes and O collinear points outside the segment.
    const int dx = dir.dx;
    const int dy = dir.dy;
    const int n1 = dy;
    const int n2 = -dx;

    auto eval = [&](int x, int y) -> int {
        return n1 * x + n2 * y;
    };

    const int s00 = 0;
    const int sN0 = eval(N, 0);
    const int s0N = eval(0, N);
    const int sNN = eval(N, N);
    const int min_s = min(min(s00, sN0), min(s0N, sNN));
    const int max_s = max(max(s00, sN0), max(s0N, sNN));
    const int range = max_s - min_s + 1;

    vector<int> counts(range, 0);
    vector<Point> starts;
    starts.reserve((N + 1) * (std::abs(dx) + std::abs(dy) + 1));

    for (int x = 0; x <= N; ++x) {
        for (int y = 0; y <= N; ++y) {
            int s = eval(x, y);
            counts[s - min_s]++;
            int px = x - dx;
            int py = y - dy;
            if (px < 0 || px > N || py < 0 || py > N) {
                starts.push_back({x, y});
            }
        }
    }

    vector<int> prefix(range + 1, 0);
    for (int i = 0; i < range; ++i) prefix[i + 1] = prefix[i] + counts[i];

    const long double p2 = p * p;
    long double sum_dir = 0.0L;
    vector<Point> line;

    for (const auto& start : starts) {
        line.clear();
        int x = start.x;
        int y = start.y;
        while (x >= 0 && x <= N && y >= 0 && y <= N) {
            line.push_back({x, y});
            x += dx;
            y += dy;
        }
        int k = static_cast<int>(line.size());
        if (k < 2) continue;

        int s0 = eval(start.x, start.y);
        int idx = s0 - min_s;
        int L = prefix[idx];
        if (L == 0) continue;
        int k_on = counts[idx];
        int R = total_points - prefix[idx] - k_on;

        long double factor = p2 * pow_q[R] * (1.0L - pow_q[L]);
        if (factor == 0.0L) continue;

        long double sum_line = 0.0L;
        for (int i = 0; i < k - 1; ++i) {
            long double wi = pow_q[i];
            int xi = line[i].x;
            int yi = line[i].y;
            for (int j = i + 1; j < k; ++j) {
                long double wj = pow_q[k - 1 - j];
                long double cross = static_cast<long double>(xi) * line[j].y
                                    - static_cast<long double>(line[j].x) * yi;
                sum_line += cross * wi * wj;
            }
        }
        sum_dir += factor * sum_line;
    }
    return sum_dir;
}

static long double compute_expected_area(int N, int threads) {
    const int total_points = (N + 1) * (N + 1);
    const long double p = 1.0L / static_cast<long double>(N + 1);
    const long double q = 1.0L - p;

    vector<long double> pow_q(total_points + 1, 1.0L);
    for (int i = 1; i <= total_points; ++i) pow_q[i] = pow_q[i - 1] * q;

    vector<Direction> dirs = build_directions(N);
    if (threads < 1) threads = 1;
    if (threads > static_cast<int>(dirs.size())) {
        threads = static_cast<int>(dirs.size());
        if (threads < 1) threads = 1;
    }

    vector<long double> partial(threads, 0.0L);
    vector<thread> pool;
    pool.reserve(threads);

    for (int t = 0; t < threads; ++t) {
        pool.emplace_back([&, t]() {
            long double local = 0.0L;
            for (size_t idx = t; idx < dirs.size(); idx += threads) {
                local += process_direction(N, dirs[idx], pow_q, p, total_points);
            }
            partial[t] = local;
        });
    }

    for (auto& th : pool) th.join();

    long double sum = 0.0L;
    for (long double v : partial) sum += v;
    return 0.5L * sum;
}

static long double round5(long double x) {
    return floor(x * 100000.0L + 0.5L) / 100000.0L;
}

static bool run_validation(int threads) {
    struct Check {
        int N;
        long double expected;
    };

    const Check checks[] = {
        {1, 0.18750L},
        {2, 0.94335L},
        {10, 55.03013L},
    };

    bool ok = true;
    for (const auto& c : checks) {
        long double got = compute_expected_area(c.N, threads);
        long double rounded = round5(got);
        if (fabsl(rounded - c.expected) > 1e-9L) {
            cerr << "Validation failed for N=" << c.N
                 << ": got " << fixed << setprecision(5) << rounded
                 << ", expected " << fixed << setprecision(5) << c.expected
                 << " (raw " << setprecision(10) << got << ")\n";
            ok = false;
        }
    }

    if (ok) {
        cerr << "Validation checkpoints passed." << '\n';
    }
    return ok;
}

int main(int argc, char** argv) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int N = 100;
    unsigned hw = thread::hardware_concurrency();
    int threads = hw ? static_cast<int>(hw) : 1;
    threads = max(1, min(threads, 8));
    bool validate = true;

    // Optional CLI: ./a.out [N] [threads] [validate(0/1)]
    if (argc >= 2) N = stoi(argv[1]);
    if (argc >= 3) threads = max(1, stoi(argv[2]));
    if (argc >= 4) validate = (stoi(argv[3]) != 0);

    if (validate && !run_validation(min(threads, 4))) {
        return 1;
    }

    long double answer = compute_expected_area(N, threads);
    cout << fixed << setprecision(5) << answer << "\n";
    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 ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


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 = run_cpp(binary=binary, src=src, root=root)
    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.ArrayList;
import java.util.List;

public class Euler514 {

    static class Direction {
        int dx, dy;

        Direction(int dx, int dy) {
            this.dx = dx;
            this.dy = dy;
        }
    }

    static class Point {
        int x, y;

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

    static int gcd(int a, int b) {
        while (b != 0) {
            int temp = a % b;
            a = b;
            b = temp;
        }
        return Math.abs(a);
    }

    static List<Direction> buildDirections(int N) {
        List<Direction> dirs = new ArrayList<>();
        for (int dx = -N; dx <= N; dx++) {
            for (int dy = -N; dy <= N; dy++) {
                if (dx == 0 && dy == 0)
                    continue;
                if (gcd(Math.abs(dx), Math.abs(dy)) != 1)
                    continue;
                dirs.add(new Direction(dx, dy));
            }
        }
        return dirs;
    }

    static double processDirection(int N, Direction dir, double[] powQ, double p, int totalPoints) {
        int n1 = dir.dy;
        int n2 = -dir.dx;

        int s00 = 0;
        int sN0 = n1 * N;
        int s0N = n2 * N;
        int sNN = n1 * N + n2 * N;

        int minS = Math.min(Math.min(s00, sN0), Math.min(s0N, sNN));
        int maxS = Math.max(Math.max(s00, sN0), Math.max(s0N, sNN));
        int range = maxS - minS + 1;

        int[] counts = new int[range];
        List<Point> starts = new ArrayList<>();

        for (int x = 0; x <= N; x++) {
            for (int y = 0; y <= N; y++) {
                int s = n1 * x + n2 * y;
                counts[s - minS]++;
                int px = x - dir.dx;
                int py = y - dir.dy;
                if (px < 0 || px > N || py < 0 || py > N) {
                    starts.add(new Point(x, y));
                }
            }
        }

        int[] prefix = new int[range + 1];
        for (int i = 0; i < range; i++) {
            prefix[i + 1] = prefix[i] + counts[i];
        }

        double p2 = p * p;
        double sumDir = 0.0;
        List<Point> line = new ArrayList<>();

        for (Point start : starts) {
            line.clear();
            int x = start.x;
            int y = start.y;
            while (x >= 0 && x <= N && y >= 0 && y <= N) {
                line.add(new Point(x, y));
                x += dir.dx;
                y += dir.dy;
            }

            int k = line.size();
            if (k < 2)
                continue;

            int s0 = n1 * start.x + n2 * start.y;
            int idx = s0 - minS;
            int L = prefix[idx];
            if (L == 0)
                continue;
            int kOn = counts[idx];
            int R = totalPoints - prefix[idx] - kOn;

            double factor = p2 * powQ[R] * (1.0 - powQ[L]);
            if (factor == 0.0)
                continue;

            double sumLine = 0.0;
            for (int i = 0; i < k - 1; i++) {
                double wi = powQ[i];
                int xi = line.get(i).x;
                int yi = line.get(i).y;
                for (int j = i + 1; j < k; j++) {
                    double wj = powQ[k - 1 - j];
                    double cross = (double) xi * line.get(j).y - (double) line.get(j).x * yi;
                    sumLine += cross * wi * wj;
                }
            }
            sumDir += factor * sumLine;
        }
        return sumDir;
    }

    static double computeExpectedArea(int N) {
        int totalPoints = (N + 1) * (N + 1);
        double p = 1.0 / (double) (N + 1);
        double q = 1.0 - p;

        double[] powQ = new double[totalPoints + 1];
        powQ[0] = 1.0;
        for (int i = 1; i <= totalPoints; i++)
            powQ[i] = powQ[i - 1] * q;

        List<Direction> dirs = buildDirections(N);

        // Multithreaded evaluation for speed
        double sumTotal = dirs.parallelStream()
                .mapToDouble(d -> processDirection(N, d, powQ, p, totalPoints))
                .sum();

        return 0.5 * sumTotal;
    }

    public static void main(String[] args) {
        double ans = computeExpectedArea(100);
        System.out.printf(java.util.Locale.US, "%.5f\n", ans);
    }
}