M-8: Nice Parameters Case Study

Our final case study has to do with answering a question about number theory (which is disguised as a 3D geometry problem). A trirectangular tetrahedron is a tetrahedron ABCD so that all three of the faces meeting at A are perpendicular. It is the 3D equivalent of a right triangle: the definition is equivalent to “all of the triangles ABC, ACD, and ADB are right-angled triangles with the right angle at A.” Our question is,

Are there trirectangular tetrahedra with all side lengths integral?

In your own work, you might come across similar situations when you want to make sure that a question you’re creating has “good” parameters so that the numbers all come out nice for the student at the end.

Phase 0: Brute Force

Let’s use the Pythagorean theorem and the type function in Maple to do a brute force search. What we mean about the Pythagorean theorem is the following: we want that the lengths of i := |AB|, j := |AC|, and k := |AD| are integers, and also that all of $|BC| = \sqrt{i^2+j^2}, |BD| = \sqrt{i^2+k^2}, |CD| = \sqrt{j^2+k^2}$ are integers. So let’s simply test all possible values of i, j and k (up to some limit) to see which ones make this true.

We need to use the type procedure to check for integrality: type(x, integer) returns true or false depending on whether x is integral. E.g.,

type(sqrt(9), integer); # returns true

It makes our code a bit easier to read if we put aside the square/sqrt/typechecking stuff into a single procedure. Define isCompatible(x, y) to be true if $\sqrt{x^2+y^2}$ is an integer:

isCompatible := proc(x, y):
   return type(sqrt(x*x+y*y), integer);
end:

Then, our basic search for triples where $i, j, k \le n$ where n is some maximum search limit we specify can be accomplished with:

findTetras0 := proc(n) local i, j, k:
   for i from 1 to n do
      for j from 1 to n do
         for k from 1 to n do
            if isCompatible(i, j) and isCompatible(j, k) 
             and isCompatible(i, k) then
                print(i, j, k):
            end:
         end:
      end:
   end:
end:

This very basic version of the code takes about $n^3$ operations and therefore its running time is approximately some constant times $n^3$. For my computer, findTetras0(50) takes about 2 seconds, and findTetras0(100) takes about 16 seconds. This does not produce any output, unfortunately. But we shouldn’t give up hope yet.

Phase 1: Be Polite

Thinking about our code, can we reduce the amount of work that it does? Yes! Notice that if i and j are not compatible, then there is no need to perform the innermost for k loop, since none of the pairs inside of that innermost loop could result in an (i, j, k) of the desired form. So let’s not ask Maple to do computations that we are sure will reach dead ends. Hence we edit the code to:

findTetras1 := proc(n) local i, j, k:
   for i from 1 to n do  
      for j from 1 to n do
         if isCompatible(i, j) then 
            for k from 1 to n do
               if isCompatible(j, k) and isCompatible(i, k) then
                  print(i, j, k):
               end:
            end:
         end:
      end:
   end:
end:

This is much faster: on the same machine, findTetras1(100) takes less than 10 seconds, and findTetras1(500) takes less than a minute. (It’s harder to precisely estimate the running time since it depends on how many compatible (i, j) pairs are found.)

More importantly, this version of the code does find some results! However, you’ll notice that it prints out multiple versions of each triple, for example it prints out both 44, 117, 240 and also 117, 44, 240 etc. So let’s only make the code check triples with $i \le j \le k$ by setting the loop lower bounds to for j from i… and for k from j… After we do that, the code will give us data that looks like this:

The next few phases won’t give dramatic improvements in running speed, but will make the code more modular and put us in a place to get a very fast algorithm at the end once we use a mathematical idea. Along the way we’ll introduce a few more detailed ideas on writing code that runs more quickly.

One thing we can notice immediately is that there are still some “redundancy” in our output in the sense that some entries are multiples of one another: twice 44, 117, 240 is the same as 88, 234, 480. It is easy to convince yourself that whenever (i, j, k) satisfy the properties we’re looking for, so does (ti, tj, tk) for all integers t ≥ 1. So let’s only print out primitive triples, where there is no common factor to i, j, and k. Maple has a gcd function and while it only takes in two arguments, you can see that gcd(gcd(i, j), k) is the common gcd of all three numbers. So in our next version, we’ll only print out triples where this value is 1.

We’ll say that the rest of this lesson gets fairly technical. It would already be quite impressive, given that this is just a beginner course, if you were able to write your own code that worked without big redundancies (i.e., getting as far as our Phase 1 here). But if you are interested in knowing more, continue reading below.

Phase 2: Pre-processing

There is a nice idiom in Maple which lets you return a list or set when you don’t know in advance what length it will have. The idea is to initialize a value to {}, the empty set. From then on, you will replace that set by its union with a singleton set. For example,

mySet := {};        
for i from 1 to 10 do
   mySet := mySet union {i};
end;
mySet;

sets mySet equal to the set {1, 2, …, 10}. Try this yourself and play around with it to see if it works like you expect.

The next major idea that our program evolution will use is the following: by doing some pre-processing, we might be able to save repeated work later on. In particular, let’s pre-compute a list of all of the compatible pairs. Specifically, for each i, let’s pre-compute all of the other values which are compatible with i. Having such a list could save us time since instead of searching through all j and k, we’ll only have to look at those j which are compatible with i, and those k that are compatible with i as well as j. By using the option remember concept introduced earlier, this can be written in the following way:

compatibleWith := proc(k, n) option remember:
   # returns a set of all i s.t. 1<=i<=n and i*i+k*k is a square
   local i, result:
   result := {};  # the empty set
   for i from 1 to n do
      if isCompatible(i, k) then result := result union {i}: end:
   end:
   return result:
end:

For example, compatibleWith(40, 100) will return {9, 30, 42, 75, 96} since these are all of the integers ≤ 100 which, together with 40, can form the legs of a Pythagorean triangle. Now we rewrite our code so that instead of searching all possible j and k, it only searches the compatible ones:

findTetras2 := proc(n)
   local i, jIndex, kIndex, cwi, j, k, numcwi:
   for i from 1 to n do
      # cwi means "compatible with i" 
      cwi := compatibleWith(i, n):
      numcwi := nops(cwi);
      for jIndex from 1 to numcwi do
         j := cwi[jIndex]:
         if (j>=i) then
            for kIndex from 1 to numcwi do
               k := cwi[kIndex]:
               if (k>=j) and gcd(gcd(i, j), k)=1 and
                isCompatible(j, k) then
                   print(i, j, k):
               end:
            end:
         end:
      end:
   end:
end:

Notice that «set»[i] works similarly to how it does for lists, and that the inner for loop variables go over jIndex and kIndex, with j and k taking on all of the values in cwi in turn. This version of the code, unfortunately, isn’t all that faster despite that the pre-processing is really looking at the problem in a novel way.

Phase 3: Cleaning Up Our Code

Maple has another kind of for loop: for «loopVar» in «structure» works when the structure is either a list or a set. It iterates the loop variable directly over the elements of the structure… that’s exactly what we’ve been doing with the above example, but we can get rid of jIndex and kIndex. So for these edits, for example, we’d change the j loop to

for j in cwi do   # for all j that are compatible with i
   ...            # do stuff with j, no jIndex necessary

Another thing which might be annoying is the fact that there are so many lines of input. Let’s make our procedure stop printing things and instead return a set of all of the working {i, j, k} triples that we found. As before, this means initializing result := {} and taking unions.

Additionally, the for k loop could possibly be implemented in a more symmetric way. The current code ensures that i and k are compatible by forcing k to be drawn from cwi (the precomputed list compatibleWith(i, n) of integers compatible with i), while it is explicitly testing the compatibility of j and k using isCompatible. Instead, let’s use the former precomputation method for both i and j. We basically want to restrict k to the intersection of compatibleWith(i, n) and compatibleWith(j, n). Using intersect, this gives the following evolved version of our code:

findTetras3 := proc(n) local i, cwi, j, k, result:
   result := {};
   for i from 1 to n do
      # cwi means "compatible with i"
      cwi := compatibleWith(i, n):
      for j in cwi do
         if (j>=i) then
            for k in cwi intersect compatibleWith(j, n) do
               if (k>=j) and gcd(gcd(i, j), k)=1 then
                  result := result union {{i, j, k}}:
               end:
            end:
         end:
      end:
   end:
   return result;
end:

The speed is about the same as before, but it’s a little easier on the eyes and has simpler logic.

Phase 4: Number Theory and `Array`s

The main thing that will get us to a blazingly fast version of the code is the following fact about Pythagorean triples:

Every Pythagorean triple can be expressed in the form $\{k(2ab), k(a^2-b^2), k(a^2+b^2)\}$ for some positive integers k, a, and b where a and b are relatively prime. The hypotenuse is $k(a^2+b^2)$ and the other two are the legs.

This parameterization, called Euclid’s formula, is very useful because it means that instead of searching for compatible pairs and missing a lot of the time, we can generate all of them and only get the hits.

We want our code to search through all values of a, b, and k, and then update our “compatibility sets” for i and j where $\{i, j\} = \{2kab, k(a^2-b^2)\}$. But how can we store the data if we search in this way rather than with for i and for j? We want to iteratively fill out a data structure whose [i]th element is the set of all numbers which are compatible with i, similar to before.

We can do this in principle with lists, except that we want to handle more than 100 elements so we’ll use a table (see the last lesson). This is what building the compatibility table looks like; note that it uses a hybrid for–while loop (see ?for):

compatibilityTable := proc(n)
   # returns a table indexed from 1 to n, where 
   # the item at index i is the set of all j<=n such that
   # i*i+j*j is a square
   local result, a, b, k, P, Q, i;
   result := table():  
   for i from 1 to n do
      result[i] := {}:
   end:  
   for b from 1 to n do
      for a from b+1 to n while max(2*a*b, a*a-b*b)<=n do
         if gcd(a, b)=1 then 
            for k from 1 to floor(n/max(2*a*b, a*a-b*b)) do
               P := k*2*a*b;
               Q := k*(a*a-b*b);
               result[P] := result[P] union {Q};
               result[Q] := result[Q] union {P};
            end:
         end:
      end:
   end:
   return result;
end:

Note that we use the floor function to get the search to stop at an integral value of k. Finally, the main procedure is the following — it’s the same as before except it uses compatibilityTable instead of compatibleWith.

findTetras4 := proc(n)
   local i, j, k, CT, result:
   result := {};
   CT := compatibilityTable(n):
   for i from 1 to n do
      for j in CT[i] do
         if (j>=i) then
            for k in CT[i] intersect CT[j] do
               if (k>=j) and gcd(gcd(i, j), k)=1 then
                  result := result union {{i, j, k}}:
               end:
            end:
         end:
      end:
   end:
   return result;
end:

This version is ridiculously fast!

If you want, you can download an entire .mw file with all of these commands.

In summary,

start with a simple version of working code before making it complicated,
optimize in small steps,
keep your code idiomatic,
see if you can use math to help you,
and use data structures if you really need to.

We should also caution that it’s not always worthwhile to try to squeeze out every last drop of optimization from your routines. For example in a recursive function, maybe you can either treat the n=2 case like all of the rest, or alternatively come up with a hand-coded optimized version for that special case, by moving some if and for statements around and copying most of the body of your procedure. Avoid the temptation to do the latter! In general, trying to include tiny improvements that save a single step of wasted computation at the cost of making your code bloated is dangerous: it is not likely to help that much, but is likely to make your code more complex and error-prone. (Perhaps the only exception is inside of a deeply nested loop, but even then you should try to rigorously compare your “improvements” with the old version to see if they really help.)

Another thing to keep in mind is that debugging code is both an art and a science. Your code seems like it’s running forever? Try it on small cases and see if it’s totally broken, or just slow. Are some variables taking on unexpected values? Use print statements to inspect things more deeply. If you are extremely hardcore, check out ?debugger to interactively debug your code. Keep your code human-readable by using proper indentation and re-read the code frequently to see that you typed what you really meant to type.

Good luck coding with Maple and have fun!