Problem 390: Triangles with Non Rational Sides and Integral Area

Project Euler Problem 390 Solution

EulerSolve provides an optimized solution for Project Euler Problem 390, Triangles with Non Rational Sides and Integral Area, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The repository solutions work with the triangle family whose side lengths are $$\sqrt{b^2+1},\qquad \sqrt{c^2+1},\qquad \sqrt{b^2+c^2},$$ where $b$ and $c$ are positive integers and the brute-force validator counts only unordered pairs $b \le c$. If $A$ denotes the area, the goal is $$S(L)=\sum A,$$ summed over all triangles in this family whose area is an integer and satisfies $A \le L$. Mathematical Approach Step 1: Reduce the Geometry to an Integer Square Test Let $$x=\sqrt{b^2+1},\qquad y=\sqrt{c^2+1},\qquad z=\sqrt{b^2+c^2}.$$ Applying Heron's identity in the form $$16A^2=2x^2y^2+2y^2z^2+2z^2x^2-x^4-y^4-z^4$$ and substituting $x^2=b^2+1$, $y^2=c^2+1$, $z^2=b^2+c^2$ gives $$16A^2=4\left(b^2c^2+b^2+c^2\right).$$ Therefore $$A=\frac{1}{2}\sqrt{b^2c^2+b^2+c^2}.$$ This is exactly the relation used by the C++ brute-force checkpoint. Writing $$m^2=b^2c^2+b^2+c^2,$$ we have $A=m/2$. So the area is integral precisely when $m$ is an even integer square root of that expression. Step 2: Parity Forces Even Parameters Because $A$ is an integer, $m=2A$ is even, so $m^2 \equiv 0 \pmod 4$. But a square is only $0$ or $1$ modulo $4$. If either $b$ or $c$ were odd, then $$b^2c^2+b^2+c^2 \equiv 1 \text{ or } 3 \pmod 4,$$ which is impossible. Hence every valid solution has even $b$ and even $c$....

Detailed mathematical approach

Problem Summary

The repository solutions work with the triangle family whose side lengths are

$$\sqrt{b^2+1},\qquad \sqrt{c^2+1},\qquad \sqrt{b^2+c^2},$$

where $b$ and $c$ are positive integers and the brute-force validator counts only unordered pairs $b \le c$. If $A$ denotes the area, the goal is

$$S(L)=\sum A,$$

summed over all triangles in this family whose area is an integer and satisfies $A \le L$.

Mathematical Approach

Step 1: Reduce the Geometry to an Integer Square Test

Let

$$x=\sqrt{b^2+1},\qquad y=\sqrt{c^2+1},\qquad z=\sqrt{b^2+c^2}.$$

Applying Heron's identity in the form

$$16A^2=2x^2y^2+2y^2z^2+2z^2x^2-x^4-y^4-z^4$$

and substituting $x^2=b^2+1$, $y^2=c^2+1$, $z^2=b^2+c^2$ gives

$$16A^2=4\left(b^2c^2+b^2+c^2\right).$$

Therefore

$$A=\frac{1}{2}\sqrt{b^2c^2+b^2+c^2}.$$

This is exactly the relation used by the C++ brute-force checkpoint. Writing

$$m^2=b^2c^2+b^2+c^2,$$

we have $A=m/2$. So the area is integral precisely when $m$ is an even integer square root of that expression.

Step 2: Parity Forces Even Parameters

Because $A$ is an integer, $m=2A$ is even, so $m^2 \equiv 0 \pmod 4$. But a square is only $0$ or $1$ modulo $4$. If either $b$ or $c$ were odd, then

$$b^2c^2+b^2+c^2 \equiv 1 \text{ or } 3 \pmod 4,$$

which is impossible. Hence every valid solution has even $b$ and even $c$. We may therefore write

$$b=2p,\qquad c=2q,$$

with integers $p,q \ge 1$. Substituting into the area formula gives

$$A^2=4p^2q^2+p^2+q^2.$$

This is the equation actually encoded by the fast state generator.

Step 3: A Pell-Type Equation for Fixed $p$

