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:

- Find the first N primes
- 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.
- 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.)
- 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.
- 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.
- 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…
- 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!

- Find all the triples that have at least one prime.
- 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:
- A = 2*m*n
- B=m*m-n*n
- C=m*m+n*n

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

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

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.