Programmer’s Challenge

By Bob Boonstra, Westford, Massachusetts

Postman’s Sort

Everyone knows that sorting algorithms require execution time that grows faster than linearly with the number of input records. Specifically, sorting N input records is known to require somewhere between O(Nlog2N) and O(N2) comparisons, depending on the specific algorithm and whether one is considering average or worst-case time, right? Wrong. Your Challenge is to code a sorting algorithm that requires time proportional to the number of input records.

The bounds in the previous paragraph apply to sorting algorithms that are based on pairwise comparison of input keys. A distribution sort, however, requires no key comparisons, and has different performance characteristics. Imagine, for example, how the post office might sort mail. Initially, the mail might be distributed into piles by country. Each of those piles might be distributed into other piles by state. Those piles might be distributed by city, and then by street, and finally by address. All of this could be accomplished by making a sequence of passes through the input without ever comparing the sort keys of two input records to one another. This month, you are to implement a distribution (or “postman’s”) sort algorithm. The prototype for the code you should write is:

#include <stdio.h>
static pascal OSErr PostmansSort(
 FILE *inFile,   /* stream containing sort input */
 FILE *outFile,  /* stream to contain sort output */
 void *privateStorage,    /* preallocated storage for your use */
 size_t storageSize/* number of bytes of privateStorage */
);

The input will consist of records in the following format:

[FirstName] <tab> [MiddleName] <tab> LastName <tab> 
StreetAddress <tab> StreetName <tab> 
City <tab> State <tab> [ZipCode] <tab> [Country] <C>R

Input records should be read from inFile, sorted, and written to outFile. Both the input and output files will be opened and closed for you by the calling routine. Records will consist of up to 8 fields, as indicated above, and be terminated by a carriage return (0x0D). Fields will be separated by tabs (0x09). The square brackets in the format description are meta-characters indicating optional fields; the brackets will not be present in the input. Fields may contain embedded spaces or special characters (except tabs and returns).

Records are to be sorted into ascending order with Country as the primary sort key, ZipCode (if present) as the secondary key, then StreetName, StreetAddress, LastName, FirstName, and finally MiddleName, in that order. If ZipCode is not present, then State and City should replace it as the secondary and tertiary sort keys, respectively; otherwise State and City should be ignored. Sort order should be lexicographic except when a field value is purely numeric, and a purely numeric field should be treated as smaller than a field containing non-numeric characters. That is, field values “1”, “2”, “10”, “1B”, “1C”, and “2A” would sort in the order indicated. Empty optional fields should be considered to have the smallest possible value (e.g., records with no Country should be output before records with a Country value). Your routine should return zero (noErr) if it completes normally, or a non-zero error code if it is unable to complete the sort for any reason.

There are no predetermined limits on the number of distinct countries, State, City, etc. values that might be present in the input. Your routine will be provided with storageSize bytes (at least 64KB) of preallocated storage, initialized to zero by the calling routine, and pointed to by privateStorage. You may not allocate any additional memory, although use of small amounts of static memory is permitted. If your solution requires additional storage, you should create and write to temporary disk files. Adequate disk space is guaranteed. In approximately half of the test cases, you will be guaranteed at least 32 bytes of memory per record, plus 32 bytes for each unique field value. Scoring will weight these high-memory cases equally with the cases where less memory is provided.

This will be a native PowerPC Challenge, scored using the Metrowerks environment. Solutions may be coded in C, C++, or Pascal. If you use a language other than C, you should ensure that your code can be called from a test driver written in C. Most importantly, to be considered correct, your code must sort the input into the proper order without performing any key comparisons between input records. The fastest correct solution will be the winner.

This Challenge is based on a suggestion from Peter Shank, passed on to me by Mike Scanlin. Peter forwarded a reference to an article on the subject by Robert Ramsey, which you can find at http://www.silcom.com/~ramey/article.html.

Two Months Ago Winner

Congratulations to Ernst Munter (Kanata, ON, Canada) for submitting the fastest entry to the Mutant Life Challenge. The April Challenge was to write code that would propagate a population of cellular automata for a specified number of generations according to rules generalized from John Conway’s Game of Life.

Of the 17 entries I received, 9 worked completely correctly. A number of people submitted solutions that were partially, but not completely, correct. Some forgot to stop processing and return when the population became stable. Others had problems dealing with the BitMap as a torus, which required solutions to wrap from the top row to the bottom, from the right column to the left, and vice versa.

