Project Euler

Almost Right-angled Triangles I

How many ordered triples (a,b,c) with a <= b <= c, a^2+b^2=c^2+1, a+b+c <= 25,000,000 exist?

Source sync May 21, 2026

Problem #0223

Level Level 11

Solved By 1,814

Languages C++, Python

Answer 61614848

Length 262 words

number_theorygeometrydynamic_programming

Let us call an integer sided triangle with sides $a \le b \le c$ barely acute if the sides satisfy $a^2 + b^2 = c^2 + 1$.

How many barely acute triangles are there with perimeter $\le 25\,000\,000$?

Problem 223: Almost Right-angled Triangles I

Mathematical Development

Starting from

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

we factor the difference of squares:

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

Introduce

d=c-b, \qquad e=c+b.

Then

de=(a-1)(a+1), \qquad b=\frac{e-d}{2}, \qquad c=\frac{e+d}{2},

and the perimeter becomes

a+b+c = a+e.

So for each fixed $a$ , we only need factor pairs $(d,e)$ of $a^2-1$ satisfying:

$d \le e$ ,
$d$ and $e$ have the same parity,
$e-d \ge 2a$ so that $b \ge a$ ,
$a+e \le L$ .

The arithmetic simplifies because the two factors around $a$ are almost coprime.

If $a$ is even, then $a-1$ and $a+1$ are odd and coprime, so every divisor of $a^2-1$ is a product of one divisor of $a-1$ and one divisor of $a+1$ .

If $a$ is odd, then $d$ and $e$ must both be even. Writing

d=2d', \qquad e=2e'

gives

d'e'=\frac{a-1}{2}\cdot\frac{a+1}{2},

and the two factors on the right are consecutive integers, hence coprime.

That is why a smallest-prime-factor sieve up to $L/3$ is enough: each case reduces to enumerating divisors of two coprime numbers.

Editorial

The factorization

(c-b)(c+b)=(a-1)(a+1)

turns the triangle count into a divisor problem. Once $a$ is fixed, the perimeter is just $a+e$ , so the large variable is not $c$ but the factor $e=c+b$ . That is the useful reduction: instead of searching in three dimensions, we enumerate factorizations of $a^2-1$ and recover $b$ and $c$ from one complementary divisor pair.

The reason this is fast enough is the coprimality built into $a-1$ and $a+1$ . For even $a$ , they are already coprime. For odd $a$ , the only common factor is the expected factor $2$ , and removing it leaves consecutive coprime integers. So the code never factors the full product from scratch. It sieves smallest prime factors once, generates divisors of the two smaller coprime pieces, combines them into $(d,e)$ , and checks the geometric constraints.

Pseudocode

Count the special family a = 1 separately:
    then b = c and the perimeter condition is immediate.

Build a smallest-prime-factor sieve up to L / 3.

For each even a from 2 to L / 3:
    factor a - 1 and a + 1
    generate all divisors of each factor
    for every product d = d1 * d2:
        let e be the complementary factor
        keep it only if d <= e, e - d >= 2a, and a + e <= L

For each odd a from 3 to L / 3:
    set u = (a - 1) / 2 and v = (a + 1) / 2
    factor u and v
    generate all divisors of each factor
    for every product d' = d1 * d2:
        let e' be the complementary factor
        turn them back into d = 2d' and e = 2e'
        keep the pair only if the same three inequalities hold

Return the total count.

Complexity Analysis

Time: Roughly $O(L \log L)$ in practice, dominated by divisor generation over all $a \le L/3$ .
Space: $O(L)$ for the smallest-prime-factor sieve.

Answer

\boxed{61614848}

C++ project_euler/problem_223/solution.cpp

#include <bits/stdc++.h>
using namespace std;

