A recent program contest at work challenged us to find as many Pythagorean triples as we could that contained at minimum of one prime. We had about a month to work on it. The first question was whether to write up a quick brute force attempt or to take some time to come up with a more elegant solution but have less time to run it. I of course chose the latter because it’s more fun. Here’s what I ended up doing:

1. Find the first N primes
   1. To do this, I calculate the all the primes from 1 to 10,000,000. This is done with the sieve of Eratosthenes. Those primes are stored in a bitmap. Each prime is represented by a single bit that is either 0 or 1 indicating not prime or prime.
   2. That bitmap is saved and used throughout the whole prime finding process. I use those primes to calculate the next primes in batches of 10,000,000 (again using a modified sieve of Eratosthenes.)
   3. The problem is that using one bit for each prime ends up being inefficient because there are a lot of numbers that are not prime. For a n numbers, roughly ln(n) of those numbers are prime. To save space, I implemented a compression scheme. I store the distance from one prime to the next. (Actually I store distance/2 since the distance between primes is always even (except for 2 and 3.)) This algorithm still wastes a lot of space (more on that later), but it’s about twice as good as a bitmap for space usage. Unfortunately it takes a long time to find out if a number is prime.
   4. For a quicker answer to “is X prime?”, I store bookmarks into the compressed array of primes. So I know that the number at index Y by looking it up in the array of bookmarks.
   5. The compression array is a series of UInt16. This compression uses less than half the space of a bitmap but could theoretically still be improved. A UInt16 is wasteful for many of the lower numbers. I wanted to implement a dynamic compression scheme that would take a chunk of numbers and decide how many bits to devote to the distance from one prime to the next. The number of bits per distance would vary by compression page. I never got around to this though…
   6. I calculate as many primes as I could based on the amount of RAM on the machine. After calculating the primes, I wrote the compressed array and the page bookmarks out to disk so that I could quickly restart the program without having to recalculate all the primes!
2. Find all the triples that have at least one prime.
   1. This was mostly brute force. There are two loops: M goes from 1 to sqrt of max prime and N goes from M to sqrt of max prime. For every value of M and N, a triple can be made with these formulas:
      1. A = 2*m*n
      2. B=m*m-n*n
      3. C=m*m+n*n
   2. 2 is always a factor of A so it is never prime. I just check B and C to see if they are prime. If so then I write the triple out to disk.
   3. This part took the longest but was easy to multithread. I used the new Parallel Extensions for .NET, and they were fantastically simple to use. The app automatically scales to the number of cores on the machine.

The only other thing to mention is that the command line arguments tell the app to either calculate primes or load them from the disk. The arguments also give a range of M for the app to examine for triples. This allowed it to be split over multiple machines easily. I wrote another small app that took all the output, cleaned it up, and combined it into a single file.

I ended up finding just under a billion triples that have at least one prime number. Since I was the only person who answered that question, I spent approximately 99.9% too much time on it.

You can download the app if you’re really curious to learn more and/or improve it.