Optimizing Conway's Game of Life in C++

How could I further optimize my implementation of Conway's Game of Life? And how would you critique my current strategies? I'm taking a C++ optimization class, the deadline has passed and my assignment has already been submitted. We were supposed to compile with these options: g++ me.cpp -std=c++11 -O3 -march=native -o me.

I'm trying to minimize real time after executing it like this: time ./me 0 1000 0. My program's first argument is the random seed, the second is the number of iterations, and the last is 1 or 0 meaning print or don't print. I don't want to print while timing. The random seed and number of iterations could be any value.

How it works

I iterate through living cells instead of all spaces on the board. I manually unrolled iterating through all 9 tiles encompassing each living cell. Living neighbors, including the central living cell, become incremented positively. Dead neighbors of a living cell become decremented.

That means after the encode function a value of -3 on a dead tile means there are 3 living neighbors, therefore that tile should become alive. 1 is not possible because living cells are incremented at least once. Positive 2 means no living neighbors, 3 means 1 living neighbor, 4 means 2 living neighbors, 5 means 3 living neighbors, etc.

According to the rules there are 3 circumstances where a cell can be alive, but my living_cipher array includes a 4th living state. That's because within the decode method a code value of 1 means I've already decided the cell is alive and now I should ignore that tile. Because 1 should be ignored, the counter_cipher has a zero at that entry. This method transforms codes ranging from -8 to 10, to either a 1 or a 0 meaning living or dead. Rise and repeat for the desired number of iterations.

I used macros as an optimization because constant variables were slower. I did not try an inline function.

#include <iostream>

using namespace std;

#define HEIGHT ((1 << 10) + 2)
#define WIDTH ((1 << 11) + 2)
#define MAX_Y (HEIGHT - 1)
#define MAX_X (WIDTH - 1)
#define NEXT_ROW (WIDTH - 2)
#define AREA (HEIGHT * WIDTH)

constexpr signed char offset_living_cipher[19] = {
0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0
};
constexpr signed char offset_counter_cipher[19] = {
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0
};

class Matthew_Conway{
private:
bool print;
signed char * board;
signed char ** cells, ** next_cells;
unsigned long total_cells;

void make_cells(){
unsigned long next_total_cells = 0;
cells = new signed char * [AREA];
next_cells = new signed char * [AREA];

for (unsigned long row = 1; row <= HEIGHT - 2; ++row){
for(unsigned long column = 1; column <= WIDTH - 2; ++column){
unsigned long index = row * WIDTH + column;
signed char cell = rand() % 2 == 0;
board[index] = cell;
cells[next_total_cells] = board + index;
next_total_cells += cell;
}
}
total_cells = next_total_cells;
}

public:
unsigned long iterations;

Matthew_Conway(char ** argv){
unsigned random_seed = (unsigned) atoi(argv[1]);
srand(random_seed);
iterations = (unsigned long) atoi(argv[2]);
print = (bool) atoi(argv[3]);
board = (signed char*) calloc(AREA, sizeof(signed char));
make_cells();
}

void print_board(){
if (!print)
return;
for (unsigned long row = 1; row <= HEIGHT - 2; ++row){
for (unsigned long column = 1; column <= WIDTH - 2; ++column){
signed char value = board[row * WIDTH + column];
cout << (int) value;
}
cout << endl;
}
}

void encode() {
unsigned long next_total_cells = 0;

for (unsigned long cell_i = 0; cell_i < total_cells; ++cell_i) {
signed char * cell = cells[next_total_cells] = cells[cell_i];
unsigned long index = cell - board;
unsigned long y = index / WIDTH;
unsigned long x = index % WIDTH;
if (x == 0 || x == MAX_X || y == 0 || y == MAX_Y){
continue;
}
++next_total_cells;

signed char * board_address = cell - WIDTH - 1; // Upper Left
INCREMENT;
INCREMENT;
INCREMENT;
board_address += NEXT_ROW; // Middle Left
INCREMENT;
INCREMENT;
board_address += NEXT_ROW; // Lower Left
INCREMENT;
INCREMENT;
INCREMENT;
}
total_cells = next_total_cells;
}

void decode() {
#define LIVING_CIPHER (offset_living_cipher + 8)
#define COUNTER_CIPHER (offset_counter_cipher + 8)
#define DECIPHER { \
next_total_cells += COUNTER_CIPHER[code]; \
}
signed char code;
unsigned long next_total_cells = 0;

for (unsigned long cell_i = 0; cell_i < total_cells; ++cell_i) {
signed char * cell = cells[cell_i];

signed char * board_address = cell - WIDTH - 1; // Upper Left
DECIPHER;
DECIPHER;
DECIPHER;
board_address += NEXT_ROW; // Middle Left
DECIPHER;
DECIPHER;
DECIPHER;
board_address += NEXT_ROW; // Lower Left
DECIPHER;
DECIPHER;
DECIPHER;
}
total_cells = next_total_cells;
swap(cells, next_cells);
}
};

int main(int argc, char ** argv){
if (argc != 4){
cout << "usage game-of-life <seed> <generations> <0:don't print, 1:print>" << endl;
return 1;
}
Matthew_Conway conway(argv);
conway.print_board();
for (unsigned long i = 0; i < conway.iterations; ++i){
conway.encode();
conway.decode();
}
conway.print_board();
return 0;
}

• How well does the following statement fit to this assignment and the review you are interested in? "I don't care how terrible the code is. As long as it executes fast and spits out the right answer I'm happy." I have a feeling answering this will help avoiding people giving you reviews you don't want.
– nwp
Oct 9, 2019 at 9:51
• @nwp That's true with regard to the assignment. Absolute speed at any cost except correctness of output. But I'm not sure what optimizations are never appropriate even when every bit of performance is desired. And to cloud things further, I'm not sure if there's a gray area when certain methods are reserved for the most extreme applications.
– Matt
Oct 9, 2019 at 10:45
• In addition to my answer below, profile the code, this will tell you where the bottle necks are. Keep in mind the first rule of optimization is don't optimize, the the compiler do it for you. The second rule of optimization is generally find a faster algorithm. Oct 9, 2019 at 16:44
• Having the compilation flags is nice, but which architecture are you targeting? Where will the program run for the measurements? Unless that information is given, we can at most give you some superficial answers. Oct 9, 2019 at 16:53
• @hoffmale Our assignment is graded on a Xeon chip of unknown specification, but I would prefer to optimize my programs for my personal i7-8565U since I can become most familiar it. ark.intel.com/content/www/us/en/ark/products/149091/…
– Matt
Oct 9, 2019 at 19:13

The code may be over optimized already and adding additional optimizations will only make it harder to read and maintain the code, but that said:

Missing Optimizations

• When printing prefer '\n' over std::endl because std::endl calls a system function to flush the output, '\n' just inserts a new line.
• Put the check to print outside the function print_board() so that the function isn't even called. This saves the time cost of pushing everything onto the stack and altering the flow of execution.
• Use iterators or pointers in print_board() and make_cells() rather than indexing. Compiling with -O3 may do this for you, but direct addressing is going to be faster than indexing.
• Depending on the processor and whether the instruction set includes an auto decrement and test instruction performance may be improved by counting down in for loops rather than counting up.
• Check the generated assembly code and know the instruction sets for the processor being used.
• Use the natural word size of the processor rather than specifying a word length. In the code this mean use unsigned rather than unsigned long. On most current processors you will be guaranteed a word size of at least 32 bits and probably a word size of 64 bits.
• Make member functions that don't change member values const.

void print_board() const
{
...
}


Magic Numbers

There are a few numeric constants that should be removed or changed to symbolic constants. Numeric constants are sometimes referred to as Magic Numbers because it's not clear what the values represent. There is a discussion of this on stackoverflow.

Examples would be the value 4 in main in the comparison of argc and the value 19 in these two declarations:

constexpr signed char offset_living_cipher[19] =
{
0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0
};

constexpr signed char offset_counter_cipher[19] =
{
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0
};


The number 19 in the previous two declarations is not necessary since these will be calculated by the compiler.

constexpr signed char offset_living_cipher[] =
{
0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0
};

constexpr signed char offset_counter_cipher[] =
{
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0
};


It might also be good to comment these declarations since it's not clear what is being done here.

It might be better to include <cstdlib> and use the system-defined symbolic constants EXIT_FAILURE and EXIT_SUCCESS in the return statements in main.

Don't Mix Memory Allocation Methodologies

In the constructor there is a call to the C programming memory allocation function calloc(), but in the function make_cells() there are two calls to new. There are a couple of problems here, the first is that the results of the call to calloc() are not checked. If the array allocated by calloc() isn't allocated (calloc returns NULL or nullptr) the first use of board will cause the program to terminate, however, if new fails it throws an exception. In this case since there are no try {} catch {} blocks the program will also terminate, but if there was a try {} catch {} block the program would still terminate if the calloc() failed. It would be best to be consistent with the memory allocation, and it would probably be better to stick with new rather than calloc().

Avoid Using Namespace std

If you are coding professionally you probably should get out of the habit of using the using namespace std; directive. The code will more clearly define where cout and other identifiers are coming from (std::cin, std::cout). As you start using namespaces in your code, it is better to identify where each function comes from because there may be function name collisions from different namespaces. The object cout you may override within your own classes. This Stack Overflow question discusses this in more detail.

• Using '\n' instead of std::endl will probably not change a thing, as '\n' usually flushes the ostream as well. To avoid unnecessary flushing, you must use '\n', but also specify std::cout.sync_with_stdio(false). For reference, see en.cppreference.com/w/cpp/io/manip/endl , Notes, §2 Oct 11, 2019 at 15:51

Coding for performance doesn't mean that you should write unmaintainable code. Remember that developer time is valuable, too: any time you spend tracking down avoidable bugs is time you could have spent tweaking program performance. Some particular things that can be improved without impacting performance:

• Include the headers for standard identifiers we use (std::swap requires <utility>; std::atoi and std::calloc require <cstdlib>).
• Avoid using namespace std;
• Prefer properly-typed constexpr values instead of preprocessor macros (I don't believe your claim that this made your code slower - I get the exact same assembly code with that change). These values should be (private) static members, not globals.
• Prefer C++ memory allocation to <cstdlib> allocation (malloc() and family). That makes the error checking much simpler (exceptions rather than return-value checking).
• Use a std::vector rather than raw new[]; that stops us leaking memory.
• Keep the argument parsing out of the core logic, and make it more robust (e.g. use std::stoul() for converting to unsigned long).
• Initialise members in the constructor's initializer list.
• If you really need short-term macros, then #undef them after use.

Once we have maintainable code, we can work on the algorithm. We're very inefficient with storage (making for poor data locality, and thus, thrashing cache unnecessarily). We have two arrays of pointers, each the same dimension as board, but much bigger on platforms where sizeof (unsigned char *) is greater than 1. The pointers all point into board, so we could store indexes instead; better still, store separate x and y coordinates, so we don't have to do divisions.

Given that our task is to minimize real time, we should look to parallelize as much as we can. On multi-core processors, Game of Life naturally parallelizes with each thread writing its own part of the next state, and with all threads sharing read access to the old state. SIMD parallelisation within each thread is possible, too, if your processor has support (e.g. NEON or AVX).

Exploiting sparseness is a nice trick, though not as good as Hashlife (it's in a class of its own). Hashlife aside, there are simple approaches that work well on actual computers, not by being algorithmically clever but by being computer-friendly. Compared to a a completely plain approach, yours is definitely better. But it does not really tap into the potential of a modern CPU.

For example on my PC (using a 4770K, so not very new any more), your benchmark takes around 2.85 seconds, but if I use some really simple no-tricks AVX2 implementation (so not using sparseness or anything else clever) of Game of Life, that takes around 0.38 seconds. That's around 7.5x faster, even though a lot more work is done. But work isn't time, and SIMD is great at crunching through a lot of work in a short time especially if you're working with 8-bit types.

In a benchmark with a much sparser grid, the SIMD approach wouldn't do so well. Even in this benchmark, simulating more steps decreases the speedup, going down to ~3.7x for 100000 steps.

AVX2 is an instruction set extension available in many x86 CPUs today, for example Haswell and newer Intel processors (including i7-8565U, but the Mystery Xeon CPU in the grading server may or may not support AVX2 depending on how old it is), and AMD Excavator and Zen. The instructions are exposed in C++ as intrinsic functions if you use a compiler that supports them. AVX2 includes instructions that apply an operation to 32 bytes at once, for example in the code below I used _mm256_add_epi8 (aka vpaddb) several times, each performing 32 additions at once, to sum up the number of live neighbours for a row of 32 cells simultaneously.

Example code:

void step() {
#define INDEX(a, b) ((a) + (b) * WIDTH)
for (size_t y = 0; y < HEIGHT; y++) {
if (y == 0 || y == MAX_Y) {
for (size_t x = 0; x < WIDTH; x++) {
next_board[INDEX(x, y)] = board[INDEX(x, y)];
}
}
else {
next_board[INDEX(0, y)] = board[INDEX(0, y)];
next_board[INDEX(MAX_X, y)] = board[INDEX(MAX_X, y)];
size_t x = 1;
for (; x < WIDTH - 33; x += 32) {
__m256i n = _mm256_loadu_si256((__m256i*)&board[INDEX(x - 1, y - 1)]);
__m256i is3 = _mm256_cmpeq_epi8(n, _mm256_set1_epi8(3));
__m256i is2 = _mm256_cmpeq_epi8(n, _mm256_set1_epi8(2));
__m256i isNextLive = _mm256_or_si256(is3, _mm256_and_si256(is2, cellItself));
isNextLive = _mm256_abs_epi8(isNextLive);
_mm256_storeu_si256((__m256i*)&next_board[INDEX(x, y)], isNextLive);
}
for (; x < WIDTH - 1; x++) {
int n = board[INDEX(x - 1, y - 1)];
n += board[INDEX(x, y - 1)];
n += board[INDEX(x + 1, y - 1)];
n += board[INDEX(x - 1, y)];
n += board[INDEX(x + 1, y)];
n += board[INDEX(x - 1, y + 1)];
n += board[INDEX(x, y + 1)];
n += board[INDEX(x + 1, y + 1)];
next_board[INDEX(x, y)] = n == 3 || (n == 2 && board[INDEX(x, y)]);
}
}
}

swap(board, next_board);
}


There are more clever uses of AVX2, for example using vpshufb to implement the automaton rule.

An other popular approach is bit-packing the cells, and then using bitwise operations to calculate the next state, relying on the bit-parallel nature of bitwise operations to speed up the computation. This page discusses that trick and some related history.

• forgive me for being naive, but could you explain AVX2. Oct 9, 2019 at 19:28
• @pacmaninbw in general or just what I used here? Oct 9, 2019 at 19:32
• In general and what you used here, add it to the answer. Oct 9, 2019 at 19:34