Problem 283: Integer Sided Triangles with Integral Area/perimeter Ratio

Project Euler Problem 283 Solution

EulerSolve provides an optimized solution for Project Euler Problem 283, Integer Sided Triangles with Integral Area/perimeter Ratio, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let a triangle with integer sides $a,b,c$ have area $A$ and perimeter $P=a+b+c$. We want all triangles for which $$\frac{A}{P}=k$$ is an integer. For each $k$, let $p(k)$ be the sum of the corresponding perimeters. The program computes $$\sum_{k=1}^{1000} p(k),$$ but the final numeric answer is intentionally omitted here. Mathematical Approach 1) Replace side lengths by semiperimeter gaps. Let $$s=\frac{a+b+c}{2},\qquad x=s-a,\qquad y=s-b,\qquad z=s-c.$$ Then $x,y,z$ are positive integers, and conversely $$a=y+z,\qquad b=z+x,\qquad c=x+y,\qquad s=x+y+z.$$ So every integer triangle corresponds to one positive integer triple $(x,y,z)$. If we sort the sides as $a \ge b \ge c$, then $x \le y \le z$, which is exactly the ordering used in the code to avoid duplicates. 2) Heron's formula gives a Diophantine equation. Heron's formula is $$A^2=s(s-a)(s-b)(s-c)=sxyz.$$ The condition $A/P=k$ means $$A=k(a+b+c)=2ks.$$ Squaring and cancelling one factor of $s$ gives $$xyz=4k^2s=4k^2(x+y+z).$$ If we define $$n=(2k)^2=4k^2,$$ the triangle problem becomes the pure integer equation $$xyz=n(x+y+z),\qquad x \le y \le z.$$ The perimeter of the recovered triangle is $$P=2s=2(x+y+z),$$ which is why the program adds $2(x+y+z)$ for each valid triple. 3) Why this is a bijection. Starting from a triangle, we get one positive triple $(x,y,z)$....

Detailed mathematical approach

Problem Summary

Let a triangle with integer sides $a,b,c$ have area $A$ and perimeter $P=a+b+c$. We want all triangles for which

$$\frac{A}{P}=k$$

is an integer. For each $k$, let $p(k)$ be the sum of the corresponding perimeters. The program computes

$$\sum_{k=1}^{1000} p(k),$$

but the final numeric answer is intentionally omitted here.

Mathematical Approach

1) Replace side lengths by semiperimeter gaps. Let

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

Then $x,y,z$ are positive integers, and conversely

$$a=y+z,\qquad b=z+x,\qquad c=x+y,\qquad s=x+y+z.$$

So every integer triangle corresponds to one positive integer triple $(x,y,z)$. If we sort the sides as $a \ge b \ge c$, then $x \le y \le z$, which is exactly the ordering used in the code to avoid duplicates.

2) Heron's formula gives a Diophantine equation. Heron's formula is

$$A^2=s(s-a)(s-b)(s-c)=sxyz.$$

The condition $A/P=k$ means

$$A=k(a+b+c)=2ks.$$

Squaring and cancelling one factor of $s$ gives

$$xyz=4k^2s=4k^2(x+y+z).$$

If we define

$$n=(2k)^2=4k^2,$$

the triangle problem becomes the pure integer equation

$$xyz=n(x+y+z),\qquad x \le y \le z.$$

The perimeter of the recovered triangle is

$$P=2s=2(x+y+z),$$

which is why the program adds $2(x+y+z)$ for each valid triple.

3) Why this is a bijection. Starting from a triangle, we get one positive triple $(x,y,z)$. Starting from a positive solution of $xyz=n(x+y+z)$, the inverse map

$$a=y+z,\qquad b=z+x,\qquad c=x+y$$

automatically satisfies the triangle inequalities because $x,y,z>0$. So no valid triangle is lost and no extra object is counted.

4) Solve one variable in terms of the other two. Rearranging the main equation gives

$$z=\frac{n(x+y)}{xy-n}.$$

Therefore $xy>n$ is necessary, and once $x$ and $y$ are fixed there is at most one possible $z$. The brute-force checkpoint routine in the code uses exactly this formula.

5) The key factorization for the fast algorithm. For the production algorithm, fixing $k$ and $x$ is enough. Starting from

$$xyz=n(x+y+z),$$

one checks that

