Problem 299: Three Similar Triangles

Project Euler Problem 299 Solution

EulerSolve provides an optimized solution for Project Euler Problem 299, Three Similar Triangles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The problem starts with $$A(a,0),\qquad B(b,0),\qquad C(0,a),\qquad D(0,d),$$ where $0<a<b$ and $0<a<d$. We seek integer triplets $(a,b,d)$ for which there exists an integer point $P$ on $AC$ such that the three triangles $ABP$, $CDP$, and $BDP$ are all similar. The quantity to count is the number of such triplets with $$b+d<N.$$ Mathematical Approach 1. Why $a=c$ and why two cases appear Project Euler already states that similarity is possible only when $a=c$, so the segment $AC$ lies on the line $$x+y=a.$$ Because $OA$ and $OC$ are symmetric, the angles at $A$ and $C$ are both $135^\circ$. Matching the three similar triangles leaves two possible angle correspondences: Incenter case. The lines $PB$ and $PD$ are angle bisectors of the right triangle $OBD$, so $P$ is the incenter of $\triangle OBD$. Parallel case. The line $AC$ is parallel to $BD$. These two families are disjoint, and the code counts them separately. 2. Incenter case: reduction to a Pythagorean triple If $P$ is the incenter of the right triangle with legs $b$ and $d$, then its coordinates must be $$P=(i,i),$$ where $i$ is the inradius. Since $P\in AC$, we have $$a=2i.$$ For a right triangle, the inradius is $$i=\frac{b+d-\sqrt{b^2+d^2}}{2}.$$ Therefore $a$ is integral exactly when $\sqrt{b^2+d^2}$ is integral....

Detailed mathematical approach

Problem Summary

The problem starts with

$$A(a,0),\qquad B(b,0),\qquad C(0,a),\qquad D(0,d),$$

where $0<a<b$ and $0<a<d$. We seek integer triplets $(a,b,d)$ for which there exists an integer point $P$ on $AC$ such that the three triangles $ABP$, $CDP$, and $BDP$ are all similar. The quantity to count is the number of such triplets with

$$b+d<N.$$

Mathematical Approach

1. Why $a=c$ and why two cases appear

Project Euler already states that similarity is possible only when $a=c$, so the segment $AC$ lies on the line

$$x+y=a.$$

Because $OA$ and $OC$ are symmetric, the angles at $A$ and $C$ are both $135^\circ$. Matching the three similar triangles leaves two possible angle correspondences:

Incenter case. The lines $PB$ and $PD$ are angle bisectors of the right triangle $OBD$, so $P$ is the incenter of $\triangle OBD$.

Parallel case. The line $AC$ is parallel to $BD$.

These two families are disjoint, and the code counts them separately.

2. Incenter case: reduction to a Pythagorean triple

If $P$ is the incenter of the right triangle with legs $b$ and $d$, then its coordinates must be

$$P=(i,i),$$

where $i$ is the inradius. Since $P\in AC$, we have

$$a=2i.$$

For a right triangle, the inradius is

$$i=\frac{b+d-\sqrt{b^2+d^2}}{2}.$$

Therefore $a$ is integral exactly when $\sqrt{b^2+d^2}$ is integral. So the incenter family is equivalent to counting right triangles with integer legs $b,d$ and hypotenuse

$$z=\sqrt{b^2+d^2},$$

under the limit $b+d<N$.

The sample $(a,b,d)=(2,3,4)$ comes from the primitive Pythagorean triple

$$3^2+4^2=5^2,$$

because

$$a=b+d-z=3+4-5=2,$$

and then $P=(1,1)$ lies on $x+y=2$.

3. Why the code uses $s_1$ and $s_2$

Write

$$u=b-a,\qquad v=d-a.$$

Then

$$b=a+u,\qquad d=a+v,\qquad z=b+d-a=a+u+v.$$

Substituting $z^2=b^2+d^2$ gives

$$a^2=2uv.$$

So the primitive solutions come in two orientations:

$$u=g m^2,\qquad v=2g n^2,\qquad a=2gmn,$$

$$u=2g m^2,\qquad v=g n^2,\qquad a=2gmn,$$

with $\gcd(m,n)=1$. After reconstructing $b$ and $d$, the perimeter bound becomes

$$b+d=g(m^2+4mn+2n^2)$$

in the first orientation, and

$$b+d=g(2m^2+4mn+n^2)$$

in the second. These are exactly the denominators

$$s_1=m^2+4mn+2n^2,\qquad s_2=2m^2+4mn+n^2$$

used by the code. For fixed primitive $(m,n)$, all scaled copies are obtained by multiplying by $g$, so the number of valid scales is

$$\left\lfloor\frac{N-1}{s_1}\right\rfloor \quad \text{or} \quad \left\lfloor\frac{N-1}{s_2}\right\rfloor.$$

The parity filters in the implementation are the primitive-normalization rules that prevent duplicate descriptions.

4. Parallel case: reduction to $Q^2+2f^2=a^2$

In the second similarity pattern we have

$$AC\parallel BD,$$

so the slope condition forces

$$d=b.$$

Let

$$f=b-a>0.$$

Then $D=(0,b)$, $B=(b,0)$, and the circumcenter of $\triangle BDP$ is $X=(b,b)$. Hence $P$ must lie on the circle

$$ (x-b)^2+(y-b)^2=b^2,$$

and also on the line $AC$, namely

$$y=a-x.$$

Substituting gives

$$x=\frac{a\pm \sqrt{a^2-2f^2}}{2}.$$

Therefore an integer point $P$ exists exactly when

$$Q^2=a^2-2f^2,$$

that is,

$$Q^2+2f^2=a^2.$$

A small example is

$$a=3,\qquad f=2,\qquad b=d=5,$$

because $1^2+2\cdot 2^2=3^2$, giving the triplet $(3,5,5)$.

5. Why the code uses $s_3$

Primitive solutions of

$$Q^2+2f^2=a^2$$

are parameterized by coprime integers $p,q$ with odd $p$:

$$Q=p^2-2q^2,\qquad f=2pq,\qquad a=p^2+2q^2.$$

After scaling by $g$, we get

$$b=d=a+f=g(p^2+2pq+2q^2).$$

Hence the bound $b+d<N$ becomes

$$2g(p^2+2pq+2q^2)<N,$$

so for each primitive pair the number of scales is

$$\left\lfloor\frac{N-1}{2(p^2+2pq+2q^2)}\right\rfloor.$$

This is exactly the denominator

$$s_3=2(p^2+2pq+2q^2)$$

used in the implementation.

6. Final counting formula and checkpoints

The two families are disjoint, so the total is

$$T(N)=\sum \left\lfloor\frac{N-1}{s_1}\right\rfloor+\sum \left\lfloor\frac{N-1}{s_2}\right\rfloor+\sum \left\lfloor\frac{N-1}{s_3}\right\rfloor.$$

The code validates itself with the published checkpoints

$$T(100)=92,\qquad T(100000)=320471.$$

How the Code Works

count_family_x_eq_y(...) enumerates the incenter family via the two forms $s_1,s_2$. count_family_u_eq_v(...) enumerates the parallel family via $s_3$. For each primitive parameter pair, the contribution is the number of scale factors $g$ allowed by $b+d<N$. The final answer is the sum of both family totals.

Complexity Analysis

The loops over $(m,n)$ and $(p,q)$ break as soon as the corresponding denominator reaches $N$. So the practical runtime is far below a naive scan over all triples $(a,b,d)$. Memory usage is $O(1)$ beyond a few counters, and the C++ / Python versions optionally split independent parameter ranges across threads or processes.

Problem 299 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <exception>
#include <iostream>
#include <limits>
#include <numeric>
#include <stdexcept>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

struct Options {
    u64 limit = 100000000ULL;
    bool run_checkpoints = true;
    unsigned requested_threads = 0U;
};

