Problem 404: Crisscross Ellipses

Project Euler Problem 404 Solution

EulerSolve provides an optimized solution for Project Euler Problem 404, Crisscross Ellipses, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The implementations count integer triples $(a,b,c)$ with $$F(N)=\#\left\{(a,b,c)\in\mathbb{Z}_{>0}^3:\ a\le N,\ a<b<c<2a,\ c^2=\frac{4a^2b^2}{5b^2-4a^2}\right\}.$$ A direct search over $(a,b,c)$ is far too slow for the actual limit. The key reduction is to convert the Diophantine condition into a rational point on a conic, classify primitive solutions once, and then count every scaled multiple in one floor division. Mathematical Approach Step 1: Rewrite the Constraint as a Conic Start from $$c^2=\frac{4a^2b^2}{5b^2-4a^2}.$$ After rearranging and dividing by $b^2c^2$, we obtain $$\left(\frac{2a}{b}\right)^2+\left(\frac{2a}{c}\right)^2=5.$$ So if we define $$X=\frac{2a}{b},\qquad Y=\frac{2a}{c},$$ then $(X,Y)$ is a rational point on the conic $X^2+Y^2=5$. The ordering $a<b<c<2a$ translates into $$1<Y<X<2.$$ This is the geometric core of the fast solution: the original triple-counting problem becomes a count of rational points in one canonical region of the conic. Step 2: Pass to Primitive Integer Data Write the rational point in lowest terms as $$X=\frac{U}{W},\qquad Y=\frac{V}{W},\qquad \gcd(U,V,W)=1.$$ Then $$U^2+V^2=5W^2.$$ Because $1<Y<X<2$, the reduced integers satisfy $$W<V<U<2W.$$ Any common divisor of two of $U,V,W$ would also divide the third, so primitiveness implies pairwise coprimality....

Detailed mathematical approach

Problem Summary

The implementations count integer triples $(a,b,c)$ with

$$F(N)=\#\left\{(a,b,c)\in\mathbb{Z}_{>0}^3:\ a\le N,\ a<b<c<2a,\ c^2=\frac{4a^2b^2}{5b^2-4a^2}\right\}.$$

A direct search over $(a,b,c)$ is far too slow for the actual limit. The key reduction is to convert the Diophantine condition into a rational point on a conic, classify primitive solutions once, and then count every scaled multiple in one floor division.

Mathematical Approach

Step 1: Rewrite the Constraint as a Conic

Start from

$$c^2=\frac{4a^2b^2}{5b^2-4a^2}.$$

After rearranging and dividing by $b^2c^2$, we obtain

$$\left(\frac{2a}{b}\right)^2+\left(\frac{2a}{c}\right)^2=5.$$

So if we define

$$X=\frac{2a}{b},\qquad Y=\frac{2a}{c},$$

then $(X,Y)$ is a rational point on the conic $X^2+Y^2=5$. The ordering $a<b<c<2a$ translates into

$$1<Y<X<2.$$

This is the geometric core of the fast solution: the original triple-counting problem becomes a count of rational points in one canonical region of the conic.

Step 2: Pass to Primitive Integer Data

Write the rational point in lowest terms as

$$X=\frac{U}{W},\qquad Y=\frac{V}{W},\qquad \gcd(U,V,W)=1.$$

Then

$$U^2+V^2=5W^2.$$

Because $1<Y<X<2$, the reduced integers satisfy

$$W<V<U<2W.$$

Any common divisor of two of $U,V,W$ would also divide the third, so primitiveness implies pairwise coprimality. Also $U$ and $V$ cannot both be odd, hence exactly one of them is even and $UV/2$ is an integer.

Step 3: Recover All Integer Triples from One Primitive Solution

From $X=2a/b=U/W$ and $Y=2a/c=V/W$, we get

$$b=\frac{2aW}{U},\qquad c=\frac{2aW}{V}.$$

Since $\gcd(U,W)=\gcd(V,W)=1$, integrality forces $U\mid 2a$ and $V\mid 2a$. Because $\gcd(U,V)=1$, we must have

$$UV\mid 2a.$$

Therefore every solution in this primitive family is obtained by writing

$$a=k\frac{UV}{2},\qquad b=kVW,\qquad c=kUW,\qquad k\ge 1.$$

The primitive family generator is thus

$$a_0=\frac{UV}{2},\qquad b_0=VW,\qquad c_0=UW,$$

and its contribution to $F(N)$ is exactly

$$\left\lfloor\frac{N}{a_0}\right\rfloor=\left\lfloor\frac{2N}{UV}\right\rfloor.$$

Step 4: Parametrize Primitive Conic Points

The implementations enumerate coprime pairs $m>n>m/2$ and define

