Project Euler

Almost Right-angled Triangles I

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

Source sync Apr 19, 2026

Problem #0223

Level Level 11

Solved By 1,814

Languages C++, Python

Answer 61614848

Length 332 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 Foundation

Theorem 1 (Difference-of-Squares Factorization). The equation $a^2 + b^2 = c^2 + 1$ is equivalent to $(c - b)(c + b) = (a - 1)(a + 1)$ .

Proof. $a^2 + b^2 = c^2 + 1 \iff c^2 - b^2 = a^2 - 1 \iff (c-b)(c+b) = (a-1)(a+1)$ . $\square$

Lemma 1 (Change of Variables). Setting $d = c - b$ and $e = c + b$ , we have $d \cdot e = (a-1)(a+1)$ , $b = (e - d)/2$ , $c = (e + d)/2$ , and the perimeter is $a + e$ .

Proof. From $d = c - b$ and $e = c + b$ : $b = (e - d)/2$ , $c = (e + d)/2$ . The perimeter is $a + b + c = a + e$ . $\square$

Theorem 2 (Case $a = 1$ ). When $a = 1$ , the equation becomes $b = c$ , contributing $(L - 1)/2$ solutions where $L = 25{,}000{,}000$ .

Proof. $a = 1$ gives $(a-1)(a+1) = 0$ , so $d = 0$ and $b = c$ . The constraint $1 + 2b \leq L$ gives $b \leq (L-1)/2 = 12{,}499{,}999$ . $\square$

Theorem 3 (Coprime Factorization, Even $a$ ). For even $a \geq 2$ , $n = (a-1)(a+1)$ is odd, $\gcd(a-1, a+1) = 1$ , and every divisor of $n$ factors uniquely as a product of a divisor of $(a-1)$ and a divisor of $(a+1)$ .

Proof. Since $a$ is even, $a - 1$ and $a + 1$ are both odd, so $n$ is odd. Since $(a+1) - (a-1) = 2$ and both are odd, $\gcd(a-1, a+1) = 1$ . By the Chinese Remainder Theorem, the divisors of $n = (a-1)(a+1)$ biject with pairs of divisors of $a-1$ and $a+1$ . $\square$

Theorem 4 (Coprime Factorization, Odd $a$ ). For odd $a \geq 3$ , $4 \mid n$ and writing $d = 2d'$ , $e = 2e'$ gives $d' e' = \frac{a-1}{2} \cdot \frac{a+1}{2}$ with $\gcd\!\left(\frac{a-1}{2}, \frac{a+1}{2}\right) = 1$ .

Proof. For odd $a$ , both $a - 1$ and $a + 1$ are even, so $4 \mid (a-1)(a+1)$ . For $b$ and $c$ to be integers, $d$ and $e$ must have the same parity. Since $n$ is even, they must both be even. Write $d = 2d'$ , $e = 2e'$ :

d'e' = \frac{n}{4} = \frac{a-1}{2} \cdot \frac{a+1}{2}.

Since $\frac{a+1}{2} - \frac{a-1}{2} = 1$ , these factors are coprime. $\square$

Lemma 2 (Constraints). For a valid triple, the divisor pair $(d, e)$ must satisfy: (i) $d \leq e$ , (ii) $e - d \geq 2a$ (so that $b \geq a$ ), and (iii) $a + e \leq L$ .

Proof. (i) ensures $b \leq c$ . (ii) follows from $b = (e-d)/2 \geq a$ . (iii) is the perimeter bound. $\square$

Editorial

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, de = (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. We else.

Pseudocode

if a is even
else

Complexity Analysis

Time: $O(L \log L)$ . For each $a$ , we enumerate divisors of two coprime factors. The average number of divisors is $O(\log n)$ , giving total work $O\!\left(\sum_{a=2}^{L/3} \tau((a-1)(a+1))\right) = O(L \log L)$ .
Space: $O(L/3)$ 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()