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

1. Problem page: Project Euler 748
2. Pythagorean parametrization: Wikipedia - Pythagorean triple
3. Composition of sums of two squares: Wikipedia - Brahmagupta-Fibonacci identity
4. Gaussian integer viewpoint: Wikipedia - Gaussian integer
5. Background on Diophantine equations: Wikipedia - Diophantine equation