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 every integer $n$, positive, zero or negative, the problem combines three expressions in $x$, $y$ and $z$ into

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

A triple $(x,y,z)$ is a golden triple of order $k$ if each of $x$, $y$ and $z$ is a rational number of the form $a/b$ with $0\lt a\lt b\le k$, and $f_n(x,y,z)=0$ for at least one integer $n$. For every golden triple of order $35$ we form the sum $s=x+y+z$, add up the distinct values of $s$, write the total in lowest terms as $t=u/v$, and must find $u+v$. The statement gives no checkpoint values. ^[1]

The whole solution in one sentence: the long expression factors as $(x+y+z)(x^n+y^n-z^n)$, so a golden triple is a solution of $x^n+y^n=z^n$ in fractions between $0$ and $1$; Fermat's Last Theorem leaves only the exponents $1$, $-1$, $2$ and $-2$, each of which fixes $z$ once $x$ and $y$ are chosen, so we try every pair $(x,y)$ of allowed fractions with four candidates for $z$ and add the distinct sums exactly.

Mathematical Approach

We simplify $f_n$, reduce the exponents to four with Fermat's Last Theorem, and turn each of them into a formula for $z$.

1. The allowed numbers

Writing a number $a/b$ with $0\lt a\lt b\le k$ in lowest terms can only lower its denominator, so the allowed numbers are exactly the fractions in lowest terms strictly between $0$ and $1$ whose denominator is at most $k$; we call this set $S_k$. For $k=5$ it consists of the nine numbers $\tfrac15$, $\tfrac14$, $\tfrac13$, $\tfrac25$, $\tfrac12$, $\tfrac35$, $\tfrac23$, $\tfrac34$ and $\tfrac45$. These are the terms of the Farey sequence of order $k$, taken between $0$ and $1$, without its end points $0/1$ and $1/1$. ^[2] For a denominator $b\ge2$ the possible numerators are the $a$ with $1\le a\lt b$ and $\gcd(a,b)=1$, and Euler's totient function $\varphi(b)$ counts exactly these numbers, ^[3] so

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

For $k=5$ this is $1+2+2+4=9$, and for $k=35$ it is $383$. The sum grows like $3k^2/\pi^2$, ^[2] which predicts about $372$ for $k=35$.

2. A first look: the exponent $n=1$

For $n=1$ the three pieces are short: $x^2+y^2-z^2$, then $(xy+yz+zx)(1+1-1)$, and $xyz\,(x^{-1}+y^{-1}-z^{-1})$${}=yz+zx-xy$. Subtracting the third piece from the sum of the first two cancels $yz$ and $zx$ and leaves

$$\begin{aligned}f_1&=x^2+2xy+y^2-z^2\\&=(x+y+z)(x+y-z).\end{aligned}$$

Since $x+y+z\gt0$, the equation $f_1=0$ simply says $z=x+y$; for example $(\tfrac13,\tfrac13,\tfrac23)$ is a golden triple of order $3$.

3. The factorisation for every exponent

The general case rests on one observation: the products $xy+yz+zx$ and $xyz$ in $f_n$ are, together with $x+y+z$, the coefficients of a cubic whose roots are $x$, $y$ and $z$. Write $p=x+y+z$, $q=xy+yz+zx$ and $r=xyz$. Multiplying out gives

$$\begin{aligned}&(T-x)(T-y)(T-z)\\&\quad=T^3-pT^2+qT-r,\end{aligned}$$

the cubic case of Vieta's formulas. ^[4] Setting $T=x$ makes the left side zero, so $x^3=px^2-qx+r$. Multiplying by $x^{n-2}$, which is allowed for every integer $n$ because $x\ne0$, gives

$$x^{n+1}+q\,x^{n-1}-r\,x^{n-2}=p\,x^{n},$$

and the same identity holds with $y$ or $z$ in place of $x$. Now $f_{1,n}$ collects the terms $w^{n+1}$, $f_{2,n}$ the terms $q\,w^{n-1}$ and $f_{3,n}$ the terms $r\,w^{n-2}$ for $w=x,y,z$, with the sign $+$ for $x$ and $y$ and $-$ for $z$. So $f_n$ is the left side of the identity for $x$, plus that for $y$, minus that for $z$, and the right sides combine in the same way:

$$\begin{aligned}&f_n(x,y,z)\\&\quad=(x+y+z)\,(x^n+y^n-z^n).\end{aligned}$$

For $n=1$ this is Section 2, and at $x=y=z=1$ both sides equal $3$ for every $n$. Since the factor $x+y+z$ is positive, a triple of positive numbers makes $f_n$ vanish exactly when

$$x^n+y^n=z^n.$$

4. Fermat's Last Theorem leaves four exponents

The problem allows every integer $n$, but almost all of them are impossible:

$n=0$: the equation reads $1+1=1$.
$n\ge3$: suppose that positive rationals satisfy $x^n+y^n=z^n$, and let $d$ be a common denominator of $x$, $y$ and $z$. Multiplying by $d^n$ gives positive integers $a=dx$, $b=dy$, $c=dz$ with $a^n+b^n=c^n$. Fermat's Last Theorem, stated by Pierre de Fermat around 1637 and proved by Andrew Wiles in 1994, says that no such integers exist. ^[5]
$n\le-3$: with $m=-n\ge3$ we have $x^n=(1/x)^m$, so the equation says $(1/x)^m+(1/y)^m=(1/z)^m$ for the positive rationals $1/x$, $1/y$, $1/z$, which the previous case excludes. ^[5]

Only $n=1$, $-1$, $2$ and $-2$ remain, and all four occur, for example in $(\tfrac13,\tfrac13,\tfrac23)$, $(\tfrac23,\tfrac23,\tfrac13)$, $(\tfrac12,\tfrac23,\tfrac56)$ and $(\tfrac13,\tfrac14,\tfrac15)$, in this order.

5. Each exponent fixes $z$

Once $x$ and $y$ are chosen, each of the four equations has exactly one positive real solution $z$, because $w\mapsto w^n$ is strictly increasing on the positive numbers for $n\gt0$ and strictly decreasing for $n\lt0$:

$n=1$: $z=x+y$.
$n=-1$: $1/z=1/x+1/y$, so $z=xy/(x+y)$. This is the optic equation, and $z$ is half the harmonic mean of $x$ and $y$. ^[6] For $x=\tfrac12$ and $y=\tfrac13$ it gives $z=\tfrac16\big/\tfrac56=\tfrac15$.
$n=2$: $z=\sqrt{x^2+y^2}$, so $x$, $y$ and $z$ are the sides of a right triangle, and when $z$ is rational, multiplying them by a common denominator gives a Pythagorean triple. ^[7] The triple $(\tfrac12,\tfrac23,\tfrac56)$ is $(3,4,5)$ divided by $6$.
$n=-2$: $1/z^2=1/x^2+1/y^2$, so $z=xy/\sqrt{x^2+y^2}$, the altitude to the hypotenuse of the right triangle with legs $x$ and $y$; this relation is the inverse Pythagorean theorem. ^[8] For $x=\tfrac13$ and $y=\tfrac14$ we get $1/z^2=9+16=25$, so $z=\tfrac15$.

The negative exponents give $z\lt\min(x,y)$, for instance $xy/(x+y)\lt xy/y=x$, so their candidates are always below $1$. The positive exponents give $z\gt\max(x,y)$, and $z$ may reach $1$: $(\tfrac12,\tfrac12)$ gives $z=1$ for $n=1$ and $(\tfrac35,\tfrac45)$ gives $z=1$ for $n=2$, and neither counts, because $a\lt b$ forces $z\lt1$. A golden triple is thus determined by the pair $(x,y)$ and the exponent, so it suffices to scan the ordered pairs, each with at most four candidates.

6. Rational square roots

For $n=\pm2$ the candidate is a square root, and it counts only if it is rational. A positive fraction $A/B$ in lowest terms is the square of a rational number exactly when $A$ and $B$ are both perfect squares. If $A=a^2$ and $B=b^2$, then $A/B=(a/b)^2$. Conversely, let $A/B=(a/b)^2$ with $\gcd(a,b)=1$, so that $Ab^2=Ba^2$. Then $B$ divides $Ab^2$ and is coprime to $A$, so $B$ divides $b^2$ by the general form of Euclid's lemma; ^[9] and $b^2$ divides $Ba^2$ and is coprime to $a^2$, so $b^2$ divides $B$. Hence $B=b^2$, and then $A=a^2$.

So the programs reduce $z^2$, take integer square roots of its numerator and its denominator, and accept the candidate only if both roots are exact; the root $a/b$ is then automatically in lowest terms. For $n=2$ and $x=y=\tfrac12$ we get $z^2=\tfrac12$, and $2$ is not a square, so there is no rational $z$. For $x=\tfrac12$ and $y=\tfrac23$ we get $z^2=\tfrac14+\tfrac49=\tfrac{25}{36}$, so $z=\tfrac56$. For $n=-2$ the same test is applied to $z^2=1/(1/x^2+1/y^2)$.

7. Which candidates count

A candidate $z$ must itself lie in $S_k$. Because every value is kept in lowest terms, this means $0\lt z\lt1$ and a reduced denominator of at most $k$, and one lookup in a set that holds the elements of $S_k$ tests both. With $k=5$, the pair $(\tfrac12,\tfrac13)$ gives $z=\tfrac56$ for $n=1$, which fails because $6\gt5$, and $z=\tfrac15$ for $n=-1$, which counts.

8. Worked example: order 5

Running through all $9^2=81$ ordered pairs of $S_5$ gives the following golden triples, listed with $x\le y$; when $x\lt y$, swapping $x$ and $y$ gives another golden triple with the same sum.

$n=1$: $(\tfrac15,\tfrac15,\tfrac25)$, $(\tfrac15,\tfrac25,\tfrac35)$, $(\tfrac15,\tfrac35,\tfrac45)$, $(\tfrac25,\tfrac25,\tfrac45)$, $(\tfrac14,\tfrac14,\tfrac12)$, $(\tfrac14,\tfrac12,\tfrac34)$ and $(\tfrac13,\tfrac13,\tfrac23)$, with the sums $\tfrac45$, $\tfrac65$, $\tfrac85$, $\tfrac85$, $1$, $\tfrac32$ and $\tfrac43$;
$n=-1$: $(\tfrac13,\tfrac12,\tfrac15)$, $(\tfrac25,\tfrac25,\tfrac15)$, $(\tfrac25,\tfrac23,\tfrac14)$, $(\tfrac12,\tfrac12,\tfrac14)$, $(\tfrac35,\tfrac34,\tfrac13)$, $(\tfrac23,\tfrac23,\tfrac13)$ and $(\tfrac45,\tfrac45,\tfrac25)$, with the sums $\tfrac{31}{30}$, $1$, $\tfrac{79}{60}$, $\tfrac54$, $\tfrac{101}{60}$, $\tfrac53$ and $2$;
$n=-2$: $(\tfrac14,\tfrac13,\tfrac15)$ and $(\tfrac12,\tfrac23,\tfrac25)$, with the sums $\tfrac{47}{60}$ and $\tfrac{47}{30}$;
$n=2$: none; the first one, $(\tfrac12,\tfrac23,\tfrac56)$, needs order $6$.

These $16$ triples give $24$ ordered golden triples, since the eight with $x\ne y$ also occur swapped. Two sums repeat: $1$ comes from $(\tfrac25,\tfrac25,\tfrac15)$ and $(\tfrac14,\tfrac14,\tfrac12)$, and $\tfrac85$ from $(\tfrac15,\tfrac35,\tfrac45)$ and $(\tfrac25,\tfrac25,\tfrac45)$. That leaves $14$ distinct sums. Counted in sixtieths, the six distinct sums for $n=1$ add up to $446$, the six new ones for $n=-1$ to $537$ and the two for $n=-2$ to $141$, so

