Problem 223: Almost Right-angled Triangles I

Project Euler Problem 223 Solution

EulerSolve provides an optimized solution for Project Euler Problem 223, Almost Right-angled Triangles I, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We want to count integer-sided triangles $(a,b,c)$ with $a \le b \le c$, perimeter $a+b+c \le P$, and $$a^2+b^2=c^2+1.$$ For the Project Euler instance, $P=25{,}000{,}000$. The equation says that the triangle is exactly one unit away from being right-angled, so a brute-force search over all triples is completely impractical. The successful idea is to fix the smallest side $a$, rewrite the equation as a factorization of $a^2-1$, and count only the divisor pairs that can still satisfy the perimeter and ordering constraints. Mathematical Approach The central observation is that once $a$ is fixed, the other two sides can be recovered from a factor pair of $a^2-1$. That turns the problem from a quadratic search in $(b,c)$ into a bounded divisor-counting problem. Fix the smallest side first Because $a \le b \le c$, every valid triangle satisfies $$3a \le a+b+c \le P.$$ Therefore $$1 \le a \le \left\lfloor \frac{P}{3}\right\rfloor.$$ The implementations loop over this range of possible smallest sides and solve the rest of the triangle from there....

Detailed mathematical approach

Problem Summary

We want to count integer-sided triangles $(a,b,c)$ with $a \le b \le c$, perimeter $a+b+c \le P$, and

$$a^2+b^2=c^2+1.$$

For the Project Euler instance, $P=25{,}000{,}000$. The equation says that the triangle is exactly one unit away from being right-angled, so a brute-force search over all triples is completely impractical.

The successful idea is to fix the smallest side $a$, rewrite the equation as a factorization of $a^2-1$, and count only the divisor pairs that can still satisfy the perimeter and ordering constraints.

Mathematical Approach

The central observation is that once $a$ is fixed, the other two sides can be recovered from a factor pair of $a^2-1$. That turns the problem from a quadratic search in $(b,c)$ into a bounded divisor-counting problem.

Fix the smallest side first

Because $a \le b \le c$, every valid triangle satisfies

$$3a \le a+b+c \le P.$$

Therefore

$$1 \le a \le \left\lfloor \frac{P}{3}\right\rfloor.$$

The implementations loop over this range of possible smallest sides and solve the rest of the triangle from there.

Rewrite the equation as a product

Start from

$$a^2+b^2=c^2+1,$$

which can be rearranged as

$$c^2-b^2=a^2-1.$$

Now define

$$u=b+c,\qquad v=c-b.$$

Then $u>0$, $v\ge 0$, and

$$uv=(b+c)(c-b)=c^2-b^2=a^2-1.$$

Conversely, once a factor pair $(u,v)$ is known, the two larger sides are

$$b=\frac{u-v}{2},\qquad c=\frac{u+v}{2}.$$

So for a fixed $a$, every admissible triangle corresponds to a divisor pair of $a^2-1$ with the parity condition

$$u\equiv v \pmod 2,$$

because $b$ and $c$ must be integers.

For $a \ge 2$, we have $a^2-1>0$, so $v>0$. Then each valid choice of the smaller factor $v$ determines $u=(a^2-1)/v$, and therefore determines $(b,c)$ uniquely.

Turn the triangle conditions into bounds on one divisor

The perimeter becomes especially simple in the new variables:

$$a+b+c=a+u.$$

Hence the perimeter bound is

$$u \le P-a.$$

Since $u=(a^2-1)/v$, this gives the lower bound

$$v \ge v_{\min}=\left\lceil\frac{a^2-1}{P-a}\right\rceil.$$

Very small divisors make $u$ too large, so they cannot contribute.

The ordering condition $a \le b$ becomes

$$a \le \frac{u-v}{2}\qquad\Longleftrightarrow\qquad u \ge v+2a.$$

Substituting $u=(a^2-1)/v$ yields

$$v^2+2av \le a^2-1.$$

Solving this quadratic inequality gives an upper bound for the smaller factor:

$$v \le v_{\max}=\left\lfloor \sqrt{2a^2-1}-a\right\rfloor.$$

Therefore the search for $v$ is completely localized:

$$v_{\min} \le v \le v_{\max},\qquad v \mid (a^2-1).$$

Only divisors of $a^2-1$ inside this interval need to be tested.

Once this upper bound is enforced, the strict triangle inequality no longer needs a separate filter. Indeed,

$$a+b-c=a-v,$$

and the bound above implies $v<a$. So every divisor that survives the bound and parity checks already yields a genuine triangle with $a \le b \le c$.

The exceptional family $a=1$

When $a=1$, the product $a^2-1$ vanishes, so the divisor parameterization above is not the right tool. The original equation becomes

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

which immediately simplifies to $b=c$. So every triangle in this family has the form

$$ (1,b,b).$$

The perimeter condition is then

$$1+2b \le P,$$

so this case contributes exactly

$$\left\lfloor \frac{P-1}{2}\right\rfloor$$

solutions and is counted separately before the main divisor loop.

Worked example

Take $a=5$. Then

$$a^2-1=24.$$

The upper bound from $a \le b$ is

$$v_{\max}=\left\lfloor \sqrt{49}-5\right\rfloor=2.$$

So only the divisors $v=1$ and $v=2$ are even worth examining.

For $v=1$, we get $u=24$, but $u-v=23$ is odd, so $b$ and $c$ would not be integers. For $v=2$, we get $u=12$, hence

$$b=\frac{12-2}{2}=5,\qquad c=\frac{12+2}{2}=7.$$

This produces the valid almost-right triangle $(5,5,7)$, and indeed

$$5^2+5^2=7^2+1.$$

The full algorithm repeats exactly this bounded divisor test for every $a \le P/3$.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical plan. They first set $a_{\max}=\lfloor P/3\rfloor$, add the separate contribution from $a=1$, and build a smallest-prime-factor table up to $a_{\max}+1$. That range is enough because the implementations never factor $a^2-1$ directly. Instead they factor $a-1$ and $a+1$, then combine those prime exponents to recover the factorization of $(a-1)(a+1)=a^2-1$.

For each fixed $a \ge 2$, the implementation computes $v_{\min}$ from the perimeter bound and an initial floating-point estimate of $v_{\max}$ from $\sqrt{2a^2-1}-a$. Because floating-point rounding can miss the true floor by one, that estimate is corrected with exact integer inequalities before any counting starts. After that, the implementation generates only those divisors of $a^2-1$ that do not exceed $v_{\max}$, and each candidate $v$ is checked against the four remaining conditions: $v \ge v_{\min}$, correct parity, $u \ge v+2a$, and $a+u \le P$.

The work splits naturally across independent ranges of $a$. The C++ and Java implementations distribute consecutive chunks of $a$-values to worker threads, while the Python implementation uses the same chunking idea with worker processes when that helps and falls back to a serial pass otherwise. Each worker needs only small temporary factor and divisor buffers, so the main shared structure is the prime-factor table.

Complexity Analysis

The smallest-prime-factor table up to $a_{\max}+1=\lfloor P/3\rfloor+1$ costs $O(P)$ time and $O(P)$ memory up to constant factors. After that, the work for one fixed $a$ is dominated by factoring $a-1$ and $a+1$ and by generating the admissible divisors of $a^2-1$.

A convenient summary of the total runtime is

$$O\!\left(P+\sum_{a\le P/3}\tau(a^2-1)\right),$$

where $\tau(n)$ is the divisor-counting function. This matches the structure of the implementations: the sieve is linear-size preprocessing, and the main loop spends its time on the divisor structure of $a^2-1$. Parallel execution improves wall-clock time but does not change the asymptotic bound. Extra memory beyond the sieve is only small thread-local or process-local working storage.

Footnotes and References

Problem page: Project Euler 223 - Almost Right-angled Triangles I
Integer triangles: Wikipedia - Integer triangle
Diophantine equations: Wikipedia - Diophantine equation
Difference of two squares: Wikipedia - Difference of two squares
Fast factorization with a smallest-prime-factor sieve: cp-algorithms - Integer factorization

Problem 223 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <thread>
#include <vector>

namespace {

constexpr int kDefaultLimit = 25'000'000;

struct Factor {
    int p;
    int e;
};

std::vector<int> build_spf(int n) {
    std::vector<int> spf(n + 1, 0);
    std::vector<int> primes;
    primes.reserve(n / 10);

    for (int i = 2; i <= n; ++i) {
        if (spf[i] == 0) {
            spf[i] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            long long v = 1LL * p * i;
            if (v > n || p > spf[i]) break;
            spf[static_cast<std::size_t>(v)] = p;
        }
    }
    if (n >= 1) spf[1] = 1;
    return spf;
}

void factorize_int(int x, const std::vector<int>& spf, std::vector<Factor>& out) {
    out.clear();
    while (x > 1) {
        int p = spf[static_cast<std::size_t>(x)];
        int e = 0;
        do {
            x /= p;
            ++e;
        } while (x > 1 && spf[static_cast<std::size_t>(x)] == p);
        out.push_back({p, e});
    }
}

void merge_factors(const std::vector<Factor>& a,
                   const std::vector<Factor>& b,
                   std::vector<Factor>& out) {
    out.clear();
    out.reserve(a.size() + b.size());

    std::size_t i = 0;
    std::size_t j = 0;
    while (i < a.size() || j < b.size()) {
        if (j == b.size() || (i < a.size() && a[i].p < b[j].p)) {
            out.push_back(a[i++]);
        } else if (i == a.size() || b[j].p < a[i].p) {
            out.push_back(b[j++]);
        } else {
            out.push_back({a[i].p, a[i].e + b[j].e});
            ++i;
            ++j;
        }
    }
}

std::uint64_t count_for_a(int a,
                          int perimeter_limit,
                          const std::vector<int>& spf,
                          std::vector<Factor>& left_factors,
                          std::vector<Factor>& right_factors,
                          std::vector<Factor>& merged_factors,
                          std::vector<std::uint64_t>& divisors) {
    const std::uint64_t n = static_cast<std::uint64_t>(a) * static_cast<std::uint64_t>(a) - 1ULL;
    if (n == 0) return 0;

    const std::uint64_t denom = static_cast<std::uint64_t>(perimeter_limit - a);
    if (denom == 0) return 0;

    const std::uint64_t v_min = (n + denom - 1) / denom;

    const long double root = std::sqrt(static_cast<long double>(2) * a * a - 1.0L);
    std::int64_t v_max = static_cast<std::int64_t>(std::floor(root - a));
    if (v_max <= 0) return 0;

    while ((static_cast<__int128>(v_max + 1) * (v_max + 1) + static_cast<__int128>(2) * a * (v_max + 1)) <=
           static_cast<__int128>(n)) {
        ++v_max;
    }
    while ((static_cast<__int128>(v_max) * v_max + static_cast<__int128>(2) * a * v_max) >
           static_cast<__int128>(n)) {
        --v_max;
    }

    if (v_max < static_cast<std::int64_t>(v_min)) return 0;

    factorize_int(a - 1, spf, left_factors);
    factorize_int(a + 1, spf, right_factors);
    merge_factors(left_factors, right_factors, merged_factors);

    divisors.clear();
    divisors.push_back(1);
    for (const Factor& f : merged_factors) {
        const std::size_t base_size = divisors.size();
        std::uint64_t mul = 1;
        for (int e = 1; e <= f.e; ++e) {
            mul *= static_cast<std::uint64_t>(f.p);
            for (std::size_t i = 0; i < base_size; ++i) {
                const std::uint64_t candidate = divisors[i] * mul;
                if (candidate <= static_cast<std::uint64_t>(v_max)) {
                    divisors.push_back(candidate);
                }
            }
        }
    }

    std::uint64_t count = 0;
    for (std::uint64_t v : divisors) {
        if (v < v_min) continue;

        const std::uint64_t u = n / v;
        if (((u - v) & 1ULL) != 0ULL) continue;
        if (u < v + static_cast<std::uint64_t>(2 * a)) continue;
        if (static_cast<std::uint64_t>(a) + u > static_cast<std::uint64_t>(perimeter_limit)) continue;

        ++count;
    }

    return count;
}

std::uint64_t count_barely_acute(int perimeter_limit, int threads) {
    if (perimeter_limit < 3) return 0;

    const int a_max = perimeter_limit / 3;
    const std::uint64_t a1_count = static_cast<std::uint64_t>((perimeter_limit - 1) / 2);
    const std::vector<int> spf = build_spf(a_max + 1);

    if (threads < 1) threads = 1;
    if (threads > a_max) threads = a_max;

    std::atomic<int> next_a{2};
    constexpr int kChunk = 2048;

    std::vector<std::uint64_t> partial(static_cast<std::size_t>(threads), 0);
    std::vector<std::thread> pool;
    pool.reserve(static_cast<std::size_t>(threads));

    for (int t = 0; t < threads; ++t) {
        pool.emplace_back([&, t]() {
            std::vector<Factor> left_factors;
            std::vector<Factor> right_factors;
            std::vector<Factor> merged_factors;
            std::vector<std::uint64_t> divisors;

            left_factors.reserve(8);
            right_factors.reserve(8);
            merged_factors.reserve(16);
            divisors.reserve(256);

            std::uint64_t local = 0;
            while (true) {
                const int start = next_a.fetch_add(kChunk, std::memory_order_relaxed);
                if (start > a_max) break;
                const int end = std::min(a_max, start + kChunk - 1);

                for (int a = start; a <= end; ++a) {
                    local += count_for_a(a,
                                         perimeter_limit,
                                         spf,
                                         left_factors,
                                         right_factors,
                                         merged_factors,
                                         divisors);
                }
            }

            partial[static_cast<std::size_t>(t)] = local;
        });
    }

    for (auto& th : pool) th.join();

    std::uint64_t total = a1_count;
    for (std::uint64_t v : partial) total += v;
    return total;
}

std::uint64_t brute_count(int perimeter_limit) {
    std::uint64_t count = 0;
    for (int a = 1; a <= perimeter_limit / 3; ++a) {
        for (int b = a; b <= (perimeter_limit - a) / 2; ++b) {
            const long long c2 = 1LL * a * a + 1LL * b * b - 1LL;
            if (c2 <= 0) continue;

            long long c = static_cast<long long>(std::sqrt(static_cast<long double>(c2)));
            while ((c + 1) * (c + 1) <= c2) ++c;
            while (c * c > c2) --c;

            if (c < b) continue;
            if (c * c != c2) continue;
            if (a + b + c > perimeter_limit) continue;
            if (a + b <= c) continue;

            ++count;
        }
    }
    return count;
}

bool validate() {
    const std::vector<int> brute_limits = {80, 120, 200, 300, 500};
    for (int limit : brute_limits) {
        const std::uint64_t slow = brute_count(limit);
        const std::uint64_t fast = count_barely_acute(limit, 1);
        if (slow != fast) {
            std::cerr << "Validation failed at perimeter=" << limit
                      << ": brute=" << slow << ", fast=" << fast << "\n";
            return false;
        }
    }

    const int thread_check_limit = 50'000;
    const std::uint64_t single = count_barely_acute(thread_check_limit, 1);

    unsigned hw = std::thread::hardware_concurrency();
    if (hw == 0) hw = 2;
    const int multi_threads = static_cast<int>(std::min<unsigned>(hw, 8));
    const std::uint64_t multi = count_barely_acute(thread_check_limit, multi_threads);
    if (single != multi) {
        std::cerr << "Thread consistency failed at perimeter=" << thread_check_limit
                  << ": single=" << single << ", multi=" << multi << "\n";
        return false;
    }

    return true;
}

}  // namespace

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

    int perimeter_limit = kDefaultLimit;
    int threads = static_cast<int>(std::thread::hardware_concurrency());
    if (threads <= 0) threads = 1;

    if (argc > 1) perimeter_limit = std::max(3, std::atoi(argv[1]));
    if (argc > 2) threads = std::max(1, std::atoi(argv[2]));

    const std::uint64_t answer = count_barely_acute(perimeter_limit, threads);
    std::cout << answer << '\n';
    return 0;
}

Python

import math
import multiprocessing

def get_spf(limit):
    spf = [0] * (limit + 1)
    primes = []
    for i in range(2, limit + 1):
        if spf[i] == 0:
            spf[i] = i
            primes.append(i)
        for p in primes:
            if p * i > limit or p > spf[i]:
                break
            spf[p * i] = p
    spf[1] = 1
    return spf

def factorize_int(x, spf):
    factors = []
    while x > 1:
        p = spf[x]
        e = 0
        while x > 1 and spf[x] == p:
            x //= p
            e += 1
        factors.append((p, e))
    return factors

def merge_factors(a, b):
    out = []
    i, j = 0, 0
    while i < len(a) or j < len(b):
        if j == len(b) or (i < len(a) and a[i][0] < b[j][0]):
            out.append(a[i])
            i += 1
        elif i == len(a) or b[j][0] < a[i][0]:
            out.append(b[j])
            j += 1
        else:
            out.append((a[i][0], a[i][1] + b[j][1]))
            i += 1
            j += 1
    return out

def get_divisors(factors, v_max):
    divisors = [1]
    for p, e in factors:
        base_size = len(divisors)
        mul = 1
        for _ in range(1, e + 1):
            mul *= p
            for i in range(base_size):
                candidate = divisors[i] * mul
                if candidate <= v_max:
                    divisors.append(candidate)
    return divisors

def worker(args):
    start_a, end_a, perimeter_limit, spf = args
    local_count = 0
    
    for a in range(start_a, end_a + 1):
        n = a * a - 1
        if n == 0: continue
        
        denom = perimeter_limit - a
        if denom <= 0: continue
        
        v_min = (n + denom - 1) // denom
        
        root = math.sqrt(2.0 * a * a - 1.0)
        v_max = int(math.floor(root - a))
        if v_max <= 0: continue
        
        while (v_max + 1) * (v_max + 1) + 2 * a * (v_max + 1) <= n:
            v_max += 1
        while v_max * v_max + 2 * a * v_max > n:
            v_max -= 1
            
        if v_max < v_min: continue
        
        left_factors = factorize_int(a - 1, spf)
        right_factors = factorize_int(a + 1, spf)
        merged_factors = merge_factors(left_factors, right_factors)
        
        divisors = get_divisors(merged_factors, v_max)
        
        for v in divisors:
            if v < v_min: continue
            u = n // v
            if (u - v) % 2 != 0: continue
            if u < v + 2 * a: continue
            if a + u > perimeter_limit: continue
            local_count += 1
            
    return local_count

def solve(perimeter_limit=25000000):
    if perimeter_limit < 3: return "0"
    
    a_max = perimeter_limit // 3
    a1_count = (perimeter_limit - 1) // 2
    
    threads = multiprocessing.cpu_count() or 1
    threads = min(threads, a_max)
    
    spf = get_spf(a_max + 1)
    
    chunk = 2048
    tasks = []
    start = 2
    while start <= a_max:
        end = min(a_max, start + chunk - 1)
        tasks.append((start, end, perimeter_limit, spf))
        start += chunk
        
    total = a1_count
    
    if threads > 1 and len(tasks) > 1:
        with multiprocessing.Pool(threads) as pool:
            results = pool.map(worker, tasks)
            total += sum(results)
    else:
        for t in tasks:
            total += worker(t)
            
    return str(total)

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

Java

import java.util.*;
import java.util.concurrent.*;

public class Euler223 {
    static class Factor {
        int p, e;

        Factor(int p, int e) {
            this.p = p;
            this.e = e;
        }
    }

