Problem 583: Heron Envelopes

Project Euler Problem 583 Solution

EulerSolve provides an optimized solution for Project Euler Problem 583, Heron Envelopes, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A Heron envelope can be modeled as a rectangle of width $2a$ and height $h$, topped by an isosceles flap whose base is also $2a$, whose equal side length is $s$, and whose altitude is $t$. Its outer perimeter is $$P=2(a+h+s).$$ The task is to compute $S(p)$, the sum of all such perimeters with $P\le p$, subject to the geometric condition $t<h$ and the integrality conditions encoded by the three reference implementations. Writing $L=p/2$ is convenient, because every valid envelope then satisfies $$a+h+s\le L.$$ Mathematical Approach The key observation is that every admissible envelope is governed by three right triangles sharing the same half-width $a$. Once those right triangles are parameterized, the problem becomes a structured matching problem rather than a brute-force search over polygons. Step 1: Convert the Geometry into Three Square Conditions The flap itself splits into two congruent right triangles, so its side length $s$ must satisfy $$a^2+t^2=s^2.$$ The rectangle has full width $2a$, so its diagonal condition is $$\left(2a\right)^2+h^2=d^2$$ for some integer $d$. The long diagonal running from a lower corner to the top vertex of the flap has horizontal offset $a$ and vertical offset $h+t$, so it must satisfy $$a^2+(h+t)^2=q^2$$ for some integer $q$....

Detailed mathematical approach

Problem Summary

A Heron envelope can be modeled as a rectangle of width $2a$ and height $h$, topped by an isosceles flap whose base is also $2a$, whose equal side length is $s$, and whose altitude is $t$. Its outer perimeter is

$$P=2(a+h+s).$$

The task is to compute $S(p)$, the sum of all such perimeters with $P\le p$, subject to the geometric condition $t<h$ and the integrality conditions encoded by the three reference implementations. Writing $L=p/2$ is convenient, because every valid envelope then satisfies

$$a+h+s\le L.$$

Mathematical Approach

The key observation is that every admissible envelope is governed by three right triangles sharing the same half-width $a$. Once those right triangles are parameterized, the problem becomes a structured matching problem rather than a brute-force search over polygons.

Step 1: Convert the Geometry into Three Square Conditions

The flap itself splits into two congruent right triangles, so its side length $s$ must satisfy

$$a^2+t^2=s^2.$$

The rectangle has full width $2a$, so its diagonal condition is

$$\left(2a\right)^2+h^2=d^2$$

for some integer $d$.

The long diagonal running from a lower corner to the top vertex of the flap has horizontal offset $a$ and vertical offset $h+t$, so it must satisfy

$$a^2+(h+t)^2=q^2$$

for some integer $q$.

Therefore a valid envelope is exactly a quadruple $(a,h,t,s)$ with

$$a^2+t^2=s^2,\qquad \left(2a\right)^2+h^2=d^2,\qquad a^2+(h+t)^2=q^2,\qquad t<h,\qquad a+h+s\le L.$$

Step 2: Enumerate All Relevant Right Triangles

Every primitive Pythagorean triple can be written as

$$x=m^2-n^2,\qquad y=2mn,\qquad z=m^2+n^2,$$

with $m>n$, $\gcd(m,n)=1$, and opposite parity. Multiplying by a positive scale $k$ gives every non-primitive triple:

$$x=k(m^2-n^2),\qquad y=k(2mn),\qquad z=k(m^2+n^2).$$

For flap candidates, one leg is $a$, the other is $t$, and the hypotenuse is $s$. Since either leg may play the role of $a$, the implementations keep both leg orders from each generated triple.

For rectangle candidates, the full width must be $2a$. So whenever a generated leg is even, that leg can be halved to obtain $a$, while the other leg becomes $h$. This is why the rectangle catalogue is built from Pythagorean triples with an even leg.

The perimeter bound gives safe search limits. Because $a+h+s\le L$, we only need flap triples with $s\le L$, and rectangle triples with $2a\le 2L$ and $h\le L$.

Step 3: Group Everything by the Shared Half-Width

Define

