Problem 914: Triangles inside Circles

Project Euler Problem 914 Solution

EulerSolve provides an optimized solution for Project Euler Problem 914, Triangles inside Circles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a positive radius $R$, let $F(R)$ denote the largest possible inradius of a primitive integer right triangle that fits strictly inside a circle of radius $R$. If the legs are $a$ and $b$, the hypotenuse is $c$, and the inradius is $r$, the challenge is to evaluate $F(10^{18})$. The geometric condition is much simpler than it first appears. A right triangle has circumradius $c/2$, so fitting strictly inside the circle is equivalent to requiring $c<2R$. The solution therefore becomes an optimization over primitive Pythagorean triples, and the implementations succeed by combining an exact formula for $r$, a monotonicity argument for fixed parameters, and an exact continuous upper bound that collapses the search interval. Mathematical Approach From the circle constraint to Euclid's parameters Every primitive Pythagorean triple can be written uniquely as $$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,$$ with integers $m>n>0$, $\gcd(m,n)=1$, and $m\not\equiv n\pmod 2$....

Detailed mathematical approach

Problem Summary

For a positive radius $R$, let $F(R)$ denote the largest possible inradius of a primitive integer right triangle that fits strictly inside a circle of radius $R$. If the legs are $a$ and $b$, the hypotenuse is $c$, and the inradius is $r$, the challenge is to evaluate $F(10^{18})$.

The geometric condition is much simpler than it first appears. A right triangle has circumradius $c/2$, so fitting strictly inside the circle is equivalent to requiring $c<2R$. The solution therefore becomes an optimization over primitive Pythagorean triples, and the implementations succeed by combining an exact formula for $r$, a monotonicity argument for fixed parameters, and an exact continuous upper bound that collapses the search interval.

Mathematical Approach

From the circle constraint to Euclid's parameters

Every primitive Pythagorean triple can be written uniquely as

$$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,$$

with integers $m>n>0$, $\gcd(m,n)=1$, and $m\not\equiv n\pmod 2$. Since the triangle is right-angled, its smallest enclosing circle has radius $c/2$, so the condition $c<2R$ becomes

$$m^2+n^2<2R.$$

Because the inequality is strict, it is convenient to set

$$T=2R-1,$$

so that admissible pairs are exactly those with

$$m^2+n^2\le T.$$

For a right triangle the inradius is

$$r=\frac{a+b-c}{2},$$

and substituting Euclid's formulas gives the remarkably simple objective

$$r=\frac{(m^2-n^2)+2mn-(m^2+n^2)}{2}=n(m-n).$$

So the entire problem becomes

$$F(R)=\max\left\{n(m-n):m>n>0,\ \gcd(m,n)=1,\ m\not\equiv n\pmod 2,\ m^2+n^2\le T\right\}.$$

Why each $n$ has a single best partner

Fix $n$. Then the score $n(m-n)$ increases strictly with $m$, so only the largest admissible $m$ can be optimal for that $n$. Ignoring parity and coprimality for a moment, the geometric ceiling is

$$m_{\max}(n)=\left\lfloor \sqrt{T-n^2}\right\rfloor.$$

If $m_{\max}(n)\le n$, then no triangle exists for that $n$. Otherwise the primitive-triple conditions are enforced in the only direction that matters: start from $m_{\max}(n)$, decrease once if $m$ has the same parity as $n$, and then keep decreasing by $2$ until $\gcd(m,n)=1$. The first pair that survives is already the best one for that $n$, because every later candidate has a smaller $m$ and therefore a smaller inradius.

This also yields the natural outer range for the search. Since every feasible pair has $m>n$, we must have

$$2n^2<m^2+n^2\le T,$$

hence

$$1\le n\le \left\lfloor \sqrt{\frac{T}{2}}\right\rfloor=\left\lfloor \sqrt{\frac{2R-1}{2}}\right\rfloor.$$

The continuous envelope that makes pruning exact

Now relax the arithmetic conditions and allow $m$ to move on the real boundary $m=\sqrt{T-n^2}$. Then every integer candidate with the same $n$ satisfies

$$r\le U(n):=n\left(\sqrt{T-n^2}-n\right).$$

This function is a true upper bound, not a heuristic estimate. Differentiating $U(n)$ shows that its unique peak occurs when the boundary value of $m$ and the chosen $n$ satisfy

$$m=(1+\sqrt 2)\,n.$$

Equivalently, the maximizing real $n$ satisfies

$$n_0^2=\frac{T}{4+2\sqrt 2}=\frac{2R-1}{4+2\sqrt 2}.$$

So the best triangle must lie near one continuous optimum rather than somewhere arbitrary in the whole interval. Once a current record is known, every $n$ with $U(n)$ no larger than that record can be discarded immediately. Because $U(n)$ rises to a single peak and then falls, the surviving region is one interval around $n_0$, which can be found accurately by bisection on the left and right sides.

Worked Example: $R=100$

Here $T=199$, so admissibility means

$$m^2+n^2\le 199.$$

The continuous peak is near

$$n_0\approx \sqrt{\frac{199}{4+2\sqrt 2}}\approx 5.40,$$

which already tells us that the optimum should occur close to $n=5$. For $n=4$,

$$m_{\max}(4)=\left\lfloor \sqrt{199-16}\right\rfloor=13.$$

The pair $(13,4)$ already has opposite parity and $\gcd(13,4)=1$, so it produces

$$a=13^2-4^2=153,\qquad b=2\cdot 13\cdot 4=104,\qquad c=13^2+4^2=185,$$

and therefore

$$r=4(13-4)=36.$$

Since $c/2=92.5<100$, the triangle really does fit inside the circle. Nearby candidates are slightly worse: $n=5$ leads to $m=12$ and $r=35$, while $n=6$ leads to $m=11$ and $r=30$. Thus $F(100)=36$.

How the Code Works

Shared search strategy

The C++, Python, and Java implementations all follow the same logic. They compute exact integer square roots so that the strict boundary $m^2+n^2<2R$ is handled as $m^2+n^2\le 2R-1$ without any floating-point off-by-one errors. They then estimate the location of the continuous peak with $n_0\approx\sqrt{(2R-1)/(4+2\sqrt 2)}$, scan a wide window around that point to obtain a strong initial record, and evaluate each tested $n$ by starting from $\lfloor\sqrt{2R-1-n^2}\rfloor$, fixing parity if necessary, and descending by $2$ until a coprime partner is found.

Exact pruning and final sweep

After the initial window, the implementations use the upper envelope $U(n)=n(\sqrt{2R-1-n^2}-n)$ to determine which $n$-values can still beat the current best answer. Because $U$ is unimodal, bisection on each side of the peak isolates the real interval where improvement is still possible. That interval is then converted back to integers and widened slightly to stay safe against rounding.

Inside the final sweep, the code applies the sharper integer bound

$$n\left(\left\lfloor\sqrt{2R-1-n^2}\right\rfloor-n\right)$$

before spending time on gcd tests. If even this value cannot improve the current record, that $n$ is skipped immediately. The pruning is therefore exact: every discarded case has been eliminated by a proved upper bound.

Complexity Analysis

The natural search variable is $n$, and its full admissible range has size $O(\sqrt R)$. For each tested $n$, the implementation performs one integer square root, possibly one parity correction, and then a downward walk through values of the correct parity until it reaches a coprime partner. So the basic structure is an $O(\sqrt R)$-sized search with very small per-candidate work in practice.

The continuous bound $U(n)$ makes the real running time much smaller than a full scan, because once a good record is found only a narrow band around the peak can remain relevant. Memory usage is $O(1)$, since the algorithm stores only a small number of integer and floating-point variables.

Footnotes and References

Problem page: https://projecteuler.net/problem=914
Pythagorean triple: Wikipedia - Pythagorean triple
Incircle and excircles of a triangle: Wikipedia - Incircle and excircles of a triangle
Thales's theorem: Wikipedia - Thales's theorem
Greatest common divisor: Wikipedia - Greatest common divisor
Integer square root: Wikipedia - Integer square root

Problem 914 source code

C++

#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>