For a fixed seed $p$, move the $q$-term to the left:

$$A^2-(4p^2+1)q^2=p^2.$$

If we define

$$D_p=4p^2+1,$$

then every admissible state satisfies the Pell-type norm equation

$$A^2-D_p q^2=p^2.$$

The trivial solution is $q=0$, $A=p$. The fast solver starts from this trivial point and moves to the next positive one on the same Pell branch.

Step 4: Derive the Linear Recurrence Used in Code

The fundamental unit for $x^2-D_p y^2=1$ used by the implementation is

$$u_p=(8p^2+1)+4p\sqrt{D_p},$$

because

$$\left(8p^2+1\right)^2-D_p(4p)^2=1.$$

Multiplying one solution $A+q\sqrt{D_p}$ by $u_p$ produces the next solution on the same branch:

$$A'+q'\sqrt{D_p}=u_p\left(A+q\sqrt{D_p}\right).$$

Comparing coefficients yields

$$q'=(8p^2+1)q+4pA,$$

$$A'=4p(4p^2+1)q+(8p^2+1)A.$$

This is exactly the transition implemented in all three solution files with

$$a=8p^2+1,\qquad b=4p,\qquad c=4p(4p^2+1),$$

$$q'=a q+b A,\qquad A'=c q+a A.$$

Step 5: First Nontrivial Area and the Seed Bound

Starting from the trivial state $(p,0,p)$, one recurrence step gives

$$q_1=4p^2,\qquad A_1=(8p^2+1)p=8p^3+p.$$

So the first actual triangle for seed $p$ is

$$b=2p,\qquad c=2q_1=8p^2,\qquad A=8p^3+p.$$

This explains the helper function seed_first_area(p) and the binary search in max_seed(limit): once $8p^3+p>L$, that seed cannot contribute any valid area $\le L$.

For $p=1$, the first solution is

$$q_1=4,\qquad A_1=9,$$

so $(b,c)=(2,8)$. The checkpoint in the C++ file confirms

$$2^2\cdot 8^2+2^2+8^2=324=18^2,\qquad A=18/2=9.$$

Step 6: Why the Solver Uses Three Branches

The equation

$$A^2=4p^2q^2+p^2+q^2$$

is symmetric in $p$ and $q$, and it only involves their squares. Therefore whenever the recurrence produces $(p,q',A')$, the states

$$(q',p,A')\qquad \text{and} \qquad (q',-p,A')$$

represent valid orientations of the same quadratic surface. The implementation therefore pushes three states:

$$(p,q',A'),\qquad (q',p,A'),\qquad (q',-p,A').$$

The first one continues the same Pell branch with fixed $p$. The two swapped states let the new value $q'$ become the fixed coordinate of later Pell branches. The negative-sign branch is essential: for example, from area $9$ the state $(4,-1,9)$ leads to the next solution $(4,15,121)$, which corresponds to $(b,c)=(8,30)$.

The repository does not store a separate formal proof of uniqueness for this tree, but it does verify the construction exhaustively against brute force for a small limit and against the statement checkpoint $S(10^6)=18018206$.

How the Code Works

The C++ file contains both a slow validator and the optimized recurrence solver. The validator iterates over unordered pairs $b \le c$, uses

$$m^2=b^2c^2+b^2+c^2$$

and checks whether $m$ is an even square root. The bounds

$$bc \le 2A \le 2L$$

justify the scan limit $b \le 2L$ and the inner bound $c \lesssim 2L/b$.

The fast solver works in half-coordinates $(p,q)$. Function sum_for_seed performs an explicit depth-first traversal with a stack of states $(p,q,A)$. Function next_state applies the Pell step above, discards branches with $A' \le 0$ or $A' > L$, and creates the three follow-up states. The C++ version uses 128-bit intermediates and distributes seeds across threads with an atomic counter; Python relies on arbitrary-precision integers, and Java uses BigInteger for the recurrence arithmetic.

Complexity Analysis

The brute-force validator is roughly

$$\sum_{b \le 2L} O\!\left(\frac{L}{b}\right)=O(L\log L),$$