$$T(a)=\left\{y\in\mathbb{Z}_{>0}:a^2+y^2\text{ is a perfect square}\right\}.$$

Then the flap condition says $t\in T(a)$, and the long-diagonal condition says $h+t\in T(a)$.

The rectangle condition is separate: $h$ must belong to the set of heights for which $\left(2a\right)^2+h^2$ is a square.

This immediately suggests the data organization used by the implementations: build all flap candidates and all rectangle candidates, then group both collections by the same value of $a$. Each group can be processed independently.

Step 4: Match Rectangle Heights with Flap Heights

Fix one half-width $a$. For every rectangle height $h$ in that group and every flap candidate $(t,s)$ in the same group, we test three conditions:

$$t<h,\qquad s\le L-a-h,\qquad h+t\in T(a).$$

The first inequality is the sensible-envelope condition. The second is just the perimeter bound rewritten from $a+h+s\le L$. The third guarantees that the long diagonal is integral.

If all three tests pass, then the envelope is valid and contributes

$$2(a+h+s)$$

to $S(p)$.

This is much cheaper than checking arbitrary polygons, because the search only combines already-valid right-triangle building blocks.

Step 5: Why Sorting Helps

Within a fixed $a$-group, the flap candidates are sorted by $t$. For fixed $a$, the flap side length is

$$s=\sqrt{a^2+t^2},$$

so increasing $t$ also increases $s$. Therefore, once $t\ge h$, no later flap can satisfy $t<h$; and once $s>L-a-h$, no later flap can satisfy the perimeter bound. Both facts allow early exits in the inner scan.

The rectangle heights are also sorted, so if $L-a-h\le 0$, no larger height can work either. These monotonicity properties are exactly what make the grouped scan practical.

Worked Example

One concrete valid envelope is obtained from

$$a=60,\qquad h=119,\qquad t=25,\qquad s=65.$$

The flap is integral because

$$60^2+25^2=65^2.$$

The rectangle diagonal is integral because

$$120^2+119^2=169^2.$$

The long diagonal is also integral because

$$60^2+(119+25)^2=60^2+144^2=156^2.$$

Since $25<119$, the flap is lower than the rectangle height as required. Its perimeter is

$$P=2(60+119+65)=488.$$

So this envelope contributes $488$ to $S(p)$ whenever $p\ge 488$.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical plan. They first set $L=p/2$, generate all flap right triangles up to hypotenuse $L$, and generate all rectangle candidates by taking Pythagorean triples whose even leg can serve as the full width $2a$.

Next, they regroup both catalogues by the shared half-width $a$. Inside each group, flap candidates are stored as altitude-side pairs $(t,s)$, and rectangle candidates are stored by their heights $h$. A lookup structure containing all admissible values of $t$ for that $a$ lets the implementation test the long-diagonal condition $h+t\in T(a)$ in constant expected time.

Each group is then scanned in sorted order. For every rectangle height, the implementation checks flap candidates until either $t\ge h$ or $s>L-a-h$. Every successful match adds $2(a+h+s)$ to the running total. The C++ version additionally parallelizes these independent $a$-groups for large limits, while the Python and Java versions keep the same logic in a single thread.

Complexity Analysis

Let $T$ be the number of generated flap candidates and $R$ the number of generated rectangle candidates after the Euclidean enumeration and scaling steps. Building those candidate lists is $O(T+R)$ in output size, up to the arithmetic cost of the gcd filters and loop bounds in the triple generator.

Grouping and sorting cost

$$O(T\log T+R\log R)$$

overall. The matching phase is best described groupwise: if $M$ is the total number of flap-rectangle pairs actually inspected before the early breaks trigger, then the scan costs $O(M)$, and each long-diagonal test is $O(1)$ on average thanks to the lookup structure.

Memory usage is linear in the stored candidates, plus the auxiliary membership structure used for the $h+t$ test. In practice, the method is efficient because it never searches arbitrary envelopes directly; it only combines prevalidated right triangles with the same half-width.

Footnotes and References

Problem page: Project Euler 583 - Heron Envelopes
Pythagorean triples: Wikipedia - Pythagorean triple
Euclid's formula: Wikipedia - Generating a Pythagorean triple
Isosceles triangles: Wikipedia - Isosceles triangle
Heron's formula: Wikipedia - Heron's formula