$$(xy-n)(xz-n)=x^2yz-nx(y+z)+n^2=n(x^2+n).$$

Define

$$m=x^2+n,\qquad M=nm,\qquad u=xy-n,\qquad v=xz-n.$$

Then

$$uv=M,$$

and the original variables are recovered by

$$y=\frac{u+n}{x},\qquad z=\frac{v+n}{x}.$$

So for fixed $(k,x)$, every divisor pair $(u,v)$ of $M$ gives a candidate, and the only remaining tests are:

$$x \mid (u+n),\qquad x \mid (v+n),\qquad y \ge x,\qquad z \ge y.$$

Because $y \le z$ is equivalent to $u \le v$, the code only needs divisor pairs with $u^2 \le M$.

6) Why the bound $x \le \sqrt{3n}$ is correct. For fixed $x$, the function

$$f(y,z)=\frac{xyz}{x+y+z}$$

is increasing in both $y$ and $z$ on positive inputs. Since $y,z \ge x$, the smallest possible value is at $y=z=x$, so

$$n=f(y,z)\ge \frac{x^3}{3x}=\frac{x^2}{3}.$$

Hence

$$x \le \sqrt{3n}.$$

This is exactly the outer bound used by both the fast solver and the brute-force verifier.

7) The brute-force upper bound for $y$. In the verifier, the extra condition $z \ge y$ and the formula for $z$ imply

$$\frac{n(x+y)}{xy-n}\ge y,$$

which rearranges to

$$xy^2-2ny-nx \le 0.$$

Solving this quadratic inequality yields

$$y \le \frac{n+\sqrt{n(n+x^2)}}{x},$$

which is the exact bound used in the checkpoint routine.

Worked Examples

Example 1: $k=1$. Then $n=4$. Take $x=2$. We get

$$M=n(x^2+n)=4(4+4)=32.$$

Choose the divisor pair $(u,v)=(4,8)$. Then

$$y=\frac{4+4}{2}=4,\qquad z=\frac{8+4}{2}=6.$$

So $(x,y,z)=(2,4,6)$, which gives the triangle

$$a=y+z=10,\qquad b=z+x=8,\qquad c=x+y=6.$$

Its semiperimeter is $s=12$, its area is $\sqrt{12\cdot 2\cdot 4\cdot 6}=24$, and its perimeter is $24$. Therefore $A/P=1$, exactly as required.

Example 2: $k=2$. Then $n=16$. The triple $(x,y,z)=(6,7,8)$ satisfies

$$6\cdot 7\cdot 8 = 16(6+7+8)=336.$$

The corresponding sides are

$$a=15,\qquad b=14,\qquad c=13,$$

so we recover the classical $13$-$14$-$15$ triangle. Its perimeter is $42$, its area is $84$, and $84/42=2$.

Why the Code Is Correct

The code first chooses $k$, hence $n=(2k)^2$. It then loops over all admissible $x$, factors

$$M=n(x^2+n),$$

generates all divisors of $M$, keeps only $u$ with $u^2 \le M$, reconstructs $y$ and $z$, and accepts exactly the candidates that satisfy the divisibility and ordering constraints. By the derivation above, every accepted candidate is a valid triangle, and every valid triangle appears exactly once.

Checks and Complexity

The file contains two explicit checkpoints. The fast method and the brute-force method both give

$$\sum_{k \le 5} p(k)=140098,\qquad \sum_{k \le 10} p(k)=3781786.$$

For complexity, the expensive part is factoring $m=x^2+n$ and enumerating divisors of $M=n(x^2+n)$. The SPF sieve is built up to

$$4(2K)^2,$$

because $x^2 \le 3n$ implies $m=x^2+n \le 4n$. In practice this makes each factorization cheap, and divisor enumeration is far smaller than brute forcing all $(y,z)$.

Problem 283 source code

C++

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <thread>
#include <utility>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;
using i64 = std::int64_t;

constexpr int kDefaultLimit = 1000;

struct Options {
    int k_limit = kDefaultLimit;
    unsigned threads = std::thread::hardware_concurrency();
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg, const std::string& prefix, unsigned& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    unsigned parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10U + static_cast<unsigned>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--k-limit=", options.k_limit)) {
            continue;
        }
        if (parse_unsigned_after_prefix(arg, "--threads=", options.threads)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.k_limit <= 0) {
        return false;
    }
    if (options.threads == 0) {
        options.threads = 1;
    }
    return true;
}