$$X_0=m^2-4mn-n^2,\qquad Y_0=2(m^2+mn-n^2),\qquad W_0=m^2+n^2.$$

A direct expansion shows

$$X_0^2+Y_0^2=5W_0^2.$$

So $(|X_0|,Y_0,W_0)$ is an integer solution of the same quadratic form. Let

$$g=\gcd(|X_0|,Y_0,W_0).$$

The reduced primitive data is

$$U=\frac{\max(|X_0|,Y_0)}{g},\qquad V=\frac{\min(|X_0|,Y_0)}{g},\qquad W=\frac{W_0}{g}.$$

Cases with $5\mid g$ are discarded by the implementations because they are non-canonical $5$-lifts of smaller primitive families and would otherwise be counted twice.

Step 5: Why the Canonical Inequalities Are Automatic

The condition $n>m/2$ is exactly what puts the parameterization into the correct region. Indeed,

$$|X_0|-W_0=2m(2n-m)>0,$$

$$Y_0-W_0=(m-n)(m+3n)>0,$$

$$2W_0-|X_0|=3m^2+n^2-4mn>0,$$

$$2W_0-Y_0=2n(2n-m)>0.$$

Hence both $|X_0|$ and $Y_0$ lie strictly between $W_0$ and $2W_0$. After dividing by $g$ and sorting the two numerators, we obtain

$$W<V<U<2W,$$

which is exactly the reduced form corresponding to $a<b<c<2a$.

Step 6: Why the Common Divisor Is So Small

Suppose an odd prime $p$ divides all of $|X_0|,Y_0,W_0$. Then $p$ divides

$$2W_0-Y_0=2n(2n-m),\qquad W_0+X_0=2m(m-2n).$$

Because $p$ is odd and $\gcd(m,n)=1$, the only possible shared cause is

$$m\equiv 2n \pmod p.$$

Substituting this into $W_0=m^2+n^2$ gives $5n^2\equiv 0\pmod p$, so $p=5$. Thus the only odd common prime factor is $5$. For the factor $2$, if $m$ and $n$ are both odd then $|X_0|,Y_0,W_0$ are even but not divisible by $4$; otherwise at least one of them is odd. Therefore

$$g\in\{1,2,5,10\},$$

and once the duplicate $5$-lifts are skipped, the surviving cases have $g=1$ or $g=2$.

Worked Example

Take $(m,n)=(3,2)$. Then

$$X_0=3^2-4\cdot 3\cdot 2-2^2=-19,\qquad Y_0=2(3^2+3\cdot 2-2^2)=22,\qquad W_0=3^2+2^2=13.$$

Here $g=1$, so after ordering we get

$$U=22,\qquad V=19,\qquad W=13.$$

The primitive triple generated by this family is

$$a_0=\frac{22\cdot 19}{2}=209,\qquad b_0=19\cdot 13=247,\qquad c_0=22\cdot 13=286.$$

Thus the whole family is

$$ (a,b,c)=k(209,247,286),\qquad k\ge 1,$$

and it contributes $\lfloor N/209\rfloor$ solutions. In particular, this already explains why $F(100)=0$: the first primitive family starts above $100$.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematics. They first compute a safe upper bound for $m$, then iterate every coprime pair $m>n>m/2$. For each pair they build the quadratic forms above, divide out the common factor, reject the duplicate $5$-lift cases, sort the two reduced numerators into the canonical order, and finally add

$$\left\lfloor\frac{N}{UV/2}\right\rfloor.$$

The C++ and Java implementations parallelize the outer loop over $m$, while the Python implementation uses the same arithmetic sequentially.

Complexity Analysis

Let $M$ be the largest $m$ that can still contribute. Using $g\le 2$ after the $5$-lift filter, we have

$$a_0=\frac{UV}{2}=\frac{|X_0|Y_0}{2g^2}\ge \frac{|X_0|Y_0}{8}.$$

Write $r=n/m$, so $r\in(1/2,1)$. Then

$$|X_0|=m^2(4r+r^2-1),\qquad Y_0=2m^2(1+r-r^2),$$

hence

$$a_0\ge \frac{m^4}{4}(4r+r^2-1)(1+r-r^2).$$

The factor on the right is increasing on $(1/2,1)$, so its minimum occurs at $r=1/2$ and equals $25/16$. Therefore

$$a_0\ge \frac{25}{64}m^4>\frac14 m^4.$$

So no value $m>(4N)^{1/4}$ can contribute, which matches the fourth-root cutoff used by the implementations. The enumeration therefore visits $O(M^2)=O(N^{1/2})$ candidate pairs, with a constant amount of arithmetic and one gcd test per pair. Memory usage is $O(1)$ apart from thread-local accumulators.

References