Problem 583 source code

C++

#include <pthread.h>
#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <unistd.h>
#include <vector>

using u32 = std::uint32_t;
using u64 = std::uint64_t;

struct Tri {
    int a;
    int t;
    int s;
};

struct Rect {
    int a;
    int h;
};

struct GroupRange {
    int a;
    std::size_t tri_begin;
    std::size_t tri_end;
    std::size_t rect_begin;
    std::size_t rect_end;
};

static std::vector<Tri> generate_triangles(int L) {
    std::vector<Tri> tris;

    const int mmax = (int)std::sqrt((long double)L) + 2;
    for (int m = 2; m <= mmax; ++m) {
        for (int n = 1; n < m; ++n) {
            if (((m - n) & 1) == 0) continue;
            if (std::gcd(m, n) != 1) continue;

            const int a0 = m * m - n * n;
            const int b0 = 2 * m * n;
            const int c0 = m * m + n * n;
            if (c0 > L) continue;

            const int kmax = L / c0;
            for (int k = 1; k <= kmax; ++k) {
                const int a = k * a0;
                const int t = k * b0;
                const int s = k * c0;
                tris.push_back(Tri{a, t, s});
                tris.push_back(Tri{t, a, s});
            }
        }
    }

    return tris;
}

static std::vector<Rect> generate_rectangles(int L) {
    const int max_leg = 2 * L;
    std::vector<Rect> rects;

    const int mmax = (int)std::sqrt((long double)max_leg) + 2;
    for (int m = 2; m <= mmax; ++m) {
        for (int n = 1; n < m; ++n) {
            if (((m - n) & 1) == 0) continue;
            if (std::gcd(m, n) != 1) continue;

            const int a0 = m * m - n * n;
            const int b0 = 2 * m * n;
            const int max0 = std::max(a0, b0);
            if (max0 > max_leg) continue;

            const int kmax = max_leg / max0;
            for (int k = 1; k <= kmax; ++k) {
                const int x = k * a0;
                const int y = k * b0;

                if ((x & 1) == 0) {
                    const int a = x / 2;
                    const int h = y;
                    if (a > 0 && a <= L && h > 0 && h <= L) rects.push_back(Rect{a, h});
                }
                if ((y & 1) == 0) {
                    const int a = y / 2;
                    const int h = x;
                    if (a > 0 && a <= L && h > 0 && h <= L) rects.push_back(Rect{a, h});
                }
            }
        }
    }

    return rects;
}

static u64 compute_group_sum(const GroupRange& g,
                             const std::vector<Tri>& tris,
                             const std::vector<Rect>& rects,
                             const int L,
                             std::vector<u32>& mark,
                             u32& stamp) {
    if (++stamp == 0U) {
        std::fill(mark.begin(), mark.end(), 0U);
        stamp = 1U;
    }
    for (std::size_t t = g.tri_begin; t < g.tri_end; ++t) {
        const int tv = tris[t].t;
        if (tv >= 0 && tv <= L) {
            mark[static_cast<std::size_t>(tv)] = stamp;
        }
    }

    const int a = g.a;
    u64 sum = 0ULL;
    for (std::size_t r = g.rect_begin; r < g.rect_end; ++r) {
        const int h = rects[r].h;
        if (h <= 0) continue;
        const int limit_s = L - a - h;
        if (limit_s <= 0) break;

        for (std::size_t t = g.tri_begin; t < g.tri_end; ++t) {
            const int tv = tris[t].t;
            if (tv >= h) break;
            const int s = tris[t].s;
            if (s > limit_s) break;

            const int u = h + tv;
            if (u > L) continue;
            if (mark[static_cast<std::size_t>(u)] != stamp) continue;

            sum += 2ULL * (u64)(a + h + s);
        }
    }
    return sum;
}

static int detect_thread_count(std::size_t work_items) {
    long cores = ::sysconf(_SC_NPROCESSORS_ONLN);
    int threads = (cores > 0) ? static_cast<int>(cores) : 4;
    if (threads < 1) threads = 1;
    if (work_items == 0) return 1;
    if ((std::size_t)threads > work_items) threads = static_cast<int>(work_items);
    return threads;
}

