At The Bakery

After the topical excitement of the last couple of posts, let’s look at an all-time great – Leslie Lamport’s Bakery Algorithm (and of course this is still topical; Lamport is the most recent winner of the Turing Award).

The problem is mutual exclusion without mutual exclusion primitives. Usually, it’s described in the context of a shared memory system (and that is what we will implement here), but will work equally well in a message-passing system with only local state (each thread or process only needs to write to its own part of the store).

For further details, and Lamport’s later thoughts see http://research.microsoft.com/en-us/um/people/lamport/pubs/pubs.html#bakery: “For a couple of years after my discovery of the bakery algorithm, everything I learned about concurrency came from studying it.” – and since Lamport understands more about concurrency than just about anyone on the planet, it’s maybe worth spending some time looking at it ourselves.

I’m not going to attempt to prove the algorithm correct, I’ll leave that to Lamport, but the crucial idea seems to me to be that a thread reading a particular value from another thread is a synchronization signal from that thread – here, reading a false value for the entering variable is a signal that the other thread isn’t in the process of deciding on it’s own number, therefore it is safe for the reading process to proceed.

Implementing on a real multiprocessor system, we find that use of memory barriers or synchronization primitives is essential – the algorithm requires that reads and writes are serialized in the sense that once a value is written, other processes won’t see an earlier value (or earlier values of other variables). This doesn’t conflict with what Lamport says about not requiring low-level atomicity – we can allow reads and writes to happen simultaneously, with the possibility of a read returning a bogus value – and in fact we can simulate this in the program by writing a random value just before a process selects its real ticket number, but once a write has completed, all processes should see the new value.

Another essential feature is the volatile flag – as many have pointed out, this isn’t enough by itself for correct thread synchronization, but for shared memory systems, prevents the compiler from making invalid assumptions about consistency of reads from shared variables.

A final point – correctness requires that ticket numbers can increase without bound, this is hard to arrange in practice, so we just assert if they grow too large (this rarely happens in reality, unless we get carried away with our randomization).

#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <assert.h>

// Compile with: g++ -Wall -O3 bakery.cpp -pthread -o bakery

static const int NTHREADS = 4;

// Some features to play with
//#define NCHECK     // Disable crucial check
//#define NSYNC      // Disable memory barrier
//#define NVOLATILE  // Disable volatile

#if defined NVOLATILE
#define VOLATILE
#else
#define VOLATILE volatile
#endif

VOLATILE bool entering[NTHREADS];
VOLATILE unsigned number[NTHREADS];
VOLATILE int count = 0;
VOLATILE int total = 0;

unsigned getmax(int n)
{
  unsigned max = 0;
  for (int i = 0; i < n; i++) {
    if (number[i] > max) max = number[i];
  }
  return max;
}

bool check(int i, int j)
{
  return number[j] < number[i] || 
         (number[j] == number[i] && j < i);
}

inline void synchronize()
{
#if !defined NSYNC
    // gcc builtin full memory barrier
  __sync_synchronize();
#endif
}

void lock(int i) {
  entering[i] = true;
  synchronize();
  // Simulate non-atomic write
  number[i] = rand();
  synchronize();
  number[i] = 1 + getmax(NTHREADS);
  assert(number[i] > 0);
  entering[i] = false;
  synchronize();
  for (int j = 0; j < NTHREADS; j++) {
    // Wait until thread j receives its number:
#if !defined NCHECK
    while (entering[j]) { /* nothing */ }
#endif
    // At this point, we have read a false value for 
    // "entering[j]", therefore any number picked by j
    // later will takes our choice into account, any value
    // chosen earlier (and so might be less than ours) 
    // will be visible to us in the following test.

    // Wait until all threads with smaller numbers or with 
    // the same number, but with higher priority, finish 
    // their work:
    while ((number[j] != 0) && check(i,j)) { /* nothing */ }
  }
}

void unlock(int i)
{
  number[i] = 0;
}

