Problem 496: Incenter and Circumcenter of Triangle

Project Euler Problem 496 Solution

EulerSolve provides an optimized solution for Project Euler Problem 496, Incenter and Circumcenter of Triangle, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The geometry of the incenter and circumcenter can be reduced to a primitive arithmetic parameterization. After that reduction, the quantity computed by the implementations is $$F(L)=\sum_{\substack{1\le r<s<2r\\ \gcd(r,s)=1\\ t\ge 1,\ rst\le L}} rst.$$ Here $(r,s)$ describes a primitive configuration and $t$ is the scaling factor. A direct scan over all admissible triples $(r,s,t)$ is too slow for $L=10^9$, so the solution turns the sum into divisor-based range queries and batches equal floor values. Mathematical Approach The arithmetic form above is the key step. Once the geometry has been compressed into coprime integer parameters, the rest of the problem becomes a careful summation problem in elementary number theory. Step 1: Describe the admissible domain The reduced parameterization gives the conditions $$1\le r,\qquad r+1\le s\le 2r-1,\qquad \gcd(r,s)=1,\qquad rst\le L.$$ For fixed $r$ and $s$, the scaling factor $t$ can range from $1$ up to $\left\lfloor L/(rs)\right\rfloor$. Also, if $r>\sqrt{L}$, then even the smallest allowed value $s=r+1$ gives $rs>r^2>L$, so no triple can contribute....

Detailed mathematical approach

Problem Summary

The geometry of the incenter and circumcenter can be reduced to a primitive arithmetic parameterization. After that reduction, the quantity computed by the implementations is

$$F(L)=\sum_{\substack{1\le r<s<2r\\ \gcd(r,s)=1\\ t\ge 1,\ rst\le L}} rst.$$

Here $(r,s)$ describes a primitive configuration and $t$ is the scaling factor. A direct scan over all admissible triples $(r,s,t)$ is too slow for $L=10^9$, so the solution turns the sum into divisor-based range queries and batches equal floor values.

Mathematical Approach

The arithmetic form above is the key step. Once the geometry has been compressed into coprime integer parameters, the rest of the problem becomes a careful summation problem in elementary number theory.

Step 1: Describe the admissible domain

The reduced parameterization gives the conditions

$$1\le r,\qquad r+1\le s\le 2r-1,\qquad \gcd(r,s)=1,\qquad rst\le L.$$

For fixed $r$ and $s$, the scaling factor $t$ can range from $1$ up to $\left\lfloor L/(rs)\right\rfloor$. Also, if $r>\sqrt{L}$, then even the smallest allowed value $s=r+1$ gives $rs>r^2>L$, so no triple can contribute. Therefore the outer loop only needs

$$1\le r\le R=\left\lfloor\sqrt{L}\right\rfloor.$$

For each such $r$, the valid interval for $s$ is

$$r+1\le s\le U_r=\min\!\left(2r-1,\left\lfloor\frac{L}{r}\right\rfloor\right).$$

Step 2: Remove the scaling loop with triangular numbers

For one fixed coprime pair $(r,s)$, the total contribution of all admissible scales is

$$\sum_{t=1}^{\lfloor L/(rs)\rfloor} rst=rs\sum_{t=1}^{\lfloor L/(rs)\rfloor} t.$$

Introduce the triangular-number function

$$T(n)=\frac{n(n+1)}{2}.$$

Then the whole problem becomes

$$F(L)=\sum_{r=1}^{R} r\sum_{\substack{r+1\le s\le U_r\\ \gcd(r,s)=1}} s\,T\!\left(\left\lfloor\frac{L}{rs}\right\rfloor\right).$$

This is already a major simplification: the implementations never iterate over $t$ explicitly.

Step 3: Encode coprimality by Möbius inversion

For a fixed $r$, define the weighted coprime prefix sum

$$C_r(x)=\sum_{\substack{1\le s\le x\\ \gcd(r,s)=1}} s.$$

Using the standard identity

$$\mathbf{1}_{\gcd(r,s)=1}=\sum_{d\mid \gcd(r,s)} \mu(d),$$

we obtain

$$\begin{aligned} C_r(x) &=\sum_{1\le s\le x} s\sum_{d\mid \gcd(r,s)}\mu(d)\\ &=\sum_{d\mid r}\mu(d)\sum_{\substack{1\le s\le x\\ d\mid s}} s\\ &=\sum_{d\mid r}\mu(d)\,d\,T\!\left(\left\lfloor\frac{x}{d}\right\rfloor\right). \end{aligned}$$

Only squarefree divisors matter, because $\mu(d)=0$ for every non-squarefree $d$. This is why the implementations precompute, for each $r$, the squarefree divisors together with their Möbius signs.

Step 4: Turn interval sums into prefix differences

Once $C_r(x)$ is available, any weighted coprime range sum follows immediately:

$$\sum_{\substack{a\le s\le b\\ \gcd(r,s)=1}} s=C_r(b)-C_r(a-1).$$

So the remaining challenge is not coprimality anymore. It is the changing floor value inside

$$T\!\left(\left\lfloor\frac{L}{rs}\right\rfloor\right).$$

Step 5: Batch all equal floor values

Fix $r$ and suppose the current left endpoint in the $s$-interval is $a$. Set

$$m=\left\lfloor\frac{L}{ra}\right\rfloor.$$

As $s$ increases, the quantity $\left\lfloor L/(rs)\right\rfloor$ stays equal to $m$ on the whole block

$$a\le s\le b=\min\!\left(U_r,\left\lfloor\frac{L}{rm}\right\rfloor\right).$$

Therefore that entire block contributes

$$r\,T(m)\left(C_r(b)-C_r(a-1)\right).$$

Then the next block starts at $b+1$. This is the same floor-grouping principle used in hyperbola-style divisor summation: one maximal interval replaces many identical floor evaluations.

Worked Example: $L=15$

Here

$$R=\left\lfloor\sqrt{15}\right\rfloor=3.$$

For $r=1$, the interval $r+1\le s\le U_r$ is empty, so there is no contribution.

For $r=2$, the only admissible value is $s=3$, and $\gcd(2,3)=1$. Then

$$\left\lfloor\frac{15}{2\cdot 3}\right\rfloor=2,$$

so the contribution is

$$2\cdot 3\cdot T(2)=6\cdot 3=18.$$

For $r=3$, the admissible values are $s=4$ and $s=5$. Both are coprime to $3$, and both satisfy

$$\left\lfloor\frac{15}{3s}\right\rfloor=1.$$

Hence their contributions are

$$3\cdot 4\cdot T(1)=12,\qquad 3\cdot 5\cdot T(1)=15.$$

Adding everything gives

$$18+12+15=45,$$

which matches the checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. They first compute $R=\lfloor\sqrt{L}\rfloor$, because no larger $r$ can appear. Next they build a smallest-prime-factor sieve up to $R$, which lets them factor every $r$ quickly.

From those factorizations they generate, for each $r$, the list of squarefree divisors and the corresponding Möbius sign. A coprime prefix query $C_r(x)$ is then evaluated by summing $d\,T(\lfloor x/d\rfloor)$ over that divisor list with the correct sign.

The main loop runs over $r$. For each $r$, it scans the valid $s$-interval in maximal blocks on which $\left\lfloor L/(rs)\right\rfloor$ is constant. On each block, the implementation asks for one coprime interval sum, multiplies it by $r$ and the corresponding triangular number, and adds the result to the total.

The numeric types are chosen to avoid overflow: the Python version uses arbitrary-precision integers automatically, the Java version uses big integers for the accumulated total, and the C++ version uses 128-bit arithmetic for the same reason.

Complexity Analysis

Let $R=\lfloor\sqrt{L}\rfloor$. The smallest-prime-factor sieve costs $O(R\log\log R)$ time and $O(R)$ memory. Building the squarefree-divisor tables costs additional work proportional to the total number of stored squarefree divisors, namely

$$O\!\left(\sum_{r\le R} 2^{\omega(r)}\right),$$

where $\omega(r)$ is the number of distinct prime factors of $r$.

During the main summation, each floor-constant block for a given $r$ needs two coprime prefix evaluations, and each such evaluation iterates over the same squarefree-divisor list. If $B(r)$ denotes the number of blocks for that $r$, then the main-loop cost is

$$O\!\left(\sum_{r\le R} B(r)\,2^{\omega(r)}\right).$$

This is far smaller in practice than iterating over every admissible triple $(r,s,t)$, because the $t$-sum has been collapsed into a triangular number and the $s$-loop is processed in batches rather than one value at a time.

Footnotes and References

Problem page: https://projecteuler.net/problem=496
Möbius function: Wikipedia — Möbius function
Möbius inversion formula: Wikipedia — Möbius inversion formula
Triangular number: Wikipedia — Triangular number
Dirichlet hyperbola method: Wikipedia — Dirichlet hyperbola method

Problem 496 source code

C++

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

namespace {

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

struct Options {
    u64 L = 1'000'000'000ULL;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& out) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        out = static_cast<u64>(std::stoull(tail));
    } catch (...) {
        return false;
    }
    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_u64_after_prefix(arg, "--L=", options.L)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

u128 tri_u128(const u64 n) {
    return static_cast<u128>(n) * static_cast<u128>(n + 1ULL) / 2U;
}

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

std::vector<int> build_spf(const int n) {
    std::vector<int> spf(static_cast<std::size_t>(n + 1), 0);
    for (int i = 2; i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0) {
            spf[static_cast<std::size_t>(i)] = i;
            if (static_cast<int64_t>(i) * i <= n) {
                for (int x = i * i; x <= n; x += i) {
                    if (spf[static_cast<std::size_t>(x)] == 0) {
                        spf[static_cast<std::size_t>(x)] = i;
                    }
                }
            }
        }
    }
    spf[1] = 1;
    return spf;
}

