Problem 180: Rational Zeros of a Function of Three Variables

Project Euler Problem 180 Solution

EulerSolve provides an optimized solution for Project Euler Problem 180, Rational Zeros of a Function of Three Variables, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a positive integer $k$, the admissible numbers are the reduced proper fractions $$S_k=\left\{\frac{a}{b}:0<a<b\le k,\ \gcd(a,b)=1\right\}.$$ Problem 180 fixes $k=35$ and asks for all triples $(x,y,z)\in S_k^3$ for which the problem's function $f_n(x,y,z)$ vanishes for at least one integer $n$. For every such triple we form $s=x+y+z$, keep only distinct rational values of $s$, add them all together, reduce the total to lowest terms $u/v$, and report $u+v$. The direct interpretation suggests searching all triples and all integers $n$. The implementations do something much sharper: they factor the defining expression, prove that only four exponents can matter, and then scan pairs $(x,y)$ instead of triples $(x,y,z)$. Mathematical Approach The whole solution rests on one algebraic simplification. Once the original zero-condition is rewritten correctly, the search space collapses from an open-ended three-variable problem to four explicit formulas for $z$. The finite state space $S_k$ Every element of $S_k$ is a positive rational smaller than 1, already stored in lowest terms. For each denominator $b$, exactly $\varphi(b)$ numerators are coprime to $b$, so $$|S_k|=\sum_{b=2}^{k}\varphi(b).$$ At the target order $k=35$, this gives $383$ admissible fractions....

Detailed mathematical approach

Problem Summary

For a positive integer $k$, the admissible numbers are the reduced proper fractions

$$S_k=\left\{\frac{a}{b}:0<a<b\le k,\ \gcd(a,b)=1\right\}.$$

Problem 180 fixes $k=35$ and asks for all triples $(x,y,z)\in S_k^3$ for which the problem's function $f_n(x,y,z)$ vanishes for at least one integer $n$. For every such triple we form $s=x+y+z$, keep only distinct rational values of $s$, add them all together, reduce the total to lowest terms $u/v$, and report $u+v$.

The direct interpretation suggests searching all triples and all integers $n$. The implementations do something much sharper: they factor the defining expression, prove that only four exponents can matter, and then scan pairs $(x,y)$ instead of triples $(x,y,z)$.

Mathematical Approach

The whole solution rests on one algebraic simplification. Once the original zero-condition is rewritten correctly, the search space collapses from an open-ended three-variable problem to four explicit formulas for $z$.

The finite state space $S_k$

Every element of $S_k$ is a positive rational smaller than 1, already stored in lowest terms. For each denominator $b$, exactly $\varphi(b)$ numerators are coprime to $b$, so

$$|S_k|=\sum_{b=2}^{k}\varphi(b).$$

At the target order $k=35$, this gives $383$ admissible fractions. A naive triple search would therefore examine $383^3$ candidates before even considering different exponents, so the algebraic reduction is the decisive step.

Factor the zero-condition

The implementations encode the original function as

$$\begin{aligned} f_n(x,y,z)={}&x^{n+1}+y^{n+1}-z^{n+1}\\ &+(xy+yz+zx)(x^{n-1}+y^{n-1}-z^{n-1})\\ &-xyz(x^{n-2}+y^{n-2}-z^{n-2}). \end{aligned}$$

If the second and third lines are expanded, the mixed terms $yz\,x^{n-1}$, $xz\,y^{n-1}$, and $xy\,z^{n-1}$ cancel, and everything that remains factors cleanly:

$$f_n(x,y,z)=(x+y+z)(x^n+y^n-z^n).$$

Because $x$, $y$, and $z$ are positive rationals, we always have $x+y+z>0$. Therefore

$$f_n(x,y,z)=0 \quad\Longleftrightarrow\quad x^n+y^n=z^n.$$

This is the key invariant behind all three implementations.

Why only four exponents survive

If $n=0$, the reduced equation becomes $1+1=1$, which is impossible. If $n>2$, any positive rational solution of $x^n+y^n=z^n$ can be cleared of denominators to produce a positive integer solution, contradicting Fermat's Last Theorem. If $n<-2$, write $m=-n>2$; then

