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.

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
    $ 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

    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
    
    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 Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: