Beating the optimizer

It’s generally believed that it’s far too difficult to beat an optimizing C compiler when producing code, at least for a simple problem. In this post we will see how that’s not true and simple “hand-optimized” code can get large improvements in some cases.

Obviously that doesn’t mean that we should program using intrinsics or assembly most of the time: it still took me far more time to write this code than the time required to compile a C file with optimizations. But it means that, if we are working in a piece of time-critical code and profiling shows us a “hot spot”, it’s reasonable to consider doing such an optimization.

The problem

The problem description is very simple: given a big buffer (many megabytes), count the number of bytes equal to a given one. The prototype of such a function could be the following one:

size_t memcount(const void *s, int c, size_t n);


The baseline solution

Before doing any optimization is important to get a baseline solution and consider how optimal it is. In this case, it’s very easy to write:

size_t memcount_naive(const void *s, int c, size_t n)
{
const uint8_t *b = (const uint8_t *)s;
size_t count = 0;
for (size_t i = 0; i < n; i++)
if (b[i] == c)
count++;
return count;
}


Evaluating the baseline

We are going to do the evaluation with a relatively slow mobile Ivy Bridge core:

$grep 'processor\|$$model name$$' /proc/cpuinfo processor : 0 model name : Intel(R) Core(TM) i5-3210M CPU @ 2.50GHz processor : 1 model name : Intel(R) Core(TM) i5-3210M CPU @ 2.50GHz processor : 2 model name : Intel(R) Core(TM) i5-3210M CPU @ 2.50GHz processor : 3 model name : Intel(R) Core(TM) i5-3210M CPU @ 2.50GHz  Timing with clock_gettime(), we get the following results: $ gcc --version
gcc (Ubuntu 4.8.4-2ubuntu1~14.04) 4.8.4
Copyright (C) 2013 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
$gcc -march=native -std=c11 -DTEST_NAIVE -O3 memcount.c -o memcount && ./memcount Generated random arrays. Count of character 45 in 'a1' using 'memcount_naive' is: 408978 Elapsed time: 103946325 ns Count of character 45 in 'a1' using 'memcount_naive' is: 408978 Elapsed time: 104022735 ns Count of character 45 in 'a2' using 'memcount_naive' is: 409445 Elapsed time: 104278289 ns Count of character 45 in 'a3' using 'memcount_naive' is: 410593 Elapsed time: 104127038 ns Count of character 45 in 'a4' using 'memcount_naive' is: 104447433 Elapsed time: 104110155 ns Count of character 45 in 'a5' using 'memcount_naive' is: 104857613 Elapsed time: 104122294 ns  The measured time to count bytes with a given value in a ~100 MB array is about 100 ms, giving us a processing speed of approximately 1 GB/s. There are no substantial caching effects, as expected given that we are processing a 100 MB array. Can we do better? Memory bandwidth Let’s start by calculating the memory bandwidth. We can get the memory parameters by using the dmidecode command: $ sudo dmidecode --type 17 | grep '$$Locator$$\|$$Clock$$\|$$Size$$\|$$Type$$'
Size: 4096 MB
Locator: ChannelA-DIMM0
Bank Locator: BANK 0
Type: DDR3
Type Detail: Synchronous
Configured Clock Speed: 1333 MHz
Size: No Module Installed
Locator: ChannelA-DIMM1
Bank Locator: BANK 1
Type: Unknown
Type Detail: None
Configured Clock Speed: Unknown
Size: 4096 MB
Locator: ChannelB-DIMM0
Bank Locator: BANK 2
Type: DDR3
Type Detail: Synchronous
Configured Clock Speed: 1333 MHz
Size: No Module Installed
Locator: ChannelB-DIMM1
Bank Locator: BANK 3
Type: Unknown
Type Detail: None
Configured Clock Speed: Unknown


We have two DDR3 DIMMs in separate channels, running at “1333 MHz”. The expected bandwidth would be:

$\displaystyle 1333\,\mathrm{MT}/\mathrm{s} \cdot 128\,\mathrm{bit}/\mathrm{T} \cdot 0.125 \,\mathrm{bit}/\mathrm{byte} \approx 20\,\mathrm{GB}/\mathrm{s}$,

so we should have plenty of room for improvement. We can get a more realistic upper limit by comparing the performance to a highly optimized normal operation, such as memchr():

