Category Archives: Games and Diversions

Code for generating a random Latin square…

The other day I mentioned that generating random Latin squares was a bit more complicated than I thought, and that an algorithm by Jacobson and Matthews was the way that people typically did it. I worked up this implementation based on a couple of different descriptions of the algorithm (the original paper was behind a paywall). The code below simply generates a random DIMxDIM Latin square using their algorithm. I haven’t done extensive testing, but it appears to mostly work. I did run a test where I generated one million 4×4 matrices, and found that it generated all 576 different matrices in approximately the same propertion, so it appears to at least be roughly operational. I’ll eventually merge the code with my nascent KenKen generator.

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

/*
 * rls.c
 *
 * A straightforward implementation of the algorithm by Jacobson and 
 * Matthews, _Generating uniformly distributed random Latin square_,
 * Journal of Combinatorial Design, 1996.
 *
 * I actually couldn't find the original paper, so this is an implementation
 * based upon other descriptions of the algorithm.
 *
 * Written by Mark VandeWettering.
 *
 */


#define DIM (6)

typedef struct t_incidence_cube {
    int s[DIM][DIM][DIM] ;
} IncidenceCube ;


void
InitIncidenceCube(IncidenceCube *c)
{
    int i, j, k ;

    /* First, zero the cube */
    for (i=0; i<DIM; i++)
        for (j=0; j<DIM; j++)
            for (k=0; k<DIM; k++)
                c->s[i][j][k] = 0 ;

    /* And then dump ones in the right place... */
    for (i=0; i<DIM; i++)
        for (j=0; j<DIM; j++)
            c->s[i][j][(i+j)%DIM] = 1 ;
}

void
PrintIncidenceCube(IncidenceCube *c)
{
    int i, j, k ;

    for (i=0; i<DIM; i++) {
        for (j=0; j<DIM; j++) {
            for (k=0; k<DIM; k++) {
                if (c->s[i][j][k]) {
                    printf("%2d ", k+1) ;
                    break ;
                }
            }
        }
        printf("\n") ;
    }
    printf("\n") ;
}

void
CompactPrintIncidenceCube(IncidenceCube *c)
{
    int i, j, k ;

    for (i=0; i<DIM; i++) {
        for (j=0; j<DIM; j++) {
            for (k=0; k<DIM; k++) {
                if (c->s[i][j][k]) {
                    printf("%01d", k+1) ;
                    break ;
                }
            }
        }
    }
    printf("\n") ;
}

void
ShuffleIncidenceCube(IncidenceCube *c)
{
    int i, j ;

    for (i=0; i<DIM*DIM*DIM; i++) {
        
        // Assume we have a "proper" matrix...
       
        // Pick a random zero from the incidence cube... 
        int rx, ry, rz ;
        do {
            rx = lrand48() % DIM ;
            ry = lrand48() % DIM ;
            rz = lrand48() % DIM ;
        } while (c->s[rx][ry][rz]) ;

        int ox, oy, oz ;

        for (j=0; j<DIM; j++) {
            if (c->s[j][ry][rz] == 1)
                ox = j ;
            if (c->s[rx][j][rz] == 1)
                oy = j ;
            if (c->s[rx][ry][j] == 1)
                oz = j ;
        }

        // do the +/- 1 move...  
        // These are all square with zeros in them...
        c->s[rx][ry][rz] ++ ;
        c->s[rx][oy][oz] ++ ;
        c->s[ox][ry][oz] ++ ;
        c->s[ox][oy][rz] ++ ;

        // These all have ones, except for maybe the last one...
        c->s[rx][ry][oz] -- ;
        c->s[rx][oy][rz] -- ;
        c->s[ox][ry][rz] -- ;
        c->s[ox][oy][oz] -- ;

        while (c->s[ox][oy][oz] < 0) {

            rx = ox ;
            ry = oy ;
            rz = oz ;

            // Pick one of the two 1's that are in conflict
            if (drand48() < 0.5) {
                for (j=0; j<DIM; j++) {
                    if (c->s[j][ry][rz] == 1) {
                        ox = j ;
                    }
                }
            } else {
                for (j=DIM-1; j>=0; j--)  {
                    if (c->s[j][ry][rz] == 1) {
                        ox = j ;
                    }
                }
            }

            if (drand48() < 0.5) {
                for (j=0; j<DIM; j++) {
                    if (c->s[rx][j][rz] == 1) {
                        oy = j ;
                    }
                }
            } else {
                for (j=DIM-1; j>=0; j--)  {
                    if (c->s[rx][j][rz] == 1) {
                        oy = j ;
                    }
                }
            }

            if (drand48() < 0.5) {
                for (j=0; j<DIM; j++) {
                    if (c->s[rx][ry][j] == 1) {
                        oz = j ;
                    }
                }
            } else {
                for (j=DIM-1; j>=0; j--)  {
                    if (c->s[rx][ry][j] == 1) {
                        oz = j ;
                    }
                }
            }
          
            // do the +/- 1 move...  
            // These are all square with zeros in them...
            c->s[rx][ry][rz] ++ ;
            c->s[rx][oy][oz] ++ ;
            c->s[ox][ry][oz] ++ ;
            c->s[ox][oy][rz] ++ ;

            // These all have ones, except for maybe the last one...
            c->s[rx][ry][oz] -- ;
            c->s[rx][oy][rz] -- ;
            c->s[ox][ry][rz] -- ;
            c->s[ox][oy][oz] -- ;
        }

    }
}

int
main()
{
    IncidenceCube c ;
    int i ;

    srand48(getpid()) ;
    InitIncidenceCube(&c) ;
    ShuffleIncidenceCube(&c) ;
    PrintIncidenceCube(&c) ;
    return 1 ;
}

Atinlay Aresquares…

My last post dealt with a solver for KenKen puzzles. Once you have one of those, then the obvious thing to work on next (in your copious spare time) is a generator for KenKen puzzles. It didn’t seem too hard. You’d begin by generating a random Latin square, then divide it up into “cages”, assign an operator to them, and then see if that has a unique solution. If it doesn’t, well, try again.

But it turns out that each of these steps hides inner difficulties, as well as a certain amount of mathematical beauty. Tom and I were discussing the first bit (generating random Latin square) over lunch, and I thought I’d summarize some of our discussion.

First of all, consider the following “canonical” Latin square of order 4:

1 2 3 4
4 1 2 3
3 4 1 2 
2 3 4 1

It’s pretty easy to generate this Latin square, and it’s pretty easy to generate other Latin squares by swapping either rows or columns. But the question is, can I generate all other Latin squares from this one by doing row and column swap operations?

Tom assured me that the answer was no, which wasn’t at all apparent to me until he presented an interesting example, so I thought I would summarize some of the ideas here. My presentation will take a slightly different path than our discussion, but it ends in the same place.

Let’s consider each square to have a canonical form, where the first row and first column are in sorted order. It should be relatively obvious that you can place any matrix into this form by a series of row and column operations. For example, the matrix above can be placed in canonical form by swapping the second and fourth row.

1 2 3 4
2 3 4 1
3 4 1 2
4 1 2 3 

But let’s consider a different matrix in canonical form formed by swapping 1 and 3 in rows 2 and 4.

1 2 3 4
2 1 4 3
3 4 1 2 
4 3 2 1 

We can’t convert this matrix into the previous one by only swapping rows and columns. Thus, if we try to generate a random Latin square by just beginning with a canonical one, there are possible Latin squares which can never be generated by just swapping row and column operations. Too bad, since this is really easy to implement.

Another easy thing to implement would be to just generate all possible Latin squares, store them in a table, and select one at random. That doesn’t seem too bad, until you consider that there are 812,851,200 possible Latin Squares to consider. Sure, disks are big, but that seems excessive.

It turns out that the algorithm that most people seem to use is the from 1996 by Jacobson and Matthews, which is hard to find online as a full paper. It uses a similar idea, but instead of swapping full rows and columns, it applies a different operator that mutates a given Latin square, and in the limit produces a distribution which is the same as the random distribution. It’s a little tricky to understand, but I am going to try to have an implementation in Python later tonight.

