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

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