Pythagorean Quadruple

February 25, 2020

A Pythagorean quadruple consists of four positive integers a, b, c and d such that a ≤ b ≤ c ≤ d 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

Posted by programmingpraxis
Filed in Exercises
5 Comments »

5 Responses to “Pythagorean Quadruple”

  1. Zack said
    February 25, 2020 at 10:20 AM

    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
    February 26, 2020 at 12:07 AM

    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
    February 26, 2020 at 10:36 AM

    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
$ g++ -O3 -Wall -fsanitize=undefined pquads.cpp -o pquads
$ time ./pquads 1000 | wc -l
85490

real	0m0.651s
user 0m0.627s
sys  0m0.032s
  4. matthew said
    February 26, 2020 at 1:48 PM

    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
    March 18, 2020 at 1:13 AM

    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:

    85490

Leave a comment