bool parse_u64_after_prefix(const std::string& arg,
                            const std::string& prefix,
                            u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }

        const u64 digit = static_cast<u64>(ch - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            throw std::overflow_error("u64 parse overflow");
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const std::string& prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        throw std::overflow_error("unsigned parse overflow");
    }

    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const 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, "--limit=", options.limit)) {
            continue;
        }
        if (parse_unsigned_after_prefix(arg, "--threads=", options.requested_threads)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

unsigned pick_thread_count(const unsigned requested_threads) {
    if (requested_threads > 0U) {
        return requested_threads;
    }

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0U) {
        threads = 1U;
    }
    return threads;
}

template <typename Worker>
u64 parallel_sum_range(const int begin,
                       const int end,
                       unsigned threads,
                       const Worker& worker) {
    if (begin >= end) {
        return 0ULL;
    }

    const int work_items = end - begin;
    if (threads == 0U) {
        threads = 1U;
    }
    if (threads == 1U || work_items < 128) {
        return worker(begin, end);
    }

    threads = std::min<unsigned>(threads, static_cast<unsigned>(work_items));

    std::vector<std::thread> pool;
    pool.reserve(threads);

    std::vector<u64> partial(threads, 0ULL);

    for (unsigned tid = 0U; tid < threads; ++tid) {
        const int lo = begin + (work_items * static_cast<int>(tid)) /
                                   static_cast<int>(threads);
        const int hi = begin +
                       (work_items * static_cast<int>(tid + 1U)) /
                           static_cast<int>(threads);

        pool.emplace_back([&, tid, lo, hi]() { partial[tid] = worker(lo, hi); });
    }

    for (std::thread& th : pool) {
        th.join();
    }

    u64 total = 0ULL;
    for (const u64 x : partial) {
        total += x;
    }
    return total;
}

u64 count_family_x_eq_y(const u64 limit, unsigned threads) {
    if (limit <= 7ULL) {
        return 0ULL;
    }

    const u64 numerator = limit - 1ULL;
    const int m_max = static_cast<int>(std::sqrt(static_cast<long double>(limit))) + 3;

    auto worker = [limit, numerator](const int lo, const int hi) -> u64 {
        u64 local = 0ULL;

        for (int m = lo; m < hi; ++m) {
            const u64 m64 = static_cast<u64>(m);
            const u64 m2 = m64 * m64;

            // Orientation 1: u = g*m^2, v = 2*g*n^2 (m must be odd).
            if ((m & 1) != 0) {
                for (int n = 1;; ++n) {
                    const u64 n64 = static_cast<u64>(n);
                    const u64 s = m2 + 4ULL * m64 * n64 + 2ULL * n64 * n64;
                    if (s >= limit) {
                        break;
                    }
                    if (std::gcd(m, n) == 1) {
                        local += numerator / s;
                    }
                }
            }

            // Orientation 2: u = 2*g*m^2, v = g*n^2 (n must be odd).
            for (int n = 1;; ++n) {
                const u64 n64 = static_cast<u64>(n);
                const u64 s = 2ULL * m2 + 4ULL * m64 * n64 + n64 * n64;
                if (s >= limit) {
                    break;
                }
                if ((n & 1) != 0 && std::gcd(m, n) == 1) {
                    local += numerator / s;
                }
            }
        }

        return local;
    };

    return parallel_sum_range(1, m_max + 1, threads, worker);
}

u64 count_family_u_eq_v(const u64 limit, unsigned threads) {
    if (limit <= 13ULL) {
        return 0ULL;
    }

    const u64 numerator = limit - 1ULL;
    const int p_max = static_cast<int>(std::sqrt(static_cast<long double>(limit / 2ULL))) + 3;
    const int odd_count = (p_max + 1) / 2;

    auto worker = [limit, numerator](const int lo, const int hi) -> u64 {
        u64 local = 0ULL;

        for (int idx = lo; idx < hi; ++idx) {
            const int p = 2 * idx + 1;
            const u64 p64 = static_cast<u64>(p);
            const u64 p2 = p64 * p64;

            for (int q = 1;; ++q) {
                const u64 q64 = static_cast<u64>(q);
                const u64 s = p2 + 2ULL * p64 * q64 + 2ULL * q64 * q64;
                const u128 denom = static_cast<u128>(2ULL) * static_cast<u128>(s);
                if (denom >= static_cast<u128>(limit)) {
                    break;
                }
                if (std::gcd(p, q) == 1) {
                    local += numerator / static_cast<u64>(denom);
                }
            }
        }

        return local;
    };

    return parallel_sum_range(0, odd_count, threads, worker);
}

