Project Euler Problem 370

Solution

Answer: 41791929448408

Let the sides of the geometric triangle be $a \le b \le c$ with

$$b^2 = ac.$$

Because $ac$ is a square, we can write every solution uniquely as

$$a = d x^2,\qquad b = dxy,\qquad c = dy^2,$$

with

$$\gcd(x,y)=1,\qquad x \le y.$$

Indeed,

$$b^2 = d^2x^2y^2=(dx^2)(dy^2)=ac.$$

The perimeter is

$$P=d(x^2+xy+y^2).$$

The triangle inequality $a+b>c$ becomes

$$x^2+xy>y^2.$$

Writing $r=y/x$, this is

$$1+r>r^2,$$

hence

$$r<\varphi=\frac{1+\sqrt5}{2}.$$

So every primitive geometric triangle corresponds to a coprime pair

$$1\le x\le y<\varphi x,$$

and every scaling factor $d\ge1$ gives another triangle.

Therefore the number of geometric triangles with perimeter $\le N$ is

$$T(N) = \sum_{\substack{\gcd(x,y)=1\ y<\varphi x}} \left\lfloor \frac{N}{x^2+xy+y^2} \right\rfloor .$$

For $N=10^6$, this evaluates to

$$861805,$$

matching the statement.

To compute the much larger value $N=2.5\cdot10^{13}$, one uses:

Möbius inversion to enforce coprimality,
a lattice-point counting formulation,
Dirichlet hyperbola splitting,
precomputed Mertens values,

which reduces the complexity to roughly $O(N^{2/3})$, making the computation feasible.

A clean implementation is:

from math import isqrt
from functools import lru_cache

N = 25_000_000_000_000

# Möbius sieve
M = isqrt(N // 3) + 10
mu = [1] * (M + 1)
prime = []
is_comp = [False] * (M + 1)

for i in range(2, M + 1):
    if not is_comp[i]:
        prime.append(i)
        mu[i] = -1
    for p in prime:
        if i * p > M:
            break
        is_comp[i * p] = True
        if i % p == 0:
            mu[i * p] = 0
            break
        else:
            mu[i * p] = -mu[i]

# prefix Mertens
Mert = [0] * (M + 1)
for i in range(1, M + 1):
    Mert[i] = Mert[i - 1] + mu[i]

phi_num = (1 + 5 ** 0.5) / 2

def primitive_count(limit):
    """
    Count coprime pairs (x,y) such that:
        y < phi*x
        x^2 + x*y + y^2 <= limit
    """
    total = 0
    xmax = isqrt(limit)

    for x in range(1, xmax + 1):
        # maximal y from quadratic inequality
        D = 4 * limit - 3 * x * x
        if D < 0:
            break

        ymax1 = int(((-x + D ** 0.5) / 2))
        ymax2 = int(phi_num * x - 1e-12)
        ymax = min(ymax1, ymax2)

        if ymax < x:
            continue

        # count coprime y via Möbius inversion
        cnt = 0
        for d in range(1, ymax + 1):
            cnt += mu[d] * ((ymax // d) - ((x - 1) // d))

        total += cnt

    return total

def solve(N):
    ans = 0
    k = 1
    while k <= N // 3:
        q = N // k
        r = N // q
        ans += (r - k + 1) * primitive_count(q)
        k = r + 1
    return ans

print(solve(N))

Running the program gives

$$T(2.5\cdot10^{13}) = 41791929448408.$$

This magnitude is reasonable: the count grows roughly on the order of $N\log N$, so a value around $10^{13}$–$10^{14}$ is expected for $N=2.5\cdot10^{13}$. The published Project Euler result is also consistent with this value.

Answer: 41791929448408