Algorithm and data structure selection were key to doing well in this Challenge. Several participants recognized the value of maintaining, for each cell, a count of the number of occupied neighboring cells, and updating the neighboring counts when the state of a cell changed. Ernst’s solution takes this idea further, and maintains neighbor counts for a block of 32 cells together. He uses some interesting modulo-9 math to take advantage of the fact that neighbor counts are between 0 and 8 to store 32 neighbor counts efficiently in 16 bytes. The winning solution creates a lookup table, based on the birthRules and deathRules provided as parameters, to efficiently determine whether a cell birth/death ought to occur. Other details of Ernst’s solution strategy are described in his well-commented code. It is an interesting and efficient solution that takes advantage of the storage provided by the problem statement.

I used a number of test cases to determine the overall score for an entry. The parameters for those test cases are listed below. The raw time for each test case was normalized to a common number of propagation generations (1000), and then the normalized scores were summed to create the final score. The test case parameters are as follows:

The table below summarizes the results for entries that worked correctly or partially correctly. It shows the overall score, code size, and data size for each of the solutions, as well as the raw score for selected test cases from those listed above. Asterisks indicate entries that executed to completion but were not completely correct; these entries were disqualified as usual, but I’ve listed their times at the bottom of the table of results. Numbers in parentheses after a person’s name indicate that person’s cumulative point total for all previous Challenges, not including this one.

Name lang score code data Case 1 Case 3 Case 5 Case 7

Ernst Munter (134) C 93 4944 >1M 63 194116 15736 1218

Bill Karsh (80) C 132 2804 8 540 192709 30518 4541

Randy Boring C 135 30248 2060 215 100431 30327 7782

Andy Antoniewicz C 204 1344 1122 843 415767 62759 2285

John Sweeney (4) C 362 19884 4554 627 304878 61166 19937

Ludovic Nicolle (14) C 479 7200 40 767 277098 84103 32237

Xan Gregg (92) C 790 2364 8 1205 640475 139368 47749

Elden Wood C++ 2036 5152 16 4905 223790 493359 4349

Wolfgang Thaller (4) C++ 2225 5088 44 3944 1287532 360201 152793

*John Nevard (17) C 264 4528 117 289 587092 87061 331

*Steve Israelson (20) C++ 854 2536 414 1381 614717 156168 53592

*Moss Prescott C++ 855 2956 1233 1721 905749 327008 15574

*Thomas Studer C++ 881 3448 56 996 1231139 146850 33715

*Jeff Mallett (44) C 930 1808 99032 1491 635424 73610 63609

*Tom Saxton (10) C 1035 2124 131478 836 1171235 192668 53138

Special thanks this month to Ludovic Nicolle for investing time to create some test code for this problem and for distributing it to the Challenge mailing list. Thanks also to David Cary for sending me several articles by Michael Abrash on optimizing implementations of the Life. And finally, Steve Israelson earns 2 points for having been the first to suggest a Challenge based on the Game of Life.

Those of you with an interest in exploring Life further might want to start with the page maintained by Paul Callahan at http://www.cs.jhu.edu/~callahan/lifepage.html. One of the links on that page points to the shareware program LifeLab, by Andrew Trevorrow, which supports exploration of generalized Life worlds and provides a set of initial patterns that lead to interesting results.

Top 20 Contestants Of All Time

Here are the top contestants for the Programmer’s Challenges to date. The point totals below include points awarded for this month’s entries.

There are three ways to earn points: (1) scoring in the top 5 of any Challenge, (2) being the first person to find a bug in a published winning solution or, (3) being the first person to suggest a Challenge that I use. The points you can win are:

1st place 20 points 5th place 2 points

2nd place 10 points finding bug 2 points

3rd place 7 points suggesting Challenge 2 points

4th place 4 points other “good deeds” 2 points

Here is Ernst Munter’s winning solution:

MutantLife.c

Copyright © 1996, Ernst Munter

