So my last article on prime numbers led me to ask the question "How many prime numbers are there between 2 and 2³²". I know the quick answer to this question is "a lot". 2³² is over four billion—and that's a fair bit of space to check. It would be easy to make a for-loop and brute force the entire number range, but this is the age of multi-core CPUs—and single threaded is so 2000. Besides, my computer hasn't broken a sweat in awhile.

I decided to do a straight C implementation. I've been doing C and C++ at work, and C++ has been irritating me lately because I tried to do a "quick and dirty" implementation without too much thought into the classes I'd need. That's a recipe for disaster, and how you end up with things like microsoft's windows ME. So vanilla C it was.

Threads and semaphores I used the POSIX standards, and I've tried to demonstrate how to make a dispatcher foreground task with one or more worker threads. The key to this is the dispatch semaphore. Normally, semaphores are used to for mutually exclusive lock, or to send signals between threads. However, they can also be used as a count. Think of them like a library that has a set number (say 5) of books. That means 5 people can check that book out, but if a 6th person wants to check that book out, they will have to wait for one of the first 5 to return their copy. The semaphore does this by keeping a count. When that count reaches zero, the next task that tries to pend the semaphore will have to wait. So the dispatcher is a loop that starts by pending the semaphore. When it gets the semaphore, it will create a task. The task will do some work, then post the semaphore. The semaphore count determines how many tasks can run at the same time. If we limit the number of tasks that can run at once to the number of CPU cores, then we can utilize the CPU to 100%.

First, the dispatcher skeleton code:

Then, the skeleton of the worker thread:

Now just launching a thread to check one number is a bit of a waste—it takes some overhead to start and stop the task. So we give each thread some range of numbers to check. A pthread can be passed parameters, and we pass the starting number, along with the amount of numbers to check. The parameter block is also used for storage of the results. So here is the full worker thread:

The dispatcher has to accumulate the results. Our could try using a global variable, but you can run into issues with non-atomic access. Adding 1 might seem a basic operation, but (at least in a RISC system) it takes 3 instructions to add one to a variable in memory: load the data to a register, add the register by 1, and store the register. You can not (without a semaphore) guarantee that those instructions will not be interrupted by an other thread. So, the accumulation happens in just one place—the dispatch loop.

The dispatcher will assign all the work, but it is finished once all the work has been assigned. We still need to wait for the work to finish. So the last part of the process looks like this:

Now we can print the results to the screen, and we're done.

I knew this code would tax the CPU rather hard, so I wanted to measure how long the process would take. "time.h" is a C standard unit. Unfortunately, "time_spec", which has a high-resolution method for checking the time, it's part of the C99 standard. But we can get time in the resolution of seconds, and that's good enough. To do this, we simply mark the time when the process begins, and mark it again at the end. The difference (for which there is a function to compute) is the total number of seconds that has elapsed.

The answer: There are 203,280,221 primes between 2 and 4,294,967,295, and it took 6,054 seconds (1 hour, 41 minutes) on a quad-core clocked at 2.5 GHz.

The full source code is here.

So I mentioned yesterday that I wrote a simple program to check to see if some numbers in a list were prime or not. I didn't post it because I wanted to spend a little more time with it—giving it an article on it's own.

Finding prime numbers goes back a long time. The easiest way to see if a number is prime is to simply check and see if the number is divisible by all the prime numbers up to the square root of the number in question.

A number (call it x) is prime if there are no two number (call them a and b) such that a * b = x. All non-prime numbers can be expressed as the product of two or more prime numbers. For example, 125 can be made from 5 * 25, but 25 can be made from 5 * 5. So 5 * 5 * 5 = 125, and represents the most factored version of 125. This is true of any number. Since it takes at least two prime numbers to create a factor, we only need to check the primes up to √x (or the square root) of the number. This is because the if x = a * a, then √x = a. If x = a * b, and b is greater a, then a must be less then √x. Thus, we only need to check primes up to √x.

This makes for a nice brute-force test method, and it's pretty simple to implement. If we're going to test a number, we need all the primes up to the square root of that number. For what I needed, the number was capped at 32-bits—or 2³². √2³² = (2³²)^1/2 = 2¹⁶ = 65536. So I need a list of all the prime numbers up to 2¹⁶. We can make the list fairly quick. Have a look at this function:

We have an array to hold our prime numbers. We set the first entry in the array to two—the first prime number (1 is prime, but who cares). From this, we can test to check the rest of the numbers. We loop from 3 to 2¹⁶. For each number, we see if any of the primes already in the list will divide evenly into the number under question. If any of them do, they are not prime.

It happens that there are 6542 prime numbers between 2 and 2¹⁶. Knowing that, we can cap the array size of our lookup table at this value—the number of primes in this range will never change.

We now have enough prime numbers to check any number up to 2³². Unlike the function that generates the prime lookup table, we can throw in one more speed improvement: we only need to check up to the square root of the number in question. So here are the guts of the prime test function:

Note that can add two short cuts. First, checking up to the root of the number like we discussed. And second, we test to see if the number is even—2 is the only even number that is prime.

That's it. Now we have a simple to implement prime number test that works on numbers up to 2³². Here is the full source code to a command line version of the test. Compile it, pass it a number (or several), and the program will tell you if the number is prime or not.

(thanks Erica for the correction)

At work on Friday I came across an interesting math problem. I needed a pseudo random bit stream, and I needed the generator function to be fast. The typical random number function is a Linear congruential generator, most often f(x) = (1103515245 * x + 12345) mod 2³²-1. The modulus remainder is implied—the arithmetic is simply allowed to overflow. Multiplies are usually expensive in terms of CPU cycles. A faster generator is a linear feedback shift register (LFSR), which are comprised of shifts and XORes. I've used LFSR a number of times when I needed random data, without looking too much at the "random" stream of data it produced.

Using a 32-bit Galois LFSR with the polynomial x³² + x³¹ + x²⁹ + x + 1, and a seed of 0x12345678 produced this output:

If you haven't noticed, the pattern "scrolls" digits to the right. This is what one would expect—this is a "shift" register after all. Each number is different, and the period between any two number (except for zero) is 2³². However, the pattern between each successive number is fairly obvious. This is because of the polynomial. The shift register only introduces bits 31, 29, 1 and 0 when clocked. Otherwise, output is simply shifted to the right by one. I've arranged the data in 4 columns because each digit (nibble) is 4 bits, so after 4 bits, the nibble is completely shifted right by one.

While this pattern probably would work fine for what I was trying to do, I thought there must be a better solution. A polynomial with more taps should produce a better random stream. But how exactly does one go about generating a polynomial for a shift register?

What is need is a irreducible primitive polynomial. One can put any polynomial (i.e. any number as the polynomial term) into a LFSR they like, but it may not (probably won't) produce a maximum length shift register. That is, a 32-bit word has 2³² possibilities. With the correct polynomial, the output of a 32-bit LFSR will produce 2³² unique outputs before repeating. Without the correct polynomial, one of two things will happen. The output will happen after less then 2³² iterations. Or the output will never repeat and some smaller sub-sequence will repeat.

The polynomial of a LFSR is XORed in with the shifted data. My thought was, if I wanted a better random pattern, we'd better XOR more bits. And not just more bits, but spread out more or less evenly. This should change the stream of data quite a bit.

So how does one create an irreducible primitive polynomial? Good question. I still don't have a complete solution, but I do have part of the solution. I found a program written by a Scott Duplichan that generates primitive polynomials—although not necessarily irreducible. It's designed to generate huge primitive polynomials—on the order of thousands of bits. However, it also would do 32-bit. What was nice about this program was that I could tell it the "weight" I wanted. That is, how many bits I wanted set. The number has to be odd, so I chose to generate polynomials with 15 bits.

There are a lot of primitive polynomials with a weight of 15—millions. And most of them did not have bit patterns that were uniformly distributed. For example, one of the first polynomials produced was 0x00041FFF. While there are 15-bits set, there are pretty much all on one side. The program didn't have a way to tell it to generate polynomials that had uniform distribution, but it did have the ability to find a set number of polynomials. That was good enough for me. I set the program off to find ten million primitive polynomials, and went off to do something else. Just over an hour latter (4152 seconds actually), I had 10 million primitive polynomials with 15-bits set.

The next task was to find polynomials with a more uniform bit distribution. To do this, I wiped up a filter program to only save polynomials that had, at most, 3-bits of the same bits in a row. This significantly reduced the number of polynomials. But I still did not have a list of irreducible primitive polynomials. I started reading about how to test for this, but it was taking too much time. It would be easier to brute force this problem. I wanted a polynomial that was maximum length—that is, produced 2³² unique values before repeating. Any polynomial that did this was what I looking for.

It only takes a few seconds to brute-force check all 2³² possibilities for a LFSR. With the bit distribution restrictions, there were only 9,441 of the original ten million polynomials to check. So, I launched a brute-force attack, and went home for the weekend.

When I came back this morning, I had a list of 40 irreducible primitive polynomials. I'm not sure how long it took to complete this search, but that was irrelevant—I had my data.

Forty was still more polynomials then I needed. So I decided to reduce the list a little further. I just wanted prime polynomials. This is an other brute-force test. Any number is prime if all none of the prime number up to the square root of the number divide evenly into it. My list was reduced from 40 to 10. I picked the first number from the list, but any would have worked.

Again, my 32-bit Galois LFSR, this time with the polynomial 0x19253292B, and a seed of 0x12345678 produced this output:

And this looks much better. But just because you can't see a pattern doesn't mean one doesn't exist. In fact, one obviously exists—the data is from an equation. But without the equation, how good does the random data look? The quickest way to find out if data is random is to try and compress it. If it compresses, it isn't very random. For example, with 256 MBs of data from the original polynomial, we gets a stream that 7-zip could compress 54% of it's original size. With the new polynomial, 7-zip was only able to compress it to 99.964% of it's original size.

Compression isn't the best judge of a file's randomness. I knew there were better algorithms out for generating an entropy number. Doing a quick search, I found this program originally written by John Walker—one of the co-founders of AutoCAD. It dates back to 1985, but the nice thing about math (unlike computer hardware) is it still works fine after 25 years.

The program produces the chi-square distribution, which is given as percentage. When this percentage is close (or near) to 0% or %100, the data stream is said to not be very random—even 95% or 5% are suspect. The chi-square distribution of the first data set is 0.01%, or not even close to random. The chi-square distribution for the new polynomial is 38.58%. For comparison, the chi-square distribution of some output of the XTEA encryption algorithm in a feedback loop had chi-square distribution of 26.5%. Encryption algorithms desire to produce highly random output.

So our new polynomial of 0x9253292B produces a good psudo-random data stream when used in a 32-bit LFSR.