void *threadfun(void *arg) 
{
  int i = *(int*)(arg);
  while (true) {
    lock(i);
    total++;
    if (total % 1000000 == 0) fprintf(stderr,"%c", 'a'+i);
    assert(count==0); // Check we have exclusive access
    count++;
    // It's not clear that these synchs are unnecessary,
    // but nothing seems to break if I remove them.
    //synchronize();
    count--;
    //synchronize();
    unlock(i);
    // non-critical section...
  }
  return NULL;
}

int main()
{
  pthread_t t[NTHREADS];
  int n[NTHREADS];
  for (int i = 0; i < NTHREADS; i++) {
    n[i] = i;
    pthread_create(&t[i], NULL, threadfun, (void*)&n[i]);
  }
  for (int i = 0; i < NTHREADS; i++) {
    pthread_join(t[i], NULL);
  }
}

Loopless Ruler

This function (based on Algorithm L for loopless Gray code generation in TAOCP 7.2.1.1):

void ruler(int n)
{
  int f[n+1];
  for (int j = 0; j <= n; j++) {
    f[j] = j;
  }
  while(true) {
    visit(f,n+1);
    int j = f[0];
    f[0] = 0;
    if (j == n) break;
    f[j] = f[j+1];
    f[j+1] = j+1;
  }
}

enables us to compute the “ruler function”, ρ for 1, 2, 3…:

0 1 0 2 0 1 0 3 0 1 0 2 0 1 0...

It’s a “loopless” algorithm and uses f, an array of “focus pointers” to keep track of what’s going on.

The ruler function ρ(n) is essentially the number of trailing 0’s in the binary representation of n.

This is what’s printed by ruler(4):

0 1 2 3 4 
1 1 2 3 4 
0 2 2 3 4 
2 1 2 3 4 
0 1 3 3 4 
1 1 3 3 4 
0 3 2 3 4 
3 1 2 3 4 
0 1 2 4 4 
1 1 2 4 4 
0 2 2 4 4 
2 1 2 4 4 
0 1 4 3 4 
1 1 4 3 4 
0 4 2 3 4 
4 1 2 3 4 

The left hand side is the ruler function of course, what are the other values?

It’s clearer if we print f in reverse and match up with binary representations:

00000
43210
00001
43211
00010
43220
00011
43212
00100
43310
00101
43311
00110
43230
00111
43213
01000
44210
01001
44211
01010
44220
01011
44212
01100
43410
01101
43411
01110
43240
01111
43214

If n[j] is the jth binary digit of n, then f[j] is just j, unless n[j] is the first of a block of 1’s, in which case f[j] is the index of the first 0 after that block.

For example, 01101, has block starting at 0 and ending at 1, so f[0] = 1, and another block starting at 2, ending at 4, so f[2] = 4; for other j, f[j] = j, so the focus pointers are 43411.

If n is of the form:

j   k
01110111

with f[0] = k, f[k] = k, f[k+1] = j

0n incrementing we get:

j   k
01111000

Now f[0] = 0, f[k] = j and f[k+1] = k+1, other values of f are the same

If k is 0, then we go from:

j   k
01110

with f[0] = 0, f[1] = j

to:

j   k
01111

with f[0] = 0, f[1] = 1

Either way, we can see that the inner loop of ruler is correct.

Ringing the Changes

Apart from a brief visit many years ago to Macclesfield church with my friend Jo, I have never taken part in the ancient English art of bell-ringing, but permutation generation has always seemed a fascinating topic.

In TAOCP 7.2.1.2, Knuth has Algorithm P (“plain changes”) for generating all permutations using adjacent interchanges:

void taocp1(int *a, int n)
{
  int c[n+1];
  int o[n+1];
  for (int j = 1; j <= n; j++) {
    c[j] = 0; o[j] = 1;
  }
 P2:
  visit(a,n);
 P3:
  int j = n; int s = 0;
 P4:
  int q = c[j] + o[j];
  if (q < 0) goto P7;
  if (q == j) goto P6;
 P5:
  swap(a[j-c[j]+s],a[j-q+s]);
  c[j] = q;
  goto P2;
 P6:
  if (j == 1) return;
  s = s+1;
 P7:
  o[j] = -o[j];
  j = j-1;
  goto P4;
}