$gcc -march=native -std=c11 -DTEST_MEMCHR -O3 memcount.c -o memcount && ./memcount Generated random arrays. Timing 'memchr()' to get character 13' in 'a5'. Result: (nil) Elapsed time: 6475968 ns Timing 'memchr()' to get character 13' in 'a5'. Result: (nil) Elapsed time: 6919187 ns  Now we know that bandwidths of ~14 GB/s are achievable in this system, so there should be plenty of room for improving performance. Using SSE intrinsics My solution (that is probably far from optimal) is the following one: size_t memcount_sse(const void *s, int c, size_t n) { size_t nb = n / 16; size_t count = 0; __m128i ct = _mm_set1_epi32(c * ((~0UL) / 0xff)); __m128i z = _mm_set1_epi32(0); __m128i acr = _mm_set1_epi32(0); for (size_t i = 0; i < nb; i++) { __m128i b = _mm_lddqu_si128((const __m128i *)s + i); __m128i cr = _mm_cmpeq_epi8 (ct, b); acr = _mm_add_epi8(acr, cr); if (i % 0xff == 0xfe) { acr = _mm_sub_epi8(z, acr); __m128i sacr = _mm_sad_epu8(acr, z); count += _mm_extract_epi64(sacr, 0); count += _mm_extract_epi64(sacr, 1); acr = _mm_set1_epi32(0); } } acr = _mm_sub_epi8(z, acr); __m128i sacr = _mm_sad_epu8(acr, z); count += _mm_extract_epi64(sacr, 0); count += _mm_extract_epi64(sacr, 1); for (size_t i = nb * 16; i < n; i++) if (((const uint8_t *)s)[i] == c) count++; return count; }  Testing it we get substantially better times: $ gcc -march=native -std=c11 -DTEST_SSE -O3 memcount.c -o memcount && ./memcount
Generated random arrays.
Count of character 45 in 'a1' using 'memcount_sse' is: 408978
Elapsed time: 10908773 ns
Count of character 45 in 'a2' using 'memcount_sse' is: 409445
Elapsed time: 11195531 ns
Count of character 45 in 'a3' using 'memcount_sse' is: 410593
Elapsed time: 12147576 ns
Count of character 45 in 'a4' using 'memcount_sse' is: 104447433
Elapsed time: 10930570 ns
Count of character 45 in 'a5' using 'memcount_sse' is: 104857613
Elapsed time: 10856319 ns


The new speeds are between 8 and 9 GB/s, almost an order of magnitude faster than the “naive” code.

Explanation

As this is my first program using intrinsics, it’s probably far from optimal. But, as it can seem pretty impenetrable, let’s explain it step by step.

Initialization