struct WorkerTask {
    const std::vector<Tri>* tris = nullptr;
    const std::vector<Rect>* rects = nullptr;
    const std::vector<GroupRange>* groups = nullptr;
    int L = 0;
    std::atomic<std::size_t>* next_idx = nullptr;
    std::vector<u32> mark;
    u32 stamp = 0U;
    u64 partial = 0ULL;
};

static void* worker_entry(void* raw) {
    auto* task = static_cast<WorkerTask*>(raw);
    u64 local = 0ULL;
    if (task->mark.empty()) {
        task->mark.resize(static_cast<std::size_t>(task->L + 1), 0U);
        task->stamp = 0U;
    }
    while (true) {
        const std::size_t idx = task->next_idx->fetch_add(1, std::memory_order_relaxed);
        if (idx >= task->groups->size()) break;
        local += compute_group_sum((*task->groups)[idx], *task->tris, *task->rects, task->L,
                                   task->mark, task->stamp);
    }
    task->partial = local;
    return nullptr;
}

static u64 S(u64 p) {
    const int L = (int)(p / 2);
    const auto tris_raw = generate_triangles(L);
    const auto rects_raw = generate_rectangles(L);

    std::vector<u32> tri_count(static_cast<std::size_t>(L + 1), 0U);
    std::vector<u32> rect_count(static_cast<std::size_t>(L + 1), 0U);
    for (const auto& t : tris_raw) {
        if (t.a >= 1 && t.a <= L) {
            ++tri_count[static_cast<std::size_t>(t.a)];
        }
    }
    for (const auto& r : rects_raw) {
        if (r.a >= 1 && r.a <= L) {
            ++rect_count[static_cast<std::size_t>(r.a)];
        }
    }

    std::vector<u32> tri_off(static_cast<std::size_t>(L + 2), 0U);
    std::vector<u32> rect_off(static_cast<std::size_t>(L + 2), 0U);
    for (int a = 1; a <= L; ++a) {
        tri_off[static_cast<std::size_t>(a + 1)] =
            tri_off[static_cast<std::size_t>(a)] + tri_count[static_cast<std::size_t>(a)];
        rect_off[static_cast<std::size_t>(a + 1)] =
            rect_off[static_cast<std::size_t>(a)] + rect_count[static_cast<std::size_t>(a)];
    }

    std::vector<Tri> tris(static_cast<std::size_t>(tri_off[static_cast<std::size_t>(L + 1)]));
    std::vector<Rect> rects(static_cast<std::size_t>(rect_off[static_cast<std::size_t>(L + 1)]));

    std::vector<u32> tri_cursor = tri_off;
    std::vector<u32> rect_cursor = rect_off;
    for (const auto& t : tris_raw) {
        if (t.a < 1 || t.a > L) continue;
        u32& pos = tri_cursor[static_cast<std::size_t>(t.a)];
        tris[static_cast<std::size_t>(pos++)] = t;
    }
    for (const auto& r : rects_raw) {
        if (r.a < 1 || r.a > L) continue;
        u32& pos = rect_cursor[static_cast<std::size_t>(r.a)];
        rects[static_cast<std::size_t>(pos++)] = r;
    }

    std::vector<GroupRange> groups;
    groups.reserve(std::min(tris.size(), rects.size()) / 16);
    for (int a = 1; a <= L; ++a) {
        const std::size_t tb = static_cast<std::size_t>(tri_off[static_cast<std::size_t>(a)]);
        const std::size_t te = static_cast<std::size_t>(tri_off[static_cast<std::size_t>(a + 1)]);
        const std::size_t rb = static_cast<std::size_t>(rect_off[static_cast<std::size_t>(a)]);
        const std::size_t re = static_cast<std::size_t>(rect_off[static_cast<std::size_t>(a + 1)]);
        if (tb == te || rb == re) continue;
        groups.push_back(GroupRange{a, tb, te, rb, re});
    }
    if (groups.empty()) return 0ULL;

    for (const auto& g : groups) {
        std::sort(tris.begin() + static_cast<std::ptrdiff_t>(g.tri_begin),
                  tris.begin() + static_cast<std::ptrdiff_t>(g.tri_end),
                  [](const Tri& x, const Tri& y) { return x.t < y.t; });
        std::sort(rects.begin() + static_cast<std::ptrdiff_t>(g.rect_begin),
                  rects.begin() + static_cast<std::ptrdiff_t>(g.rect_end),
                  [](const Rect& x, const Rect& y) { return x.h < y.h; });
    }

    if (L <= 100'000) {
        u64 sum = 0ULL;
        std::vector<u32> mark(static_cast<std::size_t>(L + 1), 0U);
        u32 stamp = 0U;
        for (const auto& g : groups) {
            sum += compute_group_sum(g, tris, rects, L, mark, stamp);
        }
        return sum;
    }

    const int threads = detect_thread_count(groups.size());
    if (threads <= 1) {
        u64 sum = 0ULL;
        std::vector<u32> mark(static_cast<std::size_t>(L + 1), 0U);
        u32 stamp = 0U;
        for (const auto& g : groups) {
            sum += compute_group_sum(g, tris, rects, L, mark, stamp);
        }
        return sum;
    }

    std::vector<pthread_t> handles(static_cast<std::size_t>(threads));
    std::vector<WorkerTask> tasks(static_cast<std::size_t>(threads));
    std::atomic<std::size_t> next_idx{0};

    for (int t = 0; t < threads; ++t) {
        auto& task = tasks[static_cast<std::size_t>(t)];
        task.tris = &tris;
        task.rects = &rects;
        task.groups = &groups;
        task.L = L;
        task.next_idx = &next_idx;
        task.mark.clear();
        task.stamp = 0U;
        task.partial = 0ULL;
        pthread_create(&handles[static_cast<std::size_t>(t)], nullptr, worker_entry, &task);
    }

    u64 sum = 0ULL;
    for (int t = 0; t < threads; ++t) {
        pthread_join(handles[static_cast<std::size_t>(t)], nullptr);
        sum += tasks[static_cast<std::size_t>(t)].partial;
    }
    return sum;
}

