Problem 742: Minimum Area of a Convex Grid Polygon

Project Euler Problem 742 Solution

EulerSolve provides an optimized solution for Project Euler Problem 742, Minimum Area of a Convex Grid Polygon, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $A(N)$ be the minimum area of a convex lattice polygon with $N$ edges. The supplied C++, Python, and Java implementations solve the target case $N\equiv 0 \pmod{4}$ by writing $N=4Q$, constructing one monotone quarter-chain, and recovering the full polygon by symmetry. The search is therefore over ordered primitive edge directions and compact area summaries rather than over polygons themselves. Mathematical Approach The implementation organizes the polygon as four reflected copies of a quarter construction. That reduces the problem to a structured optimization on primitive lattice directions. Step 1: Reduce to One Quarter of the Polygon Write $$N=4Q.$$ The implemented strategy treats one quarter of the optimal polygon as a convex monotone chain. A steep half-chain runs from the vertical boundary toward the diagonal, a mirrored half-chain finishes that quadrant, and reflections across the coordinate axes produce the remaining three quadrants. Inside the steep half-chain, every nontrivial edge direction is represented by a primitive lattice vector $$\gcd(x,y)=1,\qquad 1\le x<y.$$ Sorting these vectors by slope $x/y$ enforces the monotone turning order needed for convexity. The program therefore works with primitive directions in increasing slope order....

Detailed mathematical approach

Problem Summary

Let $A(N)$ be the minimum area of a convex lattice polygon with $N$ edges. The supplied C++, Python, and Java implementations solve the target case $N\equiv 0 \pmod{4}$ by writing $N=4Q$, constructing one monotone quarter-chain, and recovering the full polygon by symmetry. The search is therefore over ordered primitive edge directions and compact area summaries rather than over polygons themselves.

Mathematical Approach

The implementation organizes the polygon as four reflected copies of a quarter construction. That reduces the problem to a structured optimization on primitive lattice directions.

Step 1: Reduce to One Quarter of the Polygon

Write

$$N=4Q.$$

The implemented strategy treats one quarter of the optimal polygon as a convex monotone chain. A steep half-chain runs from the vertical boundary toward the diagonal, a mirrored half-chain finishes that quadrant, and reflections across the coordinate axes produce the remaining three quadrants.

Inside the steep half-chain, every nontrivial edge direction is represented by a primitive lattice vector

$$\gcd(x,y)=1,\qquad 1\le x<y.$$

Sorting these vectors by slope $x/y$ enforces the monotone turning order needed for convexity. The program therefore works with primitive directions in increasing slope order.

Step 2: Keep Only Directions That Can Still Fit

A quarter-chain has exactly $Q$ edges, but its two extreme directions are fixed by the construction, so only $Q-2$ internal primitive directions can be chosen freely. For a primitive vector $(x,y)$, the implementations precompute

$$\Pi(x,y)=\#\left\{(u,v):1\le u\le x,\ 1\le v\le y,\ u<v,\ \gcd(u,v)=1\right\}.$$

This is the number of primitive directions in the southwest rectangle up to $(x,y)$. If

$$\Pi(x,y)>Q-2,$$

then $(x,y)$ is already too deep in that ordered set to fit into a quarter-chain with only $Q$ total edges, so it is discarded immediately. This pruning step is why the candidate list stays manageable even when $Q$ is much larger than the visible number of states in the later dynamic program.

Step 3: Compress a Partial Chain into Two Numbers

Suppose the currently chosen internal directions in one steep half-chain are

$$ (x_1,y_1),\ (x_2,y_2),\ \dots,\ (x_k,y_k), $$

already sorted by slope. The implementations summarize this entire partial chain by two quantities:

$$\Sigma_k=1+2\sum_{i=1}^{k} y_i,$$

$$\Lambda_k=2\sum_{j=1}^{k} x_j\left(1+y_j+2\sum_{i<j} y_i\right).$$

The base state is

$$ (\Sigma,\Lambda)=(1,0). $$

Appending a new primitive direction $(x,y)$ updates the summary by

$$\Sigma'=\Sigma+2y,\qquad \Lambda'=\Lambda+2x(\Sigma+y).$$

This recurrence is the core of the dynamic program. The important point is that every future contribution of the already chosen directions is mediated only through $\Sigma$ and $\Lambda$, so no finer geometric description is needed once a state has been compressed into this pair.

Step 4: Prune by the Lower Convex Envelope

For each possible chain length, many reachable states are useless. If two states have the same $\Sigma$, only the one with smaller $\Lambda$ can ever help. The implementations go further and keep only the lower convex envelope of the reachable $(\Sigma,\Lambda)$ points.

The reason is visible in the final merge formula. If the opposite half of the quarter contributes $(\Sigma_1,\Lambda_1)$, then using a state $(\Sigma_0,\Lambda_0)$ on the current side leads to

$$\mathcal{V}=\Lambda_0+\Lambda_1+(\Sigma_0+2)(\Sigma_1+2)-2.$$

For fixed $(\Sigma_1,\Lambda_1)$, this is

$$\mathcal{V}=\bigl(\Lambda_0+(\Sigma_1+2)\Sigma_0\bigr)+C,$$

which is a linear functional of $(\Sigma_0,\Lambda_0)$ with positive slope $\Sigma_1+2$. Therefore any state lying above the lower convex envelope can never be optimal for any later merge, and the cross-product test used in the implementations deletes exactly those states.

Step 5: Meet in the Middle

The algorithm stores one filtered frontier for every possible number of quarter-edges. If one side of the quarter uses $\ell$ edges and the other side uses $Q-\ell$ edges, every pair of frontier states produces the candidate area

$$\mathcal{V}=\Lambda_0+\Lambda_1+(\Sigma_0+2)(\Sigma_1+2)-2.$$

The minimum is taken over all splits

$$1\le \ell < Q$$

and over all pairs of filtered states in the corresponding two frontiers. The special case $Q=1$, which means $N=4$, is handled directly and gives area $1$.

Worked Example: $N=8$

Here $Q=2$. Since a quarter-chain has only $Q-2=0$ internal slots, no extra primitive direction can be inserted. The only stored frontier state on each side is still the base state

$$ (\Sigma,\Lambda)=(1,0). $$

The only possible split is $1+1$, so the candidate value is

$$A(8)=0+0+(1+2)(1+2)-2=9-2=7.$$

This matches the checkpoint built into the implementations.

How the Code Works

The C++, Python, and Java implementations first reject inputs that are not multiples of $4$, then set $Q=N/4$. They enumerate primitive first-octant directions, count how many smaller rectangle-contained primitive directions each one has, and keep only those that satisfy the cutoff from Step 2. The surviving directions are then sorted by slope.

Next, the implementation keeps one frontier for each possible quarter-length. Starting from the base frontier $(1,0)$, every candidate direction is appended to every existing frontier of smaller length, creating new summarized states through the recurrence

$$\Sigma'=\Sigma+2y,\qquad \Lambda'=\Lambda+2x(\Sigma+y).$$

After each batch of insertions, the frontier is merged with the previously known states of that length and filtered back down to its lower convex envelope.

After all candidates have been processed, complementary frontiers of lengths $\ell$ and $Q-\ell$ are joined. For every pair, the program evaluates

$$\Lambda_0+\Lambda_1+(\Sigma_0+2)(\Sigma_1+2)-2$$

and keeps the smallest result. The C++ implementation also checks several known small values before evaluating the target case.

Complexity Analysis

Let $E$ be the number of retained primitive directions after the rectangle cutoff, and let $H_r$ be the size of the filtered frontier for quarter-length $r$. Building the primitive table and the two-dimensional prefix-count table costs $O(Q^2)$ memory. The gcd tests and the final slope sort contribute $O(Q^2\log Q+E\log E)$ time.

The dynamic-programming phase generates roughly

$$O\left(E\sum_{r=1}^{Q-1} H_r\right)$$

state extensions, followed by linear-time hull merges on the produced frontier lists. The final meet-in-the-middle stage costs

$$O\left(\sum_{r=1}^{Q-1} H_r H_{Q-r}\right).$$

So the practical running time is controlled far more by the filtered frontier sizes than by the astronomically large number of convex lattice polygons one might try to enumerate directly. Total memory usage is

$$O\left(Q^2+\sum_{r=1}^{Q} H_r\right).$$

Footnotes and References

Problem page: Project Euler 742
Convex polygon: Wikipedia - Convex polygon
Lattice polygon: Wikipedia - Lattice polygon
Coprime integers: Wikipedia - Coprime integers
Convex hull: Wikipedia - Convex hull
Dynamic programming: Wikipedia - Dynamic programming

Problem 742 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <vector>

using i64 = std::int64_t;

struct Edge {
    int x;
    int y;
};

static bool slope_less(const Edge& a, const Edge& b) {
    return static_cast<i64>(a.x) * b.y < static_cast<i64>(b.x) * a.y;
}

static std::vector<Edge> get_edge_candidates(int num_edges) {
    std::vector<std::vector<int>> primitive(num_edges + 1, std::vector<int>(num_edges + 1, 0));
    for (int y = 1; y <= num_edges; ++y) {
        for (int x = 1; x < y; ++x) {
            if (std::gcd(x, y) == 1) primitive[y][x] = 1;
        }
    }

    std::vector<std::vector<int>> pref(num_edges + 1, std::vector<int>(num_edges + 1, 0));
    for (int y = 1; y <= num_edges; ++y) {
        int row = 0;
        for (int x = 1; x <= num_edges; ++x) {
            row += primitive[y][x];
            pref[y][x] = pref[y - 1][x] + row;
        }
    }

    std::vector<Edge> edges;
    edges.reserve(num_edges * num_edges / 2);
    for (int y = 1; y <= num_edges; ++y) {
        for (int x = 1; x < y; ++x) {
            if (!primitive[y][x]) continue;
            const int count_smaller_edges = pref[y][x];
            if (2 + count_smaller_edges > num_edges) break;
            edges.push_back({x, y});
        }
    }
    std::sort(edges.begin(), edges.end(), slope_less);
    return edges;
}

using SAPair = std::pair<int, i64>;

static std::vector<SAPair> filter_convex_hull(std::vector<SAPair> domain) {
    std::sort(domain.begin(), domain.end());
    std::vector<SAPair> out;
    out.reserve(domain.size());

    for (const auto& cur : domain) {
        const int s = cur.first;
        const i64 a = cur.second;

        if (!out.empty() && out.back().first == s) continue;

        while (out.size() >= 2) {
            const int s1 = out[out.size() - 1].first;
            const i64 a1 = out[out.size() - 1].second;
            const int s2 = out[out.size() - 2].first;
            const i64 a2 = out[out.size() - 2].second;

            __int128 lhs = static_cast<__int128>(a - a1) * static_cast<__int128>(s - s2);
            __int128 rhs = static_cast<__int128>(a - a2) * static_cast<__int128>(s - s1);
            if (lhs >= rhs) break;
            out.pop_back();
        }

        if (!out.empty() && out.back().second <= a) continue;
        out.push_back(cur);
    }
    return out;
}

static std::vector<SAPair> combine_convex_hulls(const std::vector<SAPair>& d0,
                                                const std::vector<SAPair>& d1) {
    std::vector<SAPair> merged;
    merged.reserve(d0.size() + d1.size());
    std::size_t i = 0;
    std::size_t j = 0;
    while (i < d0.size() && j < d1.size()) {
        if (d0[i] <= d1[j]) {
            merged.push_back(d0[i++]);
        } else {
            merged.push_back(d1[j++]);
        }
    }
    while (i < d0.size()) merged.push_back(d0[i++]);
    while (j < d1.size()) merged.push_back(d1[j++]);
    return filter_convex_hull(std::move(merged));
}

static std::vector<SAPair> update_sa_domain(const std::vector<SAPair>& domain, const Edge& edge) {
    std::vector<SAPair> out;
    out.reserve(domain.size());
    const int x = edge.x;
    const int y = edge.y;
    for (const auto& [s, a] : domain) {
        const int ns = s + 2 * y;
        const i64 na = a + 2LL * x * (s + y);
        out.push_back({ns, na});
    }
    return out;
}

static i64 solve(int num_edges) {
    if (num_edges < 4 || (num_edges & 3) != 0) return -1;
    const int q_edges = num_edges / 4;
    if (q_edges == 1) return 1;

    std::vector<std::vector<SAPair>> domain(q_edges + 1);
    domain[1].push_back({1, 0});

    const auto edges = get_edge_candidates(q_edges);
    for (const auto& edge : edges) {
        for (int used = q_edges - 1; used >= 1; --used) {
            if (domain[used].empty()) continue;
            auto with_new = update_sa_domain(domain[used], edge);
            domain[used + 1] = combine_convex_hulls(domain[used + 1], with_new);
        }
    }

    i64 best = (1LL << 62);
    for (int left = 1; left < q_edges; ++left) {
        const auto& d0 = domain[left];
        const auto& d1 = domain[q_edges - left];
        for (const auto& [s0, a0] : d0) {
            for (const auto& [s1, a1] : d1) {
                const i64 val = a0 + a1 + static_cast<i64>(s0 + 2) * static_cast<i64>(s1 + 2) - 2;
                if (val < best) best = val;
            }
        }
    }
    return best;
}

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    if (solve(4) != 1) {
        std::cerr << "Validation failed: A(4) != 1\n";
        return 1;
    }
    if (solve(8) != 7) {
        std::cerr << "Validation failed: A(8) != 7\n";
        return 1;
    }
    if (solve(40) != 1039) {
        std::cerr << "Validation failed: A(40) != 1039\n";
        return 1;
    }
    if (solve(100) != 17473) {
        std::cerr << "Validation failed: A(100) != 17473\n";
        return 1;
    }

    std::cout << solve(1000) << '\n';
    return 0;
}