which is already too slow for the real limit. The optimized solver is output-sensitive: if $T(L)$ is the number of generated states whose next area stays within the bound, then the runtime is $O(T(L))$ arithmetic transitions. Memory is proportional to the active DFS stack size per worker, not to a full $(b,c)$ grid. The important practical improvement is that the search follows only Pell-generated solutions instead of testing almost all integer pairs.

Footnotes and References

Problem page: https://projecteuler.net/problem=390
Heron's formula: Wikipedia — Heron's formula
Pell equation: Wikipedia — Pell's equation
Diophantine equation: Wikipedia — Diophantine equation

Problem 390 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <thread>
#include <vector>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;

struct Options {
    u64 limit = 10000000000ULL;
    unsigned int threads = std::max(1U, std::thread::hardware_concurrency());
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg,
                            const std::string& prefix,
                            u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(ch - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_uint_after_prefix(const std::string& arg,
                             const std::string& prefix,
                             unsigned int& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned int>::max())) {
        return false;
    }

    value = static_cast<unsigned int>(parsed);
    return true;
}

bool parse_arguments(const 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_u64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        if (parse_uint_after_prefix(arg, "--threads=", options.threads)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.threads == 0U) {
        options.threads = 1U;
    }

    if (options.limit > static_cast<u64>(std::numeric_limits<i64>::max())) {
        std::cerr << "--limit must be <= " << std::numeric_limits<i64>::max()
                  << " to keep intermediate signed states safe.\n";
        return false;
    }

    return true;
}

std::string to_string_u128(u128 value) {
    if (value == 0U) {
        return "0";
    }

    std::string digits;
    while (value > 0U) {
        const int digit = static_cast<int>(value % 10U);
        digits.push_back(static_cast<char>('0' + digit));
        value /= 10U;
    }

    std::reverse(digits.begin(), digits.end());
    return digits;
}

u64 isqrt_u128(const u128 n) {
    u128 lo = 0U;
    u128 hi = static_cast<u128>(std::numeric_limits<u64>::max());

    while (lo < hi) {
        const u128 mid = (lo + hi + 1U) >> 1;
        if (mid <= n / mid) {
            lo = mid;
        } else {
            hi = mid - 1U;
        }
    }

    return static_cast<u64>(lo);
}

u128 seed_first_area(const i64 p) {
    const u128 pp = static_cast<u128>(p);
    return 8U * pp * pp * pp + pp;
}

i64 max_seed(const u64 limit) {
    i64 lo = 0;
    i64 hi = 1;

    while (seed_first_area(hi) <= static_cast<u128>(limit)) {
        if (hi > std::numeric_limits<i64>::max() / 2) {
            break;
        }
        hi *= 2;
    }

    while (lo < hi) {
        const i64 mid = lo + (hi - lo + 1) / 2;
        if (seed_first_area(mid) <= static_cast<u128>(limit)) {
            lo = mid;
        } else {
            hi = mid - 1;
        }
    }

    return lo;
}

struct State {
    i64 p;
    i64 q;
    u64 area;
};

bool next_state(const State& current,
                const u64 limit,
                State& same_p,
                State& swap_branch,
                State& sign_branch) {
    const i128 p = static_cast<i128>(current.p);
    const i128 p2 = p * p;

    const i128 a = 8 * p2 + 1;
    const i128 b = 4 * p;
    const i128 c = 4 * p * (4 * p2 + 1);

    const i128 q1 = a * static_cast<i128>(current.q) + b * static_cast<i128>(current.area);
    const i128 s1 = c * static_cast<i128>(current.q) + a * static_cast<i128>(current.area);

    if (s1 <= 0 || static_cast<u128>(s1) > static_cast<u128>(limit)) {
        return false;
    }
    if (q1 <= 0 || q1 > static_cast<i128>(std::numeric_limits<i64>::max())) {
        return false;
    }

    const i64 nq = static_cast<i64>(q1);
    const u64 ns = static_cast<u64>(s1);

    same_p = State{current.p, nq, ns};
    swap_branch = State{nq, current.p, ns};
    sign_branch = State{nq, -current.p, ns};
    return true;
}