$$x^n+y^n=z^n \quad\Longleftrightarrow\quad \frac{1}{x^m}+\frac{1}{y^m}=\frac{1}{z^m},$$

which is the same kind of forbidden equation after replacing $x,y,z$ by their reciprocals.

So the only exponents that can possibly contribute are

$$n\in\{-2,-1,1,2\}.$$

The infinite search over all integers has collapsed to four explicit cases.

Turn each pair $(x,y)$ into at most four candidates for $z$

Once $x$ and $y$ are fixed, each admissible exponent gives a concrete formula:

$$n=1:\quad z=x+y,$$

$$n=-1:\quad z=\frac{xy}{x+y},$$

$$n=2:\quad z=\sqrt{x^2+y^2},$$

$$n=-2:\quad z=\frac{1}{\sqrt{1/x^2+1/y^2}}=\frac{xy}{\sqrt{x^2+y^2}}.$$

This is why the optimized solver loops over ordered pairs $(x,y)$ rather than triples. For the square-root branches, the implementations compute the reduced fraction for $z^2$ and accept it only when both numerator and denominator are perfect squares. After reduction, a membership test against $S_k$ simultaneously checks that $0<z<1$ and that the reduced denominator is at most $k$.

Worked examples and distinct sums

The arithmetic is easiest to see on small examples. In the $n=-1$ branch, taking $x=\tfrac12$ and $y=\tfrac13$ gives

$$z=\frac{xy}{x+y}=\frac{\tfrac12\cdot\tfrac13}{\tfrac12+\tfrac13}=\frac{1/6}{5/6}=\frac15,$$

so $s=\tfrac12+\tfrac13+\tfrac15=\tfrac{31}{30}$. In the $n=2$ branch,

$$\left(\frac{3}{10}\right)^2+\left(\frac{2}{5}\right)^2=\frac{9}{100}+\frac{16}{100}=\frac{25}{100}=\left(\frac12\right)^2,$$

so $(x,y,z)=\left(\tfrac{3}{10},\tfrac25,\tfrac12\right)$ is valid. For $n=-2$, choosing $x=\tfrac13$ and $y=\tfrac14$ gives

$$\frac{1}{x^2}+\frac{1}{y^2}=9+16=25=\frac{1}{(1/5)^2},$$

hence $z=\tfrac15$.

Different triples, and even different exponents, can lead to the same final value $s=x+y+z$. The set of sums is therefore deduplicated at the rational-number level, not at the triple level.

How the Code Works

Enumerating the allowed fractions

The C++, Python, and Java implementations first generate every reduced fraction in $S_{35}$. These fractions are also inserted into a hash-based membership structure so that a candidate $z$ can be accepted or rejected in constant expected time.

Exact candidate generation

The main loop runs through every ordered pair $(x,y)\in S_k^2$. For each pair, the implementation constructs the four candidates listed above. All arithmetic is exact: the Python version relies on built-in rational arithmetic, while the C++ and Java versions keep reduced numerator-denominator pairs and normalize after each addition, multiplication, division, or reciprocal. The $n=\pm2$ branches use integer square-root tests on the reduced numerator and denominator to detect whether the square root stays rational.

Membership filtering and deduplication

If a candidate $z$ belongs to $S_k$, the implementation forms $s=x+y+z$ and inserts that reduced fraction into a set of distinct sums. The scan does not try to avoid the symmetry between $(x,y)$ and $(y,x)$; it simply allows duplicates to arise and lets the set remove them automatically.

Exact final accumulation

After the pair scan finishes, the distinct $s$-values are added exactly and reduced to a single fraction $u/v$. The required output is then $u+v$. The C++ implementation additionally contains small-order validation checkpoints: it compares the factored equation with the original definition, checks the optimized pair scan against brute-force triple enumeration on small $k$, and verifies consistency between single-threaded and multithreaded execution.

Complexity Analysis

Let

