## The percolation algorithm

This page contains supplementary material that goes with the paper A fast Monte Carlo algorithm for site or bond percolation, M. E. J. Newman and R. M. Ziff, Phys. Rev. E 64, 016706 (2001).

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

For convenience, the code for the algorithm is also given below, along with a blow-by-blow description, reproduced from our paper.

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];        /* 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 231, 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] = (i+1)%N;
nn[i] = (i+N-1)%N;
nn[i] = (i+L)%N;
nn[i] = (i+N-L)%N;
if (i%L==0) nn[i] = i+L-1;
if ((i+1)%L==0) nn[i] = 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();
}
```