u128 sum_for_seed(const i64 seed, const u64 limit) {
    u128 sum = 0U;
    std::vector<State> stack;
    stack.reserve(256);
    stack.push_back(State{seed, 0, static_cast<u64>(seed)});

    while (!stack.empty()) {
        const State current = stack.back();
        stack.pop_back();

        State same_p{};
        State swap_branch{};
        State sign_branch{};
        if (!next_state(current, limit, same_p, swap_branch, sign_branch)) {
            continue;
        }

        sum += static_cast<u128>(same_p.area);

        stack.push_back(same_p);
        stack.push_back(swap_branch);
        stack.push_back(sign_branch);
    }

    return sum;
}

u128 solve(const u64 limit, const unsigned int requested_threads) {
    const i64 pmax = max_seed(limit);
    if (pmax <= 0) {
        return 0U;
    }

    const unsigned int workers = std::max(
        1U,
        std::min<unsigned int>(requested_threads, static_cast<unsigned int>(pmax)));

    std::atomic<i64> next_seed(1);
    std::vector<u128> partials(workers, 0U);
    std::vector<std::thread> threads;
    threads.reserve(workers);

    for (unsigned int t = 0; t < workers; ++t) {
        threads.emplace_back([&, t]() {
            u128 local = 0U;
            while (true) {
                const i64 seed = next_seed.fetch_add(1, std::memory_order_relaxed);
                if (seed > pmax) {
                    break;
                }
                local += sum_for_seed(seed, limit);
            }
            partials[t] = local;
        });
    }

    for (std::thread& thread : threads) {
        thread.join();
    }

    u128 total = 0U;
    for (const u128 partial : partials) {
        total += partial;
    }
    return total;
}

u128 brute_sum_unordered(const u64 limit) {
    if (limit == 0U) {
        return 0U;
    }

    u128 total = 0U;
    const u64 max_b = 2U * limit;

    for (u64 b = 1U; b <= max_b; ++b) {
        const u64 max_c = max_b / b + 1U;
        if (max_c < b) {
            continue;
        }

        for (u64 c = b; c <= max_c; ++c) {
            const u128 m2 =
                static_cast<u128>(b) * b * c * c + static_cast<u128>(b) * b + static_cast<u128>(c) * c;
            const u64 m = isqrt_u128(m2);
            if (static_cast<u128>(m) * m != m2) {
                continue;
            }
            if ((m & 1U) != 0U) {
                continue;
            }

            const u64 area = m / 2U;
            if (area <= limit) {
                total += area;
            }
        }
    }

    return total;
}