std::vector<std::vector<std::pair<int, int>>> build_squarefree_mu_divisors(const int n,
                                                                            const std::vector<int>& spf) {
    std::vector<std::vector<std::pair<int, int>>> divs(static_cast<std::size_t>(n + 1));
    divs[1].push_back({1, 1});

    for (int x = 2; x <= n; ++x) {
        int t = x;
        std::vector<int> primes;
        while (t > 1) {
            const int p = spf[static_cast<std::size_t>(t)];
            primes.push_back(p);
            while (t % p == 0) {
                t /= p;
            }
        }

        std::vector<std::pair<int, int>> cur;
        cur.push_back({1, 1});
        for (const int p : primes) {
            const std::size_t before = cur.size();
            for (std::size_t i = 0; i < before; ++i) {
                cur.push_back({cur[i].first * p, -cur[i].second});
            }
        }
        divs[static_cast<std::size_t>(x)] = std::move(cur);
    }
    return divs;
}

i128 coprime_prefix_sum(const int r, const u64 x,
                        const std::vector<std::vector<std::pair<int, int>>>& sqf_divs) {
    if (x == 0ULL) {
        return 0;
    }
    i128 out = 0;
    for (const auto& [d, mu] : sqf_divs[static_cast<std::size_t>(r)]) {
        const u64 q = x / static_cast<u64>(d);
        const u128 tri = tri_u128(q);
        const i128 term = static_cast<i128>(static_cast<u128>(d) * tri);
        out += (mu > 0 ? term : -term);
    }
    return out;
}

u128 coprime_range_sum(const int r, const u64 l, const u64 rr,
                       const std::vector<std::vector<std::pair<int, int>>>& sqf_divs) {
    if (l > rr) {
        return 0;
    }
    const i128 hi = coprime_prefix_sum(r, rr, sqf_divs);
    const i128 lo = coprime_prefix_sum(r, l - 1ULL, sqf_divs);
    return static_cast<u128>(hi - lo);
}

u128 solve(const u64 L) {
    if (L < 2ULL) {
        return 0;
    }

    const int r_max = static_cast<int>(std::sqrt(static_cast<long double>(L)));
    const std::vector<int> spf = build_spf(r_max);
    const auto sqf_divs = build_squarefree_mu_divisors(r_max, spf);

    u128 total = 0;

    for (int r = 1; r <= r_max; ++r) {
        u64 s_lo = static_cast<u64>(r) + 1ULL;
        const u64 s_hi = std::min<u64>(2ULL * static_cast<u64>(r) - 1ULL, L / static_cast<u64>(r));
        if (s_lo > s_hi) {
            continue;
        }

        while (s_lo <= s_hi) {
            const u64 m = L / (static_cast<u64>(r) * s_lo);
            const u64 s_end =
                std::min<u64>(s_hi, L / (static_cast<u64>(r) * m));

            const u128 sum_s = coprime_range_sum(r, s_lo, s_end, sqf_divs);
            const u128 add = static_cast<u128>(r) * tri_u128(m) * sum_s;
            total += add;

            s_lo = s_end + 1ULL;
        }
    }

    return total;
}

bool run_checkpoints() {
    if (solve(15ULL) != static_cast<u128>(45ULL)) {
        std::cerr << "Checkpoint failed: F(15)\n";
        return false;
    }
    return true;
}

}  // namespace

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

    const u128 ans = solve(options.L);
    std::cout << to_string_u128(ans) << '\n';
    return 0;
}

Python

import math

def tri_u128(n):
    return n * (n + 1) // 2

def build_spf(n):
    spf = [0] * (n + 1)
    for i in range(2, n + 1):
        if spf[i] == 0:
            spf[i] = i
            if i * i <= n:
                for x in range(i * i, n + 1, i):
                    if spf[x] == 0:
                        spf[x] = i
    if n >= 1: spf[1] = 1
    return spf

def build_squarefree_mu_divisors(n, spf):
    divs = [[] for _ in range(n + 1)]
    if n >= 1: divs[1].append((1, 1))
    
    for x in range(2, n + 1):
        t = x
        primes = []
        while t > 1:
            p = spf[t]
            primes.append(p)
            while t % p == 0:
                t //= p
                
        cur = [(1, 1)]
        for p in primes:
            for i in range(len(cur)):
                cur.append((cur[i][0] * p, -cur[i][1]))
        divs[x] = cur
    return divs