std::vector<int> build_spf(const int limit) {
    std::vector<int> spf(static_cast<std::size_t>(limit + 1), 0);
    if (limit >= 1) {
        spf[1] = 1;
    }
    for (int i = 2; i <= limit; ++i) {
        if (spf[static_cast<std::size_t>(i)] != 0) {
            continue;
        }
        spf[static_cast<std::size_t>(i)] = i;
        if (static_cast<i64>(i) * i > limit) {
            continue;
        }
        for (int j = i * i; j <= limit; j += i) {
            if (spf[static_cast<std::size_t>(j)] == 0) {
                spf[static_cast<std::size_t>(j)] = i;
            }
        }
    }
    return spf;
}

void factorize_int(int n, const std::vector<int>& spf, std::vector<std::pair<int, int>>& out) {
    out.clear();
    while (n > 1) {
        const int p = spf[static_cast<std::size_t>(n)];
        int cnt = 0;
        while (n % p == 0) {
            n /= p;
            ++cnt;
        }
        out.push_back({p, cnt});
    }
}

void merge_factorizations(const std::vector<std::pair<int, int>>& a,
                         const std::vector<std::pair<int, int>>& b,
                         std::vector<std::pair<int, int>>& out) {
    out.clear();
    std::size_t i = 0;
    std::size_t j = 0;

    while (i < a.size() || j < b.size()) {
        if (j == b.size() || (i < a.size() && a[i].first < b[j].first)) {
            out.push_back(a[i]);
            ++i;
        } else if (i == a.size() || b[j].first < a[i].first) {
            out.push_back(b[j]);
            ++j;
        } else {
            out.push_back({a[i].first, a[i].second + b[j].second});
            ++i;
            ++j;
        }
    }
}

void generate_divisors(const std::vector<std::pair<int, int>>& factors, std::vector<u64>& divisors) {
    divisors.clear();
    divisors.push_back(1);

    for (const auto& [p, e] : factors) {
        const std::size_t base = divisors.size();
        u64 mul = 1;
        for (int i = 1; i <= e; ++i) {
            mul *= static_cast<u64>(p);
            for (std::size_t j = 0; j < base; ++j) {
                divisors.push_back(divisors[j] * mul);
            }
        }
    }
}

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }
    std::string s;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        s.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

u128 solve_range(const int k_start, const int k_end, const std::vector<int>& spf) {
    std::vector<std::pair<int, int>> factors_n;
    std::vector<std::pair<int, int>> factors_m;
    std::vector<std::pair<int, int>> factors_all;
    std::vector<u64> divisors;

    u128 total = 0;

    for (int k = k_start; k <= k_end; ++k) {
        const u64 r = static_cast<u64>(2 * k);
        const u64 n = r * r;
        factorize_int(static_cast<int>(n), spf, factors_n);

        const u64 x_max = static_cast<u64>(std::sqrt(static_cast<long double>(3 * n)));
        for (u64 x = 1; x <= x_max; ++x) {
            const u64 m = x * x + n;
            factorize_int(static_cast<int>(m), spf, factors_m);
            merge_factorizations(factors_n, factors_m, factors_all);
            generate_divisors(factors_all, divisors);

            const u64 M = n * m;
            for (u64 u : divisors) {
                if (u * u > M) {
                    continue;
                }

                const u64 v = M / u;
                if ((u + n) % x != 0 || (v + n) % x != 0) {
                    continue;
                }

                const u64 y = (u + n) / x;
                if (y < x) {
                    continue;
                }
                const u64 z = (v + n) / x;

                total += static_cast<u128>(2) * (x + y + z);
            }
        }
    }

    return total;
}

u128 solve_fast(const int k_limit, const std::vector<int>& spf, unsigned threads) {
    if (threads <= 1 || k_limit < 100) {
        return solve_range(1, k_limit, spf);
    }

    const unsigned use_threads = std::min<unsigned>(threads, static_cast<unsigned>(k_limit));
    const int chunk = (k_limit + static_cast<int>(use_threads) - 1) / static_cast<int>(use_threads);

    std::vector<std::thread> pool;
    std::vector<u128> partial(use_threads, 0);
    pool.reserve(use_threads);

    for (unsigned t = 0; t < use_threads; ++t) {
        const int start = static_cast<int>(t) * chunk + 1;
        const int end = std::min(k_limit, start + chunk - 1);
        if (start > end) {
            continue;
        }

        pool.emplace_back([&, start, end, t]() {
            partial[t] = solve_range(start, end, spf);
        });
    }

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

    u128 total = 0;
    for (u128 part : partial) {
        total += part;
    }
    return total;
}