void run_checkpoints(const Options& options) {
    {
        constexpr u64 b = 2U;
        constexpr u64 c = 8U;
        const u128 m2 = static_cast<u128>(b) * b * c * c +
                        static_cast<u128>(b) * b +
                        static_cast<u128>(c) * c;
        const u64 m = isqrt_u128(m2);
        if (static_cast<u128>(m) * m != m2 || m != 18U) {
            throw std::runtime_error("Checkpoint failed: sample triangle should have m=18.");
        }
        if (m / 2U != 9U) {
            throw std::runtime_error("Checkpoint failed: sample triangle area should be 9.");
        }
    }

    {
        constexpr u64 small_limit = 1000U;
        const u128 brute = brute_sum_unordered(small_limit);
        const u128 fast = solve(small_limit, 1U);
        if (brute != fast) {
            throw std::runtime_error(
                "Checkpoint failed: recurrence and brute-force mismatch at N=1000.");
        }
    }

    {
        constexpr u64 statement_limit = 1000000U;
        constexpr u64 statement_value = 18018206U;
        const u128 sample = solve(statement_limit, options.threads);
        if (sample != static_cast<u128>(statement_value)) {
            throw std::runtime_error("Checkpoint failed: S(10^6) mismatch.");
        }
    }
}

}  // namespace

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

    try {
        if (options.run_checkpoints) {
            run_checkpoints(options);
        }

        const u128 answer = solve(options.limit, options.threads);
        std::cout << "S(" << options.limit << ") = " << to_string_u128(answer)
                  << '\n';
    } catch (const std::exception& ex) {
        std::cerr << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

def max_seed(limit):
    lo = 0
    hi = 1
    def seed_first_area(p):
        return 8 * p * p * p + p
        
    while seed_first_area(hi) <= limit:
        hi *= 2
        
    while lo < hi:
        mid = (lo + hi + 1) // 2
        if seed_first_area(mid) <= limit:
            lo = mid
        else:
            hi = mid - 1
    return lo

def sum_for_seed(seed, limit):
    total = 0
    stack = [(seed, 0, seed)]
    
    while stack:
        p, q, area = stack.pop()
        
        p2 = p * p
        a = 8 * p2 + 1
        b = 4 * p
        c = 4 * p * (4 * p2 + 1)
        
        q1 = a * q + b * area
        s1 = c * q + a * area
        
        if s1 > 0 and s1 <= limit:
            total += s1
            stack.append((p, q1, s1))
            stack.append((q1, p, s1))
            stack.append((q1, -p, s1))
            
    return total

def solve():
    limit = 10000000000
    pmax = max_seed(limit)
    if pmax <= 0: return "0"
    
    ans = 0
    for seed in range(1, pmax + 1):
        ans += sum_for_seed(seed, limit)
        
    return str(ans)

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

Java

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

public class Euler390 {
    static long seedFirstArea(long p) {
        return 8 * p * p * p + p;
    }

    static long maxSeed(long limit) {
        long lo = 0;
        long hi = 1;
        while (true) {
            if (hi > 2000000)
                break;
            if (seedFirstArea(hi) <= limit)
                hi *= 2;
            else
                break;
        }

        while (lo < hi) {
            long mid = lo + (hi - lo + 1) / 2;
            if (seedFirstArea(mid) <= limit) {
                lo = mid;
            } else {
                hi = mid - 1;
            }
        }
        return lo;
    }

    static class State {
        long p;
        long q;
        long area;

        State(long p, long q, long area) {
            this.p = p;
            this.q = q;
            this.area = area;
        }
    }

    static long sumForSeed(long seed, long limit) {
        long total = 0;
        List<State> stack = new ArrayList<>();
        stack.add(new State(seed, 0, seed));

        BigInteger limitBi = BigInteger.valueOf(limit);

        while (!stack.isEmpty()) {
            State current = stack.remove(stack.size() - 1);
            BigInteger p = BigInteger.valueOf(current.p);
            BigInteger q = BigInteger.valueOf(current.q);
            BigInteger area = BigInteger.valueOf(current.area);

            BigInteger p2 = p.multiply(p);
            BigInteger a = BigInteger.valueOf(8).multiply(p2).add(BigInteger.ONE);
            BigInteger b = BigInteger.valueOf(4).multiply(p);
            BigInteger c = BigInteger.valueOf(4).multiply(p).multiply(
                    BigInteger.valueOf(4).multiply(p2).add(BigInteger.ONE));

            BigInteger q1 = a.multiply(q).add(b.multiply(area));
            BigInteger s1 = c.multiply(q).add(a.multiply(area));

            if (s1.compareTo(BigInteger.ZERO) > 0 && s1.compareTo(limitBi) <= 0) {
                long s1Long = s1.longValue();
                total += s1Long;

                long q1Long = q1.longValue();
                long pLong = current.p;

                stack.add(new State(pLong, q1Long, s1Long));
                stack.add(new State(q1Long, pLong, s1Long));
                stack.add(new State(q1Long, -pLong, s1Long));
            }
        }
        return total;
    }

    static String solve() {
        long limit = 10000000000L;
        long pmax = maxSeed(limit);
        if (pmax <= 0)
            return "0";

        long ans = 0;
        for (long seed = 1; seed <= pmax; seed++) {
            ans += sumForSeed(seed, limit);
        }
        return Long.toString(ans);
    }

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