$$M=|S_k|=\sum_{b=2}^{k}\varphi(b).$$

Building the admissible fraction list and the membership set costs $O(M)$. The main solver then examines every ordered pair $(x,y)$, so the dominant work is $O(M^2)$. Each pair performs four candidate constructions, a constant number of gcd reductions, a few hash lookups, and at most two perfect-square checks.

The memory usage is $O(M+U)$, where $U$ is the number of distinct sums retained. For the target $k=35$, $M=383$, so the optimized search examines only $383^2=146{,}689$ ordered pairs instead of $383^3$ triples. The brute-force $O(M^3)$ method appears only in the C++ validation checkpoints for tiny orders.

Footnotes and References

Problem page: Project Euler 180
Rational number: Wikipedia - Rational number
Euler's totient function: Wikipedia - Euler's totient function
Fermat's Last Theorem: Wikipedia - Fermat's Last Theorem
Pythagorean triple: Wikipedia - Pythagorean triple

Problem 180 source code

C++

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

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using boost::multiprecision::cpp_int;
using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr u64 kDefaultOrder = 35;

struct Options {
    u64 order = kDefaultOrder;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

struct Fraction {
    u64 num = 0;
    u64 den = 1;

    bool operator==(const Fraction& other) const noexcept {
        return num == other.num && den == other.den;
    }
};

struct FractionHash {
    std::size_t operator()(const Fraction& f) const noexcept {
        const u64 x = f.num * 0x9E3779B185EBCA87ULL;
        const u64 y = f.den + 0xC2B2AE3D27D4EB4FULL;
        return static_cast<std::size_t>(x ^ (y + (x << 6) + (x >> 2)));
    }
};

u128 gcd_u128(u128 a, u128 b) {
    while (b != 0) {
        const u128 r = a % b;
        a = b;
        b = r;
    }
    return a;
}

u64 cast_u64_checked(u128 value) {
    if (value > static_cast<u128>(std::numeric_limits<u64>::max())) {
        std::cerr << "Overflow while converting to uint64_t.\n";
        std::exit(1);
    }
    return static_cast<u64>(value);
}

Fraction make_fraction(u64 num, u64 den) {
    if (den == 0) {
        std::cerr << "Attempted to build fraction with zero denominator.\n";
        std::exit(1);
    }
    const u64 g = std::gcd(num, den);
    return Fraction{num / g, den / g};
}

Fraction add_fraction(const Fraction& a, const Fraction& b) {
    const u64 g = std::gcd(a.den, b.den);
    u128 num = static_cast<u128>(a.num) * static_cast<u128>(b.den / g) +
               static_cast<u128>(b.num) * static_cast<u128>(a.den / g);
    u128 den = static_cast<u128>(a.den / g) * static_cast<u128>(b.den);
    const u128 h = gcd_u128(num, den);
    num /= h;
    den /= h;
    return Fraction{cast_u64_checked(num), cast_u64_checked(den)};
}

Fraction mul_fraction(const Fraction& a, const Fraction& b) {
    u128 num = static_cast<u128>(a.num) * static_cast<u128>(b.num);
    u128 den = static_cast<u128>(a.den) * static_cast<u128>(b.den);
    const u128 g = gcd_u128(num, den);
    num /= g;
    den /= g;
    return Fraction{cast_u64_checked(num), cast_u64_checked(den)};
}

Fraction reciprocal(const Fraction& f) {
    if (f.num == 0) {
        std::cerr << "Attempted reciprocal of zero.\n";
        std::exit(1);
    }
    return Fraction{f.den, f.num};
}

Fraction div_fraction(const Fraction& a, const Fraction& b) {
    return mul_fraction(a, reciprocal(b));
}

Fraction square_fraction(const Fraction& f) {
    return mul_fraction(f, f);
}

bool is_square_u64(u64 value, u64& root) {
    const long double approx = std::sqrt(static_cast<long double>(value));
    u64 r = static_cast<u64>(approx);
    while (static_cast<u128>(r) * static_cast<u128>(r) < value) {
        ++r;
    }
    while (static_cast<u128>(r) * static_cast<u128>(r) > value) {
        --r;
    }
    if (static_cast<u128>(r) * static_cast<u128>(r) == value) {
        root = r;
        return true;
    }
    return false;
}

bool sqrt_fraction(const Fraction& f, Fraction& root) {
    u64 nroot = 0;
    u64 droot = 0;
    if (!is_square_u64(f.num, nroot) || !is_square_u64(f.den, droot)) {
        return false;
    }
    root = make_fraction(nroot, droot);
    return true;
}

std::vector<Fraction> build_fraction_list(u64 order) {
    std::vector<Fraction> fractions;
    for (u64 den = 2; den <= order; ++den) {
        for (u64 num = 1; num < den; ++num) {
            if (std::gcd(num, den) == 1) {
                fractions.push_back(Fraction{num, den});
            }
        }
    }
    return fractions;
}

unsigned decide_thread_count(bool allow_multithreading,
                             unsigned requested_threads,
                             std::size_t x_count) {
    if (!allow_multithreading) {
        return 1U;
    }
    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }
    if (x_count < 128U) {
        threads = 1U;
    }
    if (threads > x_count) {
        threads = static_cast<unsigned>(x_count);
    }
    return std::max(1U, threads);
}

