Problem 748: Upside Down Diophantine Equation
View on Project EulerProject Euler Problem 748 Solution
EulerSolve provides an optimized solution for Project Euler Problem 748, Upside Down Diophantine Equation, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary We seek primitive positive integer triples \((x,y,z)\) satisfying $$\frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2},\qquad \gcd(x,y,z)=1,$$ together with the bound \(x,y,z\le N\). For every admissible triple we add \(x+y+z\), so the target quantity is $$S(N)=\sum_{\substack{(x,y,z)\text{ primitive}\\ \frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2}\\ x,y,z\le N}} (x+y+z).$$ A brute-force search over all triples up to \(N=10^{16}\) is impossible, so the solution reorganizes the equation until the search space is generated by primitive Pythagorean parameters and a sharp bound on one auxiliary variable. Mathematical Approach Step 1: Remove the common factor of \(x\) and \(y\) Write $$x=g\,a,\qquad y=g\,b,\qquad \gcd(a,b)=1.$$ Substituting into the reciprocal equation gives $$z^2(a^2+b^2)=13g^2a^2b^2.$$ Because \(\gcd(a,a^2+b^2)=\gcd(a,b^2)=1\), the factor \(a^2\) must divide \(z^2\), so \(a\mid z\). The same argument shows \(b\mid z\). Since \(\gcd(a,b)=1\), we can write $$z=abk$$ for some positive integer \(k\). After substitution, $$k^2(a^2+b^2)=13g^2.$$ The primitive condition \(\gcd(x,y,z)=1\) implies \(\gcd(g,z)=1\), hence \(\gcd(g,k)=1\)....
Detailed mathematical approach
Problem Summary
We seek primitive positive integer triples \((x,y,z)\) satisfying
$$\frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2},\qquad \gcd(x,y,z)=1,$$
together with the bound \(x,y,z\le N\). For every admissible triple we add \(x+y+z\), so the target quantity is
$$S(N)=\sum_{\substack{(x,y,z)\text{ primitive}\\ \frac{1}{x^2}+\frac{1}{y^2}=\frac{13}{z^2}\\ x,y,z\le N}} (x+y+z).$$
A brute-force search over all triples up to \(N=10^{16}\) is impossible, so the solution reorganizes the equation until the search space is generated by primitive Pythagorean parameters and a sharp bound on one auxiliary variable.
Mathematical Approach
Step 1: Remove the common factor of \(x\) and \(y\)
Write
$$x=g\,a,\qquad y=g\,b,\qquad \gcd(a,b)=1.$$
Substituting into the reciprocal equation gives
$$z^2(a^2+b^2)=13g^2a^2b^2.$$
Because \(\gcd(a,a^2+b^2)=\gcd(a,b^2)=1\), the factor \(a^2\) must divide \(z^2\), so \(a\mid z\). The same argument shows \(b\mid z\). Since \(\gcd(a,b)=1\), we can write
$$z=abk$$
for some positive integer \(k\). After substitution,
$$k^2(a^2+b^2)=13g^2.$$
The primitive condition \(\gcd(x,y,z)=1\) implies \(\gcd(g,z)=1\), hence \(\gcd(g,k)=1\). Therefore \(k^2\) divides \(13\), so the only possibility is
$$k=1.$$
Thus every primitive solution has the shape
$$x=q\,r,\qquad y=p\,r,\qquad z=pq,$$
where
$$p^2+q^2=13r^2,\qquad \gcd(p,q)=1.$$
The original reciprocal equation has been reduced to finding primitive representations of \(13r^2\) as a sum of two squares.
Step 2: Parametrize \(r^2\) by primitive Pythagorean triples
Any primitive solution of
$$u^2+v^2=r^2$$
is given by Euclid's parametrization
$$u=m^2-n^2,\qquad v=2mn,\qquad r=m^2+n^2,$$
with
$$m>n\ge 0,\qquad m-n\text{ odd},\qquad \gcd(m,n)=1.$$
Allowing \(n=0\) keeps the degenerate base case \(u=1\), \(v=0\), \(r=1\); for larger \(m\) the coprimality test automatically rejects \(n=0\). So the implementation can handle the smallest primitive solution without special-case code.
Step 3: Compose with \(13=3^2+2^2\)
The identity
$$ (A^2+B^2)(C^2+D^2)=(AC-BD)^2+(AD+BC)^2=(AC+BD)^2+(AD-BC)^2 $$
turns a representation of \(r^2\) into a representation of \(13r^2\). Using \(13=3^2+2^2\) and the pair \((u,v)\) from Step 2 gives two candidate branches:
$$p=\left|3u-2v\right|,\qquad q=\left|2u+3v\right|,$$
or
$$p=\left|3u+2v\right|,\qquad q=\left|2u-3v\right|.$$
In both cases,
$$p^2+q^2=(3^2+2^2)(u^2+v^2)=13r^2.$$
These two branches are the real-coordinate form of multiplying by \(3+2i\) or \(3-2i\) in the Gaussian integers. That viewpoint also explains why the parametrization is complete for primitive solutions.
Step 4: Return to \((x,y,z)\) and enforce primitiveness
Once a pair \((p,q)\) with \(p^2+q^2=13r^2\) is known, define
$$x=q\,r,\qquad y=p\,r,\qquad z=pq.$$
Then
$$ \frac{1}{x^2}+\frac{1}{y^2} =\frac{1}{q^2r^2}+\frac{1}{p^2r^2} =\frac{p^2+q^2}{p^2q^2r^2} =\frac{13r^2}{p^2q^2r^2} =\frac{13}{z^2}. $$
The implementations keep only positive pairs and require \(\gcd(p,q)=1\). That already forces the final triple to be primitive. Indeed, if a prime divided both \(r\) and \(p\), then from \(q^2=13r^2-p^2\) it would also divide \(q\), contradicting \(\gcd(p,q)=1\). So \(\gcd(r,p)=\gcd(r,q)=1\), hence \(\gcd(x,y,z)=1\).
To avoid counting the same solution twice with \(x\) and \(y\) swapped, the implementation orders the pair so that \(p\ge q\), which implies \(y\ge x\). It also discards the rare case where the two branches produce the same pair.
Step 5: Derive the search bound
The formulas above give
$$x^2+y^2=r^2(p^2+q^2)=13r^4.$$
If \(x\le N\) and \(y\le N\), then
$$x^2+y^2\le 2N^2,$$
so every valid triple must satisfy
$$13r^4\le 2N^2.$$
This yields the sharp search limit
$$r\le r_{\max}=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
Since \(r=m^2+n^2\), the outer parameter only needs to run up to \(\lfloor\sqrt{r_{\max}}\rfloor\).
Worked Example: \(N=100\)
For \(N=100\), the bound is
$$13r^4\le 2\cdot 100^2=20000,$$
so \(r_{\max}=6\).
The pair \((m,n)=(1,0)\) gives
$$u=1,\qquad v=0,\qquad r=1.$$
Both branches collapse to the same pair \((p,q)=(3,2)\). Therefore
$$x=2,\qquad y=3,\qquad z=6,$$
and the contribution is \(2+3+6=11\).
The pair \((m,n)=(2,1)\) gives
$$u=3,\qquad v=4,\qquad r=5.$$
The two branches produce \((p,q)=(18,1)\) and \((17,6)\). They map to
$$ (x,y,z)=(5,90,18)\quad \text{and}\quad (30,85,102). $$
Only the first triple satisfies \(x,y,z\le 100\), so its contribution is \(5+90+18=113\).
No other primitive parameter pair has \(r\le 6\). Hence
$$S(100)=11+113=124,$$
which matches the implementation's checkpoint.
How the Code Works
The C++, Python, and Java implementations follow the same pipeline. First they compute the largest admissible \(r\) by integer binary search on the inequality \(13r^4\le 2N^2\). That keeps the search exact and avoids floating-point edge cases at the upper boundary.
Next they enumerate all primitive Euclid parameters \((m,n)\) with the parity and coprimality conditions from Step 2. For each such pair, they form the two linear transforms from Step 3, reorder the resulting values so the larger one comes first, discard zero values, suppress a duplicated second branch when both branches coincide, and keep only coprime pairs.
Each surviving pair is converted into \((x,y,z)=(qr,pr,pq)\). If all three coordinates are at most \(N\), the implementation adds \(x+y+z\) to the running total. The arithmetic is done with sufficiently wide integer types so that intermediate products remain exact.
Complexity Analysis
Let
$$R=\left\lfloor\left(\frac{2N^2}{13}\right)^{1/4}\right\rfloor.$$
The parameter search visits primitive pairs \((m,n)\) with \(m^2+n^2\le R\), so the number of tested pairs is \(O(R)\). Because \(R=\Theta(\sqrt{N})\), the total running time is \(O(\sqrt{N})\) with only constant extra memory beyond the integer arithmetic used for safe products and accumulation.
Footnotes and References
- Problem page: Project Euler 748
- Pythagorean parametrization: Wikipedia - Pythagorean triple
- Composition of sums of two squares: Wikipedia - Brahmagupta-Fibonacci identity
- Gaussian integer viewpoint: Wikipedia - Gaussian integer
- Background on Diophantine equations: Wikipedia - Diophantine equation
Problem 748 source code
C++
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <numeric>
namespace {
using u64 = std::uint64_t;
using u128 = __uint128_t;
constexpr u64 kMod = 1'000'000'000ULL;
u64 isqrt_u64(const u64 n) {
const long double approx = static_cast<long double>(n);
u64 x = static_cast<u64>(std::sqrt(approx));
while ((x + 1ULL) * (x + 1ULL) <= n) {
++x;
}
while (x * x > n) {
--x;
}
return x;
}
u64 max_r(const u64 n_limit) {
const u128 rhs = 2ULL * static_cast<u128>(n_limit) * static_cast<u128>(n_limit);
u64 lo = 0ULL;
u64 hi = static_cast<u64>(std::sqrt(static_cast<long double>(n_limit))) + 10ULL;
while (13ULL * static_cast<u128>(hi) * hi * hi * hi <= rhs) {
hi *= 2ULL;
}
while (lo + 1ULL < hi) {
const u64 mid = lo + (hi - lo) / 2ULL;
const u128 lhs = 13ULL * static_cast<u128>(mid) * mid * mid * mid;
if (lhs <= rhs) {
lo = mid;
} else {
hi = mid;
}
}
return lo;
}
u128 S(const u64 n_limit) {
const u64 r_limit = max_r(n_limit);
const u64 m_limit = isqrt_u64(r_limit);
u128 total = 0;
for (u64 m = 1ULL; m <= m_limit; ++m) {
for (u64 n = 0ULL; n < m; ++n) {
if (((m - n) & 1ULL) == 0ULL) {
continue;
}
if (std::gcd(m, n) != 1ULL) {
continue;
}
const u64 r = m * m + n * n;
if (r > r_limit) {
continue;
}
const std::int64_t u = static_cast<std::int64_t>(m * m - n * n);
const std::int64_t v = static_cast<std::int64_t>(2ULL * m * n);
std::pair<u64, u64> cand[2] = {
{static_cast<u64>(std::llabs(3LL * u - 2LL * v)), static_cast<u64>(std::llabs(2LL * u + 3LL * v))},
{static_cast<u64>(std::llabs(3LL * u + 2LL * v)), static_cast<u64>(std::llabs(2LL * u - 3LL * v))},
};
for (int i = 0; i < 2; ++i) {
if (i == 1 && cand[1] == cand[0]) {
continue;
}
u64 p = cand[i].first;
u64 q = cand[i].second;
if (p == 0ULL || q == 0ULL) {
continue;
}
if (p < q) {
std::swap(p, q);
}
if (std::gcd(p, q) != 1ULL) {
continue;
}
const u128 x = static_cast<u128>(q) * static_cast<u128>(r);
const u128 y = static_cast<u128>(p) * static_cast<u128>(r);
const u128 z = static_cast<u128>(p) * static_cast<u128>(q);
if (x > n_limit || y > n_limit || z > n_limit) {
continue;
}
total += x + y + z;
}
}
}
return total;
}
u64 to_u64(const u128 value) {
return static_cast<u64>(value);
}
} // namespace
int main() {
assert(to_u64(S(100ULL)) == 124ULL);
assert(to_u64(S(1'000ULL)) == 1'470ULL);
assert(to_u64(S(100'000ULL)) == 2'340'084ULL);
const u64 answer = static_cast<u64>(S(10'000'000'000'000'000ULL) % kMod);
std::cout << std::setw(9) << std::setfill('0') << answer << '\n';
return 0;
}
Python
import math
def solve():
N = 10**16; KMOD = 10**9
def isqrt(n):
x = int(n**0.5)
while (x+1)*(x+1) <= n: x += 1
while x*x > n: x -= 1
return x
def max_r(nl):
rhs = 2*nl*nl; lo = 0; hi = int(nl**0.5)+10
while 13*hi**4 <= rhs: hi *= 2
while lo+1 < hi:
mid = (lo+hi)//2
if 13*mid**4 <= rhs: lo = mid
else: hi = mid
return lo
rl = max_r(N); ml = isqrt(rl)
total = 0
for m in range(1, ml+1):
for n in range(0, m):
if (m-n)%2 == 0: continue
if math.gcd(m, n) != 1: continue
r = m*m + n*n
if r > rl: continue
u = m*m - n*n; v = 2*m*n
cands = [
(abs(3*u-2*v), abs(2*u+3*v)),
(abs(3*u+2*v), abs(2*u-3*v)),
]
seen = set()
for p, q in cands:
if p == 0 or q == 0: continue
if p < q: p, q = q, p
if (p, q) in seen: continue
seen.add((p, q))
if math.gcd(p, q) != 1: continue
x = q*r; y = p*r; z = p*q
if x > N or y > N or z > N: continue
total += x + y + z
return str(total % KMOD).zfill(9)
if __name__ == '__main__':
print(solve())
Java
import java.math.BigInteger;
public class Euler748 {
static final long kMod = 1000000000L;
static long gcd(long a, long b) {
while (b != 0) {
long t = b;
b = a % b;
a = t;
}
return a;
}
static long isqrtU64(long n) {
if (n == 0)
return 0;
long x = (long) Math.sqrt(n);
while ((x + 1) * (x + 1) <= n)
x++;
while (x * x > n)
x--;
return x;
}
static long maxR(long nLimit) {
BigInteger bigN = BigInteger.valueOf(nLimit);
BigInteger rhs = BigInteger.valueOf(2).multiply(bigN).multiply(bigN);
long lo = 0;
long hi = (long) Math.sqrt(nLimit) + 10;
BigInteger thirteen = BigInteger.valueOf(13);
while (true) {
BigInteger bigHi = BigInteger.valueOf(hi);
BigInteger lhs = thirteen.multiply(bigHi).multiply(bigHi).multiply(bigHi).multiply(bigHi);
if (lhs.compareTo(rhs) <= 0) {
hi *= 2;
} else {
break;
}
}
while (lo + 1 < hi) {
long mid = lo + (hi - lo) / 2;
BigInteger bigMid = BigInteger.valueOf(mid);
BigInteger lhs = thirteen.multiply(bigMid).multiply(bigMid).multiply(bigMid).multiply(bigMid);
if (lhs.compareTo(rhs) <= 0) {
lo = mid;
} else {
hi = mid;
}
}
return lo;
}
static BigInteger S(long nLimit) {
long rLimit = maxR(nLimit);
long mLimit = isqrtU64(rLimit);
BigInteger total = BigInteger.ZERO;
BigInteger bigNLimit = BigInteger.valueOf(nLimit);
for (long m = 1; m <= mLimit; ++m) {
for (long n = 0; n < m; ++n) {
if (((m - n) & 1) == 0)
continue;
if (gcd(m, n) != 1)
continue;
long r = m * m + n * n;
if (r > rLimit)
continue;
long u = m * m - n * n;
long v = 2 * m * n;
long[][] cand = {
{ Math.abs(3 * u - 2 * v), Math.abs(2 * u + 3 * v) },
{ Math.abs(3 * u + 2 * v), Math.abs(2 * u - 3 * v) }
};
for (int i = 0; i < 2; ++i) {
if (i == 1 && cand[1][0] == cand[0][0] && cand[1][1] == cand[0][1]) {
continue;
}
long p = cand[i][0];
long q = cand[i][1];
if (p == 0 || q == 0)
continue;
if (p < q) {
long temp = p;
p = q;
q = temp;
}
if (gcd(p, q) != 1)
continue;
BigInteger x = BigInteger.valueOf(q).multiply(BigInteger.valueOf(r));
BigInteger y = BigInteger.valueOf(p).multiply(BigInteger.valueOf(r));
BigInteger z = BigInteger.valueOf(p).multiply(BigInteger.valueOf(q));
if (x.compareTo(bigNLimit) > 0 || y.compareTo(bigNLimit) > 0 || z.compareTo(bigNLimit) > 0) {
continue;
}
total = total.add(x).add(y).add(z);
}
}
}
return total;
}
public static String solve() {
BigInteger ans = S(10000000000000000L).remainder(BigInteger.valueOf(kMod));
return String.format("%09d", ans.longValue());
}
public static void main(String[] args) {
System.out.println(solve());
}
}