Detailed mathematical approach

Problem Summary

We want the number of primitive Pythagorean triples \((a,b,c)\) with hypotenuse bound

$c \le N.$

Euclid's parametrization says that every primitive triple appears uniquely in the form

$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,$

with

$m>n>0,\qquad \gcd(m,n)=1,\qquad m\not\equiv n \pmod{2}.$

So the problem is equivalent to counting parameter pairs \((m,n)\) satisfying those three conditions together with \(m^2+n^2\le N\).

Mathematical Approach

Let \(P(X)\) denote the number of primitive parameter pairs \((m,n)\) with \(m^2+n^2\le X\), \(m>n>0\), coprime, and of opposite parity. The desired answer is \(P(N)\). The implementations compute it in two layers: first a fast raw counter without the gcd condition, then a recurrence that removes non-primitive pairs.

Step 1: Convert Primitive Triples into Parameter Pairs

The classical parametrization gives a bijection between primitive triples and pairs \((m,n)\) such that

$m>n>0,\qquad \gcd(m,n)=1,\qquad m\not\equiv n \pmod{2}.$

The hypotenuse becomes

$c=m^2+n^2,$

so counting primitive triples with \(c\le N\) is exactly the same as counting those admissible pairs inside the quarter-circle

$m^2+n^2\le N.$

This change of viewpoint is the foundation of the whole solution: geometry gives the bound, and arithmetic gives the primitivity conditions.

Step 2: Count Opposite-Parity Pairs Before Enforcing Coprimality

Define \(R(X)\) as the number of pairs \((m,n)\) with

$m>n>0,\qquad m^2+n^2\le X,\qquad m\not\equiv n \pmod{2},$

but without the condition \(\gcd(m,n)=1\).

For a fixed \(m\), the largest allowed \(n\) is

$n_{\max}(m)=\min\left(m-1,\left\lfloor\sqrt{X-m^2}\right\rfloor\right).$

The row contribution is therefore

$r_X(m)=\begin{cases} \left\lfloor \dfrac{n_{\max}(m)}{2}\right\rfloor, & m\text{ odd},\\[6pt] \left\lfloor \dfrac{n_{\max}(m)+1}{2}\right\rfloor, & m\text{ even}. \end{cases}$

Hence

$R(X)=\sum_{m\ge 1} r_X(m).$

The code evaluates this quickly by separating a full prefix of rows. If

$m^2+(m-1)^2\le X,$

then the entire interval \(1\le n\le m-1\) is inside the circle. Solving this inequality gives

$m\le \frac{1+\sqrt{2X-1}}{2}.$

Let

$m_0=\left\lfloor\frac{1+\sqrt{2X-1}}{2}\right\rfloor,\qquad K=\left\lfloor\frac{m_0}{2}\right\rfloor.$

The first \(2K\) rows contribute exactly

$1+3+5+\cdots+(2K-1)=K^2,$

because row \(2j-1\) contributes \(j-1\) and row \(2j\) contributes \(j\). After that prefix, the remaining upper bound \(\lfloor\sqrt{X-m^2}\rfloor\) decreases monotonically, so one moving square-root pointer is enough for the tail.

Step 3: Remove Non-Primitive Pairs by Their Odd Common Divisor

Now restore the coprimality condition. Opposite parity already rules out a factor \(2\), so any common divisor of \(m\) and \(n\) must be odd.

If \((m,n)\) has odd gcd \(d\), then

$m=d\,m',\qquad n=d\,n',$

and the reduced pair \((m',n')\) is primitive. The circle bound becomes

$m'^2+n'^2\le \left\lfloor \frac{X}{d^2}\right\rfloor.$

Therefore every raw pair is obtained uniquely from a primitive pair and an odd scaling factor, which gives

$R(X)=\sum_{\substack{d\ge 1\\ d\equiv 1 \pmod{2}}} P\left(\left\lfloor \frac{X}{d^2}\right\rfloor\right).$

Rearranging,

$P(X)=R(X)-\sum_{\substack{d\ge 3\\ d\equiv 1 \pmod{2}}} P\left(\left\lfloor \frac{X}{d^2}\right\rfloor\right).$

This is the recurrence used directly by the implementations. It is equivalent to odd-only Möbius inversion, but the code never needs an explicit table of Möbius values.

Step 4: Split the Recurrence at the Cube Root and Group Equal Quotients

Let

$L=\left\lfloor N^{1/3}\right\rfloor.$