$$\begin{aligned}t&=\frac{446+537+141}{60}\\&=\frac{1124}{60}=\frac{281}{15},\end{aligned}$$

and $u+v=296$ for order $5$. The same procedure gives $t=3$, $\tfrac{27}{4}$, $\tfrac{423}{20}$ and $\tfrac{1735}{42}$ for the orders $3$, $4$, $6$ and $7$.

9. Distinct sums and exact arithmetic

Each value of $s$ counts once, so the programs collect the sums in a set of reduced fractions, where swapped pairs and different triples with equal sums collapse automatically. At order $35$, the $13320$ ordered golden triples give only $3677$ distinct sums, $193$ of which arise from more than one exponent. Exact fractions are essential: rounded floating-point values of one sum, reached along two different routes, could differ and be counted twice.

How the Code Works

The C++ program ^[10] keeps every number as a Fraction of two unsigned 64-bit integers num and den in lowest terms, so equal numbers have equal pairs and FractionHash can store them in a std::unordered_set. add_fraction and mul_fraction compute in unsigned 128-bit integers, reduce with gcd_u128 and return through cast_u64_checked, which stops the program on overflow. is_square_u64 takes a long double square root as a first guess and corrects it by exact integer comparisons, and sqrt_fraction applies it to the numerator and the denominator (Section 6).

build_fraction_list returns the $383$ elements of $S_{35}$ (Section 1), and collect_sums_fast stores them in the set allowed and calls process_pair for every ordered pair $(x,y)$. process_pair forms the candidates of Section 5 for $n=1,-1,2,-2$: $z_1=x+y$, then $xy/z_1$, then the roots of $x^2+y^2$ and of the reciprocal of $1/x^2+1/y^2$. add_sum_if_member looks $z$ up in allowed (Section 7) and inserts $x+y+z$ into the set sums. sum_distinct_s_values adds the distinct sums with a numerator and a denominator of Boost's arbitrary-precision type cpp_int, reducing after every step, and main prints $u+v$. The option --order=K selects another order.

Unless --skip-checkpoints is given, run_validation_checkpoints runs first. At order $7$, collect_sums_bruteforce_original evaluates the three-piece formula of the statement in exact big rationals (RationalBig, pow_rational) for $n=-2,-1,1,2$ on all $17^3=4913$ triples, and its set of distinct sums must equal that of collect_sums_bruteforce_derived, which tests the four equations of Section 5. This is a finite consistency check: the factorisation for every triple and every exponent rests on the proof of Section 3, and the exclusion of all other exponents on Section 4. At order $12$ the pair scan must match the derived brute force over $45^3=91125$ triples, and a last check repeats the scan at the target order.

The Python function solve builds the list fracs and the set allowed, runs the four branches for every ordered pair, collects the distinct sums in sums_set and adds them with sum. It uses fractions.Fraction, which keeps lowest terms with integers of unlimited size, and math.isqrt for the square-root test. ^[11] The Java program stores fractions as long[] pairs and uses gcd64, addFrac, mulFrac, recipFrac and divFrac in plain 64-bit arithmetic; no numerator, denominator or intermediate product of these operations exceeds $2750125$. sqrtFrac uses Math.sqrt, which is correctly rounded and hence exact for perfect squares of this size, ^[11] and also tries the numerator root plus one. Its sets hold the key fracHash, the numerator times $1000000007$ plus the denominator; every key stays below $10^{14}$ and is unique because every denominator is far smaller than $1000000007$. A second pass over all pairs recovers the fraction behind each key of sums, and BigInteger adds them. The Python and Java programs fix the order at $35$ and have no checkpoints.

Complexity and Verification

With $M=|S_k|\approx3k^2/\pi^2$ allowed fractions, the scan visits $M^2$ ordered pairs. Each pair needs a fixed number of fraction operations, two rational-square tests with up to four integer square roots, and at most four lookups in allowed of expected constant time. Counting each operation on these small numbers as one step, the time grows like $M^2$, that is like $k^4$; the Java program scans the pairs twice. The memory holds $S_k$ and the distinct sums. For $k=35$ the scan visits $383^2=146689$ pairs instead of the $383^3=56181887$ triples of a direct search, and on the computer used for this article each program finished in under two seconds, the C++ checkpoints included. The denominator of a sum divides the least common multiple of the three denominators, which is at most the product of the distinct ones and so at most $35\cdot34\cdot33=39270$, and the golden triple $(\tfrac2{33},\tfrac2{35},\tfrac1{34})$ of the exponent $-1$ attains this. The common denominator of thousands of sums can be far larger, so all three programs do the final addition with arbitrary-precision integers.

The following checks were made for this article and are not part of the published programs:

Computer algebra confirms the identity of Section 3 symbolically for every $n$ from $-6$ to $7$, and exact arithmetic confirms it at $2000$ random rational points.
For every order $k\le10$, a brute force evaluates the three-piece formula for all triples of $S_k$ and all exponents with $|n|\le20$. It finds exactly the golden triples of the pair scan, and each of them satisfies exactly one exponent, always one of $\pm1$ and $\pm2$.
A separately written program tests all $56181887$ triples of $S_{35}$ by cross-multiplication in integers, without square roots and without the formulas for $z$. It finds the same $13320$ ordered golden triples, $6289$ for $n=1$, $6487$ for $n=-1$, $302$ for $n=2$ and $242$ for $n=-2$, none of them with two exponents, and the same $3677$ distinct sums.
The C++, Python and Java programs print the same value, which agrees with the total of that separate program, and the C++ program passes its checkpoints and prints $4$, $31$, $296$, $443$ and $1777$ for the orders $3$ to $7$.

Footnotes and References

The references below supply background facts; the factorisation, the reduction to four exponents, the square-root test and all counts are derived or computed here.

Project Euler 180 — Golden Triplets. The official statement defines f_n = f_{1,n} + f_{2,n} − f_{3,n} for every integer n and calls (x, y, z) a golden triple of order k if x, y and z are rationals a/b with 0 < a < b ≤ k and f_n(x, y, z) = 0 for some integer n; it asks for u + v, where u/v is the reduced sum of the distinct values of x + y + z for order 35.
Wikipedia — Farey sequence. The Farey sequence of order n consists of the completely reduced fractions with denominators at most n, and some authors restrict it to the interval between 0 and 1; that version has length 1 + φ(1) + … + φ(n), which grows like 3n²/π².
Wikipedia — Euler's totient function. φ(n) counts the positive integers up to n that are relatively prime to n.
Wikipedia — Vieta's formulas. A monic polynomial is the product of the factors x − rᵢ over its roots; for a cubic the coefficients are, up to sign, r₁ + r₂ + r₃, r₁r₂ + r₁r₃ + r₂r₃ and r₁r₂r₃.
Wikipedia — Fermat's Last Theorem. No positive integers satisfy aⁿ + bⁿ = cⁿ for an integer n > 2; Fermat stated this around 1637, and Andrew Wiles released the first proof in 1994. The article also gives the equivalent form for rational numbers and shows that exponents n < −2 lead back to the positive case.
Wikipedia — Optic equation. In 1/a + 1/b = 1/c, the number c is half the harmonic mean of a and b.
Wikipedia — Pythagorean triple. Positive integers a, b, c with a² + b² = c², such as (3, 4, 5); every multiple of a triple is again a triple.
Wikipedia — Inverse Pythagorean theorem. In a right triangle with legs a and b, hypotenuse c and altitude h to the hypotenuse, 1/a² + 1/b² = 1/h² and h = ab/c.
Wikipedia — Euclid's lemma. In its general form: if an integer n divides a product ab and is coprime to a, then n divides b.
C++ — Euler180.cpp. The linked immutable revision, with process_pair, sum_distinct_s_values and run_validation_checkpoints, is the implementation record for this article.
Python — fractions: a Fraction keeps its numerator and denominator in lowest terms. Python — math.integer: isqrt(n) is the greatest integer a with a² ≤ n, and math.isqrt is the same function. Java — Math.sqrt: the result is the correctly rounded square root, the double closest to the exact value.

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));
    }
}