int main() {
    const u64 s1e4 = S(10000);
    if (s1e4 != 884680ULL) {
        std::cerr << "Validation failed: S(1e4) got " << s1e4 << "\n";
        return 1;
    }

    std::cout << S(10000000ULL) << "\n";
    return 0;
}

Python

import math

def solve():
    P = 10000000
    L = P // 2

    def gen_triples(limit):
        tris = []
        mmax = int(limit**0.5) + 2
        for m in range(2, mmax+1):
            for nn in range(1, m):
                if (m - nn) % 2 == 0: continue
                if math.gcd(m, nn) != 1: continue
                a0 = m*m - nn*nn; b0 = 2*m*nn; c0 = m*m + nn*nn
                if c0 > limit: continue
                kmax = limit // c0
                for k in range(1, kmax+1):
                    a, t, s = k*a0, k*b0, k*c0
                    tris.append((a, t, s)); tris.append((t, a, s))
        return tris

    def gen_rects(limit):
        max_leg = 2*limit; rects = []
        mmax = int(max_leg**0.5) + 2
        for m in range(2, mmax+1):
            for nn in range(1, m):
                if (m-nn) % 2 == 0: continue
                if math.gcd(m, nn) != 1: continue
                a0 = m*m - nn*nn; b0 = 2*m*nn
                max0 = max(a0, b0)
                if max0 > max_leg: continue
                kmax = max_leg // max0
                for k in range(1, kmax+1):
                    x, y = k*a0, k*b0
                    if x % 2 == 0:
                        a, h = x//2, y
                        if 0 < a <= limit and 0 < h <= limit: rects.append((a, h))
                    if y % 2 == 0:
                        a, h = y//2, x
                        if 0 < a <= limit and 0 < h <= limit: rects.append((a, h))
        return rects

    tris = gen_triples(L)
    rects = gen_rects(L)

    # Group by 'a'
    from collections import defaultdict
    tri_by_a = defaultdict(list); rect_by_a = defaultdict(list)
    for a, t, s in tris: tri_by_a[a].append((t, s))
    for a, h in rects: rect_by_a[a].append(h)

    total = 0
    for a in range(1, L+1):
        if a not in tri_by_a or a not in rect_by_a: continue
        ts = sorted(tri_by_a[a], key=lambda x: x[0])
        rs = sorted(rect_by_a[a])
        t_set = set(t for t, s in ts)
        for h in rs:
            limit_s = L - a - h
            if limit_s <= 0: break
            for t, s in ts:
                if t >= h: break
                if s > limit_s: break
                u = h + t
                if u > L: continue
                if u in t_set:
                    total += 2 * (a + h + s)

    return str(total)

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