int main(){
    // a^2 + b^2 = c^2 + 1, a <= b <= c, a+b+c <= L
    // (c-b)(c+b) = (a-1)(a+1), c-b = d, c+b = e, d*e = (a-1)(a+1)
    // b = (e-d)/2, c = (e+d)/2, perimeter = a + e
    // Need: d,e same parity, d <= e, b >= a <=> e - d >= 2a, a + e <= L

    // Case a=1: b=c, count = (L-1)/2
    // Case a>=2: d*e = a^2-1, d <= e, same parity, e >= d+2a, a+e <= L

    // Approach: iterate d from 1 to ~sqrt(L^2) and for each d, iterate
    // over values of a such that d | (a-1)(a+1) and constraints hold.

    // Actually, better: iterate d and e directly.
    // d*e = a^2 - 1 = (a-1)(a+1). Given d and e, a = sqrt(d*e + 1).
    // So iterate d, e with d <= e, d*e + 1 is a perfect square, same parity.

    // But this is still too many pairs...

    // Best approach: iterate d from 1 to some bound.
    // For each d, iterate a such that d | (a^2-1).
    // Given d | (a^2-1), i.e., a^2 ≡ 1 (mod d), i.e., a ≡ ±1 (mod p^k) for each prime power in d.
    // Then e = (a^2-1)/d. Conditions: e >= d (i.e., a^2-1 >= d^2), same parity, e-d >= 2a, a+e <= L.

    // For each d, the valid a values form arithmetic progressions mod d.
    // We can find them by CRT. But implementing CRT for each d is complex.

    // Simplest fast approach: iterate (d, a) where d | (a-1) or d | (a+1) or
    // d shares factors with both.

    // Actually, let's factor differently. Write a-1 = s, a+1 = s+2.
    // d*e = s*(s+2). For each factorization d*e = s*(s+2) with d <= e:
    // We can write d = gcd(d, s) * gcd(d, s+2) roughly...

    // Let me try a different decomposition. Let g = gcd(s, s+2) = gcd(s, 2).
    // If s even (a odd): g = 2, s = 2u, s+2 = 2(u+1), product = 4u(u+1).
    //   d and e must be both even. d = 2d', e = 2e', d'*e' = u(u+1).
    //   u and u+1 are coprime. So divisors of u(u+1) are products of
    //   a divisor of u and a divisor of u+1.
    //   For each divisor d1 of u and d2 of (u+1): d' = d1*d2, e' = (u/d1)*((u+1)/d2).
    //   Need d' <= e'.

    // If s odd (a even): g = 1, product = s*(s+2), gcd(s, s+2) = 1.
    //   d and e must both be odd (since product is odd).
    //   Divisors of s*(s+2) are products of divisor of s and divisor of s+2.
    //   For each d1|s, d2|(s+2): d = d1*d2, e = (s/d1)*((s+2)/d2). Need d <= e.

    // So the approach is: iterate s (= a-1) and for each s, factor s and s+2,
    // then enumerate divisor pairs.

    // But factoring each s is expensive for s up to 8M... unless we use a sieve.

    const long long L = 25000000LL;
    long long count = 0;

    // Case a = 1
    count += (L - 1) / 2;

    // Sieve: for each n up to L, store its smallest prime factor
    // a ranges from 2 to L/3, so s = a-1 ranges from 1 to L/3-1 ~ 8.3M
    // s+2 ranges up to L/3+1. We need to factor s and s+2.
    // Max value to factor: L/3 + 1 ~ 8.3M

    const int SMAX = L / 3 + 2;
    vector<int> spf(SMAX + 1, 0); // smallest prime factor
    for(int i = 2; i <= SMAX; i++){
        if(spf[i] == 0){
            for(int j = i; j <= SMAX; j += i){
                if(spf[j] == 0) spf[j] = i;
            }
        }
    }

    auto get_divisors = [&](long long n) -> vector<long long> {
        if(n <= 0) return {};
        if(n == 1) return {1};
        vector<long long> divs = {1};
        long long tmp = n;
        while(tmp > 1){
            int p = spf[tmp];
            int cnt = 0;
            while(tmp % p == 0){ tmp /= p; cnt++; }
            int sz = divs.size();
            long long pk = 1;
            for(int i = 0; i < cnt; i++){
                pk *= p;
                for(int j = 0; j < sz; j++){
                    divs.push_back(divs[j] * pk);
                }
            }
        }
        return divs;
    };

    // Case a even (s = a-1 odd): s*(s+2) is odd. d, e both odd.
    // gcd(s, s+2) = 1. d = d1*d2 where d1|s, d2|(s+2).
    for(long long a = 2; a <= L / 3; a += 2){
        long long s = a - 1; // odd
        long long s2 = a + 1; // odd
        auto divs_s = get_divisors(s);
        auto divs_s2 = get_divisors(s2);
        for(long long d1 : divs_s){
            for(long long d2 : divs_s2){
                long long d = d1 * d2;
                long long e = (s / d1) * (s2 / d2);
                if(d > e) continue;
                // b = (e-d)/2, need b >= a => e - d >= 2a
                if(e - d < 2*a) continue;
                // perimeter = a + e <= L
                if(a + e > L) continue;
                count++;
            }
        }
    }

    // Case a odd (s = a-1 even): 4 | s*(s+2). d, e both even.
    // s = 2u, s+2 = 2(u+1), d = 2d', e = 2e', d'*e' = u*(u+1).
    // gcd(u, u+1) = 1. d' = d1*d2 where d1|u, d2|(u+1).
    for(long long a = 3; a <= L / 3; a += 2){
        long long u = (a - 1) / 2;
        long long u1 = u + 1;
        auto divs_u = get_divisors(u);
        auto divs_u1 = get_divisors(u1);
        for(long long d1 : divs_u){
            for(long long d2 : divs_u1){
                long long dp = d1 * d2;
                long long ep = (u / d1) * (u1 / d2);
                if(dp > ep) continue;
                long long d = 2 * dp;
                long long e = 2 * ep;
                long long b = (e - d) / 2;
                if(b < a) continue;
                if(a + e > L) continue;
                count++;
            }
        }
    }

    cout << count << endl;
    return 0;
}

