Problem 681: Maximal Area

Project Euler Problem 681 Solution

EulerSolve provides an optimized solution for Project Euler Problem 681, Maximal Area, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each positive integer $m$, let $P(m)$ denote the sum of the perimeters of all integer-sided quadrilaterals whose maximal possible area is exactly $m$. The target quantity is $$SP(n)=\sum_{m=1}^{n} P(m).$$ The implementations verify the approach with the checkpoints $SP(10)=186$ and $SP(100)=23238$, then evaluate the full case $n=10^6$. Mathematical Approach The geometry becomes manageable after rewriting the maximal-area condition as a factorization problem for $m^2$. Step 1: Start from the maximal-area formula For side lengths $a,b,c,d$, write the semiperimeter as $$s=\frac{a+b+c+d}{2}.$$ Among all quadrilaterals with these side lengths, the area is maximized when the quadrilateral is cyclic....

Detailed mathematical approach

Problem Summary

For each positive integer $m$, let $P(m)$ denote the sum of the perimeters of all integer-sided quadrilaterals whose maximal possible area is exactly $m$. The target quantity is

$$SP(n)=\sum_{m=1}^{n} P(m).$$

The implementations verify the approach with the checkpoints $SP(10)=186$ and $SP(100)=23238$, then evaluate the full case $n=10^6$.

Mathematical Approach

The geometry becomes manageable after rewriting the maximal-area condition as a factorization problem for $m^2$.

Step 1: Start from the maximal-area formula

For side lengths $a,b,c,d$, write the semiperimeter as

$$s=\frac{a+b+c+d}{2}.$$

Among all quadrilaterals with these side lengths, the area is maximized when the quadrilateral is cyclic. Brahmagupta's formula then gives

$$A_{\max}^2=(s-a)(s-b)(s-c)(s-d).$$

Problem 681 asks us to study the cases where $A_{\max}=m$, so the defining equation is

$$m^2=(s-a)(s-b)(s-c)(s-d).$$

Step 2: Introduce a symmetric substitution

Set

$$x=s-a,\qquad y=s-b,\qquad z=s-c,\qquad w=s-d.$$

Then the area condition becomes

$$xyzw=m^2.$$

The side lengths are recovered by solving the linear system:

$$a=\frac{-x+y+z+w}{2},\qquad b=\frac{x-y+z+w}{2},$$

$$c=\frac{x+y-z+w}{2},\qquad d=\frac{x+y+z-w}{2}.$$

The perimeter is especially simple:

$$a+b+c+d=x+y+z+w.$$

So once a valid quadruple $(x,y,z,w)$ is known, its contribution to $P(m)$ is exactly the sum $x+y+z+w$.

Step 3: Characterize the valid quadruples

Permuting the side lengths only permutes $x,y,z,w$, so we sort them and count only

$$x\ge y\ge z\ge w\ge 1.$$

Under this ordering, $a$ is the smallest recovered side because it subtracts the largest term $x$. Therefore positivity of all four sides is equivalent to the single inequality

$$a>0\iff x<y+z+w.$$

Integrality is controlled by parity. Since each recovered side is half of an integer expression, it is enough to require

$$x+y+z+w\equiv 0 \pmod{2}.$$

Hence valid quadrilaterals with maximal area $m$ are in bijection with ordered positive integer quadruples satisfying

$$xyzw=m^2,\qquad x\ge y\ge z\ge w,\qquad x<y+z+w,\qquad x+y+z+w\equiv 0\pmod{2}.$$

Step 4: Reduce the search to divisors of $m^2$

$$m=\prod_i p_i^{e_i},$$

then

$$m^2=\prod_i p_i^{2e_i}.$$

Every candidate value $w,z,y,x$ must therefore be a divisor of $m^2$. The implementation generates the full divisor list of $m^2$, sorts it, chooses $w\le z\le y$, and determines the last factor by

$$x=\frac{m^2}{wzy}.$$

This automatically enforces the product constraint, and the sorted divisor list makes monotone early stopping possible: once a partial product is too large, all later divisors are too large as well.

Step 5: Worked example for $m=4$

Here $m^2=16$. The ordered factor quadruples $(x,y,z,w)$ with product $16$ are

$$\begin{aligned} &(16,1,1,1),\\ &(8,2,1,1),\\ &(4,4,1,1),\\ &(4,2,2,1),\\ &(2,2,2,2). \end{aligned}$$

The first two fail $x<y+z+w$. The fourth has odd perimeter $4+2+2+1=9$, so it does not produce integer side lengths. The remaining two are valid:

$$\begin{aligned} (4,4,1,1)&\longrightarrow (a,b,c,d)=(1,1,4,4),\quad x+y+z+w=10,\\ (2,2,2,2)&\longrightarrow (a,b,c,d)=(2,2,2,2),\quad x+y+z+w=8. \end{aligned}$$

Therefore

$$P(4)=10+8=18.$$

Step 6: Sum over all maximal areas

For each fixed $m$, the program adds the perimeters of all valid quadruples to obtain $P(m)$. The final answer is then

$$SP(n)=\sum_{m=1}^{n}P(m).$$

Because different values of $m$ are independent, the outer summation is a natural place for parallel execution.

How the Code Works

The C++, Python, and Java implementations begin by preparing prime-factor information for all integers up to $n$. For a given $m$, the implementation factors $m$, doubles the exponents to describe $m^2$, and recursively generates every divisor of $m^2$.

After sorting those divisors, the implementation scans candidates in the order $w\le z\le y$. Whenever a partial product already forces the missing factor to fall below the required order, the loop stops early. If the remaining quotient is divisible by the chosen factors, that quotient becomes $x$, and the code applies exactly the two mathematical filters derived above: the positivity inequality $x<y+z+w$ and the parity condition $x+y+z+w\equiv 0\pmod{2}$.

Each surviving quadruple contributes its perimeter $x+y+z+w$ to $P(m)$. The C++ and Java implementations distribute distinct values of $m$ across worker threads, while the Python implementation reuses the compiled C++ computation path and returns the parsed final result.

Complexity Analysis

Let $D(m)=\tau(m^2)$ be the number of divisors of $m^2$. Building the prime-factor table up to $n$ uses $O(n)$ memory and $O(n\log\log n)$ time in the usual sieve analysis, with the C++ version using a linear smallest-prime-factor construction. For one fixed $m$, generating all divisors costs $O(D(m))$. The subsequent ordered scan over divisor triples has worst-case $O(D(m)^3)$ behavior, but the monotone break conditions and divisibility tests cut away most combinations in practice. Extra working memory beyond the sieve is dominated by the divisor list for one value of $m$.

Footnotes and References

Problem page: https://projecteuler.net/problem=681
Brahmagupta's formula: Wikipedia — Brahmagupta's formula
Cyclic quadrilateral: Wikipedia — Cyclic quadrilateral
Divisor function: Wikipedia — Divisor function

Problem 681 source code

C++

#include <algorithm>
#include <cassert>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

std::vector<int> build_spf(int n) {
    std::vector<int> spf(static_cast<std::size_t>(n + 1), 0);
    std::vector<int> primes;
    primes.reserve(static_cast<std::size_t>(n / 10));

    for (int i = 2; i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0) {
            spf[static_cast<std::size_t>(i)] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            const u64 v = static_cast<u64>(p) * static_cast<u64>(i);
            if (v > static_cast<u64>(n) || p > spf[static_cast<std::size_t>(i)]) {
                break;
            }
            spf[static_cast<std::size_t>(v)] = p;
        }
    }

    return spf;
}

std::vector<std::pair<u64, int>> factor_square(int m, const std::vector<int>& spf) {
    std::vector<std::pair<u64, int>> fac;
    int x = m;
    while (x > 1) {
        const int p = spf[static_cast<std::size_t>(x)];
        int e = 0;
        do {
            x /= p;
            ++e;
        } while (x % p == 0);
        fac.push_back({static_cast<u64>(p), 2 * e});
    }
    return fac;
}

void build_divisors_rec(int idx, u64 cur, const std::vector<std::pair<u64, int>>& fac,
                        std::vector<u64>& out) {
    if (idx == static_cast<int>(fac.size())) {
        out.push_back(cur);
        return;
    }

    const u64 p = fac[static_cast<std::size_t>(idx)].first;
    const int e = fac[static_cast<std::size_t>(idx)].second;
    u64 mul = 1ULL;
    for (int i = 0; i <= e; ++i) {
        build_divisors_rec(idx + 1, cur * mul, fac, out);
        mul *= p;
    }
}

u128 solve_single_m(int m, const std::vector<int>& spf, std::vector<u64>& divs) {
    const u64 n2 = static_cast<u64>(m) * static_cast<u64>(m);

    const auto fac = factor_square(m, spf);
    divs.clear();
    divs.reserve(256);
    build_divisors_rec(0, 1ULL, fac, divs);
    std::sort(divs.begin(), divs.end());

    const int L = static_cast<int>(divs.size());
    u128 ans = 0;

    for (int i = 0; i < L; ++i) {
        const u64 w = divs[static_cast<std::size_t>(i)];

        for (int j = i; j < L; ++j) {
            const u64 z = divs[static_cast<std::size_t>(j)];
            if (w > n2 / z) {
                break;
            }
            const u64 wz = w * z;
            if (n2 % wz != 0ULL) {
                continue;
            }

            const u64 rem = n2 / wz;
            if (z > rem / z) {
                break;
            }

            for (int k = j; k < L; ++k) {
                const u64 y = divs[static_cast<std::size_t>(k)];
                if (y > rem / y) {
                    break;
                }
                if (rem % y != 0ULL) {
                    continue;
                }

                const u64 x = rem / y;
                if (x < y) {
                    continue;
                }
                if (x >= y + z + w) {
                    continue;
                }
                const u64 p = x + y + z + w;
                if ((p & 1ULL) != 0ULL) {
                    continue;
                }

                ans += static_cast<u128>(p);
            }
        }
    }

    return ans;
}