Problem page: https://projecteuler.net/problem=404
Rational parametrization of conics: Wikipedia — Conic section
Euclidean algorithm and gcd arguments: Wikipedia — Euclidean algorithm
Primitive quadratic-form parametrizations and related ideas: Wikipedia — Pythagorean triple

Problem 404 source code

C++

#include <algorithm>
#include <array>
#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 i128 = __int128_t;
using u128 = __uint128_t;

struct Options {
    u64 limit = 100000000000000000ULL;
    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 overflow");
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const std::string& prefix,
                                 unsigned& 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;
        }
        parsed = parsed * 10ULL + static_cast<u64>(ch - '0');
        if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
            throw std::overflow_error("unsigned 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, "--n=", options.limit)) {
            continue;
        }
        if (parse_unsigned_after_prefix(arg, "--threads=", options.requested_threads)) {
            continue;
        }

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

    return true;
}

u64 abs_i128(const i128 value) {
    return static_cast<u64>(value < 0 ? -value : value);
}

u64 isqrt_u128(const u128 x) {
    if (x == 0U) {
        return 0U;
    }

    long double approx = std::sqrt(static_cast<long double>(x));
    u64 r = static_cast<u64>(approx);

    while (static_cast<u128>(r + 1U) * static_cast<u128>(r + 1U) <= x) {
        ++r;
    }
    while (static_cast<u128>(r) * static_cast<u128>(r) > x) {
        --r;
    }

    return r;
}

u64 fourth_root_bound(const u64 n) {
    const u128 target = static_cast<u128>(4U) * static_cast<u128>(n);
    long double approx = std::pow(static_cast<long double>(target), 0.25L);
    u64 m = static_cast<u64>(approx) + 4U;

    auto pow4 = [](const u64 x) -> u128 {
        const u128 y = static_cast<u128>(x);
        return y * y * y * y;
    };

    while (m > 0U && pow4(m) > target) {
        --m;
    }
    while (pow4(m + 1U) <= target) {
        ++m;
    }

    return m + 2U;
}

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

    unsigned hw = std::thread::hardware_concurrency();
    if (hw == 0U) {
        hw = 4U;
    }
    return hw;
}

u64 count_triplets_fast(const u64 limit, unsigned thread_count) {
    if (limit == 0U) {
        return 0U;
    }

    const u64 m_max = fourth_root_bound(limit);
    if (m_max < 2U) {
        return 0U;
    }

    thread_count = std::max(1U, std::min(thread_count, static_cast<unsigned>(m_max - 1U)));

    std::vector<u128> partial(thread_count, 0U);
    std::vector<std::thread> workers;
    workers.reserve(thread_count);

    for (unsigned tid = 0U; tid < thread_count; ++tid) {
        workers.emplace_back([=, &partial]() {
            u128 local = 0U;

            for (u64 m = 2U + static_cast<u64>(tid); m <= m_max; m += static_cast<u64>(thread_count)) {
                const u64 mm = m * m;
                const u64 n_start = m / 2U + 1U;

                for (u64 n = n_start; n < m; ++n) {
                    if (std::gcd(m, n) != 1U) {
                        continue;
                    }

                    const u64 nn = n * n;

                    // Parameterization of rational points on u^2+v^2=5.
                    const i128 x_raw = static_cast<i128>(mm) - static_cast<i128>(4U) *
                                       static_cast<i128>(m) * static_cast<i128>(n) -
                                       static_cast<i128>(nn);
                    const u64 x_abs = abs_i128(x_raw);
                    const u64 y_abs = 2U * (mm + m * n - nn);
                    const u64 w = mm + nn;

                    u64 g = std::gcd(x_abs, y_abs);
                    g = std::gcd(g, w);

                    // g divisible by 5 corresponds to the 5x-lift duplicate representation.
                    if (g % 5U == 0U) {
                        continue;
                    }

                    u64 u = x_abs / g;
                    u64 v = y_abs / g;
                    if (u < v) {
                        std::swap(u, v);
                    }

                    // n > m/2 guarantees the canonical region W < V < U < 2W.
                    const u128 a_primitive = (static_cast<u128>(u) * static_cast<u128>(v)) / 2U;
                    if (a_primitive > static_cast<u128>(limit)) {
                        continue;
                    }

                    local += static_cast<u128>(limit) / a_primitive;
                }
            }

            partial[tid] = local;
        });
    }

    for (std::thread& worker : workers) {
        worker.join();
    }

    u128 total = 0U;
    for (const u128 value : partial) {
        total += value;
    }

    return static_cast<u64>(total);
}

u64 count_triplets_bruteforce(const u64 limit) {
    u64 total = 0U;

    for (u64 a = 1U; a <= limit; ++a) {
        const u128 a2 = static_cast<u128>(a) * static_cast<u128>(a);
        const u128 b2_max = (static_cast<u128>(8U) * a2 - 1U) / 5U;
        const u64 b_max = isqrt_u128(b2_max);

        for (u64 b = a + 1U; b <= b_max; ++b) {
            const u128 b2 = static_cast<u128>(b) * static_cast<u128>(b);
            const u128 den = static_cast<u128>(5U) * b2 - static_cast<u128>(4U) * a2;
            if (den == 0U) {
                continue;
            }

            const u128 num = static_cast<u128>(4U) * a2 * b2;
            if (num % den != 0U) {
                continue;
            }

            const u128 c2 = num / den;
            const u64 c = isqrt_u128(c2);
            if (static_cast<u128>(c) * static_cast<u128>(c) != c2) {
                continue;
            }
            if (!(b < c && c < 2U * a)) {
                continue;
            }

            ++total;
        }
    }

    return total;
}

void run_checkpoints() {
    struct Checkpoint {
        u64 n;
        u64 expected;
    };

    constexpr std::array<Checkpoint, 2U> known{{
        {100U, 0U},
        {10000U, 106U},
    }};

    for (const Checkpoint& checkpoint : known) {
        const u64 got = count_triplets_fast(checkpoint.n, 1U);
        if (got != checkpoint.expected) {
            throw std::runtime_error("Checkpoint failed at N=" + std::to_string(checkpoint.n));
        }
    }

    constexpr u64 brute_n = 5000U;
    const u64 fast_small = count_triplets_fast(brute_n, 1U);
    const u64 brute_small = count_triplets_bruteforce(brute_n);
    if (fast_small != brute_small) {
        throw std::runtime_error("Fast/bruteforce mismatch at N=" + std::to_string(brute_n));
    }
}

}  // namespace

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

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

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

Python

import math

def solve():
    limit = 100000000000000000
    # Fourth root bound of 4*limit
    def fourth_root_bound(n):
        target = 4 * n
        m = int(target ** 0.25) + 4
        while m > 0 and m**4 > target: m -= 1
        while (m+1)**4 <= target: m += 1
        return m + 2

    m_max = fourth_root_bound(limit)
    total = 0
    for m in range(2, m_max + 1):
        mm = m * m
        n_start = m // 2 + 1
        for n in range(n_start, m):
            if math.gcd(m, n) != 1: continue
            nn = n * n
            x_raw = mm - 4 * m * n - nn
            x_abs = abs(x_raw)
            y_abs = 2 * (mm + m * n - nn)
            w = mm + nn
            g = math.gcd(x_abs, math.gcd(y_abs, w))
            if g % 5 == 0: continue
            u = x_abs // g
            v = y_abs // g
            if u < v: u, v = v, u
            a_prim = (u * v) // 2
            if a_prim > limit: continue
            total += limit // a_prim
    return str(total)

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

Java

import java.util.stream.LongStream;

public class Euler404 {
    static long gcd(long a, long b) {
        while (b != 0) {
            long temp = b;
            b = a % b;
            a = temp;
        }
        return Math.abs(a);
    }

    static long fourthRootBound(long limit) {
        long target = 4 * limit;
        long m = (long) Math.pow(target, 0.25) + 4;

        while (m > 0) {
            long p = m * m;
            if (p * p <= target && p * p > 0) {
                break;
            }
            m--;
        }
        while (true) {
            long p = (m + 1) * (m + 1);
            if (p * p <= target && p * p > 0) {
                m++;
            } else {
                break;
            }
        }
        return m + 2;
    }

    static String solve() {
        long limit = 100000000000000000L;
        long mMax = fourthRootBound(limit);

        long total = LongStream.rangeClosed(2, mMax)
                .parallel()
                .map(m -> {
                    long local = 0;
                    long mm = m * m;
                    long nStart = m / 2 + 1;

                    for (long n = nStart; n < m; n++) {
                        if (gcd(m, n) != 1)
                            continue;

                        long nn = n * n;
                        long xRaw = mm - 4 * m * n - nn;
                        long xAbs = Math.abs(xRaw);
                        long yAbs = 2 * (mm + m * n - nn);
                        long w = mm + nn;

                        long g = gcd(xAbs, yAbs);
                        g = gcd(g, w);

                        if (g % 5 == 0)
                            continue;

                        long u = xAbs / g;
                        long v = yAbs / g;
                        if (u < v) {
                            long tmp = u;
                            u = v;
                            v = tmp;
                        }

                        long aPrimitive = (u * v) / 2;
                        if (aPrimitive > limit)
                            continue;

                        local += limit / aPrimitive;
                    }
                    return local;
                })
                .sum();

        return Long.toString(total);
    }

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