Primes, parallelism and caches Part I

Our first excursion is going to be into using the parallel processing capabilities of a modern PC. To have an interesting problem to work with, I will make programs to count the number of prime numbers less than a specified big number, like 100,000,000,000 (100 billions, or 10¹¹). Which, by the way, is 4,118,054,813.¹

The basis for this algorithm is Eratosthenes’ Sieve, which was a brilliant idea thought out many thousands years ago.²

Imagine a sieve that can sift numbers. This sieve will let all non-prime numbers fall out, and keep the primes, those numbers that can only be divided by 1 or by itself. If you drop a number in there, the sieve can determine if it’s a prime or not, and if you drop a whole range of numbers, you can count how many primes there are in that range.

To prepare the sieve, the prime numbers are used to punch out the holes in the sieve, punching out all multiples of the prime number. If you start with the first prime you know – 2 – you can look for the next prime, the one after 2, that was not punched out. And carry on in that manner until you have a sieve that is large enough for your purpose.

The algorithm is shown here in C++ code, I will explain it step by step.

itype numPrimes(itype maxNum)
{
	auto numbers = new char[maxNum + 1] {0};
	numbers[0] = numbers[1] = 1;

	itype sqrtMaxNum = static_cast<itype>(std::sqrt(maxNum));

	itype p = 2;
	do
	{
		for (itype i = p * p; i <= maxNum; i += p)
			numbers[i] = 1;

		p += 1;
		while (p <= sqrtMaxNum && numbers[p] == 1)
			p++;

	} while (p <= sqrtMaxNum);

	itype result = 0;
	for (itype i = 2; i <= maxNum; i++)
		if (numbers[i] == 0)
			result++;

	delete [] numbers;
	return result;
}

Line 3 and 4 creates the model of the sieve as an array of bytes, exactly large enough to represent all the numbers, I want to look at. Since 0 and 1 are always non-prime, I didn’t really have to include them in the array, but doing it this way, the number under investigation will always have itself as the index into the array, which makes calculations easier, and I only waste 2 bytes. This is a common trick in programming.

I am using bytes to represent if there is a hole or not at a number position. This could actually have been presented by a 0 or a 1 – one bit. But accessing and manipulating computer memory at bit level is rather complicated and time consuming, so I have chosen the smallest easily and speedily accessible data element – the 8 bit byte. A modern computer has lots of memory, which is the main reason I don’t want to investigate if working at bit level could be faster. For this blog, other principles are more important.

I have chosen a 0 to represent no hole (yet) and a 1 for a hole. 0 has an advantage in C++ (inherited from C) when initializing data, so I can write a very compact expression, which will make sure the whole array is initialized to 0 in all positions, and let the compiler determine how to do that in the best way.

So after line 4, the numbers array will look like this:

Line 6 calculates the square root of the largest number we want to investigate, why this is relevant will follow in a few paragraphs. For the example in the images, the maxNum is 25 and the square root of that is 5.

In line 8 I select the first prime – 2- and I will then punch out holes in the sieve in lines 11 and 12 for all multiples of 2. I start with 4 (the p*p as the starting point I will explain in a minute). 2 itself should of course not be punched out. The principle looks like this:

That’s how the array will look after the first time lines 11 and 12 have been executed. This finishes the handling of the prime 2.

In lines 14-16 I now look for the next prime. I look for the next number that has not been punched out, and very quickly I arrive at 3, as you can see.

In line 18 I check if this new prime is relevant – must be less or equal to the above mentioned square root value, and since 3 is less than 5, I repeat the loop, and execute lines 11 and 12 again:

It is now obvious why I can start with 3*3 – 9 – as you can see all smaller multiples of 3 have already been taken care of, in this case 6.

Again I will look for the next prime in lines 14-16. After 2 steps I reach 5, and execute lines 11-12 again:

This time I can start with 5*5 – 25 – and the next multiple of 5 is 30, which is outside maxNum of 25, so I only punch a hole for 25. As you can see, all smaller multiples: 10, 15, 20 have been taken care of by the previous steps with 2 and 3.

The next found prime would be 7, but both the conditions in line 15 and in line 18 will tell us to stop punching more holes in the sieve. Why is that so?

According to this implementation of the algorithm, the first multiple of 7 to punch out should be 7*7 – 49 – and that number is already higher than maxNum, so there is no reason to continue.

Another way to assure yourself about this stop limit is to realize that if some number higher or equal to the square root of a number divides the same number, the quotient must be less or equal to the square root, and then we will already have found this as a divisor.

All that is left is to count the number of primes found, and this is done in lines 20-23, counting all the positions in the array that have not been punched out. For 25, we get that there is 9 primes less or equal to 25, which is the correct answer.

Below I have timed a program running the above algorithm on my AMD Ryzen 9 5950X based PC. I have compiled both with Microsoft’s C++ compiler, and the MinGW64 gcc based compiler, with full optimization. The times in columns 3 and 4 are in seconds.

Since I am working with big numbers, I must use the 64 bits versions of the compiler, to be able to handle a numbers array larger than 4G. The modern compilers will give me a really big array from the heap with no problems, as you can see, I am working with up to 10¹¹ numbers. (Previous versions of MSVC only allowed 4Gbyte arrays even in 64 bit mode, but this seems to have been solved.)

As numbers grow, this algorithm seems to behave fairly economically, going up a factor of 10 in numbers only gives me sligthly more than a factor 10 in time. A primitive algorithm would scale with something like 100 – it’s complexity would be what is called O(n²), as it would try to divide all numbers with all lower numbers to check for primality. This algorithm uses less than the squareroot of n of those numbers, so I would expect a complexity of O(n^1.5) which is better than O(n²), but the numbers indicate something even better. The distribution of prime numbers apparently also helps a lot.

There are probably many more ways to try to optimize this algorithm, eg. ignoring all even numbers, since they cannot be primes. But this is not the aim of this post.

But what happened at 10¹¹? I let the algorithm run overnight, and its still hadn’t finished the next morning!

My PC has 32GBytes of RAM and most of these are accessible for a 64 bit program. But an array with 10¹¹ entries is more than 93 GBytes, so how can this run without crashing? With a 64 bit pointer and index, I can in principle handle a lot more than the 32GBytes, because of a mechanism called virtual memory. The operating system and the hardware will cooperate to make it seem as if I have all the desired memory, and it will use the RAM combined with the system hard disk (whether a spinning or a solid state kind) for this. The CPU can only work to and from RAM, so when another part of this address space is needed, it will have to be brought in from the system disk, and in the worst case, another page of RAM will have to be written to disk before the RAM is available, to save the data stored in that page.

This situation is called a page fault, and when a lot of these are happening all of the time, the PC will use more time to administrate the virtual memory scheme than to calculate my sieve. The following picture shows a screen dump from the performance application that tells us, that our PrimeSieve.exe program is generating 27,536 new page faults every second, and has already generated more than 2 billions of these.

This behavior for a system with virtual memory is called “threshing”.

I had expected that, so now it is fully demonstrated that my understanding of virtual memory is correct.

Given enough time, the program would eventually have finished its job. But not in a very effective way. This algorithm has very bad behavior using virtual memory, since for every new prime to use for punching out multiples, all of the big array has to be accessed in order, and since there is not enough room in RAM, a lot of page faults will happen, before it starts over with the next prime. This characteristic is called “low locality” and that is bad news for an algorithm on a computer with overcommitted memory allocation.

To contrast the situation, here is a snap shot of the screen when running the program looking just at 10¹⁰ numbers. Since I have so much RAM, normally more then 24 GBytes are free, and the 9.7 GBytes needed can easily be in there, meaning that very few page faults will arise when running the program.

Where did C++ take us this time?

I personally am amazed by how a modern computer can sift through many billions (GBytes) of data in such a short time. It was fun to make my computer work on such a big amount of data, at the same time demonstrating how a simple algorithm could get the correct results. The results aren’t news to the world, but that was not the purpose this time on my web blog.

I got a good understanding of Erathostenes’ Sieve, and this understanding will hopefully help me overcome the limitations I also found: raw power and raw size, even on my gigantic home computer, could only take me so far in counting primes.

I could demonstrate the problems with virtual memory: threshing can arise, if your algorithm is not localized.

The square root optimization trick is the basis for an idea to solve this: if I can look at the sieve in smaller sections, I only have to keep the primes less than the square root of my desired maximum number at all times. Thus I should be able to punch holes in one section at a time, keeping these small enough to stay in physical RAM, count the found primes, delete the section and go on to the next. It seems to me that this scheme could make me count up to the maximum of numbers a 64 bit computer can work on: 2⁶⁴ without problems – except I don’t now how long that would take! Mind bending! I only would need to keep the primes less than 2³² which must amount to less than 4GBytes, giving me plenty of RAM to look at the remaining numbers section by section.

This whole line of thought started when I wanted to investigate parallel computations. If you look at the screen dumps from the performance monitor you can see that the algorithm is only using 3% of the CPU power. My PC has 16 cores, each capable of running 2 threads of computations at the same time! It’s like having the power of 32 PC’s, incidentally running at 1000 times the speed of the original PC from 1981. Let’s put them to work! And with the section method, it seems possible!

That will be the step taken in the next post! Look for it!

About the code

The code can be downloaded from github. Feel free to use it however you want.

It is organized as a Visual Studio 2022 solution, but there is also a build.sh batch file to build it with gcc. In both cases a 64 bit compiler should be used, in VS select the Release configuration.

The main algorithm is encapsulated in the numPrimes() function at the bottom of the PrimeSieve.cpp file. The main() function takes care of setting everything up, including reading a parameter from the command line and timing the execution.

To parse the command line parameters I use a header only library by Matthias S. Benkmann which can be found here: https://optionparser.sourceforge.net/. He calls it the Lean Mean C++ Option Parser. It works fine, I would have liked a bit more automation when writing the help messages.

This time I have strived to be completely C++17 compatible, so the program can be compiled on several platforms. It won’t be like that every time.

1. See this article on Wikipedia on the number of primes. ↩︎
2. See all about Eratosthenes’ Sieve on Wikipedia here. ↩︎