u64 SP(int n) {
    const std::vector<int> spf = build_spf(n);

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) {
        threads = 8;
    }

    std::atomic<int> next_m{1};
    std::vector<u128> partial(threads, 0);
    std::vector<std::thread> workers;
    workers.reserve(threads);

    for (unsigned tid = 0; tid < threads; ++tid) {
        workers.emplace_back([&, tid]() {
            std::vector<u64> divs;
            u128 local = 0;
            while (true) {
                const int m = next_m.fetch_add(1, std::memory_order_relaxed);
                if (m > n) {
                    break;
                }
                local += solve_single_m(m, spf, divs);
            }
            partial[tid] = local;
        });
    }

    for (auto& t : workers) {
        t.join();
    }

    u128 ans = 0;
    for (u128 x : partial) {
        ans += x;
    }
    return static_cast<u64>(ans);
}

}  // namespace

int main() {
    assert(SP(10) == 186ULL);
    assert(SP(100) == 23'238ULL);

    std::cout << SP(1'000'000) << "\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.Arrays;
import java.util.List;
import java.util.stream.LongStream;

public class Euler681 {

    static int[] buildSpf(int n) {
        int[] spf = new int[n + 1];
        for (int i = 2; i * i <= n; ++i) {
            if (spf[i] == 0) {
                for (int j = i * i; j <= n; j += i) {
                    if (spf[j] == 0)
                        spf[j] = i;
                }
            }
        }
        for (int i = 2; i <= n; ++i) {
            if (spf[i] == 0)
                spf[i] = i;
        }
        return spf;
    }

    static class Pair {
        long p;
        int e;

        Pair(long p, int e) {
            this.p = p;
            this.e = e;
        }
    }

    static List<Pair> factorSquare(int m, int[] spf) {
        List<Pair> fac = new ArrayList<>();
        int x = m;
        while (x > 1) {
            int p = spf[x];
            int e = 0;
            do {
                x /= p;
                e++;
            } while (x % p == 0);
            fac.add(new Pair(p, 2 * e));
        }
        return fac;
    }

    static void buildDivs(int idx, long cur, List<Pair> fac, long[] out, int[] count) {
        if (idx == fac.size()) {
            out[count[0]++] = cur;
            return;
        }
        long p = fac.get(idx).p;
        int e = fac.get(idx).e;
        long mul = 1;
        for (int i = 0; i <= e; ++i) {
            buildDivs(idx + 1, cur * mul, fac, out, count);
            mul *= p;
        }
    }

    static long solveSingleM(int m, int[] spf) {
        long n2 = (long) m * m;
        List<Pair> fac = factorSquare(m, spf);

        int totalDivs = 1;
        for (Pair pair : fac)
            totalDivs *= (pair.e + 1);

        long[] divs = new long[totalDivs];
        int[] count = { 0 };
        buildDivs(0, 1L, fac, divs, count);
        Arrays.sort(divs);

        int L = totalDivs;
        long ans = 0;

        for (int i = 0; i < L; ++i) {
            long w = divs[i];
            for (int j = i; j < L; ++j) {
                long z = divs[j];
                if (w > n2 / z)
                    break;
                long wz = w * z;
                if (n2 % wz != 0)
                    continue;

                long rem = n2 / wz;
                if (z > rem / z)
                    break;

                for (int k = j; k < L; ++k) {
                    long y = divs[k];
                    if (y > rem / y)
                        break;
                    if (rem % y != 0)
                        continue;

                    long x = rem / y;
                    if (x < y)
                        continue;
                    if (x >= y + z + w)
                        continue;

                    long p = x + y + z + w;
                    if (p % 2 != 0)
                        continue;

                    ans += p;
                }
            }
        }
        return ans;
    }

    public static String solve() {
        int n = 1000000;
        int[] spf = buildSpf(n);

        long ans = LongStream.rangeClosed(1, n)
                .parallel()
                .map(m -> solveSingleM((int) m, spf))
                .sum();

        return Long.toString(ans);
    }

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