It works just fine, but I find it a little hard to follow and let’s face it, we just don’t write code like that anymore, with all those gotos, 1-based arrays, etc.

0-basing the arrays isn’t too much of a problem, for eliminating the gotos we could take the functional approach:

void taocp2(int *a, int n)
{
  int c[n];
  int o[n];
  for (int j = 0; j < n; j++) {
    c[j] = 0; o[j] = 1;
  }
  P2(a,c,o,n);
}

void P2(int *a, int *c, int *o, int n) {
  visit(a,n);
  P3(a,c,o,n);
}

void P3(int *a, int *c, int *o, int n) {
  P4(a,c,o,n,n-1,0);
}

void P4(int *a, int *c, int *o, int n, int j, int s) {
  int q = c[j] + o[j];
  if (q < 0) P7(a,c,o,n,j,s);
  else if (q == j+1) P6(a,c,o,n,j,s);
  else P5(a,c,o,n,j,s,q);
}

void P5(int *a, int *c, int *o, int n, int j, int s, int q) {
  swap(a[j-c[j]+s],a[j-q+s]);
  c[j] = q;
  P2(a,c,o,n);
}

void P6(int *a, int *c, int *o, int n, int j, int s) {
  if (j != 1) {
    P7(a,c,o,n,j,s+1);
  }
}

void P7(int *a, int *c, int *o, int n, int j, int s) {
  o[j] = -o[j];
  P4(a,c,o,n,j-1,s);
}

but somehow that seems missing the point (I don’t think that the worst thing about that code is that it still uses modifiable state).

There are actually two component parts to the algorithm: generating a “mixed radix reflected Gray code” which defines the progressive inversions of the array elements, and going from the inversion changes to the actual indexes of the swapped objects.

For generating the Gray codes, with a bit of rewriting, 0-basing our arrays etc. we get:

void gray(int n)
{
  int j, c[n], o[n];
  for (int j = 0; j < n; j++) {
    c[j] = 0; o[j] = 1;
  }
  do {
    visit(c,n);
    for (j = n-1; j > 0; j--) {
      int q = c[j] + o[j];
      if (q < 0) {
        o[j] = -o[j];
      } else if (q == j+1) {
        o[j] = -o[j];
      } else {
        c[j] = q;
        break;
      }
    } 
  } while (j != 0);
}

Each time around, find the highest element whose inversion count can be changed in the desired direction. For elements whose inversion count can’t be changed, change direction. If no element can be changed we are done.

The next step is to calculate the position of element j – but this is just j less the number of elements less than j that appear after j (ie. the number of inversions, c[j]), plus the number of elements greater than j that appear before j – but this is just the number of elements we have been passed over in the “if (q == j+1)” step above, so we can now add in the rest of algorithm P:

void taocp3(int *a, int n)
{
  int c[n];
  int o[n];
  for (int j = 0; j < n; j++) {
    c[j] = 0; o[j] = 1;
  }
  int j;
  do {
    visit(a,n);
    int s = 0;
    for (j = n-1; j > 0; j--) {
      int q = c[j] + o[j];
      if (q < 0) {
        o[j] = -o[j];
      } else if (q == j+1) {
        o[j] = -o[j];
        s++;          // This element will be before element j
      } else {
        swap(a[j-c[j]+s],a[j-q+s]);
        c[j] = q;
        break;
      }
    } 
  } while (j != 0);
}

Knuth’s algorithm essentially Dijkstra’s:

Click to access EWD553.PDF

which explains it very lucidly: there is a 1-1 correspondence between inversion counts & permutations – and a Gray enumeration of inversion gives us a sequence of permutations where only 1 element at a time changes its inversion count, and only by 1 or -1, which can only be if we exchange adjacent elements.

Knuth also gives a “loopless” algorithm for generating reflected Gray sequences, and we could use a table of element positions to construct the permutations from this:

void gray2(int *m, int n)
{
  int a[n],f[n+1],o[n];
  for (int j = 0; j < n; j++) {
    a[j] = 0; f[j] = j; o[j] = 1;
  }
  f[n] = n;
  while (true) {
    visit(a,n);
    int j = f[0];
    f[0] = 0;
    if (j == n) break;
    a[j] += o[j];
    if (a[j] == 0 || a[j] == m[j]-1) {
      o[j] = -o[j];
      f[j] = f[j+1];
      f[j+1] = j+1;
    }
  }
}

The “classic” Johnson-Trotter algorithm for plain changes involves consideration of “mobile” elements: each element has a current direction and it is “mobile” if if it greater than the next adjacent element (if there is one) in that direction. The algorithm proceeds by finding the highest mobile element and moving it accordingly:

void jt(int *a, int n)
{
  int o[n]; // Direction of element i, 1 = left, 0 = right
  int c[n]; // Position of element i
  for (int j = 0; j < n; j++) {
    o[j] = 1; c[j] = j;
  }
  int j; // Will be set to the highest mobile element
  do {
    visit(a,n);
    // Search high to low. 0 is never mobile
    for (j = n-1; j > 0; j--) {
      int p = c[j]; // Position of element j
      if (o[j]) { // Going left
        int k = a[p-1];
        if (p > 0 && k < j) {
          // Swap and adjust positions
          a[p] = k; c[k] = p;
          a[p-1] = j; c[j] = p-1; 
          break;
        }
      } else { // Going right
        int k = a[p+1];
        if (p < n-1 && k < j) {
          a[p] = k; c[k] = p;
          a[p+1] = j; c[j] = p+1; 
          break;
        }
      }
      o[j] = !o[j]; // Change direction & continue
    }
  } while (j > 0);
}

Note that we can check elements from high to low directly & stop as soon as a mobile element is found – there is no need to check every element.

Johnson-Trotter and Algorithm P are essentially doing the same thing: when moving the mobile element, we are changing its inversion count in the same way as we do directly in P.

Performance-wise, the two algorithms are very similar though (somewhat surprisingly to me) the Johnson-Trotter version seems a little faster.

Either way, my laptop generates the 3628800 permutations of 10 elements in 0.03 seconds, outputting them (even to /dev/null) takes 4.5 seconds. It takes about 3 seconds to run through the 479001600 permutations of 12 elements (I haven’t tried outputting them). We can streamline the computation further by taking advantage of the fact that most of the time we are just moving the highest element in the same direction, but that can wait for another day.

Finally, here is a function that gives the next permutation in plain changes without maintaining any state:

bool nextperm(int *a, int n)
{
  int c[n], o[n];
  for (int i = 0; i < n; i++) {
    c[i] = o[i] = 0;
  }
  // Count inversions
  for (int i = 0; i < n; i++) {
    for (int j = i+1; j < n; j++) {
      if (a[j] < a[i]) c[a[i]]++;
    }
  }
  // Find directions
  for (int i = 1, k = 0; i < n-1; i++) {
    k += c[i];
    o[i+1] = k%2;
  }
  // Find mobile element and move it.
  for (int j = n-1, s = 0; j > 0; j--) {
    if (!o[j]) {
      if (c[j] < j) {
	swap(a[j-c[j]+s],a[j-c[j]+s-1]);
	return true;
      }
      s++;
    } else {
      if (c[j] > 0) {
	swap(a[j-c[j]+s],a[j-c[j]+s+1]);
	return true;
      }
    }
  }
  return false;
}

void stateless(int *a, int n)
{
  do {
    visit(a,n);
  } while (nextperm(a,n));
}

We use the Knuth/Dijkstra trick to avoid building a position table. I suspect that n-squared inversion counting can be improved, but we still come in at under second for all permutations of 10, though that’s about 30 times slower than the optimized version.

Krazy Kompilation

Let’s take a simple factorial program:

int fact(int n)
{
  if (n == 0) return 1;
  else return n * fact(n-1);
}

