### February 25, 2020

A Pythagorean quadruple consists of four positive integers a, b, c and d such that abcd and a² + b² + c² = d². For instance, (2 3 6 7) is a Pythagorean quadruple because 2² + 3² + 6² = 4 + 9 + 36 = 49 = 7².

Your task is to write a program that counts the Pythagorian quadruples with a, b, c less than or equal to some given N, and compute the number of Pythagorian quadruples with N = 1000. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

### 5 Responses to “Pythagorean Quadruple”

1. Zack said

Interesting problem. Here is my take on it using Julia: https://pastebin.com/eYYR51sm

For N = 1000 I found 85490 quadruples (allowing for multiple instances of the same number in the quadruple set).

Cheers

2. matthew said

Here’s a C++ solution that runs in O(n²) time, using O(n²) space (assuming hash table access time is constant):

```#include <vector>
#include <utility>
#include <iostream>
#include <cstdio>

// Find all "pythagorean quadruples" (a,b,c,d) with
// 1 <= a <= b <= c <= N and
// a^2 + b^2 + c^2 = d^2

// Meet-in-the-middle: build table of all possible ways
// of making a*a+b*b, then search all possible values
// of d*d-c*c

int main(int argc, char *argv[]) {
int n = strtoul(argv,0,0);
std::unordered_map<int,std::vector<std::pair<int,int>>> m;
for (int a = 1; a <= n; a++) {
for (int b = a; b <= n; b++) {
m[a*a+b*b].push_back(std::make_pair(a,b));
}
}
// sqrt(3c^2) = sqrt(3)*c = 2c, near enough
for (int c = 1; c <= n; c++) {
for (int d = c+1; d <= 2*c; d++) {
for (auto p: m[d*d-c*c]) {
auto a = p.first, b = p.second;
if (c >= p.second) {
std::cout << a << b << c << d << "\n";
}
}
}
}
}
```
```\$ g++ -O3 -Wall pquads.cpp -o pquads
\$ time ./pquads 1000 | wc -l
85490

real	0m0.468s
user	0m0.443s
sys	0m0.032s
```
3. matthew said

Just in case overflow in all those squares is a problem, I thought I’d see what gcc can do to check things at runtime:

```\$ g++ -O3 -Wall -fsanitize=signed-integer-overflow pquads.cpp -o pquads
\$ time ./pquads 1000 | wc -l
85490

real 0m0.447s
user 0m0.426s
sys  0m0.029s
\$ time ./pquads 1000 | wc -l
85490

real	0m0.651s
user 0m0.627s
sys  0m0.032s
```
4. matthew said

Here’s a slightly improved version (that prints spaces between the values for a,b,c,d):

```int main(int argc, char *argv[]) {
int n = strtoul(argv,0,0);
std::unordered_map<int,std::vector<std::pair<int,int>>> m;
for (int a = 1; a <= n; a++) {
for (int b = a; b <= n; b++) {
m[a*a+b*b].emplace_back(a,b);
}
}
for (int c = 1; c <= n; c++) {
for (int d = c+1; d <= 2*c; d++) {
for (auto p: m[d*d-c*c]) {
auto a = p.first, b = p.second;
if (b <= c) std::cout << a << " " << b << " " << c << " " << d << "\n";
}
}
}
}
```
5. Daniel said

Here’s a solution in Python.

```from collections import defaultdict
from itertools import count

bc_lookup = defaultdict(list)
for b in range(1, N + 1):
for c in range(b, N + 1):
bc_lookup[b * b + c * c].append((b, c))
output = 0
for a in range(1, N + 1):
for d in count(a):
if d * d - a * a > 2 * N * N:
break
for b, c in bc_lookup.get(d * d - a * a, []):
if a <= b:
output += 1
return output

```85490