Here’s a report that provides some details and the skeleton of an implementation in Java.

KenKen puzzle solver…

Lately, my lunch hours have been spent working on the NYT Crossword with my lunch companion Tom. While I find that the Thursday crosswords are often beyond my ability to do in the allotted time, between the two of us, more often than not we manage to plow through them. Slowly over time, we’ve begun to solve them slightly quicker, so I’ve branched out to doing the KenKen puzzles.

For those of you who haven’t seen these before, you can learn about them on their official website here. The 4×4 ones are usually pretty easy, but the last couple days have done some 6×6 puzzles have been shockingly hard. So much so, that today I got to this state on today’s puzzle:

IMG_1128

Yes, I was using pen. And I made a mistake, although I can’t see what I was thinking when I did. And because I made a mistake, I couldn’t spot it, and I was hung up. The “4” in first column was wrong. I knew it could not be a 6 or 1, or a 3, or a 4, but for some reason I thought that the final column couldn’t be 4 and 3, but..

Never mind, it doesn’t matter what the mistake was.

By the end up twenty minutes I was annoyed, knew I had made a mistake, but decided to give into the urge to write a program to solve it. How long could it take?

Well, the answer was about 27 minutes.

This program doesn’t have any interface or flexibility, everything is hard coded. Basically it uses backtracking and fills in each square in order from bottom left to top right, and as the top most, right most square of each outlined area is finished, it checks to make sure that the sums/differences/whatever for the filled in region sum to the right values. I didn’t know how I was going to encode those constraints, so I just ended up writing a function that computes the proper constraint. This is not really the right way to do it: it’s error prone and requires recompilation of the code to change the puzzle. But it does work fairly well, once you type the expressions in correctly.

Once I got those straightened out, it blats out the solution in less than 1 millisecond on my desktop.

+---+---+---+---+---+---+
| 6 | 1 | 5 | 2 | 4 | 3 |
+---+---+---+---+---+---+
| 2 | 6 | 1 | 3 | 5 | 4 |
+---+---+---+---+---+---+
| 1 | 4 | 6 | 5 | 3 | 2 |
+---+---+---+---+---+---+
| 3 | 2 | 4 | 6 | 1 | 5 |
+---+---+---+---+---+---+
| 5 | 3 | 2 | 4 | 6 | 1 |
+---+---+---+---+---+---+
| 4 | 5 | 3 | 1 | 2 | 6 |
+---+---+---+---+---+---+

1 solutions found

It points out that the 4 should have been a 2 (IDIOT!) and then it all works out.

The code is 178 lines long, but is not particularly pretty/short/adaptable. Some day, I’ll have to tidy it up and make a better version.

But in case people find this stuff interesting, here it is. Remember: 27 minutes to write.

include <stdio.h>
#include <stdlib.h>
#include <string.h>

/*
 * kenken.c
 *
 * A program inspired by a particular tough 6x6 puzzle that appeared
 * on Thursday, May 11, 2016, which Tom Duff and I could not appear to
 * get the answer to during our lunch hour.  I thought that it was possible
 * that I could write a program to solve the puzzle in less than an hour.
 *
 * Hence, this code.
 *
 */

int total = 0 ;

#define DIM     (6)

typedef struct t_puzzle_state {
    int sq[DIM][DIM] ;
    int r[DIM] ;            // bitflag for values left in rows...
    int c[DIM] ;            // bitflag for values available in columns
} PuzzleState ;


void
PrintPuzzleState(PuzzleState *p)
{
    int i, j ;

    printf("+") ;
    for(j=0; j<DIM; j++)
        printf("---+") ;
    printf("\n") ;

    for (i=0; i<DIM; i++) {
        printf("|") ;
        for (j=0; j<DIM; j++) {

            if (p->sq[DIM-1-i][j])
                printf(" %1d |", p->sq[DIM-1-i][j]) ;
            else
                printf("   |") ;
        }
        printf("\n+") ;
        for(j=0; j<DIM; j++)
            printf("---+") ;
        printf("\n") ;
    }
    printf("\n") ;
}