and look at what gcc compiles it to (using Matthew Godbolt’s excellent GCC Explorer):

With -O1, a nice clean version of the recursive function:

fact(int):
	pushq	%rbx
	movl	%edi, %ebx
	movl	$1, %eax
	testl	%edi, %edi
	je	.L2
	leal	-1(%rdi), %edi
	call	fact(int)
	imull	%ebx, %eax
.L2:
	popq	%rbx
	ret

With -O2, that’s impressive, it’s done some program transformation and now we have a loop:

fact(int):
	testl	%edi, %edi
	movl	$1, %eax
	je	.L4
.L3:
	imull	%edi, %eax
	subl	$1, %edi
	jne	.L3
	rep
	ret
.L4:
	rep
	ret

With -O3, whaaah!:

fact(int):
	testl	%edi, %edi
	je	.L7
	movl	%edi, %edx
	movl	%edi, %esi
	shrl	$2, %edx
	leal	0(,%rdx,4), %ecx
	testl	%ecx, %ecx
	je	.L8
	cmpl	$6, %edi
	jbe	.L8
	leal	-1(%rdi), %eax
	movl	%edi, -24(%rsp)
	movdqa	.LC1(%rip), %xmm4
	movl	%eax, -20(%rsp)
	leal	-2(%rdi), %eax
	movd	-20(%rsp), %xmm2
	movl	%eax, -16(%rsp)
	leal	-3(%rdi), %eax
	movd	-16(%rsp), %xmm1
	movl	%eax, -12(%rsp)
	xorl	%eax, %eax
	movd	-12(%rsp), %xmm0
	punpckldq	%xmm0, %xmm1
	movd	-24(%rsp), %xmm0
	punpckldq	%xmm2, %xmm0
	punpcklqdq	%xmm1, %xmm0
	movdqa	.LC0(%rip), %xmm1
.L4:
	movdqa	%xmm1, %xmm3
	psrldq	$4, %xmm1
	addl	$1, %eax
	movdqa	%xmm0, %xmm2
	cmpl	%eax, %edx
	pmuludq	%xmm0, %xmm3
	paddd	%xmm4, %xmm0
	psrldq	$4, %xmm2
	pmuludq	%xmm1, %xmm2
	pshufd	$8, %xmm3, %xmm1
	pshufd	$8, %xmm2, %xmm2
	punpckldq	%xmm2, %xmm1
	ja	.L4
	movdqa	%xmm1, %xmm2
	subl	%ecx, %edi
	cmpl	%ecx, %esi
	psrldq	$8, %xmm2
	movdqa	%xmm2, %xmm0
	psrldq	$4, %xmm2
	pmuludq	%xmm1, %xmm0
	psrldq	$4, %xmm1
	pshufd	$8, %xmm0, %xmm0
	pmuludq	%xmm2, %xmm1
	pshufd	$8, %xmm1, %xmm1
	punpckldq	%xmm1, %xmm0
	movdqa	%xmm0, %xmm1
	psrldq	$4, %xmm1
	movdqa	%xmm1, %xmm2
	psrldq	$4, %xmm1
	pmuludq	%xmm0, %xmm2
	psrldq	$4, %xmm0
	pmuludq	%xmm0, %xmm1
	pshufd	$8, %xmm2, %xmm0
	pshufd	$8, %xmm1, %xmm1
	punpckldq	%xmm1, %xmm0
	movd	%xmm0, -24(%rsp)
	movl	-24(%rsp), %eax
	je	.L13
.L9:
	imull	%edi, %eax
	subl	$1, %edi
	jne	.L9
	rep
	ret
.L7:
	movl	$1, %eax
	ret
.L8:
	movl	$1, %eax
	jmp	.L9
.L13:
	ret

Well, OK, it’s unwound the loop and vectorized with SSE instructions though a more precise analysis will have to wait for a future post. Impressive optimization from GCC though and pretty groovy considering the biggest factorial we can fit in a 32-bit int is fact(15)!