Python project_euler/problem_223/solution.py

"""
Problem 223: Almost Right-angled Triangles I

Count ordered triples (a, b, c) with a <= b <= c,
a^2 + b^2 = c^2 + 1, and a + b + c <= 25,000,000.

Key: (c-b)(c+b) = (a-1)(a+1). Let d = c-b, e = c+b, d*e = (a-1)(a+1).
b = (e-d)/2, c = (e+d)/2, perimeter = a + e.
Need d,e same parity, d <= e, b >= a (e-d >= 2a), a+e <= L.

For a even: n = a^2-1 is odd, so d,e both odd. n = (a-1)(a+1), gcd=1.
For a odd: n = a^2-1 divisible by 4. d,e both even. n/4 = ((a-1)/2)*((a+1)/2), gcd=1.

Factor using coprime factorization to enumerate divisors efficiently.
"""

import math

def get_divisors_via_spf(n, spf):
    """Get all divisors of n using smallest prime factor array."""
    if n == 0:
        return []
    if n == 1:
        return [1]
    divs = [1]
    tmp = n
    while tmp > 1:
        p = spf[tmp]
        exp = 0
        while tmp % p == 0:
            tmp //= p
            exp += 1
        new_divs = []
        pk = 1
        for e in range(exp):
            pk *= p
            for d in divs:
                new_divs.append(d * pk)
        divs.extend(new_divs)
    return divs

def solve():
    L = 25_000_000
    count = 0

    # Case a = 1: b = c, perimeter = 1 + 2b <= L
    count += (L - 1) // 2

    # Sieve smallest prime factor up to L/3 + 2
    SMAX = L // 3 + 2
    spf = list(range(SMAX + 1))  # spf[i] = i initially
    for i in range(2, int(math.isqrt(SMAX)) + 1):
        if spf[i] == i:  # i is prime
            for j in range(i * i, SMAX + 1, i):
                if spf[j] == j:
                    spf[j] = i

    def get_divisors(n):
        if n <= 1:
            return [1] if n == 1 else []
        divs = [1]
        tmp = n
        while tmp > 1:
            p = spf[tmp] if tmp <= SMAX else tmp  # for large tmp, it's prime
            exp = 0
            while tmp % p == 0:
                tmp //= p
                exp += 1
            sz = len(divs)
            pk = 1
            for e in range(exp):
                pk *= p
                for j in range(sz):
                    divs.append(divs[j] * pk)
        return divs

    # Case a even (a = 2, 4, 6, ...):
    # n = (a-1)*(a+1) is odd, gcd(a-1, a+1) = 1.
    # Divisors of n = divisors of (a-1) * divisors of (a+1).
    for a in range(2, L // 3 + 1, 2):
        s = a - 1  # odd
        s2 = a + 1  # odd
        divs_s = get_divisors(s)
        divs_s2 = get_divisors(s2)
        for d1 in divs_s:
            for d2 in divs_s2:
                d = d1 * d2
                e = (s // d1) * (s2 // d2)
                if d > e:
                    continue
                if e - d < 2 * a:
                    continue
                if a + e > L:
                    continue
                count += 1

    # Case a odd (a = 3, 5, 7, ...):
    # n = (a-1)*(a+1) divisible by 4. d,e both even.
    # u = (a-1)//2, u1 = (a+1)//2, gcd(u, u1) = 1.
    # d' * e' = u * u1, d = 2d', e = 2e'.
    for a in range(3, L // 3 + 1, 2):
        u = (a - 1) // 2
        u1 = (a + 1) // 2
        divs_u = get_divisors(u)
        divs_u1 = get_divisors(u1)
        for d1 in divs_u:
            for d2 in divs_u1:
                dp = d1 * d2
                ep = (u // d1) * (u1 // d2)
                if dp > ep:
                    continue
                d = 2 * dp
                e = 2 * ep
                if e - d < 2 * a:
                    continue
                if a + e > L:
                    continue
                count += 1

    print(count)

if __name__ == "__main__":
    solve()