u64 count_triplets(const u64 limit, unsigned threads) {
    if (limit <= 0ULL) {
        return 0ULL;
    }

    if (limit < 1000000ULL) {
        threads = 1U;
    }

    const u64 family_x_eq_y = count_family_x_eq_y(limit, threads);
    const u64 family_u_eq_v = count_family_u_eq_v(limit, threads);
    return family_x_eq_y + family_u_eq_v;
}

void run_checkpoints(const unsigned threads) {
    struct Checkpoint {
        u64 limit;
        u64 expected;
    };

    const std::vector<Checkpoint> checkpoints = {
        {100ULL, 92ULL},
        {100000ULL, 320471ULL},
    };

    for (const Checkpoint& checkpoint : checkpoints) {
        const u64 got = count_triplets(checkpoint.limit, threads);
        if (got != checkpoint.expected) {
            throw std::runtime_error(
                "Checkpoint failed for limit=" + std::to_string(checkpoint.limit) +
                ": expected " + std::to_string(checkpoint.expected) +
                ", got " + std::to_string(got));
        }
    }
}

} // namespace

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

        const unsigned threads = pick_thread_count(options.requested_threads);

        if (options.run_checkpoints) {
            run_checkpoints(threads);
        }

        const u64 answer = count_triplets(options.limit, threads);
        std::cout << answer << '\n';
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

import math
from concurrent.futures import ProcessPoolExecutor
import multiprocessing

def count_x_eq_y_worker(lo, hi, limit, numerator):
    local = 0
    for m in range(lo, hi):
        m2 = m * m
        if m % 2 != 0:
            n = 1
            while True:
                s = m2 + 4 * m * n + 2 * n * n
                if s >= limit: break
                if math.gcd(m, n) == 1:
                    local += numerator // s
                n += 1
                
        n = 1
        while True:
            s = 2 * m2 + 4 * m * n + n * n
            if s >= limit: break
            if n % 2 != 0 and math.gcd(m, n) == 1:
                local += numerator // s
            n += 1
            
    return local

def count_family_x_eq_y(limit, threads):
    if limit <= 7: return 0
    numerator = limit - 1
    m_max = math.isqrt(limit) + 3
    
    if threads <= 1:
        return count_x_eq_y_worker(1, m_max + 1, limit, numerator)
        
    chunk = (m_max + threads) // threads
    ranges = [(max(1, i * chunk), min(m_max + 1, (i + 1) * chunk)) for i in range(threads)]
    
    total = 0
    with ProcessPoolExecutor(max_workers=threads) as executor:
        futures = [executor.submit(count_x_eq_y_worker, r[0], r[1], limit, numerator) for r in ranges]
        for f in futures:
            total += f.result()
            
    return total

def count_u_eq_v_worker(lo, hi, limit, numerator):
    local = 0
    for idx in range(lo, hi):
        p = 2 * idx + 1
        p2 = p * p
        q = 1
        while True:
            s = p2 + 2 * p * q + 2 * q * q
            denom = 2 * s
            if denom >= limit: break
            if math.gcd(p, q) == 1:
                local += numerator // denom
            q += 1
            
    return local