void
InitializePuzzleState(PuzzleState *p)
{
    int i, j ;

    memset((void *) p, 0, sizeof(*p)) ;

    for (i=0; i<DIM; i++) {
        for (j=1; j<=DIM; j++) {
            p->r[i] |= (1<<j) ;
            p->c[i] |= (1<<j) ;
        }
    }
}

void
HardcodePuzzleState(PuzzleState *p, int r, int c, int val)
{
    p->sq[r][c] = val ;
    p->r[r] &= ~(1<<val) ;
    p->c[c] &= ~(1<<val) ;
}

#define SQUARE(i, j)    ((i)*DIM+(j))
#define ABS(x)          ((x)<0?(-(x)):(x))

int
ConstraintSatisfied(PuzzleState *p, int i, int j)
{
    switch (SQUARE(i,j)) {
    case SQUARE(0,1):
        return ABS(p->sq[0][0] - p->sq[0][1]) == 1 ;
    case SQUARE(0,2):
        return p->sq[0][2] == 3 ;
    case SQUARE(1, 3):
        return (p->sq[0][3]+p->sq[1][3]) == 5 ;
    case SQUARE(1, 4):
        return (p->sq[0][4] * p->sq[0][5] * p->sq[1][4]) == 72 ;
    case SQUARE(2, 2):
        if (p->sq[2][2] == 2 * p->sq[1][2])
            return 1 ;
        if (p->sq[1][2] == 2 * p->sq[2][2])
            return 1 ;
        return 0 ;
    case SQUARE(2, 4):
        return ABS(p->sq[2][4] - p->sq[2][3]) == 5 ;
    case SQUARE(3, 1):
        if (p->sq[3][1] == 2 * p->sq[2][1])
           return 1 ;
        if (p->sq[2][1] == 2 * p->sq[3][1])
           return 1 ;
        return 0 ;
    case SQUARE(3, 5):
        return (p->sq[3][5]+p->sq[2][5]+p->sq[1][5]) == 8 ;
    case SQUARE(4, 2):
        return ABS(p->sq[4][2] - p->sq[4][1]) == 5 ;
    case SQUARE(4, 3):
        return p->sq[4][3] == 3 ;
    case SQUARE(4, 4):
        return (p->sq[4][4]+p->sq[3][2]+p->sq[3][3]+p->sq[3][4]) == 19 ;
    case SQUARE(5, 1):
        return (p->sq[4][0]*p->sq[5][0]*p->sq[5][1]) == 12 ;
    case SQUARE(5, 4):
        return (p->sq[5][2]+p->sq[5][3]+p->sq[5][4]) == 11 ;
    case SQUARE(5, 5):
        return (p->sq[5][5]+p->sq[4][5]) == 7 ;
    default:
        return 1 ;
    }

    return 1 ;
}

void
SolvePuzzleState(PuzzleState *p, int i, int j)
{
    int v ;

    if (i >= DIM) {
        PrintPuzzleState(p) ;
        total ++ ;
        return ;
    }
        
    int pval = p->r[i] & p->c[j] ;

    for (v=1; v<=DIM; v++) {
        if ((pval & (1<<v)) == 0)
            continue ;
        // Try to place v in sq[i][j] ;
        p->sq[i][j] = v ;

        if (ConstraintSatisfied(p, i, j)) {
            p->r[i] &= ~(1<<v) ;
            p->c[j] &= ~(1<<v) ;

            int ni, nj ;
            nj = j + 1 ;
            if (nj < DIM) {
                ni = i ;
            } else {
                ni = i + 1 ;
                nj = 0 ;
            }

            SolvePuzzleState(p, ni, nj) ;

            p->r[i] |= (1<<v) ;
            p->c[j] |= (1<<v) ;
        } 
    }
    p->sq[i][j] = 0 ;
}

main()
{
    PuzzleState p ;
    InitializePuzzleState(&p) ;
    // HardcodePuzzleState(&p, 0, 2, 3) ;
    // HardcodePuzzleState(&p, 4, 3, 3) ;

    SolvePuzzleState(&p, 0, 0) ;
    fprintf(stderr, "%d solutions found\n", total) ;
}

Combinations…

This is a math/geeky/computer post of something relatively simple, but it arose in the wild in a program that I wrote several years ago, and when I saw it again, it confused me, and I spent a bit of time thinking about it, and I thought I might just write it up. If math/geeky/computer stuff bores you, perhaps you should move on.

If I’ve still got you, today’s post is about combinations, a topic that arises in lots and lots of computer programs. For instance, the program that spawned this post was this program for computing the number of possible positions in a game of checkers that I did back in 2009. It uses a function called “comb” which takes two arguments a and b, which computes the number of different ways of selecting b objects from a collection of a objects. Here’s the implementation from that post:

def fact(n):
    if n == 0:
        return 1 
    else:
        return n*fact(n-1)

def comb(a, b):
        return fact(a) / (fact(b)*fact(a-b))

If you’ve done a lot of programming, this probably seems fairly obvious to you. This is one of those formulas that I don’t have to think about: I learned it so long ago that I don’t even remember how I originally learned it. It might have been when I first learned about Pascal’s Triangle back in grade school. I probably learned the formula to compute these coefficients using the factorial function, and I’ve probably used it dozens and dozens of times in various contexts. But really, when I thought about it, the factorial function doesn’t seem to me to be completely obvious when staring at the triangle. When I first learned about Pascal’s triangle, the rules were something like this. First, write down a at the top, then write a bunch of rows underneath it, each row containing one more number so they form a triangle. The first and last number in each row needs to be a one. The other entries are the sum of the numbers above and to the left and above and to the right. If you think of the number of path starting at the root, and reaching the desired new node, it’s clear that it’s the sum of the number of paths reaching the node above and to the left, and above and to the right.

Pascal's Triangle

Pascal’s Triangle

Which brings me to a different chunk of code. In the posting above, I had mentioned that I had written a program to compute the number of possible checkers positions before, but had lost the code, causing me to write the code above. Butr recently I found that original code. But this code didn’t contain the factorial/comb implementation above. It had code that looked like this.

def c(a, b):
    if b == 0 or b == a:
        return 1
    r = c(a-1, b-1) + c(a-1, b)
    return r

At first, this seemed less familiar to me. But if you think about Pascal’s triangle, this is exactly what I described above. If you think of trying to compute the bth number in the ath row, it’s just the sum of the numbers in the row above, to the left and right.

And this works! If you compute c(52, 5), you can verify that there are 2598960 ways to deal a 5 card poker hand from a single deck. In a way, this code seems more elegant to me: what was with all that multiplication and division stuff up above? But there is a problem: performance.

A moment’s thought should reveal why that should be the case. This code only returns 1, or the sum of two other terms. That means to come up with the 2598960, it bottoms out that recursion… 2598960 times. Modern computers can count that high, even by ones, even in python rather quickly. It took about 1.4s to compute. But to compute c(52, 6) takes about 7 times longer (because the number is roughly 7x larger). Computing c(52,13) is out of the question with that code.

But my original program computed the number of combinations this way, and was reasonably efficient. What gives?

Well, I lied. The code in my npos program wasn’t exactly like that. Here’s the complete code:

#!/usr/bin/env python
#
# some programs to compute the number of legal positions in checkers
#

import locale

locale.setlocale(locale.LC_ALL, "")

mem = {}

def c(a, b):
    if b == 0 or b == a:
	return 1 
    if mem.has_key((a, b)):
	return mem[(a, b)]
    r = c(a-1, b-1) + c(a-1, b)
    mem[(a, b)] = r 
    return r

def num(b, bk, w, wk, f):
    r = c(4,f)*c(24,b-f)*c(28-(b-f),w)*c(32-b-w,bk)*c(32-b-w-bk,wk)
    return r

def np(n):
    total = 0 
    for b in range(min(n, 12)+1):
	for bk in range(min(n, 12)-b+1):
	    for w in range(min(n-b-bk, 12)+1):
		for wk in range(min(n-b-bk, 12)-w+1):
		    for f in range(min(b, 4)+1):
			total += num(b, bk, w, wk, f)
    return total - 1 

tot = {}

for n in range(25):
	tot[n] = np(n)

maxnew, maxtot = 0, 2**64-1
for n in range(1, 25):
	new = tot[n] - tot[n-1]
	newstr = locale.format("%d", new, True)
	totstr = locale.format("%d", tot[n], True)
	maxnew = max(maxnew, len(newstr))
	maxtot = max(maxtot, len(totstr))

print " N%s" % "possible positions".center(maxnew)
print (maxnew+3)* '-'
for n in range(1, 25):
	new = tot[n] - tot[n-1]
	newstr = locale.format("%d", new, True)
	print "%2d %s" % (n, newstr.rjust(maxnew))
print (maxnew+3)* '-'
totstr = locale.format("%d", tot[24], True)

Can you spot the difference? The c function is memoized. The vast majority of invocations of the function have been evaluated once before. So, this code remembers the value of each invocation, and before computing, looks up the result to see if we’ve computed it before. An obvious way to implement this lookup would be as a two dimensional array, but Python’s dictionary types are quick and easy. When the program is run, it stores stores 306 values. With this added, the entire program, which counts 500,995,484,682,338,672,639 positions executes in just about 1 second on my laptop. And while it does rely on bignum arithmetic, it doesn’t require multiplication or division: just addition works fine.

I like it.

Some Musings on Crossword Puzzles…

puzzleWhen I was growing up, my Grandma Busch used to spend time each morning doing crossword puzzles. Off and on through my life, probably in imitation of her, I’ve enjoyed sitting down and doing them as well. The Oakland Tribune (the paper I see most often) carries two different puzzles per day: the Daily Commuter, and also the New York Times Crossword puzzle (an older one, not the current one in the NYT). I find the Daily Commuter to be pretty easy: I’m not good enough that I can just write out the answers mostly in order, I usually alternate horizontal and vertical, and in ten or twelve minutes, I usually power through them.

The NYT is a different story though. Early in the week, they are pretty simple. The Monday puzzles are just as easy as the Daily Commuters. But as the week goes on, my time to complete these puzzles goes up. By the time Thursday or Friday rolls around, I don’t often finish them without some judicious web searching. I find those puzzles often rely on what seems to me fairly obscure references to music, literature or geography. It actually makes them somewhat less fun for me. I seldom do the Sunday Crosswords for this reason: there is just something rather crude and unaesthetic to these puzzles.

Ultimately, there is a kind of beauty in the very best crosswords. The clues are hard, but not too hard. The themes are clever. The word choices are good. I like it when there are a few words that I’ve never seen before. I find it especially interesting when these are short words. Last week I encountered the word ORT (a morsel left after a meal) which I had never seen before. Today, I encountered the word TUTEE, which seems to me a rather odd word, but also one I’d never seen before. And of course you occasionally get to use words that seem cool, but aren’t commonly used. Today, one of the answers was a great word: PALAVER. It just doesn’t get used often enough.

The thing that wasn’t obvious to me (and may not be obvious to the casual puzzler) is that there are definitely “good” and “bad” puzzles. Bad puzzles just seem like they were assembled by machine, with no more pattern than the average Sudoku (I’m not very good at Sudoku, so maybe I’m just not educated enough to see beauty in them). But I never thought that there was actual criticism of crosswords on the net…

Enter Rex Harper does the NY Times Crossword Puzzle. Each day, he bats out (probably faster than I find a pencil) the NYT crossword, blogs the answers, and then, somewhat interestingly, presents his critique of the quality of the puzzle. For instance, today’s Tribune republished the NYT from November 12, 2013, which Rex blogged about here.. I must admit, I didn’t think the puzzle was that amazing, I got through it in about 22 minutes, with a question about the intersection of the clue for 55 Down (When the balcony scene occurs in “Romeo and Juliet”) and 67 Across (Lacking depth… or like 17- 23- 37- 48- and 58-Across). My recollection was that the balcony scene had to occur in Act II, but the spacing left only a single space. Was it really in Act I? I found that hard to believe. And the across clue, what was that about? I had discovered that all the clues had a pair of “Ds” to introduce them, which was the last letter of the four letter answer. INID? That didn’t seem right. I scratched my head. I finally mused upon the answer that it really was Act 2, but written “ACT2” making the cross clue “IN2D”. I didn’t like that, without a greater theme that suggests numbers are legitimate in the puzzle, this seemed arbitrary. Annoying. Rex apparently thought so too. He says…

Nobody says “In 2D” except, perhaps, ironically. “I’m gonna go see ‘My Dinner With Andre’ in glorious 2D!” So the revealer is contrived and awkward. Also, a bra-size revealer is going to be soooo much better here. You’d have to ditch DOUBLE DUTCH, but surely there are other DD answers out there for you to use. Friend of mine suggested DDS would work as a revealer as well (interesting twist on the dental degree, which is a common crossword answer). Not a big fan of this revealer, as I’ve now made overly clear.

I suspect that when I get better at these, I’ll be an even greater critic.

Addendum: The NYT Crossword App for the iPhone/iPad is great, especially on the iPad (harder to type on the iPhone, at least for my meaty fists).

Crazy programming experiment of the evening…

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!

#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 ;
}

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.

Donald Michie, Alan Turing, Martin Gardner, and Tic Tac Toe

As anyone who reads my blog with any regularity will tell you, I like to read and learn new things. The problem with being self taught and also easily distracted means that you often learn a great deal, but don’t always perceive the connections and scope of what you are learning. I found another example today while surfing.

Years ago, I remember reading one of Martin Gardner’s Mathematical Games columns (from March, 1962, in case you want to look it up) where he described an interesting machine for playing tic-tac-toe. It was made entirely out of matchboxes, each one of which had a tic tac toe position on the top. Inside was a collection of colored beads. Each color specified a possible legal move for the position on top. The idea was that you’d play a game by drawing these beads from the appropriate box, and making the appropriate move. At the end of the game, you’d remove the bead from the last box that sent you along the losing path. Eventually, all the losing moves get removed, and the machine plays perfect tic-tac-toe. Gardner showed how this same idea could be used to create a matchbox computer to play hexapawn, a simple game played with six pawns on a 3×3 board.

I really haven’t given it much thought since then. Many of you have probably read this article in one of the collections of Gardner’s columns.

But today, I was surfing through links and reread some of the details. I found that the machine was called MENACE (Matchbox Educable Naughts and Crosses Engine) and was invented in 1960 by a gentleman named Donald Michie. And it turns out that he’s a pretty interesting guy.

He was a colleague and friend of Alan Turing, and worked with him at Bletchley Park. Apparently Michie, Turing and Jack Good were all involved in the British code breaking efforts, and in the creation of Collosus, the first digital programmable computer which was used to crack the German “Tunny” teleprinter code. (Good and Michie were apparently two of the authors of the General Report on Tunny, a report on the cracking of the code which has only in recent years become declassified). None of this work could have been known by Martin Gardner at the time of this publication. Of course, this was also true of Turing’s work as well.

Turing made a huge impact in several related disciplines: in mathematical logic and computation, in his wartime efforts in code breaking, and in his role in creating some of the first digital computers. Turing also became interested in mathematics in biology, writing about the chemical foundations of morphogenesis and predicting oscillatory chemical reactions. Michie received a doctorate in biology from Oxford, but returned to have a profound and lasting influence on artificial intelligence. Oh, and that modest little paper on Tic Tac Toe? One of the first instances of reinforcement learning.

Very cool, to discover that the little bit of reading you did as a teen, which seemed like an insignificant game at the time, actually has profound threads which stretch out into lots of different areas.

Donald Michie’s Home Page
Trial and Error, the paper describing MENACE

Bit banger

Anyone who has seen my projects on the Atari 2600 might reasonably conclude that I have a thing for retro computing. The saying goes “it is no virtue to do with more, what can be done with less” and I can’t think of someone whose projects have embodied that more than demo coder Linus Akesson (lft). His latest demo runs on the tiniest of tiny Atmel processors, the ATTINY15. It has just 32 bytes of memory (just the registers, really) and room for 512 instructions. Yet, lft made a cool demo that runs on it, generating both video and sound. Oh, and he’s only running it at 1.6 Mhz (yes, 1.6 Megahertz). Very cool!