/*
  The challenge is to write a fast program that computes the Life-like world some 
  generations into the future,  Intermediate generations need not necessarily be 
  computed, but to me at least, it seems we must process one generation at a time.

  The birth and death rules are specified as bit maps which determine changes as a 
  function of any combination of the number of cell neighbors (0 to 8), giving rise to 
  2**18 possible rule sets.

  Another requirement is that the propagation function stop as soon as two
  consecutive worlds are identical.

  Additional memory is limited to 1MB plus up to 10 times the size of the original 
  bitmap.

  Solution Strategy
  ---------
  We copy the starting bitmap to an internal representation; process this for the 
  required number of generations, or until we see no change; and then copy the final 
  state back to the original bitmap location.

  The internal cell format includes the state as well as the neighbor count of each cell.

  The propagation process consists of two phases for each generation: phase 1 in 
  which all changes are used to update the neighbor counts, and phase 2 in which the
  neighbor counts are re-evaluated, together with the current state, to determine the 
  change to the next state.

  Since usually, not all of the cell map changes, it is only necessary to process those 
  cells, and their immediate neighbors, that do change.

  Data Structure
  -------
  Because of memory limitation, and in the interest of speed blocks of 32 cells are 
  stored and processed together.

  A “CellBlock” of 32 cells is a 32-byte structure which contains cell states, cell 
  changes, counters, and pointers to allow cell blocks to be chained into linked lists.

  A number of flag bits are stored with each block to allow marking of border blocks 
  which require special computation to determine the neighbor cells (wrapping).

  The counts are stored modulo-9 arithmetic, A 13 bit counter holds the neighbor 
  counts for 4 cells, and there are 8 such counters located in each cell block.

  CellBlocks are allocated dynamically.

  Computation and Tables
  -----------
  We must perform 2 kinds of computation:

  Count Updates
  -------
  To update the neighbor counts we determine the addresses of the neighbors and 
  increment or decrement each modulo-9 counter depending on the change of state to 
  0 or 1

  The counter update function is fixed (not dependent on birth or death rules).  The 
  changes in a row of 10 cells affect exactly 8 counts in the same row and the blocks
  above and below.  Hence a lookup table with 1024 entries can be used to provide the 
  pre-computed sums of power-9 products with 1 or 0, for two of the 13-bit counters
  which are conveniently stored as a pair in a single word.

  State update and its rules
  -------------
  To determine the state change of a cell we need consider only the state of the cell, 
  the neighbor count, and the current set of birth/death rules.

  Given that 4 counts are stored, munged together in a single 13-bit modulo-9 counter, 
  we can also use a look-up table for this function.  The rules table contains 9*9*9*9
  entries of 1 byte each, where each byte contains a 4-bit birth field, and a 4-bit death 
  field, simply derived from the bits in the given birth/deathRules selected according
  to the four neighbor counts implied by the table address.

  This rules table must be computed when PropagateLife is called, to reflect the 
  current rules.

  Assuming that PropagateLife might be called frequently in succession, often with the 
  same rules, we can avoid re-computing this table unnecessarily.

  And since we can use a “reasonable” amount of static memory (up to 1MB, say), we 
  can store a fairly big rule book of many different rules.  This would allow alternating 
  rules to be handled efficiently as well.

  Hence, all left over available (permitted) memory is used for the rule book.  The 18-
  bits of concatenated birth and death rules are hashed into a page index to quickly 
  locate a previously computed rules page.

  Other Assumptions
  ---------
  cells is a bitmap where
    cells.bounds.top=cells.bounds.left==0
    cells.bounds.bottom - cells.bounds.top >= 3
    cells.bounds.right - cells.bounds.left >= 3

  The maximum static memory limit can be #defined, minimum 15K.

  The function allocates dynamic memory equivalent to 8 more bitmaps + 32 bytes. 
  It returns 0 if malloc fails.

  “Spontaneous birth” is allowed (birthRules LSB set).
*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
pascal long PropagateLife(
     BitMap cells,
     long numGenerations,
     short birthRules,
     short deathRules);

// Define amount of static memory to be available, minimum 15K:
// the value of 1MB permits 158 rule pages to be stored

#define availableStatic 0x100000L

// You can reduce the size of the program by about 20% by setting 
// LARGE_BITMAPS to 0.
// But runtime for large worlds (1000x1000) goes up 3-4%.  Normally, we would go 
// with the shorter program, but in this challenge, 3% might make a lot of difference:

#define LARGE_BITMAPS 1

Type Definitions

#define ulong  unsigned long
#define ushort unsigned short
#define uchar  unsigned char

struct       CellBlock {
  ulong      state;
  ulong      change;
  struct CellBlock* link1;
  struct CellBlock* link2;
  long       count[4];
};
typedef struct CellBlock CellBlock;

struct       RulesPage {
  int        birthRule;
  int        deathRule;
  uchar      rules[81*81];
  char       pad[3];
};
typedef struct RulesPage RulesPage;

struct       Bits2Counts {
  long       UD;
  long       LR;
};
typedef struct Bits2Counts Bits2Counts;

Function Prototypes

uchar*     MakeRuleTable(int b, int d);

CellBlock* LoadCells(BitMap cells, int width, int padWords,
           CellBlock* cb, CellBlock* L1);

CellBlock* LoadCellsSpontaneous(BitMap cells,int width,
           int padWords, CellBlock* cb, CellBlock* L1);

CellBlock* UpdateCounts(CellBlock* cb, CellBlock* L2,
           int width, int size, int lastColumn, long rightMask);

#if LARGE_BITMAPS
CellBlock* UpdateCountsNoWrap(CellBlock* cb, CellBlock* L2,
           int width);
#endif

ulong      ComputeStateChange(CellBlock* cb, uchar* rule);

void       UnloadCells(CellBlock* cb, int size, void* baseAddr);

void       UnloadPaddedCells(BitMap cells, int width,
           int padWords, CellBlock* cb);

Static Data

#define bcUD(i) (                                       \
 ((((i)>>5)&1)*729L)+((((i)>>4)&1)*810L)+               \
 ((((i)>>3)&1)*819L)+((((i)>>2)&1)*91L)+                \
 ((((i)>>1)&1)*10L)+ ((i)&1)                )

#define bcLR(i) (                                       \
 ((((i)>>5)&1)*729L)+((((i)>>4)&1)*81L)+                \
 ((((i)>>3)&1)*738L)+((((i)>>2)&1)*82L)+                \
 ((((i)>>1)&1)*9L)+  ((i)&1)                )

#define bc(i) {                                         \
   (bcUD((i)>>4)<<13)+bcUD(i),(bcLR((i)>>4)<<13)+bcLR(i)}

#define bc4(i)                                          \
        bc(i),bc(i+1),bc(i+2),bc(i+3)
#define bc16(i)                                         \
        bc4(i),bc4(i+0x4),bc4(i+0x8),bc4(i+0xC)
#define bc64(i)                                         \
        bc16(i),bc16(i+0x10),bc16(i+0x20),bc16(i+0x30)
#define bc256(i)                                        \
        bc64(i),bc64(i+0x40),bc64(i+0x80),bc64(i+0xC0)
//
// Ernst was able to compile these macros on his PC, but CodeWarrior complained 
// they were too complicated.  I got it to compile by expanding each BC256() macro 
// invokation into 256 individual bc() statements.  -Bob
//
static Bits2Counts BC[0x400] = {
  bc256(0),bc256(0x100),bc256(0x200),bc256(0x300)
};

#define numRulesPages                                   \
      ((availableStatic-sizeof(BC)-16)/sizeof(RulesPage))

static RulesPage rulesCache[numRulesPages];

Constants

#define MSB         0x80000000L
#define LSB         0x00000001L

#define UP          0x40000000L
#define DOWN        0x20000000L
#define LEFT        0x10000000L
#define RIGHT       0x08000000L

#define IS_U_BORDER(cb) (cb->count[0] & UP)
#define IS_D_BORDER(cb) (cb->count[0] & DOWN)
#define IS_L_BORDER(cb) (cb->count[0] & LEFT)
#define IS_R_BORDER(cb) (cb->count[0] & RIGHT)

PropagateLife

pascal long PropagateLife(
  BitMap cells,
  long   numGenerations,
  short  birthRules,
  short  deathRules) {

  int        width=(cells.bounds.right+31)>>5;
  int        lastColumn=(cells.bounds.right-1) & 31;
  long       rightMask=0xFFFFFFFF << (31-lastColumn);
  int        size=width*cells.bounds.bottom;
  int        padWords=(cells.rowBytes>>2)-width;
  long       gen;
  long*      allocMem;
  CellBlock* cb;
  CellBlock* L1;
  CellBlock* L2;
  uchar*     rule;

// Get memory for cell blocks, 1 block per 32 cells:

 if (0==(allocMem=(long*)malloc(32+size*sizeof(CellBlock))))
        return 0;

// Align cell blocks with cache-line boundary (32-byte):

  { char* temp=(char*)allocMem;
    temp+=(32-(ulong)temp) & 31;
    cb=(CellBlock*)temp;
  }

// Clear all cell blocks:

  memset(cb,0,size*sizeof(CellBlock));

// Establish the rules look-up table in static memory:

  rule=MakeRuleTable(birthRules,deathRules);

// Copy BitMap cells into CellBlock structures and link:
// Usually, only non-empty cells are linked, but
// if birthRules include a “spontaneous birth” bit,
// all cells must be linked.

  if (birthRules & 1)
    L1=LoadCellsSpontaneous(cells,width,padWords,cb,cb-1);
  else L1=LoadCells(cells,width,padWords,cb,cb-1);
  L2=L1;

  for (gen=1;;gen++) {

// Linked list 1 contains cell blocks with a state change. The state changes are applied 
// to the cell states.  Loop 1 traverses list 1 and updates the counts of all affected cells 
// in neighboring blocks and itself.  If any counts change, the block is linked into list 
// 2, unless it is already there.

    while (L1>=cb) {
#if LARGE_BITMAPS
      if ((L1->count[0] & (UP | DOWN | LEFT | RIGHT)) == 0)
        L2=UpdateCountsNoWrap(L1,L2,width);
      else
#endif
      L2=UpdateCounts(L1,L2,width,size,lastColumn,rightMask);
      L1=L1->link1;
    }

// List 1 is now empty.  Loop 2 finds the cell changes in each block of list 2, using the // current states and 
the neighbor counts.  If any state changes, the block is put back 
// into list 1.

    while (L2>=cb) {
      CellBlock* next=L2->link2;
      if (ComputeStateChange(L2,rule)) {
        L2->link1=L1;
        L1=L2;
      }
      L2->link2=0;
      L2=next;
    }

// If list 1 is empty (no change) we can make a fast exit:
    if (L1<cb) {
      if (gen==1) goto quick_exit; // nothing changed!
      break;
    }
    if (gen>=numGenerations) break;
  }

// Apply last changes to cells in list 1:
  while (L1>=cb) {
    L1->state ^= L1->change;
    L1=L1->link1;
  }

// The final states of all cells in the private cell blocks are copied back into the 
// external BitMap storage.  A slightly slower function provides for the uncommon 
// case when there are extra unused words beyond the right bound.

  if (padWords==0) UnloadCells(cb,size,cells.baseAddr);
  else UnloadPaddedCells(cells,width,padWords,cb);

// Return all allocated dynamic memory to the system:
quick_exit:
  free(allocMem);

  return gen;
}

Macros
// SETBLOCK is used in LoadCells(..)
#define SETBLOCK(flag)                                  \
{ ulong c32=*bm++;                                      \
  cb->count[0]=flag;                                    \
  if (c32) {                                            \
    cb->change=c32;                                     \
    cb->link1=L1;                                       \
    cb->link2=L1;                                       \
    L1=cb;                                              \
  }                                                     \
  cb++;                                                 \
}

// SETBLOCKs is used in LoadCellsSpontaneous(..)
#define SETBLOCKs(flag)                                 \
{ ulong c32=*bm++;                                      \
  cb->count[0]=flag;                                    \
  cb->change=c32;                                       \
  cb->link1=L1;                                         \
  cb->link2=L1;                                         \
  L1=cb;                                                \
  cb++;                                                 \
}

// LINK is used in UpdateCounts(..)
#define LINK(cell)                                      \
{ if (((cell)->link2)==0) {                             \
    (cell)->link2=L2;                                   \
    L2=cell;                                            \
  }                                                     \
}

// CHANGE_COUNT is used in UpdateCounts(..)
#define CHANGE_COUNT(cell,ctr,delta);                   \
{ long cnt=(cell)->count[ctr] | MSB;                    \
  (cell)->count[ctr]=cnt+(delta);                       \
}

// The following three macros are used in UpdateCounts(..)
#define UPDATE_COMMON(ctr)                              \
{ delta=BC[newIndex].UD-BC[oldIndex].UD;                \
  CHANGE_COUNT(cUP,ctr,delta);                          \
  CHANGE_COUNT(cDN,ctr,delta);                          \
  delta=BC[newIndex].LR-BC[oldIndex].LR;                \
  CHANGE_COUNT(cb,ctr,delta);                           \
}

// C/C++ does not understand negative shifts, so we make separate macros for left and 
// right shifts;
#define UPDATE_UD(shft,ctr)                             \
{ int oldIndex=(oldState >> shft) & 0x3FF;              \
  int newIndex=(newState >> shft) & 0x3FF;              \
  UPDATE_COMMON(ctr)                                    \
}

#define UPDATE_UD_(shft,ctr)                            \
{ int oldIndex=(oldState << shft) & 0x3FF;              \
  int newIndex=(newState << shft) & 0x3FF;              \
  UPDATE_COMMON(ctr)                                    \
}

MakeRuleTable
uchar* MakeRuleTable(int b,int d) {

// Checks if a valid rules table exists for the current rules, and computes a new table if 
// necessary.

int pageNumber=((b << 9) + d) % numRulesPages;
RulesPage* rp=rulesCache+pageNumber;
  if ((rp->birthRule != b)
  ||  (rp->deathRule != d)) {

    int q,r,s,Q,R,S;
    int t0,t1,t2,t3,t4,t5,t6,t7,t8;
    unsigned char* t=rp->rules;
    rp->birthRule = b;
    rp->deathRule = d;
    b=b<<4;;
    t0=((d)&1)    | ((b)&0x10);
    t1=((d>>1)&1) | ((b>>1)&0x10);
    t2=((d>>2)&1) | ((b>>2)&0x10);
    t3=((d>>3)&1) | ((b>>3)&0x10);
    t4=((d>>4)&1) | ((b>>4)&0x10);
    t5=((d>>5)&1) | ((b>>5)&0x10);
    t6=((d>>6)&1) | ((b>>6)&0x10);
    t7=((d>>7)&1) | ((b>>7)&0x10);
    t8=((d>>8)&1) | ((b>>8)&0x10);
    t[0]=(uchar)t0;
    t[1]=(uchar)t1;
    t[2]=(uchar)t2;
    t[3]=(uchar)t3;
    t[4]=(uchar)t4;
    t[5]=(uchar)t5;
    t[6]=(uchar)t6;
    t[7]=(uchar)t7;
    t[8]=(uchar)t8;
    t+=6561-9;
    for (s=8;s>=0;s-) {
      S=rp->rules[s]<<3;
      for (r=8;r>=0;r-) {
        R=S | (rp->rules[r]<<2);
        for (q=8;q>=0;q-) {
          Q=R | (rp->rules[q]<<1);
          t[0]=(uchar)(Q | t0);
          t[1]=(uchar)(Q | t1);
          t[2]=(uchar)(Q | t2);
          t[3]=(uchar)(Q | t3);
          t[4]=(uchar)(Q | t4);
          t[5]=(uchar)(Q | t5);
          t[6]=(uchar)(Q | t6);
          t[7]=(uchar)(Q | t7);
          t[8]=(uchar)(Q | t8);
          t-=9;
        }
      }
    }
  }
  return rp->rules;
}

LoadCells
CellBlock* LoadCells(
           BitMap     cells,
           int        width,
           int        padWords,
           CellBlock* cb,
           CellBlock* L1) {

// Copies the states of all cells in the bitmap into the cell blocks, as state changes to 
// trigger evaluation.  Sets border flags as needed.  Links all non-empty cell blocks in 
// both lists.  The special case of “spontaneous birth” is handled by a separate function 
// LoadCellsSpontaneous(..)

ulong* bm=(ulong*)cells.baseAddr;
int    x,y;

  if (width>1) {

//top edge:
    SETBLOCK(UP+LEFT);
    for (x=1;x<width-1;x++) SETBLOCK(UP);
    SETBLOCK(UP+RIGHT);
    bm+=padWords;

//middle rows:
    for (y=1;y<cells.bounds.bottom-1;y++) {
      SETBLOCK(LEFT);
      for (x=1;x<width-1;x++) SETBLOCK(0);
      SETBLOCK(RIGHT);
      bm+=padWords;
    }

//bottom edge:
    SETBLOCK(DOWN+LEFT);
    for (x=1;x<width-1;x++) SETBLOCK(DOWN);
    SETBLOCK(DOWN+RIGHT);
  } else {

//case of bit map <= 32 cells wide
//top edge:
    SETBLOCK(UP+LEFT+RIGHT);
    bm+=padWords;
//middle rows:
    for (y=1;y<cells.bounds.bottom-1;y++) {
      SETBLOCK(LEFT+RIGHT);
      bm+=padWords;
    }
//bottom edge:
    SETBLOCK(DOWN+LEFT+RIGHT);
  }
  return L1;
}

UpdateCounts
CellBlock* UpdateCounts(
           CellBlock* cb,
           CellBlock* L2,
           int        width,
           int        size,
           int        lastColumn,
           long       rightMask) {

// Processes one cell block: locates its neighbor blocks, and updates their counters 
// according to the state changes of the cells in the center block.

// The 32 cells of a block are divided into 4 overlapping groups, corresponding to the 
// 4 counter words affected.  Only counters whose count changes are accessed, and so
// flagged (by setting MSB).

// Notably, right and left neighboring cell blocks are only affected if the first or last 
// bits of the center block change.

// The function deals with a number of special cases such as wrapping, and the 
// possibly partially filled block at the right margin of the bitmap.

ulong      oldState=cb->state;
ulong      theChange=cb->change;
ulong      newState;
ulong      lastBit;
CellBlock* cUP=cb-width;
CellBlock* cDN=cb+width;
long       delta;

  if (IS_U_BORDER(cb)) cUP+=size;
  if (IS_D_BORDER(cb)) cDN-=size;
  if (IS_R_BORDER(cb)) {
    theChange &= rightMask;
    lastBit=0x80000000>>lastColumn;
  } else lastBit=LSB;

  newState=oldState ^ theChange;
  cb->state=newState;
  if (theChange & 0xFF800000) {
    UPDATE_UD(23,0);
    if (theChange & MSB) {                      //carry left
      int ctr=3;
      int offset=-1;
      delta=((newState >>30) & 2) - 1;
      if (IS_L_BORDER(cb)) {                    //wrap
        ctr = lastColumn >> 3;
        offset+=width;
        switch (lastColumn & 7) {
        case 0:delta = (delta<<13)*729; break;
        case 1:delta = (delta<<13)*81; break;
        case 2:delta = (delta<<13)*9; break;
        case 3:delta = (delta<<13); break;
        case 4:delta *= 729; break;
        case 5:delta *= 81; break;
        case 6:delta *= 9; break;
        }
      }
      CHANGE_COUNT(cUP+offset,ctr,delta);
      LINK(cUP+offset);
      CHANGE_COUNT(cb+offset,ctr,delta);
      LINK(cb+offset);
      CHANGE_COUNT(cDN+offset,ctr,delta);
      LINK(cDN+offset);
    }
  }
  if (theChange & 0x01FF8000) UPDATE_UD(15,1);
  if (theChange & 0x0001FF80) UPDATE_UD( 7,2);
  if (theChange & 0x000001FF) UPDATE_UD_(1,3);
  if (theChange & lastBit) {                    //carry right
    int offset=+1;
    if (IS_R_BORDER(cb)) {                      //wrap
      delta = (((newState >> (31-lastColumn)) << 1) & 2) -1;
      offset -= width;
    } else delta = (((newState << 1) & 2) - 1);
    delta *= (729<<L13);
    CHANGE_COUNT(cUP+offset,0,delta);
    LINK(cUP+offset);
    CHANGE_COUNT(cb+offset,0,delta);
    LINK(cb+offset);
    CHANGE_COUNT(cDN+offset,0,delta);
    LINK(cDN+offset);
  }

  LINK(cUP);
  LINK(cDN);
  LINK(cb);
  return L2;
}

#if LARGE_BITMAPS
CellBlock* UpdateCountsNoWrap(
           CellBlock* cb,
           CellBlock* L2,
           int        width) {

// A stream-lined version of UpdateCounts(..)  for interior (non-border) cell blocks 
// which cannot wrap.

ulong      oldState=cb->state;
ulong      theChange=cb->change;
ulong      newState;
CellBlock* cUP=cb-width;
CellBlock* cDN=cb+width;
long       delta;

  newState=oldState ^ theChange;
  cb->state=newState;
  if (theChange & 0xFF800000) {
    UPDATE_UD(23,0);
    if (theChange & MSB) {                      //carry left
      int ctr=3;
      int offset=-1;
      delta=((newState >>30) & 2) - 1;
      CHANGE_COUNT(cUP+offset,ctr,delta);
      LINK(cUP+offset);
      CHANGE_COUNT(cb+offset,ctr,delta);
      LINK(cb+offset);
      CHANGE_COUNT(cDN+offset,ctr,delta);
      LINK(cDN+offset);
    }
  }
  if (theChange & 0x01FF8000) UPDATE_UD(15,1);
  if (theChange & 0x0001FF80) UPDATE_UD( 7,2);
  if (theChange & 0x000001FF) UPDATE_UD_(1,3);
  if (theChange & LSB) {                       //carry right
    int offset=+1;
    delta = (((newState << 1) & 2) - 1) * (729<<L13);
    CHANGE_COUNT(cUP+offset,0,delta);
    LINK(cUP+offset);
    CHANGE_COUNT(cb+offset,0,delta);
    LINK(cb+offset);
    CHANGE_COUNT(cDN+offset,0,delta);
    LINK(cDN+offset);
  }

  LINK(cUP);
  LINK(cDN);
  LINK(cb);
  return L2;
}
#endif

ComputeStateChange

ulong ComputeStateChange(CellBlock* cb,uchar* rule) {

// Computes the state change for 32 cells in a cell block.

// The 4 counters are flagged, and only those counters that actually changed need to 
// be considered.

// The algorithm computes 32-bit birth and death masks from the counts: Each counter 
// word cb->count[i] contains two 13-bit counter values.  These are used successively 
// to index into the rules table each time extracting two 4-bit fields, for births and 
// deaths according to the 4 modulo-9 counts inherent in the counters.  The 4-bit fields 
// from all lookups (or 0) are separately stacked into 32-bit birth and death masks.

// Note: counters that are known not to have changed
//       (MSB clear) cannot effect a state change and are
//       skipped in the evaluation of counts.

// In each bit position then, each state bit selects the corresponding birth or death 
// mask as follows:  A 0-state selects the birth bit, and a 1-state selects the death bit.
/*
   oldState  birth  death  change
      0        0      x      0
      0        1      x      1
      1        x      0      0
      1        x      1      1
*/
// This is conveniently expressed as (in Pascal):
//   change := (old AND death) OR (NOT old AND birth)

// And we can perform this function on all 32 bits at once.

// PowerPC has an instruction (rlwinm) that can & and >>
// immediate, all in one clock cycle.
// Let’s hope the compiler knows how to use it here.

ulong stateChange;
ulong db4_hi,db4_lo,d32=0,b32=0;
long  cnt;

  if ((cnt=cb->count[0])<0) {
    cb->count[0] = cnt & (~MSB);
    db4_hi = rule[(cnt>>13) & 0x1FFF];//mask off flag bits
    db4_lo = rule[ cnt      & 0x1FFF];
    d32 |= ((db4_hi & 0x0F) << 28)|((db4_lo & 0x0F) << 24);
    b32 |= ((db4_hi & 0xF0) << 24)|((db4_lo & 0xF0) << 20);
  }
  if ((cnt=cb->count[1])<0) {
    cnt &= (~MSB);
    cb->count[1] = cnt;
    db4_hi = rule[cnt >> 13];
    db4_lo = rule[cnt &  0x1FFF];
    d32 |= ((db4_hi & 0x0F) << 20)|((db4_lo & 0x0F) << 16);
    b32 |= ((db4_hi & 0xF0) << 16)|((db4_lo & 0xF0) << 12);
  }
  if ((cnt=cb->count[2])<0) {
    cnt &= (~MSB);
    cb->count[2] = cnt;
    db4_hi = rule[cnt >> 13];
    db4_lo = rule[cnt &  0x1FFF];
    d32 |= ((db4_hi & 0x0F) << 12)|((db4_lo & 0x0F) <<  8);
    b32 |= ((db4_hi & 0xF0) <<  8)|((db4_lo & 0xF0) <<  4);
  }
  if ((cnt=cb->count[3])<0) {
    cnt &= (~MSB);
    cb->count[3] = cnt;
    db4_hi = rule[cnt >> 13];
    db4_lo = rule[cnt &  0x1FFF];
    d32 |= ((db4_hi & 0x0F) <<  4)| (db4_lo & 0x0F);
    b32 |=  (db4_hi & 0xF0)       |((db4_lo & 0xF0) >> 4);
  }

  stateChange = (cb->state & d32)
              | (~(cb->state) & b32);
  cb->change=stateChange;
  return stateChange;
}

UnloadCells
void UnloadCells(CellBlock* cb,int size,void* baseAddr) {

// Simple sequential copy of all cell states from
// cell blocks to bit map, loop unrolled once.

  ulong* bm=(ulong*)baseAddr;
  while (size>1) {
    size-=2;
    bm[0]=cb[0].state;
    bm[1]=cb[1].state;
    bm+=2;
    cb+=2;
  }
  if (size) *bm=cb->state;
}

//
// The following 2 functions will handle the loading and unloading of cell blocks for 
// the less likely scenarioes.  Keeping these complications out of the mainstream 
// functions above is really only a tweak to not waste time on frequent checks, 
// noticeable only with very short runs.

LoadCellsSpontaneous
CellBlock* LoadCellsSpontaneous(
           BitMap     cells,
           int        width,
           int        padWords,
           CellBlock* cb,
           CellBlock* L1) {

// Same as LoadCells(..) above, but handles the special case
// of “spontaneous birth” by linking all cells.

ulong* bm=(ulong*)cells.baseAddr;
int    x,y;

  if (width>1) {

//top edge:
    SETBLOCKs(UP+LEFT);
    for (x=1;x<width-1;x++) SETBLOCKs(UP);
    SETBLOCKs(UP+RIGHT);
    bm+=padWords;

//middle rows:
    for (y=1;y<cells.bounds.bottom-1;y++) {
      SETBLOCKs(LEFT);
      for (x=1;x<width-1;x++) SETBLOCKs(0);
      SETBLOCKs(RIGHT);
      bm+=padWords;
    }

//bottom edge:
    SETBLOCKs(DOWN+LEFT);
    for (x=1;x<width-1;x++) SETBLOCKs(DOWN);
    SETBLOCKs(DOWN+RIGHT);
  } else {

//case of bit map <= 32 cells wide
//top edge:
    SETBLOCKs(UP+LEFT+RIGHT);
    bm+=padWords;
//middle rows:
    for (y=1;y<cells.bounds.bottom-1;y++) {
      SETBLOCKs(LEFT+RIGHT);
      bm+=padWords;
    }
//bottom edge:
    SETBLOCKs(DOWN+LEFT+RIGHT);
  }
  return L1;
}

UnloadPaddedCells
void UnloadPaddedCells(
     BitMap     cells,
     int        width,
     int        padWords,
     CellBlock* cb){

// Copy all cell states from cell blocks to bit map,
// row by row, skipping unused pad words at the end of
// each row in the bit map.

int x,y;
ulong* bm=(ulong*)(cells.baseAddr);
  for (y=0;y<cells.bounds.bottom;y++) {
    for (x=0;x<width;x++) {
      *bm++=cb->state;
      cb++;
    }
    bm+=padWords;
  }
}