def count_family_u_eq_v(limit, threads):
    if limit <= 13: return 0
    numerator = limit - 1
    p_max = math.isqrt(limit // 2) + 3
    odd_count = (p_max + 1) // 2
    
    if threads <= 1:
        return count_u_eq_v_worker(0, odd_count, limit, numerator)
        
    chunk = (odd_count + threads - 1) // threads
    ranges = [(i * chunk, min(odd_count, (i + 1) * chunk)) for i in range(threads)]
    
    total = 0
    with ProcessPoolExecutor(max_workers=threads) as executor:
        futures = [executor.submit(count_u_eq_v_worker, r[0], r[1], limit, numerator) for r in ranges]
        for f in futures:
            total += f.result()
            
    return total

def count_triplets(limit):
    threads = multiprocessing.cpu_count()
    if threads == 0: threads = 1
    
    if limit < 1000000:
        threads = 1
        
    total_x = count_family_x_eq_y(limit, threads)
    total_u = count_family_u_eq_v(limit, threads)
    return total_x + total_u

def solve(limit=100000000):
    return str(count_triplets(limit))

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

Java

import java.util.ArrayList;
import java.util.List;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;

public class Euler299 {
    static int gcd(int a, int b) {
        if (b == 0)
            return a;
        return gcd(b, a % b);
    }

    static long countXEqYWorker(int lo, int hi, long limit, long numerator) {
        long local = 0;
        for (int m = lo; m < hi; ++m) {
            long m2 = (long) m * m;

            if ((m & 1) != 0) {
                for (int n = 1;; ++n) {
                    long s = m2 + 4L * m * n + 2L * n * n;
                    if (s >= limit)
                        break;
                    if (gcd(m, n) == 1) {
                        local += numerator / s;
                    }
                }
            }

            for (int n = 1;; ++n) {
                long s = 2L * m2 + 4L * m * n + (long) n * n;
                if (s >= limit)
                    break;
                if ((n & 1) != 0 && gcd(m, n) == 1) {
                    local += numerator / s;
                }
            }
        }
        return local;
    }

    static long countFamilyXEqY(long limit, int threads) {
        if (limit <= 7)
            return 0;
        long numerator = limit - 1;
        int mMax = (int) Math.sqrt(limit) + 3;

        if (threads <= 1) {
            return countXEqYWorker(1, mMax + 1, limit, numerator);
        }

        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<Long>> futures = new ArrayList<>();
        int chunk = (mMax + threads) / threads;

        for (int i = 0; i < threads; ++i) {
            final int lo = Math.max(1, i * chunk);
            final int hi = Math.min(mMax + 1, (i + 1) * chunk);
            if (lo >= hi)
                continue;
            futures.add(executor.submit(() -> countXEqYWorker(lo, hi, limit, numerator)));
        }

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

    static long countUEqVWorker(int lo, int hi, long limit, long numerator) {
        long local = 0;
        for (int idx = lo; idx < hi; ++idx) {
            long p = 2L * idx + 1;
            long p2 = p * p;
            for (long q = 1;; ++q) {
                long s = p2 + 2L * p * q + 2L * q * q;
                long denom = 2L * s;
                if (denom >= limit)
                    break;
                if (gcd((int) p, (int) q) == 1) {
                    local += numerator / denom;
                }
            }
        }
        return local;
    }

    static long countFamilyUEqV(long limit, int threads) {
        if (limit <= 13)
            return 0;
        long numerator = limit - 1;
        int pMax = (int) Math.sqrt(limit / 2) + 3;
        int oddCount = (pMax + 1) / 2;

        if (threads <= 1) {
            return countUEqVWorker(0, oddCount, limit, numerator);
        }

        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<Long>> futures = new ArrayList<>();
        int chunk = (oddCount + threads - 1) / threads;

        for (int i = 0; i < threads; ++i) {
            final int lo = i * chunk;
            final int hi = Math.min(oddCount, (i + 1) * chunk);
            if (lo >= hi)
                continue;
            futures.add(executor.submit(() -> countUEqVWorker(lo, hi, limit, numerator)));
        }

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

    public static String solve() {
        long limit = 100000000L;
        int threads = Runtime.getRuntime().availableProcessors();
        if (threads == 0)
            threads = 1;

        long totalX = countFamilyXEqY(limit, threads);
        long totalU = countFamilyUEqV(limit, threads);

        return String.valueOf(totalX + totalU);
    }

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