The percolation algorithm

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).

Description of the code

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³¹, 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