Problem 482: The Incenter of a Triangle
View on Project EulerProject Euler Problem 482 Solution
EulerSolve provides an optimized solution for Project Euler Problem 482, The Incenter of a Triangle, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary We consider integer-sided triangles \(ABC\) whose perimeter \(p\) does not exceed a bound \(P\). Let \(I\) be the incenter, and define $$L=p+IA+IB+IC.$$ The task is to sum \(L\) over exactly those triangles for which the three distances from the incenter to the vertices are integers. A brute-force scan over all side triples is far too expensive for the real limit, so the solution rewrites the triangle in coordinates that expose the arithmetic structure of the incenter directly. Mathematical Approach The key observation is that the incenter conditions become much simpler after switching from side lengths to semiperimeter coordinates. Step 1: Replace the Sides by Semiperimeter Offsets Let $$s=\frac{a+b+c}{2},\qquad x=s-a,\qquad y=s-b,\qquad z=s-c.$$ Then \(x,y,z>0\), and the original sides can be recovered from $$a=y+z,\qquad b=x+z,\qquad c=x+y.$$ Because \(s=x+y+z\), the perimeter is $$p=2s=2(x+y+z).$$ Using Heron's formula \(\Delta^2=s(s-a)(s-b)(s-c)\) together with \(\Delta=rs\), we obtain $$r^2=\frac{\Delta^2}{s^2}=\frac{xyz}{x+y+z}.$$ So every admissible triangle corresponds to positive integers \(x,y,z,r\) satisfying that identity....
Detailed mathematical approach
Problem Summary
We consider integer-sided triangles \(ABC\) whose perimeter \(p\) does not exceed a bound \(P\). Let \(I\) be the incenter, and define
$$L=p+IA+IB+IC.$$
The task is to sum \(L\) over exactly those triangles for which the three distances from the incenter to the vertices are integers. A brute-force scan over all side triples is far too expensive for the real limit, so the solution rewrites the triangle in coordinates that expose the arithmetic structure of the incenter directly.
Mathematical Approach
The key observation is that the incenter conditions become much simpler after switching from side lengths to semiperimeter coordinates.
Step 1: Replace the Sides by Semiperimeter Offsets
Let
$$s=\frac{a+b+c}{2},\qquad x=s-a,\qquad y=s-b,\qquad z=s-c.$$
Then \(x,y,z>0\), and the original sides can be recovered from
$$a=y+z,\qquad b=x+z,\qquad c=x+y.$$
Because \(s=x+y+z\), the perimeter is
$$p=2s=2(x+y+z).$$
Using Heron's formula \(\Delta^2=s(s-a)(s-b)(s-c)\) together with \(\Delta=rs\), we obtain
$$r^2=\frac{\Delta^2}{s^2}=\frac{xyz}{x+y+z}.$$
So every admissible triangle corresponds to positive integers \(x,y,z,r\) satisfying that identity.
Step 2: Encode the Integer Vertex Distances
The tangency point on side \(BC\) forms a right triangle with legs \(r\) and \(x\), so
$$IA^2=r^2+x^2,\qquad IB^2=r^2+y^2,\qquad IC^2=r^2+z^2.$$
Therefore each of \(x,y,z\) must participate in a Pythagorean relation with the same inradius \(r\). For a fixed \(r\), write
$$t^2-x^2=r^2,$$
or equivalently
$$\left(t-x\right)\left(t+x\right)=r^2.$$
If we set \(d=t-x\) and \(e=t+x\), then \(de=r^2\), \(d\le e\), and \(d\) and \(e\) must have the same parity. This gives
$$x=\frac{e-d}{2},\qquad t=\frac{e+d}{2},\qquad e=\frac{r^2}{d}.$$
So enumerating divisors \(d\mid r^2\) generates all possible integer pairs \((x,t)\) with \(x^2+r^2=t^2\). The same list serves for \(y\) and \(z\) as well.
Step 3: Recover the Third Coordinate from Two Choices
Once two entries \(x\) and \(y\) are chosen from that list, the inradius formula determines \(z\). Starting from
$$r^2=\frac{xyz}{x+y+z},$$
we solve for \(z\):
$$z=\frac{r^2(x+y)}{xy-r^2}.$$
This immediately explains the filters used by the implementation: the denominator must be positive, so
$$xy>r^2,$$
and the numerator must be divisible by that denominator so that \(z\) is an integer. Finally, \(z\) must itself belong to the previously generated Pythagorean list; otherwise \(IC\) would not be integral.
Step 4: Avoid Duplicates and Enforce the Perimeter Bound
The variables \(x,y,z\) are symmetric, so we may impose
$$x\le y\le z$$
to count each triangle once. The perimeter condition becomes
$$2(x+y+z)\le P,$$
or equivalently
$$z\le \frac{P}{2}-x-y.$$
For each valid triple the contribution is
$$L=2(x+y+z)+\sqrt{x^2+r^2}+\sqrt{y^2+r^2}+\sqrt{z^2+r^2}.$$
There is also a global bound on \(r\). Among triangles with fixed perimeter, the equilateral triangle maximizes the inradius, so
$$r\le \frac{P\sqrt{3}}{18}.$$
That turns the search into a finite loop over
$$1\le r\le \left\lfloor\frac{P\sqrt{3}}{18}\right\rfloor.$$
Step 5: Worked Example
Take \(r=21\), so \(r^2=441\). Factor pairs of \(441\) with the same parity give several Pythagorean options. Two useful ones are
$$d=7,\ e=63\quad\Rightarrow\quad x=\frac{63-7}{2}=28,\qquad \frac{63+7}{2}=35,$$
and
$$d=3,\ e=147\quad\Rightarrow\quad z=\frac{147-3}{2}=72,\qquad \frac{147+3}{2}=75.$$
Now choose \(x=y=28\). Then
$$z=\frac{441(28+28)}{28\cdot 28-441}=\frac{441\cdot 56}{343}=72.$$
Because \(72\) is already in the same Pythagorean list, all three vertex distances are integral:
$$\sqrt{21^2+28^2}=35,\qquad \sqrt{21^2+72^2}=75.$$
The corresponding sides are
$$a=y+z=100,\qquad b=x+z=100,\qquad c=x+y=56,$$
so the perimeter is \(256\) and the contribution is
$$L=256+35+35+75=401.$$
How the Code Works
The C++, Python, and Java implementations all follow the same pipeline. First they compute the upper limit \(\left\lfloor P\sqrt{3}/18\right\rfloor\) for the inradius and build a smallest-prime-factor sieve up to that value. The sieve makes it cheap to factor each candidate \(r\).
For a fixed inradius, the implementation doubles the exponents in the factorization of \(r\) to enumerate every divisor of \(r^2\). Each divisor pair produces one candidate \((x,t)\) with \(x^2+r^2=t^2\), and these candidates are stored in sorted order by \(x\).
It then loops over ordered pairs \(x\le y\), computes the forced value of \(z\), and rejects the pair unless the divisibility condition, the positivity condition, and the perimeter bound all hold. A binary search checks whether that \(z\) already appears in the candidate list, which is equivalent to checking that \(z^2+r^2\) is a square. Every surviving triple contributes the perimeter plus the three already-determined vertex distances to the running total.
The larger implementations can split the inradius interval into chunks and sum those chunks independently. They also use small direct checks and the known checkpoint at \(P=1000\), namely \(3619\), before evaluating the full bound.
Complexity Analysis
Let
$$R=\left\lfloor\frac{P\sqrt{3}}{18}\right\rfloor.$$
Building the smallest-prime-factor sieve costs \(O(R\log\log R)\) time and \(O(R)\) memory. For a fixed \(r\), suppose \(m(r)\) candidate values survive the divisor-to-leg conversion. Factoring \(r\) and generating divisors is essentially proportional to the number of divisors of \(r^2\), sorting costs \(O(m(r)\log m(r))\), and the pair search costs \(O(m(r)^2\log m(r))\) because each pair uses one binary search for \(z\).
Thus the total running time is
$$O\left(R\log\log R+\sum_{r=1}^{R} m(r)^2\log m(r)\right),$$
with \(O(R+m_{\max})\) memory. In practice \(m(r)\) is small, so this arithmetic reformulation is dramatically faster than scanning all integer side triples.
Footnotes and References
- Problem page: https://projecteuler.net/problem=482
- Heron's formula: Wikipedia - Heron's formula
- Incenter and incircle formulas: Wikipedia - Incircle and excircles of a triangle
- Pythagorean triples and difference of squares: Wikipedia - Pythagorean triple
Problem 482 source code
C++
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>
namespace {
using int64 = long long;
using u128 = unsigned __int128;
struct Leg {
int64 x;
int64 t;
};
int max_inradius_for_perimeter(int64 P) {
long double bound = static_cast<long double>(P) * std::sqrt(3.0L) / 18.0L;
int r_max = static_cast<int>(std::floor(bound + 1e-12L));
return std::max(0, r_max);
}
std::vector<int> build_spf(int limit) {
std::vector<int> spf(limit + 1, 0);
if (limit >= 1) spf[1] = 1;
for (int i = 2; i <= limit; ++i) {
if (spf[i] == 0) {
spf[i] = i;
if ((int64)i * i <= limit) {
for (int64 j = (int64)i * i; j <= limit; j += i) {
if (spf[(int)j] == 0) spf[(int)j] = i;
}
}
}
}
return spf;
}
void factorize(int n, const std::vector<int>& spf, std::vector<std::pair<int, int>>& out) {
out.clear();
while (n > 1) {
int p = spf[n];
int cnt = 0;
while (n % p == 0) {
n /= p;
++cnt;
}
out.push_back({p, cnt});
}
}
std::vector<Leg> build_legs(int r, int64 r2, int64 max_x, const std::vector<int>& spf) {
std::vector<std::pair<int, int>> factors;
factorize(r, spf, factors);
std::vector<int64> divisors;
divisors.reserve(64);
divisors.push_back(1);
for (const auto& f : factors) {
int p = f.first;
int exp = 2 * f.second;
int64 p_pow = 1;
const std::size_t base = divisors.size();
for (int e = 1; e <= exp; ++e) {
p_pow *= p;
for (std::size_t i = 0; i < base; ++i) {
divisors.push_back(divisors[i] * p_pow);
}
}
}
std::vector<Leg> legs;
legs.reserve(divisors.size());
for (int64 d : divisors) {
if (d > r) continue;
int64 d2 = r2 / d;
if (((d + d2) & 1LL) != 0) continue;
int64 x = (d2 - d) / 2;
if (x <= 0 || x > max_x) continue;
int64 t = (d2 + d) / 2;
legs.push_back({x, t});
}
std::sort(legs.begin(), legs.end(), [](const Leg& a, const Leg& b) { return a.x < b.x; });
return legs;
}
u128 compute_sum_range(int64 P, int r_start, int r_end, const std::vector<int>& spf) {
const int64 P_half = P / 2;
u128 total = 0;
for (int r = r_start; r <= r_end; ++r) {
int64 r2 = (int64)r * r;
auto legs = build_legs(r, r2, P_half, spf);
if (legs.size() < 2) continue;
std::vector<int64> xs;
std::vector<int64> ts;
xs.reserve(legs.size());
ts.reserve(legs.size());
for (const auto& leg : legs) {
xs.push_back(leg.x);
ts.push_back(leg.t);
}
const std::size_t n = xs.size();
for (std::size_t i = 0; i < n; ++i) {
int64 x = xs[i];
int64 tx = ts[i];
for (std::size_t j = i; j < n; ++j) {
int64 y = xs[j];
int64 ty = ts[j];
int64 denom = x * y - r2;
if (denom <= 0) continue;
u128 numer = (u128)r2 * (x + y);
if (numer % denom != 0) continue;
u128 z_u = numer / denom;
int64 max_z = P_half - x - y;
if (max_z <= 0) continue;
if (z_u > (u128)max_z) continue;
int64 z = (int64)z_u;
if (z < y) continue;
auto it = std::lower_bound(xs.begin(), xs.end(), z);
if (it == xs.end() || *it != z) continue;
std::size_t k = static_cast<std::size_t>(it - xs.begin());
int64 tz = ts[k];
int64 p = 2 * (x + y + z);
int64 L = p + tx + ty + tz;
total += (u128)L;
}
}
}
return total;
}
u128 compute_sum(int64 P, const std::vector<int>& spf, int r_max, unsigned threads) {
if (r_max <= 0) return 0;
if (threads == 0) threads = 1;
if (threads == 1 || r_max < 200) {
return compute_sum_range(P, 1, r_max, spf);
}
int chunk = (r_max + (int)threads - 1) / (int)threads;
std::vector<std::thread> pool;
std::vector<u128> partial(threads, 0);
pool.reserve(threads);
for (unsigned t = 0; t < threads; ++t) {
int start = (int)t * chunk + 1;
int end = std::min(r_max, start + chunk - 1);
if (start > end) continue;
pool.emplace_back([&, start, end, t]() {
partial[t] = compute_sum_range(P, start, end, spf);
});
}
for (auto& th : pool) th.join();
u128 total = 0;
for (const auto& v : partial) total += v;
return total;
}
std::string to_string_u128(u128 value) {
if (value == 0) return "0";
std::string out;
while (value > 0) {
int digit = static_cast<int>(value % 10);
out.push_back(static_cast<char>('0' + digit));
value /= 10;
}
std::reverse(out.begin(), out.end());
return out;
}
int64 brute_sum(int P) {
int64 total = 0;
for (int a = 1; a <= P / 3; ++a) {
for (int b = a; b <= (P - a) / 2; ++b) {
int max_c = std::min(P - a - b, a + b - 1);
if (max_c < b) continue;
for (int c = b; c <= max_c; ++c) {
int p = a + b + c;
if (p > P) break;
double s = 0.5 * p;
double x = s - a;
double y = s - b;
double z = s - c;
double area2 = s * x * y * z;
if (area2 <= 0.0) continue;
double r = std::sqrt(area2) / s;
double IA = std::sqrt(r * r + x * x);
double IB = std::sqrt(r * r + y * y);
double IC = std::sqrt(r * r + z * z);
auto round_int = [](double v) { return std::llround(v); };
long long ia = round_int(IA);
long long ib = round_int(IB);
long long ic = round_int(IC);
if (std::fabs(IA - ia) > 1e-9) continue;
if (std::fabs(IB - ib) > 1e-9) continue;
if (std::fabs(IC - ic) > 1e-9) continue;
total += p + ia + ib + ic;
}
}
}
return total;
}
} // namespace
int main(int argc, char** argv) {
int64 P = 10000000LL;
if (argc >= 2) {
P = std::stoll(argv[1]);
}
unsigned threads = std::thread::hardware_concurrency();
if (argc >= 3) {
threads = static_cast<unsigned>(std::stoul(argv[2]));
}
int r_max = max_inradius_for_perimeter(P);
auto spf = build_spf(r_max);
if (P <= 200) {
u128 fast = compute_sum(P, spf, r_max, 1);
int64 slow = brute_sum(static_cast<int>(P));
if (fast != static_cast<u128>(slow)) {
std::cerr << "Validation failed for P=" << P << ": fast="
<< to_string_u128(fast) << ", brute=" << slow << "\n";
return 1;
}
std::cout << to_string_u128(fast) << "\n";
return 0;
}
if (P >= 1000) {
int r_max_check = max_inradius_for_perimeter(1000);
u128 check = compute_sum(1000, spf, r_max_check, 1);
if (check != 3619) {
std::cerr << "Validation failed for P=1000: got "
<< to_string_u128(check) << ", expected 3619\n";
return 1;
}
}
u128 answer = compute_sum(P, spf, r_max, threads);
std::cout << to_string_u128(answer) << "\n";
return 0;
}
Python
import math
import multiprocessing
import bisect
def max_inradius_for_perimeter(P):
bound = P * math.sqrt(3.0) / 18.0
r_max = int(math.floor(bound + 1e-12))
return max(0, r_max)
def build_spf(limit):
spf = [0] * (limit + 1)
if limit >= 1:
spf[1] = 1
for i in range(2, limit + 1):
if spf[i] == 0:
spf[i] = i
if i * i <= limit:
for j in range(i * i, limit + 1, i):
if spf[j] == 0:
spf[j] = i
return spf
def factorize(n, spf):
factors = []
while n > 1:
p = spf[n]
cnt = 0
while n % p == 0:
n //= p
cnt += 1
factors.append((p, cnt))
return factors
def build_legs(r, r2, max_x, spf):
factors = factorize(r, spf)
divisors = [1]
for p, cnt in factors:
exp = 2 * cnt
base_len = len(divisors)
p_pow = 1
for e in range(1, exp + 1):
p_pow *= p
for i in range(base_len):
divisors.append(divisors[i] * p_pow)
legs = []
for d in divisors:
if d > r:
continue
d2 = r2 // d
if (d + d2) % 2 != 0:
continue
x = (d2 - d) // 2
if x <= 0 or x > max_x:
continue
t = (d2 + d) // 2
legs.append((x, t))
legs.sort(key=lambda leg: leg[0])
return legs
def init_worker(shared_P, shared_spf):
global P_half, s_spf, P
P = shared_P
P_half = P // 2
s_spf = shared_spf
def worker_chunk(args):
start, end = args
total = 0
for r in range(start, end + 1):
r2 = r * r
legs = build_legs(r, r2, P_half, s_spf)
if len(legs) < 2:
continue
xs = [l[0] for l in legs]
ts = [l[1] for l in legs]
n_legs = len(legs)
for i in range(n_legs):
x = xs[i]
tx = ts[i]
for j in range(i, n_legs):
y = xs[j]
ty = ts[j]
denom = x * y - r2
if denom <= 0:
continue
numer = r2 * (x + y)
if numer % denom != 0:
continue
z_u = numer // denom
max_z = P_half - x - y
if max_z <= 0 or z_u > max_z:
continue
z = z_u
if z < y:
continue
idx = bisect.bisect_left(xs, z)
if idx == n_legs or xs[idx] != z:
continue
tz = ts[idx]
p = 2 * (x + y + z)
L = p + tx + ty + tz
total += L
return total
def compute_sum(P, spf, r_max, threads=None):
if r_max <= 0:
return 0
if threads is None:
threads = multiprocessing.cpu_count()
if threads == 1 or r_max < 200:
init_worker(P, spf)
return worker_chunk((1, r_max))
chunks = []
chunk_size = (r_max + threads - 1) // threads
for i in range(1, r_max + 1, chunk_size):
chunks.append((i, min(r_max, i + chunk_size - 1)))
with multiprocessing.Pool(processes=threads, initializer=init_worker, initargs=(P, spf)) as pool:
results = pool.map(worker_chunk, chunks)
return sum(results)
def solve():
P = 10000000
r_max = max_inradius_for_perimeter(P)
spf = build_spf(r_max)
ans = compute_sum(P, spf, r_max)
return str(ans)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.concurrent.Callable;
import java.util.concurrent.ExecutionException;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;
public class Euler482 {
static class Leg {
long x;
long t;
Leg(long x, long t) {
this.x = x;
this.t = t;
}
}
private static int maxInradiusForPerimeter(long P) {
double bound = P * Math.sqrt(3.0) / 18.0;
int rMax = (int) Math.floor(bound + 1e-12);
return Math.max(0, rMax);
}
private 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) {
spf[i] = i;
if ((long) i * i <= limit) {
for (long j = (long) i * i; j <= limit; j += i) {
if (spf[(int) j] == 0)
spf[(int) j] = i;
}
}
}
}
return spf;
}
private static List<int[]> factorize(int n, int[] spf) {
List<int[]> out = new ArrayList<>();
while (n > 1) {
int p = spf[n];
int cnt = 0;
while (n % p == 0) {
n /= p;
cnt++;
}
out.add(new int[] { p, cnt });
}
return out;
}
private static List<Leg> buildLegs(int r, long r2, long maxX, int[] spf) {
List<int[]> factors = factorize(r, spf);
long[] divisors = new long[50000];
int numDivisors = 1;
divisors[0] = 1;
for (int[] f : factors) {
int p = f[0];
int exp = 2 * f[1];
long pPow = 1;
int baseSize = numDivisors;
for (int e = 1; e <= exp; ++e) {
pPow *= p;
for (int i = 0; i < baseSize; ++i) {
divisors[numDivisors++] = divisors[i] * pPow;
}
}
}
List<Leg> legs = new ArrayList<>();
for (int i = 0; i < numDivisors; i++) {
long d = divisors[i];
if (d > r)
continue;
long d2 = r2 / d;
if (((d + d2) & 1) != 0)
continue;
long x = (d2 - d) / 2;
if (x <= 0 || x > maxX)
continue;
long t = (d2 + d) / 2;
legs.add(new Leg(x, t));
}
legs.sort((a, b) -> Long.compare(a.x, b.x));
return legs;
}
private static long computeSumRange(long P, int rStart, int rEnd, int[] spf) {
long pHalf = P / 2;
long total = 0; // Using long, P=10^7 max P is 10^7, L ~ 2*10^7, total ~ sum of L over maybe 1M
// elements, easily fits in 64-bit
for (int r = rStart; r <= rEnd; ++r) {
long r2 = (long) r * r;
List<Leg> legs = buildLegs(r, r2, pHalf, spf);
if (legs.size() < 2)
continue;
int n = legs.size();
long[] xs = new long[n];
long[] ts = new long[n];
for (int i = 0; i < n; i++) {
xs[i] = legs.get(i).x;
ts[i] = legs.get(i).t;
}
for (int i = 0; i < n; ++i) {
long x = xs[i];
long tx = ts[i];
for (int j = i; j < n; ++j) {
long y = xs[j];
long ty = ts[j];
long denom = x * y - r2;
if (denom <= 0)
continue;
long numer = r2 * (x + y);
if (numer % denom != 0)
continue;
long zU = numer / denom;
long maxZ = pHalf - x - y;
if (maxZ <= 0 || zU > maxZ)
continue;
long z = zU;
if (z < y)
continue;
int idx = Arrays.binarySearch(xs, z);
if (idx < 0)
continue;
long tz = ts[idx];
long p = 2 * (x + y + z);
long L = p + tx + ty + tz;
total += L;
}
}
}
return total;
}
private static long computeSum(long P, int[] spf, int rMax) throws InterruptedException, ExecutionException {
if (rMax <= 0)
return 0;
int threads = Runtime.getRuntime().availableProcessors();
if (threads <= 0)
threads = 1;
if (rMax < 200 || threads == 1) {
return computeSumRange(P, 1, rMax, spf);
}
int chunk = (rMax + threads - 1) / threads;
ExecutorService executor = Executors.newFixedThreadPool(threads);
List<Future<Long>> futures = new ArrayList<>();
for (int t = 0; t < threads; ++t) {
int start = t * chunk + 1;
int end = Math.min(rMax, start + chunk - 1);
if (start > end)
continue;
futures.add(executor.submit(() -> computeSumRange(P, start, end, spf)));
}
long total = 0;
for (Future<Long> f : futures) {
total += f.get();
}
executor.shutdown();
return total;
}
public static void main(String[] args) throws Exception {
long P = 10000000L;
int rMax = maxInradiusForPerimeter(P);
int[] spf = buildSpf(rMax);
long answer = computeSum(P, spf, rMax);
System.out.println(answer);
}
}