Implementing cellular automata

It is convenient to represent the state of a cellular automaton at each step by a list such as {0, 0, 1, 0, 0}, where 0 corresponds to a white cell and 1 to a black cell. An initial condition consisting of n white cells with one black cell in the middle can then be obtained with the function (see below for comments on this and other Mathematica functions)

CenterList[n_Integer] := ReplacePart[Table[0, {n}], 1, Ceiling[n/2]]

For cellular automata of the kind discussed in this chapter, the rule can also be represented by a list. Thus, for example, rule 30 on page 27 corresponds to the list {0, 0, 0, 1, 1, 1, 1, 0}. (The numbering of rules is discussed on page 53.) In general, the list for a particular rule can be obtained with the function

ElementaryRule[num_Integer] := IntegerDigits[num, 2, 8]

Given a rule together with a list representing the state a of a cellular automaton at a particular step, the following simple function gives the state at the next step:

CAStep[rule_List, a_List] := rule〚8 - (RotateLeft[a] + 2 (a + 2 RotateRight[a]))〛

A list of states corresponding to evolution for t steps can then be obtained with

CAEvolveList[rule_, init_List, t_Integer] :=NestList[CAStep[rule, #]&, init, t]

Graphics of this evolution can be generated using

CAGraphics[history_List] := Graphics[Raster[1 - Reverse[history]], AspectRatio  Automatic]

And having set up the definitions above, the Mathematica input

Show[CAGraphics[CAEvolveList[ElementaryRule[30], CenterList[103], 50]]]

will generate the image:

The description just given should be adequate for most cellular automaton simulations. In some earlier versions of Mathematica a considerably faster version of the program can be created by using the definition

CAStep = Compile[{{rule, _Integer, 1}, {a, _Integer, 1}}, rule〚8 - (RotateLeft[a] + 2 (a + 2 RotateRight[a]))〛]

In addition, in Mathematica 4 and above, one can use

CAStep[rule_, a_]:=rule〚8 - ListConvolve[{1, 2, 4}, a, 2]]]〛

or directly in terms of the rule number num

Sign[BitAnd[2^ListConvolve[{1, 2, 4}, a, 2], num]]

(In versions of Mathematica subsequent to the release of this book the built-in CellularAutomaton function can be used, as discussed on page 867.) It is also possible to have CAStep call the following external C language program via MathLink—though typically with successive versions of Mathematica the speed advantage obtained will be progressively less significant:

#include "mathlink.h"

main(argc, argv)
int argc; char *argv[];
{
MLMain(argc, argv);
}

void casteps(revrule, rlen, a, n, steps)
int *revrule, rlen, *a, n, steps;
{
int i, *ap, t, tp;

for (i = 0; i < steps; i++) { a[0] = a[n-2]; /* right boundary */ a[n-1] = a[1]; /* left boundary */
t = a[0]; for (ap = a+1; ap <= a+n-2; ap++) { tp = ap[0]; ap[0] = revrule[ap[1]+2*(tp + 2*t)]; t = tp; } }

MLPutIntegerList(stdlink, a, n);
}

The linkage of this external program to the Mathematica function CAStep is achieved with the following MathLink template (note the optional third argument which allows CAStep to perform several steps of cellular automaton evolution at a time):

:Begin:
:Function: casteps
:Pattern: CAStep[rule_List, a_List, steps_Integer:1]
:Arguments: {Reverse[rule], a, steps}
:ArgumentTypes: {IntegerList, IntegerList, Integer}
:ReturnType: Manual
:End:

There are a couple of tricky issues in the C program above. First, cellular automaton rules are always defined to use the old values of neighbors in determining the new value of any particular cell. But since the C program explicitly updates values sequentially from left to right, the left-hand neighbor of a particular cell will already have been given its new value when one tries to updates the cell itself. As a result, it is necessary to store the old value of the left-hand neighbor in a temporary variable in order to make it available for updating the cell itself. (Another approach to this problem is to maintain two copies of the array of cells, and to interchange pointers to them after every step in the cellular automaton evolution.)

Another tricky point in cellular automaton programs concerns boundary conditions. Since in a practical computer one can use only a finite array of cells, one must decide how the cellular automaton rule is to be applied to the cells at each end of the array. In both the Mathematica and the C programs above, we effectively use a cyclic array, in which the left neighbor of the leftmost cell is taken to be rightmost cell, and vice versa. In the C program, this is implemented by explicitly copying the value of the leftmost cell to the rightmost position in the array, and vice versa, before updating the values in the array. (In a sense there is a bug in the program in that the update only puts new values into n-2 of the n array elements.)