void add_sum_if_member(const Fraction& x,
                       const Fraction& y,
                       const Fraction& z,
                       const std::unordered_set<Fraction, FractionHash>& allowed,
                       std::unordered_set<Fraction, FractionHash>& sums) {
    if (allowed.find(z) == allowed.end()) {
        return;
    }
    const Fraction s = add_fraction(add_fraction(x, y), z);
    sums.insert(s);
}

void process_pair(const Fraction& x,
                  const Fraction& y,
                  const std::unordered_set<Fraction, FractionHash>& allowed,
                  std::unordered_set<Fraction, FractionHash>& sums) {
    // For positive rationals, f_n(x,y,z)=0 is equivalent to x^n + y^n = z^n.
    // Non-trivial rational solutions only occur for n in {-2,-1,1,2}.

    // n = 1: x + y = z.
    const Fraction z1 = add_fraction(x, y);
    add_sum_if_member(x, y, z1, allowed, sums);

    // n = -1: 1/x + 1/y = 1/z  =>  z = xy/(x+y).
    const Fraction zminus1 = div_fraction(mul_fraction(x, y), z1);
    add_sum_if_member(x, y, zminus1, allowed, sums);

    // n = 2: x^2 + y^2 = z^2.
    const Fraction z2_sq = add_fraction(square_fraction(x), square_fraction(y));
    Fraction z2{};
    if (sqrt_fraction(z2_sq, z2)) {
        add_sum_if_member(x, y, z2, allowed, sums);
    }

    // n = -2: 1/x^2 + 1/y^2 = 1/z^2  =>  z^2 = 1 / (1/x^2 + 1/y^2).
    const Fraction inv_x_sq = square_fraction(reciprocal(x));
    const Fraction inv_y_sq = square_fraction(reciprocal(y));
    const Fraction zminus2_sq = reciprocal(add_fraction(inv_x_sq, inv_y_sq));
    Fraction zminus2{};
    if (sqrt_fraction(zminus2_sq, zminus2)) {
        add_sum_if_member(x, y, zminus2, allowed, sums);
    }
}

std::unordered_set<Fraction, FractionHash> collect_sums_fast(
    u64 order,
    bool allow_multithreading,
    unsigned requested_threads) {
    const std::vector<Fraction> fractions = build_fraction_list(order);
    const std::unordered_set<Fraction, FractionHash> allowed(fractions.begin(), fractions.end());
    const std::size_t n = fractions.size();

    const unsigned threads = decide_thread_count(allow_multithreading, requested_threads, n);
    if (threads <= 1U) {
        std::unordered_set<Fraction, FractionHash> sums;
        sums.reserve(n * 4);
        for (const Fraction& x : fractions) {
            for (const Fraction& y : fractions) {
                process_pair(x, y, allowed, sums);
            }
        }
        return sums;
    }

    std::vector<std::unordered_set<Fraction, FractionHash>> locals(threads);
    for (auto& local : locals) {
        local.reserve((n * 4) / threads + 64U);
    }

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

    for (unsigned tid = 0; tid < threads; ++tid) {
        const std::size_t begin = (n * tid) / threads;
        const std::size_t end = (n * (tid + 1U)) / threads;

        workers.emplace_back([&, tid, begin, end]() {
            auto& local = locals[tid];
            for (std::size_t i = begin; i < end; ++i) {
                const Fraction& x = fractions[i];
                for (const Fraction& y : fractions) {
                    process_pair(x, y, allowed, local);
                }
            }
        });
    }

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

    std::unordered_set<Fraction, FractionHash> sums;
    for (const auto& local : locals) {
        sums.insert(local.begin(), local.end());
    }
    return sums;
}

bool satisfies_order_from_derived_equations(const Fraction& x,
                                            const Fraction& y,
                                            const Fraction& z) {
    if (add_fraction(x, y) == z) {
        return true;
    }

    if (div_fraction(mul_fraction(x, y), add_fraction(x, y)) == z) {
        return true;
    }

    if (add_fraction(square_fraction(x), square_fraction(y)) == square_fraction(z)) {
        return true;
    }

    const Fraction inv_x_sq = square_fraction(reciprocal(x));
    const Fraction inv_y_sq = square_fraction(reciprocal(y));
    const Fraction inv_z_sq = square_fraction(reciprocal(z));
    if (add_fraction(inv_x_sq, inv_y_sq) == inv_z_sq) {
        return true;
    }

    return false;
}

struct RationalBig {
    cpp_int num = 0;
    cpp_int den = 1;

    RationalBig() = default;
    RationalBig(cpp_int n, cpp_int d) : num(std::move(n)), den(std::move(d)) {
        normalize();
    }

    static cpp_int abs_value(cpp_int v) {
        return (v < 0) ? -v : v;
    }

    static cpp_int gcd(cpp_int a, cpp_int b) {
        a = abs_value(a);
        b = abs_value(b);
        while (b != 0) {
            cpp_int r = a % b;
            a = b;
            b = r;
        }
        return a;
    }

    void normalize() {
        if (den == 0) {
            std::cerr << "Big rational with zero denominator.\n";
            std::exit(1);
        }
        if (den < 0) {
            den = -den;
            num = -num;
        }
        const cpp_int g = gcd(num, den);
        num /= g;
        den /= g;
    }
};

RationalBig operator+(const RationalBig& a, const RationalBig& b) {
    return RationalBig(a.num * b.den + b.num * a.den, a.den * b.den);
}

RationalBig operator-(const RationalBig& a, const RationalBig& b) {
    return RationalBig(a.num * b.den - b.num * a.den, a.den * b.den);
}

RationalBig operator*(const RationalBig& a, const RationalBig& b) {
    return RationalBig(a.num * b.num, a.den * b.den);
}

RationalBig operator/(const RationalBig& a, const RationalBig& b) {
    return RationalBig(a.num * b.den, a.den * b.num);
}

RationalBig pow_rational(RationalBig base, int exponent) {
    if (exponent == 0) {
        return RationalBig(1, 1);
    }
    if (exponent < 0) {
        base = RationalBig(base.den, base.num);
        exponent = -exponent;
    }

    RationalBig result(1, 1);
    int e = exponent;
    while (e > 0) {
        if ((e & 1) != 0) {
            result = result * base;
        }
        base = base * base;
        e >>= 1;
    }
    return result;
}

RationalBig to_big(const Fraction& f) {
    return RationalBig(cpp_int(f.num), cpp_int(f.den));
}

bool satisfies_order_from_original_formula(const Fraction& x,
                                           const Fraction& y,
                                           const Fraction& z) {
    static const int exponents[] = {-2, -1, 1, 2};

    const RationalBig bx = to_big(x);
    const RationalBig by = to_big(y);
    const RationalBig bz = to_big(z);

    for (const int n : exponents) {
        const RationalBig f1 = pow_rational(bx, n + 1) + pow_rational(by, n + 1) -
                               pow_rational(bz, n + 1);

        const RationalBig symmetric = bx * by + by * bz + bz * bx;
        const RationalBig f2 = symmetric * (pow_rational(bx, n - 1) +
                                            pow_rational(by, n - 1) -
                                            pow_rational(bz, n - 1));

        const RationalBig f3 = bx * by * bz * (pow_rational(bx, n - 2) +
                                               pow_rational(by, n - 2) -
                                               pow_rational(bz, n - 2));

        const RationalBig fn = f1 + f2 - f3;
        if (fn.num == 0) {
            return true;
        }
    }
    return false;
}

std::unordered_set<Fraction, FractionHash> collect_sums_bruteforce_derived(u64 order) {
    const std::vector<Fraction> fractions = build_fraction_list(order);
    std::unordered_set<Fraction, FractionHash> sums;
    sums.reserve(fractions.size() * 8);

    for (const Fraction& x : fractions) {
        for (const Fraction& y : fractions) {
            for (const Fraction& z : fractions) {
                if (satisfies_order_from_derived_equations(x, y, z)) {
                    sums.insert(add_fraction(add_fraction(x, y), z));
                }
            }
        }
    }
    return sums;
}

std::unordered_set<Fraction, FractionHash> collect_sums_bruteforce_original(u64 order) {
    const std::vector<Fraction> fractions = build_fraction_list(order);
    std::unordered_set<Fraction, FractionHash> sums;
    sums.reserve(fractions.size() * 8);

    for (const Fraction& x : fractions) {
        for (const Fraction& y : fractions) {
            for (const Fraction& z : fractions) {
                if (satisfies_order_from_original_formula(x, y, z)) {
                    sums.insert(add_fraction(add_fraction(x, y), z));
                }
            }
        }
    }
    return sums;
}

bool set_equal(const std::unordered_set<Fraction, FractionHash>& a,
               const std::unordered_set<Fraction, FractionHash>& b) {
    if (a.size() != b.size()) {
        return false;
    }
    for (const Fraction& value : a) {
        if (b.find(value) == b.end()) {
            return false;
        }
    }
    return true;
}

std::pair<cpp_int, cpp_int> sum_distinct_s_values(
    const std::unordered_set<Fraction, FractionHash>& sums) {
    cpp_int num = 0;
    cpp_int den = 1;

    for (const Fraction& f : sums) {
        num = num * f.den + den * f.num;
        den *= f.den;
        const cpp_int g = RationalBig::gcd(num, den);
        num /= g;
        den /= g;
    }
    return {num, den};
}

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

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

    u64 parsed = 0;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }
    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(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 (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u64 parsed_order = 0;
        if (parse_u64_after_prefix(arg, "--order=", parsed_order)) {
            options.order = parsed_order;
            continue;
        }

        unsigned parsed_threads = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_threads)) {
            options.requested_threads = parsed_threads;
            continue;
        }

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

    if (options.order < 2) {
        std::cerr << "--order must be at least 2.\n";
        return false;
    }

    return true;
}

bool run_validation_checkpoints(const Options& options) {
    // Checkpoint 1: direct original-formula brute force vs derived equations.
    {
        const u64 k = 7;
        const auto original = collect_sums_bruteforce_original(k);
        const auto derived = collect_sums_bruteforce_derived(k);
        if (!set_equal(original, derived)) {
            std::cerr << "Checkpoint failed: original formula and derived equations differ at "
                      << "order " << k << ".\n";
            return false;
        }
    }

    // Checkpoint 2: fast pair-scan vs brute-force triples.
    {
        const u64 k = 12;
        const auto brute = collect_sums_bruteforce_derived(k);
        const auto fast = collect_sums_fast(k, false, 1);
        if (!set_equal(brute, fast)) {
            std::cerr << "Checkpoint failed: fast solver mismatch at order " << k << ".\n";
            return false;
        }
    }

    // Checkpoint 3: single-thread vs multi-thread consistency on target order.
    {
        const auto single = collect_sums_fast(options.order, false, 1);
        const auto multi = collect_sums_fast(options.order,
                                             options.allow_multithreading,
                                             options.requested_threads);
        if (!set_equal(single, multi)) {
            std::cerr << "Checkpoint failed: thread-consistency mismatch at order "
                      << options.order << ".\n";
            return false;
        }
    }

    return true;
}

}  // namespace

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

    if (options.run_checkpoints) {
        if (!run_validation_checkpoints(options)) {
            return 1;
        }
    }

    const auto sums = collect_sums_fast(options.order,
                                        options.allow_multithreading,
                                        options.requested_threads);
    const auto [u, v] = sum_distinct_s_values(sums);
    std::cout << (u + v) << '\n';
    return 0;
}

Python

from fractions import Fraction

def solve():
    order = 35
    fracs = []
    for d in range(2, order + 1):
        for n in range(1, d):
            from math import gcd
            if gcd(n, d) == 1:
                fracs.append(Fraction(n, d))

    allowed = set(fracs)
    sums_set = set()

    for x in fracs:
        for y in fracs:
            # n=1: z = x+y
            z1 = x + y
            if z1 in allowed:
                sums_set.add(x + y + z1)
            # n=-1: z = xy/(x+y)
            s = x + y
            if s != 0:
                zm1 = x * y / s
                if zm1 in allowed:
                    sums_set.add(x + y + zm1)
            # n=2: z = sqrt(x^2+y^2)
            z2sq = x*x + y*y
            # Check if z2sq is a perfect rational square
            if z2sq > 0:
                num = z2sq.numerator
                den = z2sq.denominator
                from math import isqrt
                rn = isqrt(num)
                rd = isqrt(den)
                if rn*rn == num and rd*rd == den:
                    z2 = Fraction(rn, rd)
                    if z2 in allowed:
                        sums_set.add(x + y + z2)
            # n=-2: 1/z^2 = 1/x^2 + 1/y^2
            if x != 0 and y != 0:
                inv_sq = Fraction(1, x*x) + Fraction(1, y*y)
                if inv_sq > 0:
                    zm2sq = Fraction(1, inv_sq)
                    num = zm2sq.numerator
                    den = zm2sq.denominator
                    from math import isqrt
                    rn = isqrt(num)
                    rd = isqrt(den)
                    if rn*rn == num and rd*rd == den:
                        zm2 = Fraction(rn, rd)
                        if zm2 in allowed:
                            sums_set.add(x + y + zm2)

    total = sum(sums_set)
    return str(total.numerator + total.denominator)

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

Java

import java.math.BigInteger;
import java.util.*;

public class Euler180 {
    static BigInteger gcd(BigInteger a, BigInteger b) {
        return a.gcd(b);
    }

    static long[] makeFrac(long n, long d) {
        if (d == 0)
            return null;
        long g = gcd64(n, d);
        return new long[] { n / g, d / g };
    }

    static long gcd64(long a, long b) {
        a = Math.abs(a);
        b = Math.abs(b);
        while (b != 0) {
            long t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    static long[] addFrac(long[] a, long[] b) {
        long g = gcd64(a[1], b[1]);
        long num = a[0] * (b[1] / g) + b[0] * (a[1] / g);
        long den = (a[1] / g) * b[1];
        long h = gcd64(Math.abs(num), Math.abs(den));
        return new long[] { num / h, den / h };
    }

    static long[] mulFrac(long[] a, long[] b) {
        long num = a[0] * b[0], den = a[1] * b[1];
        long g = gcd64(Math.abs(num), Math.abs(den));
        return new long[] { num / g, den / g };
    }

    static long[] sqFrac(long[] f) {
        return mulFrac(f, f);
    }

    static long[] recipFrac(long[] f) {
        return new long[] { f[1], f[0] };
    }

    static long[] divFrac(long[] a, long[] b) {
        return mulFrac(a, recipFrac(b));
    }

    static boolean isSq64(long v) {
        if (v < 0)
            return false;
        long r = (long) Math.sqrt(v);
        return r * r == v || (r + 1) * (r + 1) == v;
    }

    static long[] sqrtFrac(long[] f) {
        long rn = (long) Math.sqrt(f[0]), rd = (long) Math.sqrt(f[1]);
        if (rn * rn == f[0] && rd * rd == f[1])
            return makeFrac(rn, rd);
        rn++;
        if (rn * rn == f[0] && rd * rd == f[1])
            return makeFrac(rn, rd);
        return null;
    }

    static long fracHash(long[] f) {
        return f[0] * 1000000007L + f[1];
    }

    public static void main(String[] args) {
        int order = 35;
        List<long[]> fracs = new ArrayList<>();
        for (long d = 2; d <= order; d++)
            for (long n = 1; n < d; n++)
                if (gcd64(n, d) == 1)
                    fracs.add(new long[] { n, d });
        Set<Long> allowed = new HashSet<>();
        for (long[] f : fracs)
            allowed.add(fracHash(f));

        Set<Long> sums = new HashSet<>();
        for (long[] x : fracs)
            for (long[] y : fracs) {
                // n=1: z = x+y
                long[] z1 = addFrac(x, y);
                if (z1[1] > 0 && allowed.contains(fracHash(z1)))
                    sums.add(fracHash(addFrac(addFrac(x, y), z1)));
                // n=-1: z = xy/(x+y)
                long[] s = addFrac(x, y);
                if (s[0] != 0) {
                    long[] zm1 = divFrac(mulFrac(x, y), s);
                    if (zm1[1] > 0 && allowed.contains(fracHash(zm1)))
                        sums.add(fracHash(addFrac(addFrac(x, y), zm1)));
                }
                // n=2: z^2 = x^2+y^2
                long[] z2sq = addFrac(sqFrac(x), sqFrac(y));
                long[] z2 = sqrtFrac(z2sq);
                if (z2 != null && allowed.contains(fracHash(z2)))
                    sums.add(fracHash(addFrac(addFrac(x, y), z2)));
                // n=-2: 1/z^2 = 1/x^2 + 1/y^2
                long[] ixsq = sqFrac(recipFrac(x)), iysq = sqFrac(recipFrac(y));
                long[] zm2sq = recipFrac(addFrac(ixsq, iysq));
                long[] zm2 = sqrtFrac(zm2sq);
                if (zm2 != null && allowed.contains(fracHash(zm2)))
                    sums.add(fracHash(addFrac(addFrac(x, y), zm2)));
            }

        // Recover fractions from hashes and sum them
        Map<Long, long[]> hashToFrac = new HashMap<>();
        for (long[] x : fracs)
            for (long[] y : fracs) {
                long[][] zCands = new long[4][];
                zCands[0] = addFrac(x, y);
                long[] s = addFrac(x, y);
                zCands[1] = s[0] != 0 ? divFrac(mulFrac(x, y), s) : null;
                long[] z2 = sqrtFrac(addFrac(sqFrac(x), sqFrac(y)));
                zCands[2] = z2;
                long[] zm2 = sqrtFrac(recipFrac(addFrac(sqFrac(recipFrac(x)), sqFrac(recipFrac(y)))));
                zCands[3] = zm2;
                for (long[] z : zCands) {
                    if (z != null && z[1] > 0 && allowed.contains(fracHash(z))) {
                        long[] sv = addFrac(addFrac(x, y), z);
                        long h = fracHash(sv);
                        if (sums.contains(h))
                            hashToFrac.put(h, sv);
                    }
                }
            }

        BigInteger num = BigInteger.ZERO, den = BigInteger.ONE;
        for (long[] f : hashToFrac.values()) {
            BigInteger fn = BigInteger.valueOf(f[0]), fd = BigInteger.valueOf(f[1]);
            num = num.multiply(fd).add(den.multiply(fn));
            den = den.multiply(fd);
            BigInteger g = num.gcd(den);
            num = num.divide(g);
            den = den.divide(g);
        }
        System.out.println(num.add(den));
    }
}