    static int[] buildSpf(int n) {
        int[] spf = new int[n + 1];
        List<Integer> primes = new ArrayList<>(n / 10);
        for (int i = 2; i <= n; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
            }
            for (int p : primes) {
                long v = (long) p * i;
                if (v > n || p > spf[i])
                    break;
                spf[(int) v] = p;
            }
        }
        if (n >= 1)
            spf[1] = 1;
        return spf;
    }

    static void factorizeInt(int x, int[] spf, List<Factor> out) {
        out.clear();
        while (x > 1) {
            int p = spf[x];
            int e = 0;
            do {
                x /= p;
                e++;
            } while (x > 1 && spf[x] == p);
            out.add(new Factor(p, e));
        }
    }

    static void mergeFactors(List<Factor> a, List<Factor> b, List<Factor> out) {
        out.clear();
        int i = 0, j = 0;
        while (i < a.size() || j < b.size()) {
            if (j == b.size() || (i < a.size() && a.get(i).p < b.get(j).p)) {
                out.add(a.get(i++));
            } else if (i == a.size() || b.get(j).p < a.get(i).p) {
                out.add(b.get(j++));
            } else {
                out.add(new Factor(a.get(i).p, a.get(i).e + b.get(j).e));
                i++;
                j++;
            }
        }
    }

    public static String solve() {
        int perimeterLimit = 25000000;
        if (perimeterLimit < 3)
            return "0";

        int aMax = perimeterLimit / 3;
        long a1Count = (perimeterLimit - 1) / 2;
        int[] spf = buildSpf(aMax + 1);

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
        threads = Math.min(threads, aMax);

        int chunk = 2048;
        List<int[]> tasks = new ArrayList<>();
        for (int start = 2; start <= aMax; start += chunk) {
            int end = Math.min(aMax, start + chunk - 1);
            tasks.add(new int[] { start, end });
        }

        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<Long>> futures = new ArrayList<>();

        for (int[] task : tasks) {
            final int startA = task[0];
            final int endA = task[1];

            futures.add(executor.submit(() -> {
                long local = 0;
                List<Factor> leftFactors = new ArrayList<>(8);
                List<Factor> rightFactors = new ArrayList<>(8);
                List<Factor> mergedFactors = new ArrayList<>(16);
                long[] divisors = new long[256];

                for (int a = startA; a <= endA; ++a) {
                    long n = (long) a * a - 1;
                    if (n == 0)
                        continue;

                    long denom = perimeterLimit - a;
                    if (denom <= 0)
                        continue;

                    long vMin = (n + denom - 1) / denom;
                    double root = Math.sqrt(2.0 * a * a - 1.0);
                    long vMax = (long) Math.floor(root - a);
                    if (vMax <= 0)
                        continue;

                    while (true) {
                        long nextV = vMax + 1;
                        // Avoid overflow using BigInteger or checking sizes if necessary, but max a ~
                        // 8.3M, a^2 ~ 6e13. n is 6e13.
                        // (vMax+1)^2 + 2*a*(vMax+1)
                        // vMax won't exceed sqrt(n) ~ 8.3M
                        // (8.3M)^2 ~ 6.9e13, fits in signed long (9e18).
                        long test = nextV * nextV + 2L * a * nextV;
                        if (test <= n)
                            vMax++;
                        else
                            break;
                    }
                    while (true) {
                        long test = vMax * vMax + 2L * a * vMax;
                        if (test > n)
                            vMax--;
                        else
                            break;
                    }

                    if (vMax < vMin)
                        continue;

                    factorizeInt(a - 1, spf, leftFactors);
                    factorizeInt(a + 1, spf, rightFactors);
                    mergeFactors(leftFactors, rightFactors, mergedFactors);

                    divisors[0] = 1;
                    int divSize = 1;

                    for (Factor f : mergedFactors) {
                        int baseSize = divSize;
                        long mul = 1;
                        for (int e = 1; e <= f.e; ++e) {
                            mul *= f.p;
                            for (int i = 0; i < baseSize; ++i) {
                                long cand = divisors[i] * mul;
                                if (cand <= vMax) {
                                    if (divSize == divisors.length) {
                                        divisors = Arrays.copyOf(divisors, divisors.length * 2);
                                    }
                                    divisors[divSize++] = cand;
                                }
                            }
                        }
                    }

                    for (int i = 0; i < divSize; ++i) {
                        long v = divisors[i];
                        if (v < vMin)
                            continue;

                        long u = n / v;
                        if (((u - v) & 1) != 0)
                            continue;
                        if (u < v + 2L * a)
                            continue;
                        if (a + u > perimeterLimit)
                            continue;
                        local++;
                    }
                }
                return local;
            }));
        }

        long total = a1Count;
        for (Future<Long> f : futures) {
            try {
                total += f.get();
            } catch (Exception e) {
            }
        }
        executor.shutdown();

        return String.valueOf(total);
    }

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