Pythagorean Quadruple
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.
Programming Praxis 
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
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"; } } } } }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:
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"; } } } }Here’s a solution in Python.
from collections import defaultdict from itertools import count def pythagorean_quadruples(N): 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 print(pythagorean_quadruples(1000))Output: