I was struck by the lunatic programming muse again today. While reading my twitter feed, I encountered a description of SmoothLife, a generalization of Conway’s classic Game of Life. Instead of being implemented on a grid of discrete binary values, SmoothLife is implemented over a continuous, real valued field. What’s kind of neat is that “gliders” exist in this world, but unlike the discrete version, they can travel at any angle. I found it really intriguing. Witness!

Pretty neat organic shapes.

In reading up on how this works, I discovered that the implementation was unusual in a way I had not considered. It dawned on me that it would be possible to implement the “summation of neighbors” part of Conway’s life using a simple convolution. As a bonus, we’d get the “wrap around” for free: the indexing would be pretty simple. I hypothesized that I could write it using the fftw library in just a few hundred lines. So I did.

Here’s the code!

[sourcecode lang=”c”]
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <unistd.h>
#include <fftw3.h>

#define SIZE (512)
#define SHIFT (18)

fftw_complex *filter ;
fftw_complex *state ;
fftw_complex *tmp ;
fftw_complex *sum ;

int
main(int argc, char *argv[])
{
fftw_plan fwd, rev, flt ;
fftw_complex *ip, *jp ;
int x, y, g ;

srand48(getpid()) ;

filter = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ;
state = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ;
tmp = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ;
sum = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ;

flt = fftw_plan_dft_2d(SIZE, SIZE,
filter, filter, FFTW_FORWARD, FFTW_ESTIMATE) ;
fwd = fftw_plan_dft_2d(SIZE, SIZE,
state, tmp, FFTW_FORWARD, FFTW_ESTIMATE) ;
rev = fftw_plan_dft_2d(SIZE, SIZE,
tmp, sum, FFTW_BACKWARD, FFTW_ESTIMATE) ;

/* initialize the state */
for (y=0, ip=state; y<SIZE; y++) {
for (x=0; x<SIZE; x++, ip++) {
*ip = (fftw_complex) (lrand48() % 2) ;
}
}

/* initialize and compute the filter */

for (y=0, ip=filter; y<SIZE; y++, ip++) {
for (x=0; x<SIZE; x++) {
*ip = 0. ;
}
}

#define IDX(x, y) (((x + SIZE) % SIZE) + ((y+SIZE) % SIZE) * SIZE)
filter[IDX(-1, -1)] = 1. ;
filter[IDX( 0, -1)] = 1. ;
filter[IDX( 1, -1)] = 1. ;
filter[IDX(-1, 0)] = 1. ;
filter[IDX( 1, 0)] = 1. ;
filter[IDX(-1, 1)] = 1. ;
filter[IDX( 0, 1)] = 1. ;
filter[IDX( 1, 1)] = 1. ;

fftw_execute(flt) ;

for (g = 0; g < 1000; g++) {
fprintf(stderr, "generation %03d\r", g) ;

fftw_execute(fwd) ;

/* convolve */
for (y=0, ip=tmp, jp=filter; y<SIZE; y++) {
for (x=0; x<SIZE; x++, ip++, jp++) {
*ip *= *jp ;
}
}

/* go back to the sums */
fftw_execute(rev) ;

for (y=0, ip=state, jp=sum; y<SIZE; y++) {
for (x=0; x<SIZE; x++, ip++, jp++) {
int s = (int) round(creal(*ip)) ;
int t = ((int) round(creal(*jp))) >> SHIFT ;
if (s)
*ip = (t == 2 || t == 3) ;
else
*ip = (t == 3) ;
}
}

/* that’s it! dump the frame! */

char fname[80] ;
sprintf(fname, "frame.%04d.pgm", g) ;
FILE *fp = fopen(fname, "wb") ;
fprintf(fp, "P5\n%d %d\n%d\n", SIZE, SIZE, 255) ;

for (y=0, ip=state; y<SIZE; y++) {
for (x=0; x<SIZE; x++, ip++) {
int s = ((int) creal(*ip)) ;
fputc(255*s, fp) ;
}
}

fclose(fp) ;
}
fprintf(stderr, "\n") ;

return 0 ;
}
[/sourcecode]

Pardon the slightly odd code, but I wrote it while watching the VP debates and the last of the Game 5 between Detroit and the Athletics. (Thanks to the Athletics for a great season!)

It’s kind of a ridiculous way to implement Conway’s Life (there are more efficient and compact implementations), but it has one interesting feature that SmoothLife (or any other discrete life with larger neighborhoods) would use: the run time doesn’t increase with increasing size of neighborhood. I’ll try to work up a SmoothLife implementation soon.

It was just a bit of thought provoking fun, like my previous attempts at implementing bignum arithmetic with the FFT.

Addendum: Each frame (512×512 resolution) takes about 48ms on my MacBook. Not really competitive with more conventional implementations, but not miserable either.