def coprime_prefix_sum(x, sqf_divs_r):
    if x == 0: return 0
    out = 0
    for d, mu in sqf_divs_r:
        q = x // d
        term = d * tri_u128(q)
        if mu > 0:
            out += term
        else:
            out -= term
    return out

def coprime_range_sum(l, rr, sqf_divs_r):
    if l > rr: return 0
    hi = coprime_prefix_sum(rr, sqf_divs_r)
    lo = coprime_prefix_sum(l - 1, sqf_divs_r)
    return hi - lo

def solve():
    L = 1000000000
    if L < 2: return "0"
    
    r_max = int(math.sqrt(L))
    spf = build_spf(r_max)
    sqf_divs = build_squarefree_mu_divisors(r_max, spf)
    
    total = 0
    for r in range(1, r_max + 1):
        s_lo = r + 1
        s_hi = min(2 * r - 1, L // r)
        if s_lo > s_hi: continue
        
        while s_lo <= s_hi:
            m = L // (r * s_lo)
            s_end = min(s_hi, L // (r * m))
            
            sum_s = coprime_range_sum(s_lo, s_end, sqf_divs[r])
            add = r * tri_u128(m) * sum_s
            total += add
            
            s_lo = s_end + 1
            
    return str(total)

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

Java

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

public class Euler496 {

    static class Pair {
        long d;
        long mu;

        Pair(long d, long mu) {
            this.d = d;
            this.mu = mu;
        }
    }

    private static BigInteger triU128(long n) {
        BigInteger bn = BigInteger.valueOf(n);
        return bn.multiply(bn.add(BigInteger.ONE)).divide(BigInteger.valueOf(2));
    }

    private static int[] buildSpf(int n) {
        int[] spf = new int[n + 1];
        for (int i = 2; i <= n; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                if ((long) i * i <= n) {
                    for (int x = i * i; x <= n; x += i) {
                        if (spf[x] == 0)
                            spf[x] = i;
                    }
                }
            }
        }
        if (n >= 1)
            spf[1] = 1;
        return spf;
    }

    private static List<List<Pair>> buildSquarefreeMuDivisors(int n, int[] spf) {
        List<List<Pair>> divs = new ArrayList<>(n + 1);
        for (int i = 0; i <= n; i++)
            divs.add(new ArrayList<>());
        if (n >= 1)
            divs.get(1).add(new Pair(1, 1));

        for (int x = 2; x <= n; ++x) {
            int t = x;
            List<Integer> primes = new ArrayList<>();
            while (t > 1) {
                int p = spf[t];
                primes.add(p);
                while (t % p == 0) {
                    t /= p;
                }
            }

            List<Pair> cur = new ArrayList<>();
            cur.add(new Pair(1, 1));
            for (int p : primes) {
                int before = cur.size();
                for (int i = 0; i < before; ++i) {
                    cur.add(new Pair(cur.get(i).d * p, -cur.get(i).mu));
                }
            }
            divs.set(x, cur);
        }
        return divs;
    }

    private static BigInteger coprimePrefixSum(long x, List<Pair> sqfDivsR) {
        if (x == 0)
            return BigInteger.ZERO;
        BigInteger out = BigInteger.ZERO;
        for (Pair p : sqfDivsR) {
            long q = x / p.d;
            BigInteger term = BigInteger.valueOf(p.d).multiply(triU128(q));
            if (p.mu > 0) {
                out = out.add(term);
            } else {
                out = out.subtract(term);
            }
        }
        return out;
    }

    private static BigInteger coprimeRangeSum(long l, long rr, List<Pair> sqfDivsR) {
        if (l > rr)
            return BigInteger.ZERO;
        BigInteger hi = coprimePrefixSum(rr, sqfDivsR);
        BigInteger lo = coprimePrefixSum(l - 1, sqfDivsR);
        return hi.subtract(lo);
    }

    public static void main(String[] args) {
        long L = 1000000000L;
        if (L < 2) {
            System.out.println(0);
            return;
        }

        int rMax = (int) Math.sqrt(L);
        int[] spf = buildSpf(rMax);
        List<List<Pair>> sqfDivs = buildSquarefreeMuDivisors(rMax, spf);

        BigInteger total = BigInteger.ZERO;

        for (int r = 1; r <= rMax; ++r) {
            long sLo = r + 1L;
            long sHi = Math.min(2L * r - 1L, L / r);
            if (sLo > sHi)
                continue;

            while (sLo <= sHi) {
                long m = L / (r * sLo);
                long sEnd = Math.min(sHi, L / (r * m));

                BigInteger sumS = coprimeRangeSum(sLo, sEnd, sqfDivs.get(r));
                BigInteger add = BigInteger.valueOf(r).multiply(triU128(m)).multiply(sumS);
                total = total.add(add);

                sLo = sEnd + 1L;
            }
        }

        System.out.println(total.toString());
    }
}