namespace {

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

u64 isqrt_u128(u128 x) {
    long double y = std::sqrt(static_cast<long double>(x));
    u64 r = static_cast<u64>(y);
    while ((u128)(r + 1) * (r + 1) <= x) {
        ++r;
    }
    while ((u128)r * r > x) {
        --r;
    }
    return r;
}

u64 best_for_n(u64 R, u64 n) {
    const u128 C = (u128)2 * R;
    const u128 nn = (u128)n * n;
    if (nn >= C) {
        return 0;
    }

    const u128 t = C - 1 - nn;
    u64 m = isqrt_u128(t);
    if (m <= n) {
        return 0;
    }

    if (((m ^ n) & 1ULL) == 0ULL) {
        --m;
    }

    while (m > n && std::gcd(m, n) != 1) {
        m -= 2;
    }

    if (m <= n) {
        return 0;
    }

    return n * (m - n);
}

u64 exhaustive_F(u64 R) {
    const u128 C = (u128)2 * R;
    const u64 n_max = isqrt_u128((C - 1) / 2);

    u64 best = 0;
    for (u64 n = 1; n <= n_max; ++n) {
        const u64 cand = best_for_n(R, n);
        if (cand > best) {
            best = cand;
        }
    }
    return best;
}

u64 fast_F(u64 R) {
    const u128 C = (u128)2 * R;
    const long double sqrt2 = std::sqrt(2.0L);
    const long double denom = 4.0L + 2.0L * sqrt2;
    const long double Cld = static_cast<long double>(C - 1);

    u64 n0 = static_cast<u64>(std::llround(std::sqrt(Cld / denom)));
    if (n0 == 0) {
        n0 = 1;
    }

    u64 best = 0;
    constexpr int kWindow = 4096;
    for (int d = -kWindow; d <= kWindow; ++d) {
        const std::int64_t ni = static_cast<std::int64_t>(n0) + d;
        if (ni <= 0) {
            continue;
        }
        const u64 cand = best_for_n(R, static_cast<u64>(ni));
        if (cand > best) {
            best = cand;
        }
    }

    const auto ub_real = [&](long double n) {
        if (n <= 0) {
            return 0.0L;
        }
        const long double t = Cld - n * n;
        if (t <= 0) {
            return 0.0L;
        }
        return n * (std::sqrt(t) - n);
    };

    const long double n0ld = static_cast<long double>(n0);
    const long double nmaxld = std::sqrt(Cld / 2.0L);

    long double left = n0ld;
    long double right = n0ld;

    if (ub_real(n0ld) > static_cast<long double>(best)) {
        long double lo = 0.0L;
        long double hi = n0ld;
        for (int it = 0; it < 140; ++it) {
            const long double mid = (lo + hi) * 0.5L;
            if (ub_real(mid) > static_cast<long double>(best)) {
                hi = mid;
            } else {
                lo = mid;
            }
        }
        left = hi;

        lo = n0ld;
        hi = nmaxld;
        for (int it = 0; it < 140; ++it) {
            const long double mid = (lo + hi) * 0.5L;
            if (ub_real(mid) > static_cast<long double>(best)) {
                lo = mid;
            } else {
                hi = mid;
            }
        }
        right = lo;
    }

    u64 L = 1;
    if (left > 16.0L) {
        L = static_cast<u64>(std::floor(left)) - 16;
    }

    const u64 n_max = isqrt_u128((C - 1) / 2);
    u64 Rn = static_cast<u64>(std::ceil(right)) + 16;
    if (Rn > n_max) {
        Rn = n_max;
    }

    for (u64 n = L; n <= Rn; ++n) {
        const u128 nn = (u128)n * n;
        if (nn >= C) {
            break;
        }
        const u64 mmax = isqrt_u128(C - 1 - nn);
        if (mmax <= n) {
            continue;
        }
        const u64 ub = n * (mmax - n);
        if (ub <= best) {
            continue;
        }

        const u64 cand = best_for_n(R, n);
        if (cand > best) {
            best = cand;
        }
    }

    return best;
}

void validate() {
    assert(exhaustive_F(100) == 36);
    assert(fast_F(100) == 36);

    for (u64 R = 2; R <= 5000; ++R) {
        assert(fast_F(R) == exhaustive_F(R));
    }

    u64 seed = 1;
    for (int i = 0; i < 200; ++i) {
        seed = seed * 6364136223846793005ULL + 1ULL;
        const u64 R = 2 + (seed % 100'000'000ULL);
        assert(fast_F(R) == exhaustive_F(R));
    }
}

}  // namespace

int main() {
    validate();
    std::cout << fast_F(1'000'000'000'000'000'000ULL) << '\n';
    return 0;
}

Python

import math

def isqrt(x):
    return math.isqrt(x)

def best_for_n(R, n):
    C = 2 * R
    nn = n * n
    if nn >= C:
        return 0
        
    t = C - 1 - nn
    m = isqrt(t)
    if m <= n:
        return 0
        
    if ((m ^ n) & 1) == 0:
        m -= 1
        
    while m > n and math.gcd(m, n) != 1:
        m -= 2
        
    if m <= n:
        return 0
        
    return n * (m - n)

def fast_F(R):
    C = 2 * R
    sqrt2 = math.sqrt(2.0)
    denom = 4.0 + 2.0 * sqrt2
    Cld = float(C - 1)
    
    n0 = round(math.sqrt(Cld / denom))
    if n0 == 0:
        n0 = 1
        
    best = 0
    kWindow = 4096
    for d in range(-kWindow, kWindow + 1):
        ni = n0 + d
        if ni <= 0:
            continue
        cand = best_for_n(R, ni)
        if cand > best:
            best = cand
            
    def ub_real(n):
        if n <= 0:
            return 0.0
        t = Cld - n * n
        if t <= 0:
            return 0.0
        return n * (math.sqrt(t) - n)
        
    n0ld = float(n0)
    nmaxld = math.sqrt(Cld / 2.0)
    
    left = n0ld
    right = n0ld
    
    if ub_real(n0ld) > float(best):
        lo = 0.0
        hi = n0ld
        for _ in range(140):
            mid = (lo + hi) * 0.5
            if ub_real(mid) > float(best):
                hi = mid
            else:
                lo = mid
        left = hi
        
        lo = n0ld
        hi = nmaxld
        for _ in range(140):
            mid = (lo + hi) * 0.5
            if ub_real(mid) > float(best):
                lo = mid
            else:
                hi = mid
        right = lo
        
    L = 1
    if left > 16.0:
        L = int(math.floor(left)) - 16
        
    n_max = isqrt((C - 1) // 2)
    Rn = int(math.ceil(right)) + 16
    if Rn > n_max:
        Rn = n_max
        
    for n in range(L, Rn + 1):
        nn = n * n
        if nn >= C:
            break
        mmax = isqrt(C - 1 - nn)
        if mmax <= n:
            continue
        ub = n * (mmax - n)
        if ub <= best:
            continue
            
        cand = best_for_n(R, n)
        if cand > best:
            best = cand
            
    return best

def solve():
    return str(fast_F(1000000000000000000))

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

Java

public class Euler914 {

    static long isqrt(long x) {
        long r = (long) Math.sqrt((double) x);
        while ((r + 1) * (r + 1) <= x) {
            ++r;
        }
        while (r * r > x) {
            --r;
        }
        return r;
    }

    static long gcd(long a, long b) {
        return b == 0 ? a : gcd(b, a % b);
    }

    static long bestForN(long R, long n) {
        long C = 2 * R;
        long nn = n * n;
        if (nn >= C) {
            return 0;
        }

        long t = C - 1 - nn;
        long m = isqrt(t);
        if (m <= n) {
            return 0;
        }

        if (((m ^ n) & 1L) == 0L) {
            --m;
        }

        while (m > n && gcd(m, n) != 1) {
            m -= 2;
        }

        if (m <= n) {
            return 0;
        }

        return n * (m - n);
    }

    static long fastF(long R) {
        long C = 2 * R;
        double sqrt2 = Math.sqrt(2.0);
        double denom = 4.0 + 2.0 * sqrt2;
        double Cld = (double) (C - 1);

        long n0 = Math.round(Math.sqrt(Cld / denom));
        if (n0 == 0) {
            n0 = 1;
        }

        long best = 0;
        int kWindow = 4096;
        for (int d = -kWindow; d <= kWindow; ++d) {
            long ni = n0 + d;
            if (ni <= 0) {
                continue;
            }
            long cand = bestForN(R, ni);
            if (cand > best) {
                best = cand;
            }
        }

        double n0ld = (double) n0;
        double nmaxld = Math.sqrt(Cld / 2.0);

        double left = n0ld;
        double right = n0ld;

        double ubN0 = n0ld <= 0 || Cld - n0ld * n0ld <= 0 ? 0.0 : n0ld * (Math.sqrt(Cld - n0ld * n0ld) - n0ld);

        if (ubN0 > (double) best) {
            double lo = 0.0;
            double hi = n0ld;
            for (int it = 0; it < 140; ++it) {
                double mid = (lo + hi) * 0.5;
                double ubMid = mid <= 0 || Cld - mid * mid <= 0 ? 0.0 : mid * (Math.sqrt(Cld - mid * mid) - mid);
                if (ubMid > (double) best) {
                    hi = mid;
                } else {
                    lo = mid;
                }
            }
            left = hi;

            lo = n0ld;
            hi = nmaxld;
            for (int it = 0; it < 140; ++it) {
                double mid = (lo + hi) * 0.5;
                double ubMid = mid <= 0 || Cld - mid * mid <= 0 ? 0.0 : mid * (Math.sqrt(Cld - mid * mid) - mid);
                if (ubMid > (double) best) {
                    lo = mid;
                } else {
                    hi = mid;
                }
            }
            right = lo;
        }

        long L = 1;
        if (left > 16.0) {
            L = (long) Math.floor(left) - 16;
        }

        long nMax = isqrt((C - 1) / 2);
        long Rn = (long) Math.ceil(right) + 16;
        if (Rn > nMax) {
            Rn = nMax;
        }

        for (long n = L; n <= Rn; ++n) {
            long nn = n * n;
            if (nn >= C) {
                break;
            }
            long mmax = isqrt(C - 1 - nn);
            if (mmax <= n) {
                continue;
            }
            long ub = n * (mmax - n);
            if (ub <= best) {
                continue;
            }

            long cand = bestForN(R, n);
            if (cand > best) {
                best = cand;
            }
        }

        return best;
    }

    public static String solve() {
        return Long.toString(fastF(1000000000000000000L));
    }

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