size_t memcount_sse(const void *s, int c, size_t n)
{
size_t nb = n / 16;
size_t count = 0;
__m128i ct = _mm_set1_epi32(c * ((~0UL) / 0xff));
__m128i z = _mm_set1_epi32(0);
__m128i acr = _mm_set1_epi32(0);


This section just initializes some variables:

• nb: we will process the input in blocks of 16 bytes (128 bits), so we use this variable to track the number of blocks to be processed.
• count: the return value, to be filled with the byte count to be computed.
• ct: the comparison target, a 128 bit value filled with copies of c, the byte value we are counting.
• z: a zero 128 bit value.
• acr: the accumulated comparison result, a 128 bit value that accumulates the comparison results in a special format, to be described later.

Only one intrinsic is used here, _mm_set1_epi32(). It loads repeated copies of a 32 bit value in a 128 bit value.

Comparing and accumulating

    for (size_t i = 0; i < nb; i++)
{
__m128i b = _mm_lddqu_si128((const __m128i *)s + i);
__m128i cr = _mm_cmpeq_epi8 (ct, b);
acr = _mm_add_epi8(acr, cr);


The loop is self explanatory: we are just iterating over 128 bit blocks. The next intrinsic, _mm_lddqu_si128() just loads a potentially unaligned 128 bit value into the 128 bit variable b.

Now comes the key section: we compare b byte per byte against the comparison target ct, using the intrinsic _mm_cmpeq_epi8(), and we store 0xff when the bytes match and 0x00 when they don’t in the comparison result cr. That result is summed byte per byte in the accumulated counts variable acr, with the intrinsic _mm_add_epi8().

As 0xff is -1 modulo 256, the accumulation works normally, but we could have an overflow after 256 iterations. To avoid that, we need to move this counts to our count variable periodically.

Definitive accumulation

        if (i % 0xff == 0xfe)
{
acr = _mm_sub_epi8(z, acr);
__m128i sacr = _mm_sad_epu8(acr, z);
count += _mm_extract_epi64(sacr, 0);
count += _mm_extract_epi64(sacr, 1);
acr = _mm_set1_epi32(0);
}


Every 255 iterations (to avoid overflow), we do the following operations:

• Change the sign of acr by subtracting it from zero to get the positive values of the counts (using the _mm_sub_epi8() intrinsic).
• Get the absolute value of each byte and sum all these values to get the total count (using the _mm_sad_epu8() intrinsic, that leaves the sum split in the low and high parts of the 128 bit variable).
• Sum the low and high parts of the 128 bit variable to count, to get the total count, using the _mm_extract_epi64() intrinsic to get those parts.
• Set the accumulated comparison result acr variable to zero.

After loop accumulation

    acr = _mm_sub_epi8(z, acr);
__m128i sacr = _mm_sad_epu8(acr, z);
count += _mm_extract_epi64(sacr, 0);
count += _mm_extract_epi64(sacr, 1);


Accumulates the value that remains in acr after exiting the main loop.

Bytes after last block

    for (size_t i = nb * 16; i < n; i++)
if (((const uint8_t *)s)[i] == c)
count++;
return count;
}


Handles the bytes that come after the last 16 byte block and returns the total count.

Conclusion

We have seen that a substantial improvement of almost 10x can be obtained over the optimizing compiler results, even in cases where the code to be optimized is quite simple. It would be interesting to see how close to the 20 GB/s we can get with more carefully optimized code…

The code used for these tests is available in Github.

Escaping variables

(Just a short post to break the “dry spell”.)

One interesting consequence of C having non-pure functions is that every function you call can be accessing your local variables if pointers to them “escape”.

If we compile the following C code with -O3

extern void f(int *a);
extern void g(void);

int main(void)
{
int a = 0;
f(0);
a++;
g();
a++;
g();
a++;
return a;
}


we get machine code where the additions are collapsed and a is not even given a memory position:

main:
push    {r3, lr}
movs    r0, #0
bl      f(int*)
bl      g()
bl      g()
movs    r0, #3
pop     {r3, pc}


But if we pass a pointer to the local variable to a single function,

extern void f(int *a);
extern void g(void);

int main(void)
{
int a = 0;
f(&a);
a++;
g();
a++;
g();
a++;
return a;
}


we get a load/store after every function call:

main:
push    {lr}
sub     sp, sp, #12
add     r0, sp, #8
movs    r3, #0
str     r3, [r0, #-4]!
bl      f(int*)
ldr     r3, [sp, #4]
adds    r3, r3, #1
str     r3, [sp, #4]
bl      g()
ldr     r3, [sp, #4]
adds    r3, r3, #1
str     r3, [sp, #4]
bl      g()
ldr     r0, [sp, #4]
adds    r0, r0, #1
add     sp, sp, #12
ldr     pc, [sp], #4


Once a pointer to a variable has escaped, the optimizer has to assume that it’s part of the global state and every function call could be modifying it.

Abusing macros for typechecking

One commonly used macro in C programming is ASIZE(), generally defined as something like this

#define ASIZE(a) (sizeof(a)/sizeof(a[0]))


and used to calculate the number of elements in an array.

The main problem with this macro, as written, is that it doesn’t distinguish between arrays and pointers. If passed a pointer, it will silently produce wrong results:

Code

#include <stdio.h>

#define ASIZE(a) (sizeof (a) / sizeof((a)[0]))

int main(void)
{
short a[3];
short *b;
int c[2];
int *d;
long long e[5][4];
char *f[4];
char (*g)[4];
(void)a; (void)b; (void)c; (void)d; (void)e; (void)f; (void)g;
printf("ASIZE() accepts pointers, producing invalid results.\n");
printf("%zu\n", ASIZE( a ));
printf("%zu\n", ASIZE( b ));
printf("%zu\n", ASIZE( c ));
printf("%zu\n", ASIZE( d ));
printf("%zu\n", ASIZE( e ));
printf("%zu\n", ASIZE( f ));
printf("%zu\n", ASIZE( g ));
return 0;
}


Output

3
2
2
1
5
4
1


By adding a new macro checking if the parameter is an array, we can define a safer ASIZE():

#define CHECK_ARRAY(a) ((void)(0&&((int (*)(__typeof__(a[0])(*)[ASIZE(a)]))NULL)(&(a))))
#define ASIZE_SAFE(a) (CHECK_ARRAY(a), ASIZE(a))


Checking this new version, we see it gets the correct results when passed arrays, but now the compilation fails when applied to pointers:

Code

#include <stdio.h>

#define ASIZE(a) (sizeof (a) / sizeof((a)[0]))

#define CHECK_ARRAY(a) ((void)(0&&((int (*)(__typeof__(a[0])(*)[ASIZE(a)]))NULL)(&(a))))

#define ASIZE_SAFE(a) (CHECK_ARRAY(a), ASIZE(a))

int main(void)
{
short a[3];
short *b;
int c[2];
int *d;
long long e[5][4];
char *f[4];
char (*g)[4];
(void)a; (void)b; (void)c; (void)d; (void)e; (void)f; (void)g;
printf("ASIZE() accepts pointers, producing invalid results.\n");
printf("%zu\n", ASIZE( a ));
printf("%zu\n", ASIZE( b ));
printf("%zu\n", ASIZE( c ));
printf("%zu\n", ASIZE( d ));
printf("%zu\n", ASIZE( e ));
printf("%zu\n", ASIZE( f ));
printf("%zu\n", ASIZE( g ));
printf("ASIZE_SAFE() only accepts arrays (try uncommenting).\n");
printf("%zu\n", ASIZE_SAFE( a ));
//printf("%zu\n", ASIZE_SAFE( b ));
printf("%zu\n", ASIZE_SAFE( c ));
//printf("%zu\n", ASIZE_SAFE( d ));
printf("%zu\n", ASIZE_SAFE( e ));
//printf("%zu\n", ASIZE_SAFE( f ));
//printf("%zu\n", ASIZE_SAFE( g ));
return 0;
}


Output

ASIZE() accepts pointers, producing invalid results.
3
2
2
1
5
4
1
ASIZE_SAFE() only accepts arrays (try uncommenting).
3
2
5


It works in a relatively straightforward way, though I have put the details in a gist to avoid spoiling them.

Universality with only two NOT gates

In a previous post we have asked how many NOTs are needed to compute an arbitrary Boolean function. In this post we will see that only two NOT gates are enough.

Building 3 NOT gates starting from 2

If we call the inputs $X$, $Y$ and $Z$, we can make a function detecting when no more than one input is active using a single NOT gate:

$\displaystyle f(X, Y, Z) = \overline{XY + YZ + XZ}$.

Detects when no more than one input is active.

By selecting only the cases where at least one input is present, adding a term to detect when all the inputs are active and using an additional NOT gate, we can detect when exactly zero or two inputs are active:

$\displaystyle g(X, Y, Z) = \overline{f(X, Y, Z)(X + Y + Z) + XYZ}$

$\displaystyle = \overline{\overline{XY + YZ + XZ}(X + Y + Z) + XYZ}$.

Detects when zero or two inputs are active.

Now we know that if $X$ is not present, we either have:

• 0 inputs present: we can check that by simultaneously ensuring that we don’t have more than one input present and that we have either zero or two inputs present, $f(X, Y, Z)\cdot g(X, Y, Z)$.
• 1 input present: we should have no more than one input present and $Y$ or $Z$ should be present, $f(X, Y, Z)\cdot(Y + Z)$.
• 2 inputs present: we can check that by simultaneously ensuring that either zero or two inputs are present and that $Y$ and $Z$ are present, $g(X, Y, Z)\cdot YZ$.

Putting all together and adding the symmetrical cases:

$\displaystyle \overline{X} = f(X, Y, Z) \cdot (Y + Z) + (f(X, Y, Z) + YZ)\cdot g(X, Y, Z)$

$\displaystyle \overline{Y} = f(X, Y, Z) \cdot (X + Z) + (f(X, Y, Z) + XZ)\cdot g(X, Y, Z)$

$\displaystyle \overline{Z} = f(X, Y, Z) \cdot (X + Y) + (f(X, Y, Z) + XY)\cdot g(X, Y, Z)$.

Computing NOT X (the other cases are symmetrical).

Checking the solution

To get independent confirmation, let’s check the truth tables with a simple Python script:

from itertools import product

for x, y, z in product((False, True), repeat=3):
f_xyz = not (x and y or x and z or y and z)
g_xyz = not (f_xyz and (x or y or z) or x and y and z)
not_x = f_xyz and (y or z) or (f_xyz or y and z) and g_xyz
not_y = f_xyz and (x or z) or (f_xyz or x and z) and g_xyz
not_z = f_xyz and (x or y) or (f_xyz or x and y) and g_xyz
assert (not x == not_x) and (not y == not_y) and (not z == not_z)


Conclusion

As this technique allows us to expand two NOT gates to three and it can be applied repeatedly, we can compute an arbitrary Boolean function with a circuit containing only two NOT gates.

In a following post we will see how I arrived to the solution by using brute force.

How many NOTs are needed?

It’s relatively easy to see that we cannot compute an arbitrary Boolean function using only AND and OR gates. For example, even the NOT function cannot be computed using only those gates (why?).

Can we build a circuit to compute an arbitrary Boolean function using a constant number of NOT gates?

Solution to the “42 code golf” problem

This was my best result:

n=1e6,m,c,d;main(){while(c+=d==42,d=0,m=--n)while(d+=m%10,m/=10);printf("%d\n",c);}


It would have been nice to find a solution under 80 bytes but, after one hour of trying, that was the best I could manage…

Hilbert matrices are positive definite

(Argument stolen from Yahoo Answers)

It can be seen by using its integral form:

$\displaystyle H_{ij} = \int_0^1 dx\,x^{i+j-2}$

Then we can express the positive definite condition as

$\displaystyle \sum_{i,j} x_i H_{ij} x_j > 0$

for $\mathbf{x} \ne 0$.

Replacing and operating:

$\displaystyle \sum_{i,j} x_i \int_0^1 du\,u^{i+j-2} x_j = \int_0^1 du \sum_{i,j} x_i u^{i+j-2} x_j$

$\displaystyle = \int_0^1 du \sum_{i,j} x_i u^{i-1} x_j u^{j-1}$

$\displaystyle = \int_0^1 du \sum_i x_i u^{i-1} \sum_j x_j u^{j-1}$

$\displaystyle = \int_0^1 du \left( \sum_i x_i u^{i-1} \right) \left( \sum_j x_j u^{j-1} \right)$

$\displaystyle = \int_0^1 du \left( \sum_i x_i u^{i-1} \right)^2$

The term inside the parentheses is just a (polynomial) function of $u$:

$\displaystyle = \int_0^1 du\,p(u)^2$

The only way the integral can be zero is if $p(u)$ is identically zero:

$\displaystyle \sum_i x_i u^{i-1} = 0 \implies \forall i: x_i = 0$.

QED

Solving the viral Singaporean math problem

The following problem has become very popular in social media:

In this blog we have solved similar problems before, but this is one that can be easily solved by hand. We only need to be careful not to confuse our knowledge state with the one of Albert & Bernard.

Albert and Bernard just became friends with Cheryl, and they want to know when her birthday is. Cheryl gives them a list of 10 possible dates.

May 15 May 16 May 19

June 17 June 18

July 14 July 16

August 14 August 15 August 17

Cheryl then tells Albert and Bernard separately the month and the day of her birthday, respectively.

We will describe a list of the possible knowledge states of Albert and Bernard after being given that information:

Albert
1. May 15 or May or May 19
2. June 17 or June 18
3. July 14 or July 16
4. August 14 or August 15 or August 17
Bernard
1. July 14 or August 14
2. May 15 or August 15
3. May 16 or July 16
4. June 17 or August 17
5. June 18
6. May 19

Albert: I don’t know when Cheryl’s birthday is, but I know that Bernard does not know, too.

We already knew that Albert wouldn’t know the day based just on being given the month, but he is giving us additional information by telling us that he knows that Bernard doesn’t know. Bernard would know the date if the day were 18 or 19, so Albert knows that those days could not be the right ones. That excludes options 1 and 2 from our knowledge of Albert knowledge:

Albert
1. July 14 or July 16
2. August 14 or August 15 or August 17

Bernard can do the same deductions we have and eliminate the options that are not possible from his state of knowledge (all the options with months different from July and August).

Bernard
1. July 14 or August 14
2. August 15
3. July 16
4. August 17

Bernard: At first, I didn’t know when Cheryl’s birthday is, but I know now.

Updating our knowledge of Bernard knowledge:

Bernard
1. August 15
2. July 16
3. August 17

As Albert also knows what we know about Bernard knowledge…

Albert
1. July 16
2. August 15 or August 17

Albert: Then I also know when Cheryl’s birthday is.

Now we know the right date:

1. July 16
1. July 16