Problem 540: Counting Primitive Pythagorean Triples
View on Project EulerProject Euler Problem 540 Solution
EulerSolve provides an optimized solution for Project Euler Problem 540, Counting Primitive Pythagorean Triples, with C++, Python, Java, and a step-by-step mathematical explanation.
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....
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
- Problem page: https://projecteuler.net/problem=540
- Pythagorean triple: Wikipedia — Pythagorean triple
- Coprime integers: Wikipedia — Coprime integers
- Möbius inversion formula: Wikipedia — Möbius inversion formula
- Möbius function: Wikipedia — Möbius function
Problem 540 source code
C++
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <vector>
using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
static inline u64 isqrt_u64(const u64 x) {
u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
while ((u128)(r + 1) * (r + 1) <= (u128)x) ++r;
while ((u128)r * r > (u128)x) --r;
return r;
}
static inline u64 icbrt_u64(const u64 x) {
u64 r = static_cast<u64>(std::cbrt(static_cast<long double>(x)));
while ((u128)(r + 1) * (r + 1) * (r + 1) <= (u128)x) ++r;
while ((u128)r * r * r > (u128)x) --r;
return r;
}
static i64 helper(const u64 L) {
i64 res = 0;
u64 m0 = static_cast<u64>(0.5L + std::sqrt(0.25L + static_cast<long double>(L - 1) / 2.0L));
u64 K = m0 / 2;
res += static_cast<i64>(K * K);
u64 m = 2 * K + 1;
u64 n = isqrt_u64((L > m * m) ? (L - m * m) : 0ULL);
while (true) {
u64 z = (L > m * m) ? (L - m * m) : 0ULL;
if (z == 0ULL) {
break;
}
while ((u128)n * n > (u128)z) {
--n;
}
res += static_cast<i64>(n / 2);
++m;
z = (L > m * m) ? (L - m * m) : 0ULL;
if (z == 0ULL) {
break;
}
while ((u128)n * n > (u128)z) {
--n;
}
res += static_cast<i64>((n + 1) / 2);
++m;
}
return res;
}
static u64 solve(const u64 N) {
const u64 L = icbrt_u64(N);
std::vector<i64> v(static_cast<std::size_t>(L + 1), 0);
std::vector<i64> bigV(static_cast<std::size_t>(L / 2 + 1), 0);
for (u64 x = 1; x <= L; ++x) {
i64 res = helper(x);
const u64 c = icbrt_u64(x);
for (u64 g = 3; g <= c; g += 2) {
res -= v[static_cast<std::size_t>(x / (g * g))];
}
i64 prev_y = static_cast<i64>((isqrt_u64(x) - 1) / 2);
const u64 z_end = c + ((x / (c * c) != c) ? 1ULL : 0ULL);
for (u64 z = 1; z < z_end; ++z) {
i64 y = static_cast<i64>((isqrt_u64(x / (z + 1)) - 1) / 2);
res -= (prev_y - y) * v[static_cast<std::size_t>(z)];
prev_y = y;
}
v[static_cast<std::size_t>(x)] = res;
}
for (i64 a = static_cast<i64>((L - 1) / 2); a >= 0; --a) {
const u64 x = static_cast<u64>(2 * a + 1);
const u64 k = N / (x * x);
i64 res = helper(k);
const u64 c = icbrt_u64(k);
for (u64 b = 1; b <= (c - 1) / 2; ++b) {
const u64 g = 2 * b + 1;
const u64 k_gg = k / (g * g);
if (k_gg <= L) {
res -= v[static_cast<std::size_t>(k_gg)];
} else {
const u64 idx = 2ULL * static_cast<u64>(a) * b + static_cast<u64>(a) + b;
res -= bigV[static_cast<std::size_t>(idx)];
}
}
i64 prev_y = static_cast<i64>((isqrt_u64(k) - 1) / 2);
const u64 z_end = c + ((k / (c * c) != c) ? 1ULL : 0ULL);
for (u64 z = 1; z < z_end; ++z) {
i64 y = static_cast<i64>((isqrt_u64(k / (z + 1)) - 1) / 2);
res -= (prev_y - y) * v[static_cast<std::size_t>(z)];
prev_y = y;
}
bigV[static_cast<std::size_t>(a)] = res;
}
return static_cast<u64>(bigV[0]);
}
int main() {
assert(solve(20ULL) == 3ULL);
assert(solve(1'000'000ULL) == 159139ULL);
std::cout << solve(3'141'592'653'589'793ULL) << '\n';
return 0;
}
Python
import math
def solve():
N = 3141592653589793
def isqrt(x):
r = int(math.isqrt(x))
while (r+1)*(r+1) <= x: r += 1
while r*r > x: r -= 1
return r
def icbrt(x):
r = round(x ** (1/3))
while (r+1)**3 <= x: r += 1
while r**3 > x: r -= 1
return r
def helper(L):
res = 0
m0sq = 0.25 + (L-1)/2.0
m0 = int(0.5 + math.sqrt(m0sq))
K = m0 // 2
res += K * K
m = 2*K + 1
n = isqrt(max(0, L - m*m))
while True:
z = max(0, L - m*m)
if z == 0: break
while n*n > z: n -= 1
res += n // 2
m += 1
z = max(0, L - m*m)
if z == 0: break
while n*n > z: n -= 1
res += (n+1) // 2
m += 1
return res
L = icbrt(N)
v = [0] * (L + 1)
for x in range(1, L + 1):
res = helper(x)
c = icbrt(x)
for g in range(3, c+1, 2):
res -= v[x // (g*g)]
prev_y = (isqrt(x) - 1) // 2
z_end = c + (1 if x // (c*c) != c else 0)
for z in range(1, z_end):
y = (isqrt(x // (z+1)) - 1) // 2
res -= (prev_y - y) * v[z]
prev_y = y
v[x] = res
bigV = [0] * (L // 2 + 1)
for a in range((L-1)//2, -1, -1):
x = 2*a + 1
k = N // (x*x)
res = helper(k)
c = icbrt(k)
for b in range(1, (c-1)//2 + 1):
g = 2*b + 1
k_gg = k // (g*g)
if k_gg <= L:
res -= v[k_gg]
else:
idx = 2*a*b + a + b
res -= bigV[idx]
prev_y = (isqrt(k) - 1) // 2
z_end = c + (1 if k // (c*c) != c else 0)
for z in range(1, z_end):
y = (isqrt(k // (z+1)) - 1) // 2
res -= (prev_y - y) * v[z]
prev_y = y
bigV[a] = res
return str(bigV[0])
if __name__ == '__main__':
print(solve())
Java
public class Euler540 {
static long isqrt(long x) {
if (x == 0)
return 0;
long r = (long) Math.sqrt((double) x);
while ((r + 1) * (r + 1) <= x && (r + 1) * (r + 1) > 0) {
r++;
}
while (r * r > x || r * r < 0) {
r--;
}
return r;
}
static long icbrt(long x) {
if (x == 0)
return 0;
long r = (long) Math.cbrt((double) x);
while (Math.pow(r + 1, 3) <= x) {
r++;
}
while (r * r * r > x || r * r * r < 0) {
r--;
}
return r;
}
static long helper(long L) {
if (L == 0)
return 0;
long res = 0;
long m0 = (long) (0.5 + Math.sqrt(0.25 + (double) (L - 1) / 2.0));
long K = m0 / 2;
res += K * K;
long m = 2 * K + 1;
long n = isqrt(L > m * m ? L - m * m : 0);
while (true) {
long z = L > m * m ? L - m * m : 0;
if (z == 0)
break;
while (n * n > z)
n--;
res += n / 2;
m++;
z = L > m * m ? L - m * m : 0;
if (z == 0)
break;
while (n * n > z)
n--;
res += (n + 1) / 2;
m++;
}
return res;
}
public static String solve() {
long N = 3141592653589793L;
long L = icbrt(N);
long[] v = new long[(int) L + 1];
long[] bigV = new long[(int) (L / 2) + 1];
for (long x = 1; x <= L; x++) {
long res = helper(x);
long c = icbrt(x);
for (long g = 3; g <= c; g += 2) {
res -= v[(int) (x / (g * g))];
}
long prev_y = (isqrt(x) - 1) / 2;
long z_end = c + ((x / (c * c) != c) ? 1 : 0);
for (long z = 1; z < z_end; z++) {
long y = (isqrt(x / (z + 1)) - 1) / 2;
res -= (prev_y - y) * v[(int) z];
prev_y = y;
}
v[(int) x] = res;
}
for (long a = (L - 1) / 2; a >= 0; a--) {
long x = 2 * a + 1;
long k = N / (x * x);
long res = helper(k);
long c = icbrt(k);
for (long b = 1; b <= (c - 1) / 2; b++) {
long g = 2 * b + 1;
long k_gg = k / (g * g);
if (k_gg <= L) {
res -= v[(int) k_gg];
} else {
int idx = (int) (2 * a * b + a + b);
res -= bigV[idx];
}
}
long prev_y = (isqrt(k) - 1) / 2;
long z_end = c + ((k / (c * c) != c) ? 1 : 0);
for (long z = 1; z < z_end; z++) {
long y = (isqrt(k / (z + 1)) - 1) / 2;
res -= (prev_y - y) * v[(int) z];
prev_y = y;
}
bigV[(int) a] = res;
}
return Long.toString(bigV[0]);
}
public static void main(String[] args) {
System.out.println(solve());
}
}