### February 25, 2020

```(define (pyquad n)
(list-of (list a b c d)
(a range 1 (+ n 1))
(b range a (+ n 1))
(c range b (+ n 1))
(s is (+ (* a a) (* b b) (* c c)))
(d is (isqrt s))
(= (* d d) s)))```

That finds the 85490 Pythagorean quadruples less than 1000 in a little bit less than a minute and a half on my machine.

The time complexity of that function is clearly O(n³) due to the three nested loops on a, b and c. It seems to me that there ought to be an O(n²) solution by first computing a² + b² in O(n²) time, storing each of the (a b) pairs keyed by the sum of its squares, and comparing to d² – c², but I couldn’t find an appropriate method to limit c and d.

You can run the program at https://ideone.com/m3D0yF.

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[1],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[1],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