More details about Bit banger, and how he did it..

Book Review: Land of LISP, by Conrad Barski

I learned to program as a teenager back in the 1980s, starting as most of a generation of future computer professionals did. I had an Atari 400 at home, and learned to program using the most popular language of the day: BASIC. There were lots of magazines like COMPUTE! and Creative Computing that included program listings for simple games that you could type in and experiment. The interaction was very hands on and immediate.

Sometimes I feel like we’ve lost this kind of “hobbyist programming”. Because programming can be lucrative, we often concentrate on the saleability of our skills, rather than the fun of them. And while the computers are more powerful, they also are more complicated. Sometimes it feels like we’ve lost ground: that to actually be a hobbyist requires that you understand too much, and work too hard.

That’s a weird way to introduce a book review, but it’s that back story that compelled me to buy Conrad Barski’s Land of Lisp. The book is subtitled Learn to Program in LISP, One Game at a Time!, and it’s delightful. The book is chock-full of comics and cartoons, humorously illustrated and fun. But beneath the fun exterior, Barski is an avid LISP enthusiast with a mission: to convince you that using LISP will change the way you program and even the way you think about programming.

I’ve heard this kind of fanatical enthusiasm before. It’s not hard to find evangelists for nearly every programming language, but in my experience, most of the more rare or esoteric languages simply don’t seem to be able to convince you to go through the hassle of learning them. For instance, I found that none of the “real world” examples of Haskell programs in Real World Haskell were the kind of programs that I would write, or were programs that I already knew how to write in other languages, where Haskell’s power was simply not evident. We probably know how to write FORTRAN programs in nearly any language you like.

But I think that Land of Lisp succeeds where other books of this sort fail. It serves as an introductory tutorial to the Common LISP language by creating a series of retro-sounding games like the text adventure game “Wizard’s Adventure” and a web based game called “Dice of Doom”. But these games are actually not just warmed over rehashes of the kind of games that we experimented with 30 years ago (and grew bored of 29 years ago): they demonstrate interesting and powerful techniques that pull them significantly above those primitive games. You’ll learn about macros and higher-order programming. You’ll make a simple webserver. You’ll learn about domain-specific languages, and how they can be used to generate and parse XML and SVG.

In short, you’ll do some of the things that programmers of this decade want to do.

I find the book to be humorous and fun, but it doesn’t pander or treat you like an idiot. While it isn’t as strong academically as The Structure and Interpretation of Computer Programs (which is also excellent, read it!) it is a lot more fun, and at least touches on many of the same topics. I am not a huge fan of Common LISP (I prefer Scheme, which I find to be easier to understand), but Common LISP is a reasonable language, and does have good implementations for a wide variety of platforms. Much of what you learn about Common LISP can be transferred to other LISP variants like Scheme or Clojure.

But really what I like about this book is just the sense of fun that it brings back to programming. Programming after all should be fun. The examples are whimsical but not trivial, and can teach you interesting things about programming.

At the risk of scaring you off, I’ll provide the following link to their music video, which gives you a hint of what you are getting yourself into if you buy this book:



I’d rate this book 4/5. Worth having on your shelf!

Land of LISP, by Conrad Barski

Addendum: Reader “Angry Bob” has suggested that the language itself should be spelled not as an acronym, but as a word, “Lisp”. It certainly seems rather common to do so, but is it wrong to spell it all in caps? LISP actually is derived from LISt Processor (or alternatively Lots of Insufferable Superfluous Parentheses), just as FORTRAN is derived from FORmula TRANslator, so capitalizing it does seem reasonable to me. Wikipedia lists both. I suppose I should bow to convention and use the more conventional spelling “Lisp”, if for no other reason that Guy Steele’s Common Lisp the Language (available at the link in HTML and PDF format, a useful reference) spells it that way. Or maybe I’ll just continue my backward ways…

Watermarking and titling with ffmpeg and other open source tools…

WordPress › Error

There has been a critical error on this website.

Learn more about troubleshooting WordPress.