Java

import java.util.*;

public class Euler583 {
    public static String solve() {
        int P = 10000000, L = P / 2;
        List<int[]> tris = genTriples(L);
        List<int[]> rects = genRects(L);
        Map<Integer, List<int[]>> triByA = new HashMap<>();
        Map<Integer, List<Integer>> rectByA = new HashMap<>();
        for (int[] t : tris)
            triByA.computeIfAbsent(t[0], k -> new ArrayList<>()).add(new int[] { t[1], t[2] });
        for (int[] r : rects)
            rectByA.computeIfAbsent(r[0], k -> new ArrayList<>()).add(r[1]);
        long total = 0;
        for (int a = 1; a <= L; a++) {
            if (!triByA.containsKey(a) || !rectByA.containsKey(a))
                continue;
            List<int[]> ts = triByA.get(a);
            ts.sort(Comparator.comparingInt(x -> x[0]));
            List<Integer> rs = rectByA.get(a);
            Collections.sort(rs);
            Set<Integer> tSet = new HashSet<>();
            for (int[] t : ts)
                tSet.add(t[0]);
            for (int h : rs) {
                int limitS = L - a - h;
                if (limitS <= 0)
                    break;
                for (int[] t : ts) {
                    if (t[0] >= h)
                        break;
                    if (t[1] > limitS)
                        break;
                    int u = h + t[0];
                    if (u > L)
                        continue;
                    if (tSet.contains(u))
                        total += 2L * (a + h + t[1]);
                }
            }
        }
        return String.valueOf(total);
    }

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

    static List<int[]> genTriples(int limit) {
        List<int[]> tris = new ArrayList<>();
        int mmax = (int) Math.sqrt(limit) + 2;
        for (int m = 2; m <= mmax; m++)
            for (int n = 1; n < m; n++) {
                if ((m - n) % 2 == 0 || gcd(m, n) != 1)
                    continue;
                int a0 = m * m - n * n, b0 = 2 * m * n, c0 = m * m + n * n;
                if (c0 > limit)
                    continue;
                for (int k = 1; k * c0 <= limit; k++) {
                    tris.add(new int[] { k * a0, k * b0, k * c0 });
                    tris.add(new int[] { k * b0, k * a0, k * c0 });
                }
            }
        return tris;
    }

    static List<int[]> genRects(int limit) {
        List<int[]> rects = new ArrayList<>();
        int maxLeg = 2 * limit, mmax = (int) Math.sqrt(maxLeg) + 2;
        for (int m = 2; m <= mmax; m++)
            for (int n = 1; n < m; n++) {
                if ((m - n) % 2 == 0 || gcd(m, n) != 1)
                    continue;
                int a0 = m * m - n * n, b0 = 2 * m * n;
                if (Math.max(a0, b0) > maxLeg)
                    continue;
                for (int k = 1; k * Math.max(a0, b0) <= maxLeg; k++) {
                    int x = k * a0, y = k * b0;
                    if (x % 2 == 0) {
                        int a = x / 2;
                        if (a > 0 && a <= limit && y > 0 && y <= limit)
                            rects.add(new int[] { a, y });
                    }
                    if (y % 2 == 0) {
                        int a = y / 2;
                        if (a > 0 && a <= limit && x > 0 && x <= limit)
                            rects.add(new int[] { a, x });
                    }
                }
            }
        return rects;
    }

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