Here is the code for the algorithm as described in the paper. It is in ANSI standard C. Available files:

`perc.c` :
| complete source code for the algorithm |

`rnglong.c` ,
`rnglong.h` :
| source code for a suitable random number generator to go with it |

`Makefile` :
| a Makefile suitable for compiling the program on Unix-type systems (including Linux). |

First we set up some constants and global variables:

#include ‹stdlib.h› #include ‹stdio.h› #define L 128 /* Linear dimension */ #define N (L*L) #define EMPTY (-N-1) int ptr[N]; /* Array of pointers */ int nn[N][4]; /* Nearest neighbors */ int order[N]; /* Occupation order */

Sites are indexed with a single signed integer label for speed, taking
values from 0 to *N*-1. Note that on computers which represent
integers in 32 bits, this program can, for this reason, only be used for
lattices of up to 2^{31}, or about 2 billion sites. While this is
adequate for most purposes, longer labels will be needed if you wish to
study larger lattices.

The array `ptr[]`

serves triple duty: for non-root occupied
sites it contains the label for the site's parent in the tree (the
"pointer"); root sites are recognized by a negative value of
`ptr[]`

, and that value is equal to minus the size of the
cluster; for unoccupied sites `ptr[]`

takes the value
`EMPTY`

.

Next we set up the array `nn[][]`

which contains a list of the
nearest neighbors of each site. Only this array need be changed in order
for the program to work with a lattice of different topology.

void boundaries() { int i; for (i=0; i‹N; i++) { nn[i][0] = (i+1)%N; nn[i][1] = (i+N-1)%N; nn[i][2] = (i+L)%N; nn[i][3] = (i+N-L)%N; if (i%L==0) nn[i][1] = i+L-1; if ((i+1)%L==0) nn[i][0] = i-L+1; } }

Now we generate the random order in which the sites will be occupied, by
randomly permuting the integers from 0 to *N*-1:

void permutation() { int i,j; int temp; for (i=0; i‹N; i++) order[i] = i; for (i=0; i‹N; i++) { j = i + (N-i)*drand(); temp = order[i]; order[i] = order[j]; order[j] = temp; } }

Here the function `drand()`

generates a random double precision
floating point number between 0 and 1. Many people will have such a
function already to hand. For those who don't, a suitable one is supplied
here and here.

We also define a function which performs the "find" operation, returning the label of the root site of a cluster, as well as path compression. The version we use is recursive, as described in the paper:

int findroot(int i) { if (ptr[i]‹0) return i; return ptr[i] = findroot(ptr[i]); }

The code to perform the actual algorithm is quite brief. It works exactly
as described in the paper. Sites are occupied in the order specified by
the array `order[]`

. The function `findroot()`

is
called to find the roots of each of the adjacent sites. If amalgamation is
needed, it is performed in a weighted fashion, smaller clusters being added
to larger (bearing in mind that the value of `ptr[]`

for the
root nodes is *minus* the size of the corresponding cluster).

void percolate() { int i,j; int s1,s2; int r1,r2; int big=0; for (i=0; i‹N; i++) ptr[i] = EMPTY; for (i=0; i‹N; i++) { r1 = s1 = order[i]; ptr[s1] = -1; for (j=0; j‹4 j++) { s2 = nn[s1][j]; if (ptr[s2]!=EMPTY) { r2 = findroot(s2); if (r2!=r1) { if (ptr[r1]>ptr[r2]) { ptr[r2] += ptr[r1]; ptr[r1] = r2; r1 = r2; } else { ptr[r1] += ptr[r2]; ptr[r2] = r1; } if (-ptr[r1]>big) big = -ptr[r1]; } } } printf("%i %i\n",i+1,big); } }

The main program is now simple:

main() { boundaries(); permutation(); percolate(); }

Last modified: January 19, 2001

Mark Newman