## Version 3 of HashLife

Updated 2008-07-27 17:27:53 by lars_h

Lars H, 2008-08-27: I stumbled upon a Wikipedia article [L1 ] about this curious little algorithm, and once I had read enough to understand how it works, I couldn't help implementing it (especially since Wikipedia claims "While an amateur programmer might be able to write a simple Life player over an afternoon, few Hashlife implementations exist"). So here it is, as a weekend fun project.

Basic principles of the algorithm:

1. Use a quadtree to store patterns. This simplifies supporting unlimited grid sizes (just use a tree with enough levels).
2. Only ever create one node for each pattern (so the tree is really a DAG). This requires introducing the hash to map pattern represented by a node to the actual node index. The hash key is formed from the indices of the four subnodes.
3. In Life, cells are affected by all their neighbours, so a given n×n grid only contains enough information to compute the next state for the middle (n–2)×(n–2) cells. If n = 2**k (for a level k node in the quadtree), then the largest size of a quadtree node that fits within the inner square is 2**(k–1)× 2**(k–1) — one level further down. Hence this is the basic mode of operation for stepping in HashLife: given a tree node, compute the one step smaller node for the next state of the middle part of the given grid.
4. Use memoizing to speed things up when the same node occurs in several places — cache (in the node itself) the next value for a node once you've computed it.
5. Make the step size for a node proportional to the side of the node! In the code, this amounts to making next::bignode advance its components two steps instead of one, but the result is an exponential speedup in the simulated time.

## Short demo

```% encode {{1 1 1} {1 0 1} {1 0 1} {} {} {} {1 1} {1 1}}
10
% decode 10
{1 1 1 0 0 0 0 0} {1 0 1 0 0 0 0 0} {1 0 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 0} {1 1 0 0 0 0 0 0} {1 1 0 0 0 0 0 0}
% chardecode 10
OOO
O O
O O

OO
OO
% grow 10
14
% chardecode 14

OOO
O O
O O

OO
OO

% step 14
137
% chardecode 137

OOO
O   O
O     O
O     O
OOO OOO

OO
OO

% set generations
8```

Since node 14 is a 16×16 grid, stepping it forward advances time 8 generations. growing larger nodes also speeds up the processing.

```% step [grow 137]
743
% step [grow 743]
2769
% set generations
56```

Larger nodes are difficult to display in character graphics. '''node2xbm" can be used to make a bitmap instead, like so:

``` % package require Tk
% image create bitmap -background yellow
image1
% pack [label .l -image image1]
% image1 configure -data [node2xbm 2769 3] ; # 3 for making a 3x3 block of pixels for each cell.```

The 1496th generation of that pattern, displayed as above and with one pixel per cell looks like:

## The code

High level commands defined below:

encode
Compute node (number) for a matrix of 0s and 1s.
decode
Given node (number), return corresponding matrix.
chardecode
Given node (number), return corresponding block of chars.
grow
Double the side of a node, by padding with zeroes. Return new node number.
side
Given node (number), return its side length.
step
Given node (number), advance it side/2 generations, return new node number. (Combination of grow and next, so new node has the same size as the old one.)
step1
Given node (number), advance it 1 generation, return new node number. (Combination of grow and next1, so new node has the same size as the old one.)
node2xbm
Given node and magnification factor (number of pixels in side of one cell), return XBM data for a corresponding bitmap.
```# HashLife -- Hashed quadtree-based Life algorithm with memoisation

# The basic data structure is a heap (implemented as a list) and a hash
# table (array) that stores the positions in the heap of particular nodes.

# Heap item contents:
#
# cell <0 or 1>
# 2node <nw> <ne> <sw> <se>
# 4node <nw> <ne> <sw> <se> <self>
# bignode <nw> <ne> <sw> <se> <self>
# cached <nw> <ne> <sw> <se> <next>
# cached1 <nw> <ne> <sw> <se> <next>
#
# Evaluating one of them in the next namespace computes the <next> node
# and mutates the heap item to a [cached] node. Evaluating one of them
# in the next1 namespace computes the <next> node and mutates the heap
# item to a [cached1] node.
#
# The <nw>, <ne>, <sw>, <se>, and <next> are positions of other nodes.
# The <self> is the position of the node itself, which is needed for
# nodes that mutate themselves into [cached] nodes.
# In a [2node], they're [cell]s. In a [4node], they're [2node]s.
# In a [bignode], they're [4node]s, [bignode]s, or [cached]s.
# In a [cached], they're [2node]s or [cached]s.
#
#

set heap {{cell 0} {cell 1}}
# Must be items 0 and 1 -- relied upon by several procs below.
array unset hash

##
# The following procedure finds the heap position (pointer) for a
# particular node, given its four corners. It will *create* the heap
# item for the node if it doesn't already exist.
##

proc find {nw ne sw se} {
variable hash
set key [list \$nw \$ne \$sw \$se]
if {![info exists hash(\$key)]} then {
variable heap
set hash(\$key) [llength \$heap]
switch -- [lindex \$heap \$nw 0] "cell" {
lappend heap [list 2node \$nw \$ne \$sw \$se]
} "2node" {
lappend heap [list 4node \$nw \$ne \$sw \$se \$hash(\$key)]
} default {
lappend heap [list bignode \$nw \$ne \$sw \$se \$hash(\$key)]
}
}
return \$hash(\$key)
}

##
# side \$node -- Compute and return the side of the \$node.
##

proc side {node} {expr {
[lindex \$::heap \$node 0] eq "cell" ? 1 :
2*[side [lindex \$::heap \$node 1]]
}}

##
# encode \$matrix
#    Encode a matrix of booleans as a node in the heap.
#    The matrix is roughly centered within the node's grid.
#
# encode \$matrix \$side \$col \$row
#    Encode the \$side-by-\$side submatrix of \$matrix whose
#    NW corner is column \$col and row \$row. This submatrix
#    need not be contained within the \$matrix. The \$side
#    must be a power of 2.
##

proc encode {matrix {side 0} {x 0} {y 0}} {
if {!\$side} then {
set min [llength \$matrix]
foreach row \$matrix {
if {[llength \$row]>\$min} then {set min [llength \$row]}
}
set side 1
while {\$side<\$min} {incr side \$side}
set x [expr {(\$min-\$side)/2}]
set y [expr {([llength \$matrix]-\$side)/2}]
}
if {\$side==1} then {
return [lindex {0 1} [string is true -strict [lindex \$matrix \$y \$x]]]
}
set half [expr {\$side/2}]
set nw [encode \$matrix \$half \$x \$y]
set ne [encode \$matrix \$half [expr {\$x+\$half}] \$y]
set sw [encode \$matrix \$half \$x [expr {\$y+\$half}]]
set se [encode \$matrix \$half [expr {\$x+\$half}] [expr {\$y+\$half}]]
return [find \$nw \$ne \$sw \$se]
}

##
# Decode a heap node as a matrix of 0s and 1s.
##

proc decode {node} {
if {[lindex \$::heap \$node 0] == "cell"} then {
return [list [list [lindex \$::heap \$node 1]]]
}
set res {}
foreach low [decode [lindex \$::heap \$node 1]]\
high [decode [lindex \$::heap \$node 2]] {
lappend res [concat \$low \$high]
}
foreach low [decode [lindex \$::heap \$node 3]]\
high [decode [lindex \$::heap \$node 4]] {
lappend res [concat \$low \$high]
}
return \$res
}

##
# Decode a heap node as a character block. Dead cells
# are spaces, live cells some char (defaults to O).
##

proc chardecode {node {cell O}} {
string map [list { } {} 0 { } 1 \$cell] [
join [decode \$node] \n
]
}

##
# Expand a heap node \$level times, into a 2**\$level by 2**\$level
# matrix of nodes.
##

proc expand {node level} {
if {\$level<=0} then {return [list [list \$node]]}
incr level -1
set res {}
foreach low [expand [lindex \$::heap \$node 1] \$level]\
high [expand [lindex \$::heap \$node 2] \$level] {
lappend res [concat \$low \$high]
}
foreach low [expand [lindex \$::heap \$node 3] \$level]\
high [expand [lindex \$::heap \$node 4] \$level] {
lappend res [concat \$low \$high]
}
return \$res
}

##
# Contract the \$side-by-\$side submatrix, whose NW corner is at
# (\$row,\$col), of the \$matrix of nodes into a single node.
##

proc contract {matrix side row col} {
if {\$side<=1} then {return [lindex \$matrix \$row \$col]}
set half [expr {\$side/2}]
set nw [contract \$matrix \$half \$row \$col]
set ne [contract \$matrix \$half \$row [expr {\$col+\$half}]]
set sw [contract \$matrix \$half [expr {\$row+\$half}] \$col]
set se [contract \$matrix \$half [expr {\$row+\$half}] [expr {\$col+\$half}]]
return [find \$nw \$ne \$sw \$se]
}

namespace eval next {
proc cached {nw ne sw se next} {return \$next}

proc cell {c} {error "Can't advance a single cell"}
proc 2node {nw ne sw se} {error "Can't advance a 2x2 node"}

proc cached1 {nw ne sw se next} {
set self [find \$nw \$ne \$sw \$se]
if {[lindex \$::heap \$nw 0] == "2node"} then {
4node \$nw \$ne \$sw \$se \$self
} else {
bignode \$nw \$ne \$sw \$se \$self
}
}
}

##
# The following procedure computes the next generation of the middle
# [2node] in a [4node]. It does not mutate the node itself.
##

proc 4node {nw ne sw se} {
set grid {}
foreach low [decode \$nw] high [decode \$ne] {
lappend grid [concat \$low \$high]
}
foreach low [decode \$sw] high [decode \$se] {
lappend grid [concat \$low \$high]
}
set call [list find]
foreach y1 {1 2} {
foreach x1 {1 2} {
set sum 0
foreach y2 {-1 0 1} {
foreach x2 {-1 0 1} {
incr sum [lindex \$grid [expr {\$y1+\$y2}] [expr {\$x1+\$x2}]]
}
}
set sum [expr {2*\$sum - [lindex \$grid \$y1 \$x1]}]
lappend call [expr {\$sum>=5 && \$sum<=7}]
}
}
eval \$call
}

##
# The following procedure computes the next generation of the middle
# [2node] in a [4node], and mutates the node into a [cached] node.
##

proc next::4node {nw ne sw se self} {
set next [::4node \$nw \$ne \$sw \$se]
lset ::heap \$self 0 cached
lset ::heap \$self 5 \$next
return \$next
}

proc next::bignode {nw ne sw se self} {
set grid {}
foreach low [expand \$nw 1] high [expand \$ne 1] {
lappend grid [concat \$low \$high]
}
foreach low [expand \$sw 1] high [expand \$se 1] {
lappend grid [concat \$low \$high]
}
set n   [contract \$grid 2 0 1]
set w   [contract \$grid 2 1 0]
set mid [contract \$grid 2 1 1]
set e   [contract \$grid 2 1 2]
set s   [contract \$grid 2 2 1]
set grid2 {}
foreach row [list [list \$nw \$n \$ne] [list \$w \$mid \$e] [list \$sw \$s \$se]] {
set row2 {}
foreach node \$row {
lappend row2 [eval [lindex \$::heap \$node]]
}
lappend grid2 \$row2
}
set call [list find]
foreach row {0 1} {
foreach col {0 1} {
set node [contract \$grid2 2 \$row \$col]
lappend call [eval [lindex \$::heap \$node]]
}
}
set next [eval \$call]
lset ::heap \$self 0 cached
lset ::heap \$self 5 \$next
return \$next
}

namespace eval next1 {
proc cached1 {nw ne sw se next} {return \$next}

proc cell {c} {error "Can't advance a single cell"}
proc 2node {nw ne sw se} {error "Can't advance a 2x2 node"}

proc cached {nw ne sw se next} {
set self [find \$nw \$ne \$sw \$se]
if {[lindex \$::heap \$nw 0] == "2node"} then {
lset ::heap \$self 0 cached1
set next
} else {
bignode \$nw \$ne \$sw \$se \$self
}
}
}

##
# The following procedure computes the next generation of the middle
# [2node] in a [4node], and mutates the node into a [cached1] node.
##

proc next1::4node {nw ne sw se self} {
set next [::4node \$nw \$ne \$sw \$se]
lset ::heap \$self 0 cached1
lset ::heap \$self 5 \$next
return \$next
}

proc next1::bignode {nw ne sw se self} {
set M [expand \$self 3]
set call [list find]
foreach {row col} {1 1  1 3   3 1  3 3} {
set node [contract \$M 4 \$row \$col]
lappend call [eval [lindex \$::heap \$node]]
}
set next [eval \$call]
lset ::heap \$self 0 cached1
lset ::heap \$self 5 \$next
return \$next
}

##
# Commands in the zeroed namespace return a node with the same size as
# the node which they are, but with all cells 0.
##

namespace eval zeroed {
proc cell {c} {return 0}

proc 2node {a b c d} {
find 0 0 0 0
}

proc cached {a b c d next} {
set sub [eval [lindex \$::heap \$a]]
find \$sub \$sub \$sub \$sub
}

}
interp alias {} zeroed::bignode {} zeroed::cached
interp alias {} zeroed::4node {} zeroed::cached
interp alias {} zeroed::cached1 {} zeroed::cached

##
# The [grow] procedure return a node whose side is twice
# that of the argument, by padding it with a frame of all
# zeroes.
##

proc grow {node} {
if {[lindex \$::heap \$node 0] == "cell"} then {
error "You can't grow a single cell"
}
foreach {type nw ne sw se} [lindex \$::heap \$node] break
set clean [namespace eval ::zeroed [lindex \$::heap \$nw]]
set call [list find]
lappend call [find \$clean \$clean \$clean \$nw]
lappend call [find \$clean \$clean \$ne \$clean]
lappend call [find \$clean \$sw \$clean \$clean]
lappend call [find \$se \$clean \$clean \$clean]
eval \$call
}

set generations 0

proc step {node} {
variable generations
set side [side \$node]
set node2 [grow \$node]
set res [namespace eval ::next [lindex \$::heap \$node2]]
incr generations [expr {\$side/2}]
return \$res
}

proc step1 {node} {
variable generations
set node2 [grow \$node]
set res [namespace eval ::next1 [lindex \$::heap \$node2]]
incr generations
return \$res
}

# A [bitmap] actually seems more appropriate than a [photo] for
# representing Life patterns, so we need some way of converting
# nodes to XBM data. Intermediate representation:
#
#   \$side \$matrix
#
# If the \$side is divisible by 8, then the \$matrix is a list of
# \$side rows, where each row is a string matching the regexp:
#   (0x[[:xdigit:]]{2}, )+
#
# If the \$side is not divisible by 8, then the \$matrix is a list
# of \$side rows, where each row is a string of \$side binary digits.

proc node2bm {node mag} {
if {\$node <= 1} then {
set row [string repeat \$node \$mag]
if {!(\$mag % 8)} then {
binary scan [binary format b* \$row] H* row
regsub -all {..} \$row {0x&, } row
}
set mat {}
for {set n 1} {\$n<=\$mag} {incr n} {lappend mat \$row}
return [list \$mag \$mat]
}
set mat {}
foreach {side leftM} [node2bm [lindex \$::heap \$node 1] \$mag] break
foreach l \$leftM r [
lindex [node2bm [lindex \$::heap \$node 2] \$mag] 1
] {lappend mat \$l\$r}
foreach l [
lindex [node2bm [lindex \$::heap \$node 3] \$mag] 1
] r [
lindex [node2bm [lindex \$::heap \$node 4] \$mag] 1
] {lappend mat \$l\$r}
if {(\$side%8) && !(\$side%4)} then {
set mat2 {}
foreach row \$mat {
binary scan [binary format b* \$row] H* row
regsub -all {..} \$row {0x&, } row
lappend mat2 \$row
}
set mat \$mat2
}
return [list [expr {2*\$side}] \$mat]
}

proc node2xbm {node mag} {
foreach {side mat} [node2bm \$node \$mag] break
if {\$side % 8} then {
set width [expr {8*(1 + \$side/8)}]
set mat2 {}
foreach row \$mat {
binary scan [binary format b* \$row] H* row
regsub -all {..} \$row {0x&, } row
lappend mat2 \$row
}
set mat \$mat2
} else {
set width \$side
}
format {
#define data_width %d
#define data_height %d
static unsigned char data_bits[] = {
%s
};
} \$width \$side [string trimright [join \$mat "\n    "] ", "]
}```