First compute \(P(x)\) for every \(1\le x\le L\). For a fixed \(x\), odd divisors \(d\le x^{1/3}\) are handled one by one in

$P(x)=R(x)-\sum_{\substack{d\ge 3\\ d\equiv 1 \pmod{2}}} P\left(\left\lfloor \frac{x}{d^2}\right\rfloor\right).$

For larger \(d\), the quotient \(\left\lfloor x/d^2\right\rfloor\) is small, and many different \(d\) produce the same quotient. Those repeated values are grouped together. For a fixed small \(z\), the number of odd \(d\) with

$\left\lfloor \frac{x}{d^2}\right\rfloor=z$

is

$\left\lfloor\frac{\sqrt{x/z}-1}{2}\right\rfloor-\left\lfloor\frac{\sqrt{x/(z+1)}-1}{2}\right\rfloor.$

That converts many repeated recursive terms into one multiplication by a multiplicity.

After the small table is known, the implementations evaluate the transformed values

$Q(t)=P\left(\left\lfloor \frac{N}{t^2}\right\rfloor\right),\qquad t\text{ odd},$

processing odd \(t\) from large to small. When the recurrence for \(Q(t)\) reaches a reduced argument at most \(L\), it is read from the small table. Otherwise it is another transformed value \(Q(tu)\) with a larger odd factor, which has already been computed. Finally, the required answer is simply

$Q(1)=P(N).$

Worked Example: \(N=50\)

The opposite-parity pairs with \(m^2+n^2\le 50\) are

$ (2,1),\ (3,2),\ (4,1),\ (4,3),\ (5,2),\ (5,4),\ (6,1),\ (6,3). $

So \(R(50)=8\).

The only possible odd common divisor larger than \(1\) is \(3\), because \(5^2>50\). Thus

$P(50)=R(50)-P\left(\left\lfloor\frac{50}{9}\right\rfloor\right)=8-P(5).$

Now \(P(5)=1\), coming from the single pair \((2,1)\). Therefore

$P(50)=8-1=7.$

The seven primitive triples are

$ (3,4,5),\ (5,12,13),\ (8,15,17),\ (7,24,25),\ (20,21,29),\ (12,35,37),\ (9,40,41). $

This small example shows exactly how the raw geometric count and the odd-gcd recurrence fit together.

How the Code Works

The C++, Python, and Java implementations begin with exact integer square-root and cube-root routines so that every cutoff such as \(\lfloor \sqrt{x}\rfloor\) and \(\lfloor N^{1/3}\rfloor\) is computed safely despite floating-point starting estimates.

They then evaluate the raw counter \(R(X)\) by combining a closed-form prefix \(K^2\) with a monotone tail scan. The tail keeps one decreasing square-root boundary and alternates the parity formula for odd and even rows, which avoids recomputing the whole range of \(n\) for each \(m\).

Next the implementation fills a first table for all primitive counts \(P(x)\) with \(x\le L\). Each entry uses the recurrence above, split into individually handled large quotients and grouped small quotients.

After that it fills a second table for transformed arguments \(\left\lfloor N/t^2\right\rfloor\) with odd \(t\), processed in descending order. Small reduced arguments are read from the first table; large reduced arguments correspond to already computed transformed values with larger odd factors.

The entry for \(t=1\) is \(P(N)\), so it is exactly the number of primitive Pythagorean triples with hypotenuse at most \(N\).

Complexity Analysis

A single evaluation of the raw counter \(R(X)\) costs \(O(\sqrt{X})\) time because the tail scan visits only the remaining boundary rows. The recurrence work around that raw counter contributes about \(O(X^{1/3})\) grouped-divisor operations.

Building the small table up to \(L=\lfloor N^{1/3}\rfloor\) costs \(O(N^{1/2})\) time. In the transformed stage, the dominant raw-count cost is

$\sum_{\substack{t\le L\\ t\equiv 1 \pmod{2}}} O\left(\sqrt{\frac{N}{t^2}}\right)=O(\sqrt{N}\log N).$

The extra grouped-divisor work is smaller than that dominant term. The memory usage is \(O(N^{1/3})\), because both stored tables have size proportional to the cube-root threshold.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=540
2. Pythagorean triple: Wikipedia — Pythagorean triple
3. Coprime integers: Wikipedia — Coprime integers
4. Möbius inversion formula: Wikipedia — Möbius inversion formula
5. Möbius function: Wikipedia — Möbius function