Searching For Hypotenuses

October 9, 2015

A hypotenuse is a square that is the sum of two squares. So our first program generates the integers less than n and tests each to see if its square is the sum of two squares:

(define (hypotenuses n)
  (let loop ((h 1) (qs (list)) (hs (list)))
    (if (= h n) (reverse hs)
      (let* ((h2 (* h h)) (qs (cons h2 qs)))
        (if (hypotenuse? qs h2)
            (loop (+ h 1) qs (cons h hs))
            (loop (+ h 1) qs hs))))))

We determine if a number is the sum of two squares by starting from the two ends of the list of squares less than the number, discarding the low end if the sum is too small, discarding the high end if the sum is too large, and reporting success if the sum is equal. This is easy because we accumulate the squares as we go:

(define (hypotenuse? qs n2)
  (let loop ((rs (reverse qs)) (qs qs))
    (cond ((or (null? qs) (null? rs)) #f)
          ((= (+ (car qs) (car rs)) n2) #t)
          ((< (+ (car qs) (car rs)) n2)
            (loop (cdr rs) qs))
          (else (loop rs (cdr qs))))))

Here is a sample:

> (time (length (hypotenuses 10000)))
(time (length (hypotenuses 10000)))
    95 collections
    6302 ms elapsed cpu time, including 0 ms collecting
    6325 ms elapsed real time, including 18 ms collecting
    803350480 bytes allocated, including 800408704 bytes reclaimed
6448

That’s disappointingly slow. A better approach generates hypotenuses directly. We’ve discussed pythagorean triples in a previous exercise; here is a function to make a list of all primitive pythagorean triples with hypotenuse less than n:

(define (pyth n)
  (let loop ((a 3) (b 4) (c 5))
    (if (< n c) (list)
      (append
        (list (if (< a b) (list a b c) (list b a c)))
        (loop (+ a (- b) (- b) c c)
              (+ a a (- b) c c)
              (+ a a (- b) (- b) c c c))
        (loop (+ a b b c c)
              (+ a a b c c)
              (+ a a b b c c c))
        (loop (+ (- a) b b c c)
              (+ (- a) (- a) b c c)
              (+ (- a) (- a) b b c c c))))))

> (pyth 100)
((3 4 5) (5 12 13) (7 24 25) (9 40 41) (11 60 61) (13 84 85)
  (48 55 73) (28 45 53) (20 21 29) (39 80 89) (36 77 85)
  (8 15 17) (33 56 65) (65 72 97) (12 35 37) (16 63 65))

That’s a good start, but it ignores triangles like (6 8 10) that are multiples of primitive triangles. Here we extract the hypotenuses, remove duplicates,calculate all the multiples of each, and removeduplicates again:

(define (hypotenuses n)
  (unique = (sort <
    (mappend (lambda (z) (range z n z))
      (unique = (sort < (map caddr (pyth n))))))))

This is agreeably quick, about a second-and-a-half to find all the hypotenuses less than a million:

> (time (length (hypotenuses 10000)))
(time (length (hypotenuses 10000)))
    no collections
    0 ms elapsed cpu time
    6 ms elapsed real time
    4402736 bytes allocated
6448
> (time (length (hypotenuses 1000000)))
(time (length (hypotenuses 1000000)))
    54 collections
    1404 ms elapsed cpu time, including 328 ms collecting
    1442 ms elapsed real time, including 358 ms collecting
    473052144 bytes allocated, including 341499648 bytes reclaimed
721371

We used range, mappend, sort and unique from the Standard Prelude. You can run the program at http://ideone.com/B8zGvC.

Advertisements

Pages: 1 2

4 Responses to “Searching For Hypotenuses”

  1. Andras said

    Scala but it gives me one more:

      def multi(tr:Tuple3[Int, Int, Int], max:Int):List[Tuple3[Int, Int, Int]]={
      	(for(m <- 1 to max/tr._3) yield (tr._1*m,tr._2*m,tr._3*m)).toList
      }                                               
      
      def pyth(a: Int, b: Int, c: Int, max: Int): List[Tuple3[Int, Int, Int]] = {
        if (c > max) List()
        else multi((a, b, c),max) ++
          pyth(-a + b * 2 + c * 2, -a * 2 + b + c * 2, -a * 2 + b * 2 + c * 3, max) ++
          pyth(a + b * 2 + c * 2, a * 2 + b + c * 2, a * 2 + b * 2 + c * 3, max) ++
          pyth(a - b * 2 + c * 2, a * 2 - b + c * 2, a * 2 - b * 2 + c * 3, max)
      }                                               
    
      pyth(3,4,5,10)                                  //> res4: List[(Int, Int, Int)] = List((3,4,5), (6,8,10))
    
      pyth(3, 4, 5, 10*1000).map(_._3).distinct.size
                                                      //> res5: Int = 6449
    
  2. Paul said

    In Python.

    def pyth_triple(a, b, c):
        """generate 3 primitive triples from the primitive triple a, b, c
           a chain can be startedfrom 3,4,5
        """
        a2, b2, c2, c3 = 2 * a, 2 * b, 2 * c, 3 * c
        return ((-a + b2 + c2, -a2 + b + c2, -a2 + b2 + c3),
                (a + b2 + c2,  a2 + b + c2,  a2 + b2 + c3),
                (a - b2 + c2,  a2 - b + c2,  a2 - b2 + c3))
    
    def gen_triples(limit_hyp):
        Q = [(3, 4, 5)]
        while Q:
            a, b, hyp = Q.pop()
            if hyp < limit_hyp:
                yield hyp
                Q.extend(pyth_triple(a, b, hyp))
    
    MAX = 10**6
    unique = set()
    for hyp in gen_triples(MAX):
        if hyp not in unique:
            unique |= set(range(hyp, MAX, hyp)) # add all multiples of hyp
    print(len(unique))
    
  3. matthew said

    Here’s a simple sieve method: all hypotenuses are of the form k(i*i + j*j) for 0 < i 0, so we can just generate all such values under the limit and mark them off in a table:

    #include <vector>
    #include <numeric>
    #include <iostream>
    #include <stdlib.h>
    
    int main(int argc, char *argv[])
    {
      int N = atoi(argv[1]);
      std::vector<bool> a(N);
      for (int i = 1; i*i < N/2; i++) {
        for (int j = i+1; ; j++) {
          int t = i*i+j*j;
          if (t >= N) break;
          for (int k = t; k < N; k += t) {
            a[k] = true;
          }
        }
      }
      int count = std::accumulate(a.begin(),a.end(),0);
      std::cout << count << "\n";
    }
    
  4. matthew said

    Comment got a bit garbled there: I meant, for k > 0, and j > i > 0.

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: