Recursive Worlds
Clifford A. Pickover
Cliff is the author of numerous books including Computers, Pattern, Chaos, and Beauty (1990), Computers and the Imagination (1991), and Mazes for the Mind (1992), all published by St. Martin's Press. He is also a researcher at IBM's Thomas J. Watson Research Center. Cliff can be contacted at cliff@watson.ibm.com.
He watched her for a long time and she knew that he was watching her and he knew that she knew he was watching her, and he knew that she knew that he knew; in a kind of regression of images that you get when two mirrors face each other and the images go on and on and on in some kind of infinity.
--Robert Pirsig, Lila
Whatever can be done once can always be repeated," begins Louise B. Young in The Mystery of Matter when describing the shapes and structures of nature. From the branching of rivers and blood vessels, to the highly convoluted surface of brains and bark, the physical world contains intricate patterns formed from simple shapes through the recursive application of dynamic procedures. Questions about the fundamental rules underlying the variety of nature have led to the search to identify, measure, and define these patterns in precise scientific terms.
Recursion is a fundamental concept in computer science, mathematics, biology, art, and even linguistics. Imagine, for instance, you're paging through a dictionary to find the definition of "recursion," finding it just below "recuperate:"
recumbent - lying down
recuperate - to recover health
recursion - look up the definition of recursion
Here, "recursion" is defined in terms of itself--a recursive definition, in other words. ("Recursion" is distinct from "iteration," which is exemplified by turning the pages of the dictionary to find a particular word, expressed in a program as: NextPage=CurrentPage+1.)
Other linguistic recursive definitions include "a wolf pack is two wolves or a wolf pack together with a wolf," or more simple constructs such as "art is art," "a dog is a dog," or even "a dog is not a dog." Sometimes the term "self-referential" is used when referring to these constructs.
The related concept of "quining," discussed in Douglas Hofstadter's Godel, Escher, Bach (Vintage, 1980), is a process of taking a group of words, and forming a self-referential sentence by preceding the original group with the same group enclosed in quotes. For example: "is a sentence fragment of seven words" is a sentence fragment of seven words.
Perhaps the most striking application of recursion occurs in the biological world where growth starts with a bud, which grows into a pipe, which then branches into two buds, each of these two buds branching in a recursive growth process. Iteration, or repeated application of these simple rules, results in a self-similar oak tree, arterial blood system, the bronchial system of lungs, or the like. The branching patterns are thought to be the result of the simplest of growth algorithms: The steps repeat the previous ones on smaller and smaller scales. Example 2
Figure 1 shows a coral-like form I computed using simple recursive branching rules at smaller and smaller size scales. In geometry, one of the most interesting consequences of recursive processes is "self-similarity." A self-similar object appears roughly the same after increasing or shrinking in size. Like a nested collection of Russian dolls within dolls, self-similar objects contain within themselves miniature copies of themselves. Look inside a turbulent stream of water: The largest eddies contain smaller ones, and these contain smaller ones still. The beautiful consequences of self-similarity are intricate, fine-grained patterns, now generally called "fractals." The term fractals was coined in The Fractal Geometry of Nature by Benoit Mandelbrot (Freeman, 1982) to encompass many of the detailed and convoluted shapes found in nature and produced by recursion in both the mathematical and natural worlds.
One of my earliest (and still most interesting) introductions to recursion and self-similarity was a Don Martin cartoon in Mad magazine. In the first frame, a man lies anesthetized on a hospital operating table. In the second, a surgeon uses a small circular saw to cut along the circumference of the man's head. Then, he removes the patient's skull cap.
Inside the head, the surgeon doesn't find a brain, but rather, another fully formed head, identical to the original, but slightly smaller. The surgeon removes this smaller head and opens its skull cap, finding yet another head. He continues to open smaller and smaller heads until the last head sits in the palm of his hand. He opens this one and finds a slip of paper that reads, "Inspected by number 47."
In computer programming, recursion often refers to the structure and functioning of programs. The common definition of a recursive program is one that calls itself, and a recursive function is one that is defined in terms of itself. (Interestingly, recursion can be removed from any recursive program using iteration.)
Perhaps the most common example of recursion in programming and in mathematics (where recursion is called a "recurrence relation") is the factorial function: X! =Xx(X-1)! for X
1, 0! =1 where X is an integer. This can be accomplished with a simple recursive program; in Example 1(a), the program calls itself in line 3. Another common example of a recurrence relation is one that defines the Fibonacci numbers. This sequence of numbers, called the Fibonacci sequence, plays important roles in mathematics and nature. These numbers are such that, after the first two, every number in the sequence equals the sum of the two previous numbers: FN= FN-1+FN-2 for N2,F0=F1=1. This defines the sequence: 1,1,2,3,5,8,13,21,34,55_.
Like the factorial function, a simple recursive program can be written to generate the Fibonacci sequence such as in Example 1(b).
In actuality, as with many recurrence relations, it's easy to compute FN using arrays in a non-recursive program like Example 1(c). This program computes the first 30 Fibonacci numbers using an array size of 30. This method of using arrays to store previous results is usually the preferred method for evaluating recurrence relations, because it allows even complex expressions to be processed in a uniform and efficient manner. Of course, in this example, you can avoid the array by retaining the last two values in two variables.
I'll now turn to a particular class of self-similar objects I call recursive lattices because they can easily be constructed using checkerboards of different sizes. The concept of repeatedly replacing copies of a pattern at different-size scales to produce interesting patterns dates back many decades, including the work of mathematicians Koch, Hilbert, and Peano. More recent work has been done by Mandelbrot and Lindenmeyer. Artists such as Escher, Vasarely, Shepard, and Kim have also experimented with recursive patterns. I'll provide computational recipes for some intriguing, yet simple to compute, designs.
To create the intricate forms, start with a collection of squares called the initiator lattice or array. You can see what these look like in the upper-left corners in Figures 2 and 3, for instance. The initial collection of squares represents one size scale. At each filled (black) square in the initial array I place a small copy of the filled array. This is the second-size scale. At each point in this new array, I place another copy of the initial pattern. This is the third-size scale. In practice, I only use three size scales for computational speed, and because an additional-size scale does not add much to the beauty of the final pattern.
In mathematical terms, begin with an SxS square array, (A), containing all 0s to which 1s, representing filled squares or sites, are added at random locations. For example:
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 1 1 1 0 0
0 0 0 1 0 0 0 0
0 0 0 1 0 0 0 0
0 1 1 1 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
Just how many patterns can you create by randomly selecting array locations and filling them with 1s? Think of the process of filling array locations in terms of cherries and wineglasses. Consider an SxS grid of beautiful crystal wineglasses. Throw M cherries at the grid. A glass is considered occupied if it contains one or more cherries. With every throw each cherry goes into one of the glasses. How many different patterns of occupied glasses can you make? (A glass with more than one cherry is considered the same as a glass with one cherry in the pattern.)
It turns out that for an SxS array and M cherries, the number of different patterns is shown in Figure 4(a). As an example of how large the number of potential patterns is, consider that 32 cherries thrown at a 9x9 grid creates more than 1022 different patterns. This is far greater than the number of stars in the Milky Way galaxy (1012) and greater than the number of atoms in a person's breath (10x1021). In fact, it is about equal to the estimated number of stars in the universe (1022).
For patterns like Figures 2 and 3, I use S=7. Smaller arrays would lead to fewer potential patterns (particularly with the added induced symmetry, discussed later), and greater values of S lead to diffuse patterns with the scaling used. In C, the process of filling the initial array can be coded as in Example 2(a). Notice that lines 6--9 introduce a fourfold--symmetrical pattern which leads to an overall symmetry in the design at several size scales. Since each of the "cherries" is symmetrically placed in each of the four quadrants of the initial array, each pattern is really defined by the smaller quadrant subpattern. Although this symmetrization decreases the number of possible patterns to that in Figure 4(b).
For even values of S, use Floor((S+1)/2), where "Floor(x)" returns the largest integer value not greater than the x parameter. I find that inducing the symmetry of the original pattern produces more aesthetically appealing designs than those produced by a random initial array.
For convenience, a program may store the values of the A array in two 1-D arrays, x and y, whose values are centered at the origin. The final value of N is the number of filled points in the symmetrized array in Example 2(b). The major computational step begins with the generation of the three size scales from these x and y arrays, the values of which store the positions of the filled (one) sites in the original structural motif. To determine the final positions of all the black squares in the final design, a scale factor 
is applied to each position and the resulting terms summed; see Figure 4(c). For the patterns in Figures 2 and 3, 1=1,2=S,3=S2. To do this in a program, see Example 2(c).
Several numerical measures of the patterns can be computed and displayed in an effort to quantify the structure of the patterns (note the graphs in the upper-right of Figures 2 and 3). If these parameters are found to correlate with the "beauty" of certain patterns, computer programs can then examine these parameters and automatically generate classes of patterns of aesthetic value.
Perhaps the most obvious parameter to compute is the fractal dimension D, which characterizes the size-scaling behavior of the pattern. This value gives an indication of the degree to which the pattern fills the plane--how the pattern "behaves" through different magnifications. Luckily, the D-value for these recursive lattice patterns is easy to compute; in general, N=SD (where N is the number of filled sites in the original symmetrical array of squares, and S is the magnification factor used, in this case the same as the size of an edge of the original array, 7). Notice as the number of filled (black) sites (N) in the initial array increases, the dimension (D) increases. In one sense, D quantifies the degree to which the porous patterns fill the plane in which they reside. If all the sites in the initial array are filled, N is 49 and therefore D=2. This makes sense: If the entire plane is filled, the dimension of the object should be 2. In a program, D is calculated from D=log(N)/log(7). In recursive words, dimensions are tangled up like a ball of twine, and all the patterns are neither one nor two dimensions, but somewhere in between.
Another quantitative measure of the patterns' spatial characteristics is p(r), the length-distribution function. The function indicates the distribution of all the interpoint distances, rij, in the pattern. Mathematically speaking, this new function p(r) is defined by Figure 4(d) where N is the number of points in the pattern. The sums are over i and j, and what's being summed is 
(rij, r)=1 if rij=r and 0 if rij
r. In simple English, to compute p(r), select a point in the pattern and compute the distances from that point to every other point in the pattern. Do this for each point in the pattern, and create a graph showing the number of different lengths as a function of each length. To interpret these graphs, the right-most point in the p(r) graph is the maximum distance found in the structure, and the left-most point is the minimum. The most common distance is the highest point of the curve.
In a program, I use a Monte Carlo approach to compute p(r), because actually computing all the interpoint distances is computationally expensive. Instead, you can randomly select pairs of points in the structure and produce the final p(r) curve after 500,000 pairs are cataloged. In C, the Monte Carlo process looks like Example 3(a). Rmax, the longest vector in the structure, or maximum linear dimension, can also be estimated by examining the last non-zero value of p(r). The variable Npts is the number of points in the final recursive lattice pattern. Obviously there is some noise in the Monte Carlo process, and I use a simple nearest-neighbor smoother to reduce the noise so the eye can concentrate on the global structures of the curve; see Example 3(b). If you run the Monte Carlo process twice using different random numbers, the graphics look almost identical.
Another parameter, the radius of gyration, Rg, can be computed from the unsmoothed length-distribution function p(r) in Figure 5. The radius of gyration quantifies the spatial extent of the structure in the plane. Small, compact patterns have small values of Rg. LaRge, extended patterns have laRge values of Rg. To get an intuitive feeling for this parameter, if all of the squares in the final pattern were to lie on the edge of a circle of radius R, then Rg=R. If the circle were stretched in one direction to form an ellipse, the Rg value would increase because the "mass" of the pattern is further from the center of mass of the entire pattern. In C, the Rg computation looks like Example 4(a).
After examining hundreds of recursive lattice patterns, I find that many people prefer high values of the fractal dimension of around 1.8. The p(r) curves for the preferred structures usually do not exhibit global features (bumps and valleys). To date, I haven't been able to determine a correlation between perceived beauty and the Rg parameter. You may wish to compute the lattice patterns and also look for correlations between perceived beauty and the fractal dimension.
As indicated earlier, most people seem to prefer patterns with a symmetrical structure over those with purely random structures. I prefer patterns with fourfold symmetry; however, I've also experimented with patterns with inversion symmetry, bilateral symmetry, and random-walk symmetry. Figure 3, for example, is a bilaterally symmetric pattern produced by making the left and right sides contain the same patterns: Aij=1, AS-1-i,j=1. Example 4(b) shows this in C. Inversion symmetry can be computed by: Aij=1, AS-1-i, S-1-j=1. Example 4(c) shows this in C.
Random walks can be used to force greater correlation between points in the initial array. Rather than selecting each site randomly with random-walk symmetry, the position of each new site is related to the previous. For example, select a point and continually subtract 1 from, or add 1 to, the initial site's x and y coordinates. Other symmetries, such as the sixfold symmetry of a snowflake, can also be used.
(a)
1 Factorial(X);
2 IF X=0 then Factorial=1
3 ELSE Factorial=X*Factorial(X-1)
4 END
(b)
1 Fibonacci(X)
2 IF N <= 1 then Fibonacci=1
3 ELSE Fibonacci=Fibonacci(N-1)+Fibonacci(N-2)
4 END
(c)
Fibonacci
F[0]=1;F[1]=1;
For i = 2 to 30
F[i]=F[i-1]+F[i-2]
END
END
1 S = 7; Sz = S-1; M = 20;
2 for (h=1; h<=M; h++) {
3 /* rand returns a value between 0 and 32767 */
4 v = ((float) rand()/32767.)*S;
5 w = ((float) rand()/32767.)*S;
6 a[v][w]=1;
7 a[Sz-v][w]=1;
8 a[Sz-v][Sz-w]=1;
9 a[v][Sz-w]=1;
10 }
(b)
N=0; Sz = S - 1;
for (i=0; i<=S; i++)
for (j=0;j<=S; j++)
if (a[i][j]==1)
{N++; x[h]=i-Sz/2; y[h]=j-Sz/2;}
(c)
for (i=1; i<=N; i++) {
for (j=1; j<=N; j++) {
for (k=1; k<=N; k++) {
X = alpha[1]*x[h]+alpha[2]*x[h]+alpha[3]*x[h]
Y = alpha[1]*y[h]+alpha[2]*y[h]+alpha[3]*y[h]
PlotSquareAt(X,Y)
}
}
}
(a)
/* initialize p array to zero */
for(i=0; i<900; i++) p&lbrk.i&rbrk.=0;
total=500000;
/* catalog 500000 vector lengths */
for(i=0; i<total; i++) {
rnd1 = ((float) rand()/32767.)*Npts;
rnd2 = ((float) rand()/32767.)*Npts;
xterm = corx[rnd1]-corx[rnd2];
yterm = cory[rnd1]-cory[rnd2];
dist = sqrt (xterm*xterm+yterm*yterm);
dist = dist+.5;
p[(int)dist]++;
}
(b)
for(i=1; i<900;i++) {
ps[i]=.25*(float)p[i-1]+.5*(float)p[i]+.25*(float)p[i+1];
}
(a)
sum=0; summ=0;
for(i=0; i<900;i++) {
summ=summ+p[i];
sum=sum+p[i]*i*i;
}
Rg=sqrt(sum/((float) summ*2)) ;
(b)
Sz = S - 1;
for (h=1; h<=40; h++) {
v = ((float) rand()/32767.)*S;
w = ((float) rand()/32767.)*S;
a[v][w]=1;
a[Sz-v][w]=1;
}
(c)
Sz = S - 1;
for (h=1;i h<=30; h++) {
v = ((float) rand()/32767.)*S;
w = ((float) rand()/32767.)*S;
a[v][w]=1;
a[Sz-v][Sz-w]=1;
}
Figure 1: Coral-like form generated from a recursive program
Figure 2: Recursive lattice design using a fourfold-symmetrical pattern
Figure 3: Recursive lattice design using a bilaterally symmetric pattern
Figure 4: (a) The number of different patterns for an SXS array; (b) symmetrization decreases the number of possible patterns; (c) determining the final positions of all the black squares in the final design; (d) p(r), the length-distribution function, is a quantitative measure of the patterns' spatial characteristics
Figure 5: The radius of gyration R sub g, can be computed from the unsmoothed length-distribution function p(r).
Copyright © 1993, Dr. Dobb's Journal