Python

import math
import functools

def gcd(a, b):
    while b:
        a, b = b, a % b
    return a

def cmp_edges(a, b):
    lhs = a[0] * b[1]
    rhs = b[0] * a[1]
    if lhs < rhs: return -1
    elif lhs > rhs: return 1
    return 0

def get_edge_candidates(num_edges):
    primitive = [[0] * (num_edges + 1) for _ in range(num_edges + 1)]
    for y in range(1, num_edges + 1):
        for x in range(1, y):
            if gcd(x, y) == 1:
                primitive[y][x] = 1
                
    pref = [[0] * (num_edges + 1) for _ in range(num_edges + 1)]
    for y in range(1, num_edges + 1):
        row = 0
        for x in range(1, num_edges + 1):
            row += primitive[y][x]
            pref[y][x] = pref[y - 1][x] + row
            
    edges = []
    for y in range(1, num_edges + 1):
        for x in range(1, y):
            if not primitive[y][x]:
                continue
            count_smaller_edges = pref[y][x]
            if 2 + count_smaller_edges > num_edges:
                break
            edges.append((x, y))
            
    edges.sort(key=functools.cmp_to_key(cmp_edges))
    return edges

def filter_convex_hull(domain):
    domain.sort()
    out = []
    
    for s, a in domain:
        if out and out[-1][0] == s:
            continue
            
        while len(out) >= 2:
            s1, a1 = out[-1]
            s2, a2 = out[-2]
            
            lhs = (a - a1) * (s - s2)
            rhs = (a - a2) * (s - s1)
            
            if lhs >= rhs:
                break
            out.pop()
            
        if out and out[-1][1] <= a:
            continue
            
        out.append((s, a))
        
    return out

def combine_convex_hulls(d0, d1):
    merged = []
    i = 0
    j = 0
    while i < len(d0) and j < len(d1):
        if d0[i] <= d1[j]:
            merged.append(d0[i])
            i += 1
        else:
            merged.append(d1[j])
            j += 1
            
    while i < len(d0):
        merged.append(d0[i])
        i += 1
    while j < len(d1):
        merged.append(d1[j])
        j += 1
        
    return filter_convex_hull(merged)

def update_sa_domain(domain, edge):
    out = []
    x, y = edge
    for s, a in domain:
        ns = s + 2 * y
        na = a + 2 * x * (s + y)
        out.append((ns, na))
    return out

def solve():
    num_edges = 1000
    if num_edges < 4 or num_edges % 4 != 0:
        return "-1"
        
    q_edges = num_edges // 4
    if q_edges == 1:
        return "1"
        
    domain = [[] for _ in range(q_edges + 1)]
    domain[1].append((1, 0))
    
    edges = get_edge_candidates(q_edges)
    
    for edge in edges:
        for used in range(q_edges - 1, 0, -1):
            if not domain[used]:
                continue
            with_new = update_sa_domain(domain[used], edge)
            domain[used + 1] = combine_convex_hulls(domain[used + 1], with_new)
            
    best = 1 << 62
    for left in range(1, q_edges):
        d0 = domain[left]
        d1 = domain[q_edges - left]
        for s0, a0 in d0:
            for s1, a1 in d1:
                val = a0 + a1 + (s0 + 2) * (s1 + 2) - 2
                if val < best:
                    best = val
                    
    return str(best)

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

public class Euler742 {

    static class Edge implements Comparable<Edge> {
        int x, y;

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

        @Override
        public int compareTo(Edge other) {
            long lhs = (long) this.x * other.y;
            long rhs = (long) other.x * this.y;
            return Long.compare(lhs, rhs);
        }
    }

    static int gcd(int a, int b) {
        while (b != 0) {
            int t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    static List<Edge> getEdgeCandidates(int numEdges) {
        int[][] primitive = new int[numEdges + 1][numEdges + 1];
        for (int y = 1; y <= numEdges; ++y) {
            for (int x = 1; x < y; ++x) {
                if (gcd(x, y) == 1) {
                    primitive[y][x] = 1;
                }
            }
        }

        int[][] pref = new int[numEdges + 1][numEdges + 1];
        for (int y = 1; y <= numEdges; ++y) {
            int row = 0;
            for (int x = 1; x <= numEdges; ++x) {
                row += primitive[y][x];
                pref[y][x] = pref[y - 1][x] + row;
            }
        }

        List<Edge> edges = new ArrayList<>();
        for (int y = 1; y <= numEdges; ++y) {
            for (int x = 1; x < y; ++x) {
                if (primitive[y][x] == 0)
                    continue;
                int countSmallerEdges = pref[y][x];
                if (2 + countSmallerEdges > numEdges)
                    break;
                edges.add(new Edge(x, y));
            }
        }
        Collections.sort(edges);
        return edges;
    }

    static class SAPair implements Comparable<SAPair> {
        int s;
        long a;

        SAPair(int s, long a) {
            this.s = s;
            this.a = a;
        }

        @Override
        public int compareTo(SAPair o) {
            if (this.s != o.s) {
                return Integer.compare(this.s, o.s);
            }
            return Long.compare(this.a, o.a);
        }
    }

    static List<SAPair> filterConvexHull(List<SAPair> domain) {
        Collections.sort(domain);
        List<SAPair> out = new ArrayList<>(domain.size());

        for (SAPair cur : domain) {
            int s = cur.s;
            long a = cur.a;

            if (!out.isEmpty() && out.get(out.size() - 1).s == s)
                continue;

            while (out.size() >= 2) {
                int s1 = out.get(out.size() - 1).s;
                long a1 = out.get(out.size() - 1).a;
                int s2 = out.get(out.size() - 2).s;
                long a2 = out.get(out.size() - 2).a;

                BigInteger diffAS1 = BigInteger.valueOf(a - a1);
                BigInteger diffSS2 = BigInteger.valueOf(s - s2);
                BigInteger diffAS2 = BigInteger.valueOf(a - a2);
                BigInteger diffSS1 = BigInteger.valueOf(s - s1);

                BigInteger lhs = diffAS1.multiply(diffSS2);
                BigInteger rhs = diffAS2.multiply(diffSS1);

                if (lhs.compareTo(rhs) >= 0)
                    break;
                out.remove(out.size() - 1);
            }

            if (!out.isEmpty() && out.get(out.size() - 1).a <= a)
                continue;
            out.add(cur);
        }
        return out;
    }

    static List<SAPair> combineConvexHulls(List<SAPair> d0, List<SAPair> d1) {
        List<SAPair> merged = new ArrayList<>(d0.size() + d1.size());
        int i = 0, j = 0;
        while (i < d0.size() && j < d1.size()) {
            if (d0.get(i).compareTo(d1.get(j)) <= 0) {
                merged.add(d0.get(i++));
            } else {
                merged.add(d1.get(j++));
            }
        }
        while (i < d0.size())
            merged.add(d0.get(i++));
        while (j < d1.size())
            merged.add(d1.get(j++));
        return filterConvexHull(merged);
    }

    static List<SAPair> updateSADomain(List<SAPair> domain, Edge edge) {
        List<SAPair> out = new ArrayList<>(domain.size());
        int x = edge.x;
        int y = edge.y;
        for (SAPair p : domain) {
            int ns = p.s + 2 * y;
            long na = p.a + 2L * x * (p.s + y);
            out.add(new SAPair(ns, na));
        }
        return out;
    }

    public static String solve() {
        int numEdges = 1000;
        if (numEdges < 4 || (numEdges % 4) != 0)
            return "-1";
        int qEdges = numEdges / 4;
        if (qEdges == 1)
            return "1";

        List<List<SAPair>> domain = new ArrayList<>(qEdges + 1);
        for (int i = 0; i <= qEdges; ++i) {
            domain.add(new ArrayList<>());
        }
        domain.get(1).add(new SAPair(1, 0L));

        List<Edge> edges = getEdgeCandidates(qEdges);

        for (Edge edge : edges) {
            for (int used = qEdges - 1; used >= 1; --used) {
                if (domain.get(used).isEmpty())
                    continue;
                List<SAPair> withNew = updateSADomain(domain.get(used), edge);
                domain.set(used + 1, combineConvexHulls(domain.get(used + 1), withNew));
            }
        }

        long best = 1L << 62;
        for (int left = 1; left < qEdges; ++left) {
            List<SAPair> d0 = domain.get(left);
            List<SAPair> d1 = domain.get(qEdges - left);
            for (SAPair p0 : d0) {
                for (SAPair p1 : d1) {
                    long val = p0.a + p1.a + (long) (p0.s + 2) * (p1.s + 2) - 2;
                    if (val < best) {
                        best = val;
                    }
                }
            }
        }

        return Long.toString(best);
    }

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