u128 solve_bruteforce_small(const int k_limit) {
    u128 total = 0;

    for (int k = 1; k <= k_limit; ++k) {
        const u64 r = static_cast<u64>(2 * k);
        const u64 n = r * r;

        const u64 x_max = static_cast<u64>(std::sqrt(static_cast<long double>(3 * n)));
        for (u64 x = 1; x <= x_max; ++x) {
            const long double disc = std::sqrt(static_cast<long double>(n) * (static_cast<long double>(n) + x * x));
            const u64 y_max = static_cast<u64>((n + disc) / x + 1e-12L);

            for (u64 y = x; y <= y_max; ++y) {
                const i64 d = static_cast<i64>(x * y) - static_cast<i64>(n);
                if (d <= 0) {
                    continue;
                }

                const u64 num = n * (x + y);
                if (num % static_cast<u64>(d) != 0) {
                    continue;
                }

                const u64 z = num / static_cast<u64>(d);
                if (z < y) {
                    continue;
                }

                if (x * y * z != n * (x + y + z)) {
                    continue;
                }

                total += static_cast<u128>(2) * (x + y + z);
            }
        }
    }

    return total;
}

bool run_checkpoints(const std::vector<int>& spf) {
    const u128 fast5 = solve_fast(5, spf, 1);
    const u128 brute5 = solve_bruteforce_small(5);
    if (fast5 != brute5 || fast5 != static_cast<u128>(140098ULL)) {
        std::cerr << "Checkpoint failed for k<=5: fast=" << to_string_u128(fast5)
                  << ", brute=" << to_string_u128(brute5) << '\n';
        return false;
    }

    const u128 fast10 = solve_fast(10, spf, 1);
    const u128 brute10 = solve_bruteforce_small(10);
    if (fast10 != brute10 || fast10 != static_cast<u128>(3781786ULL)) {
        std::cerr << "Checkpoint failed for k<=10: fast=" << to_string_u128(fast10)
                  << ", brute=" << to_string_u128(brute10) << '\n';
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    // max(x^2 + n) with n=(2k)^2 and x<=sqrt(3n) is at most 4*n
    const int n_max = 4 * (2 * options.k_limit) * (2 * options.k_limit);
    const std::vector<int> spf = build_spf(n_max);

    if (options.run_checkpoints && !run_checkpoints(spf)) {
        return 2;
    }

    const u128 answer = solve_fast(options.k_limit, spf, options.threads);
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

import math

def solve():
    k_limit = 1000
    n_max = 4 * (2 * k_limit) ** 2

    # Build smallest prime factor
    spf = list(range(n_max + 1))
    spf[0] = 0
    if n_max >= 1:
        spf[1] = 1
    for i in range(2, math.isqrt(n_max) + 1):
        if spf[i] == i:
            for j in range(i*i, n_max + 1, i):
                if spf[j] == j:
                    spf[j] = i

    def factorize(n):
        factors = []
        while n > 1:
            p = spf[n]
            e = 0
            while n % p == 0:
                n //= p
                e += 1
            factors.append((p, e))
        return factors

    def merge_factors(a, b):
        result = []
        i = j = 0
        while i < len(a) or j < len(b):
            if j == len(b) or (i < len(a) and a[i][0] < b[j][0]):
                result.append(a[i]); i += 1
            elif i == len(a) or b[j][0] < a[i][0]:
                result.append(b[j]); j += 1
            else:
                result.append((a[i][0], a[i][1] + b[j][1])); i += 1; j += 1
        return result

    def gen_divisors(factors):
        divs = [1]
        for p, e in factors:
            base = len(divs)
            mul = 1
            for _ in range(e):
                mul *= p
                for j in range(base):
                    divs.append(divs[j] * mul)
        return divs

    total = 0
    for k in range(1, k_limit + 1):
        r = 2 * k
        n = r * r
        factors_n = factorize(n)

        x_max = int(math.sqrt(3 * n))
        for x in range(1, x_max + 1):
            m = x * x + n
            factors_m = factorize(m)
            factors_all = merge_factors(factors_n, factors_m)
            divs = gen_divisors(factors_all)

            M = n * m
            for u in divs:
                if u * u > M:
                    continue
                v = M // u
                if (u + n) % x != 0 or (v + n) % x != 0:
                    continue
                y = (u + n) // x
                if y < x:
                    continue
                z = (v + n) // x
                total += 2 * (x + y + z)

    return str(total)

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

Java

import java.util.*;
import java.util.concurrent.*;

public class Euler283 {
    static class Pair {
        int first, second;

        Pair(int f, int s) {
            first = f;
            second = s;
        }
    }

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

    static void factorizeInt(int n, int[] spf, List<Pair> out) {
        out.clear();
        while (n > 1) {
            int p = spf[n];
            int cnt = 0;
            while (n % p == 0) {
                n /= p;
                ++cnt;
            }
            out.add(new Pair(p, cnt));
        }
    }

    static void mergeFactorizations(List<Pair> a, List<Pair> b, List<Pair> out) {
        out.clear();
        int i = 0, j = 0;
        while (i < a.size() || j < b.size()) {
            if (j == b.size() || (i < a.size() && a.get(i).first < b.get(j).first)) {
                out.add(a.get(i));
                ++i;
            } else if (i == a.size() || b.get(j).first < a.get(i).first) {
                out.add(b.get(j));
                ++j;
            } else {
                out.add(new Pair(a.get(i).first, a.get(i).second + b.get(j).second));
                ++i;
                ++j;
            }
        }
    }

    static void generateDivisors(List<Pair> factors, List<Long> divisors) {
        divisors.clear();
        divisors.add(1L);
        for (Pair p : factors) {
            int base = divisors.size();
            long mul = 1;
            for (int i = 1; i <= p.second; ++i) {
                mul *= p.first;
                for (int j = 0; j < base; ++j) {
                    divisors.add(divisors.get(j) * mul);
                }
            }
        }
    }

    static long solveRange(int kStart, int kEnd, int[] spf) {
        List<Pair> factorsN = new ArrayList<>();
        List<Pair> factorsM = new ArrayList<>();
        List<Pair> factorsAll = new ArrayList<>();
        List<Long> divisors = new ArrayList<>();

        long total = 0;

        for (int k = kStart; k <= kEnd; ++k) {
            long r = 2L * k;
            long n = r * r;
            factorizeInt((int) n, spf, factorsN);

            long xMax = (long) Math.sqrt(3L * n);
            for (long x = 1; x <= xMax; ++x) {
                long m = x * x + n;
                factorizeInt((int) m, spf, factorsM);
                mergeFactorizations(factorsN, factorsM, factorsAll);
                generateDivisors(factorsAll, divisors);

                long M = n * m;
                for (long u : divisors) {
                    if (u > M / u)
                        continue;
                    long v = M / u;
                    if ((u + n) % x != 0 || (v + n) % x != 0)
                        continue;

                    long y = (u + n) / x;
                    if (y < x)
                        continue;
                    long z = (v + n) / x;

                    total += 2L * (x + y + z);
                }
            }
        }
        return total;
    }

    public static String solve() {
        int kLimit = 1000;
        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());

        int nMax = 4 * (2 * kLimit) * (2 * kLimit);
        int[] spf = buildSpf(nMax);

        if (threads <= 1 || kLimit < 100) {
            long ans = solveRange(1, kLimit, spf);
            return String.valueOf(ans);
        }

        int useThreads = Math.min(threads, kLimit);
        int chunk = (kLimit + useThreads - 1) / useThreads;

        ExecutorService executor = Executors.newFixedThreadPool(useThreads);
        List<Future<Long>> futures = new ArrayList<>();

        for (int t = 0; t < useThreads; ++t) {
            final int start = t * chunk + 1;
            final int end = Math.min(kLimit, start + chunk - 1);
            if (start > end)
                continue;

            futures.add(executor.submit(() -> solveRange(start, end, spf)));
        }

        long total = 0;
        for (Future<Long> f : futures) {
            try {
                total += f.get();
            } catch (Exception e) {
                e.printStackTrace();
            }
        }
        executor.shutdown();

        return String.valueOf(total);
    }

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