Here's the program, written in ANSI C.
Code: Select all
/*
* epsrc.c
* blah 2020
*
* Exhaustive period search.
*
*/
/* cc -ansi -O3 epsrc.c -o epsrc */
#include <stdio.h>
#include <string.h> /* strcmp */
#include <stdlib.h> /* malloc */
/*
* MAP rules are defined as 64-element char arrays in little endian ordering.
* Where the bit positions of the cells in the neighbourhood are:
* 876
* 543
* 210
* indexing into the map at position bit is:
*/
#define MAPBIT(map, bit) (((map)[(bit)>>3] >> ((bit)&7))&1)
/*
* ie B0 is map[0]&1, S0 is map[2]&1, S8 is map[63]>>7, etc
*/
#define TOGGLEMAPBIT(map, bit) ((map)[(bit)>>3] ^= 1<<((bit)&7))
#define SETMAPBIT(map, bit) ((map)[(bit)>>3] |= 1<<((bit)&7))
#define RESETMAPBIT(map, bit) ((map)[(bit)>>3] &= ~(1<<((bit)&7)))
typedef short coord;
typedef unsigned char map[64];
#define OP_ON_MAP(op, a, b, dest) for(i=0;i<64;i++) (dest)[i]=((a)[i] op (b)[i])
#define MAP_ANDNOT(a, b, dest) OP_ON_MAP(&~, a, b, dest)
#define MAP_NOT(src, dest) for(i=0;i<64;i++) (dest)[i] = ~(src)[i]
#define MIN(a, b) ((a)<(b)?(a):(b))
#define MAX(a, b) ((a)>(b)?(a):(b))
coord *w, *next_w;
int pop;
int max_period = 20;
int min_period = 0;
void runstep(int par, map *rule, map *gen_to){
/* Run one step of *rule on w, writing to next_w, with parity par; par is 0
* on even generations (when the background of the universe is state 0),
* and 1 on odd generations (when the background is 1). If gen_to is nonzero
* then *gen_to has a 1 written to each bit whose corresponding transition
* occurs at least once (this may or may not include 0777 or 0000).
*
* The world is represented as a singly-linked list of sparse arrays. That
* is, more precisely, each row consists of a number n, which is the amount
* of cells in the row plus one, followed by n-1 numbers, which are the x
* coordinates of the cells, sorted (!) such that the leftmost (with the
* lowest x coordinate) is at the lowest address. The top and bottom rows
* are non-empty. The value after the end of the final row is 0. The first
* number in the array is the y coordinate of the first row, and the first
* row immediately follows. For example:
* ###
* # = 0 4 0 1 2 2 2 2 1 0
* #
*
* 3 pointers (r0, r1, r2) point to the 3 rows around and including the one
* currently being simulated, if such a row is within bounds, else they are
* NULL. They are dereferenced and decremented to get the offsets to the end
* (r0o, r1o, r2o).
*
* The row being simulated has its cells iterated through linearly. The
* neighbourhood is calculated and then used (inverted if par=1) as an index
* into *rule to determine the outcome of a cell. Since 6 cells are shared
* between the neighbourhoods of x and x+1, bitwise shifting is used and
* only the 3 cells added onto the right are fetched.
*
* 3 index counters (r0i, r1i, r2i) are used to keep track of the cells that
* are to be added to the right of the neighbourhood. If rn[rni] (n=0,1,2)
* is to the right of x (rn[rni]==x+1), the cell is added to the
* neighbourhood and rni is incremented.
*/
int y=*w;
int next_pop=0;
const coord *r0, *r1, *r2;
coord *p=next_w+1;
coord *end=p;
int r0o, r1o, r2o;
r0 = NULL;
r1 = NULL;
r2 = w+1;
for(;r2 || r0; r0=r1,r1=r2,r2?r2+=*r2:0, y++){ /* for each row */
/* create new row if the first one isn't empty (trim leading rows) */
if(p==next_w+2) p--;
coord *row_start = p++;
*row_start = 1;
/* if pointer to row is NULL, treat it as empty to avoid segfaults */
r2o=r2?*r2-1:0;
r1o=r1?*r1-1:0;
r0o=r0?*r0-1:0;
/* prevent confusion on null row terminator */
if(r2o==-1) r2o=0;
if(r1o==-1) r1o=0;
/* (r0 never points to the end within the loop body) */
if(r2 && !*r2) r2 = NULL; /* stop loading new addrs at the end */
if(! (r2o|r1o|r0o)) continue; /* all 3 rows empty */
/* calculate left and right bounds */
int left, right;
/* make sure horz. bounds are valid for at least one row */
if(r2o) left=r2[1], right=r2[r2o];
if(r1o) left=r1[1], right=r1[r1o];
if(r0o) left=r0[1], right=r0[r0o];
/* then adjust to max and min */
if(r2o) left =MIN(left, r2[1]),
right=MAX(right,r2[r2o]);
if(r1o) left =MIN(left, r1[1]),
right=MAX(right,r1[r1o]);
/* we don't need to adjust for r0; if it's nonempty it was already the
* starting value */
left--, right++;
int x, r0i=1, r1i=1, r2i=1;
int neigh=(-par)&511;
for(x = left;x<=right;x++){ /* simulate row */
int i;
if(r2i<=r2o && r2[r2i]<x+1) r2i++;
if(r2i<=r2o && r2[r2i]==x+1) neigh^=1<<0;
if(r1i<=r1o && r1[r1i]<x+1) r1i++;
if(r1i<=r1o && r1[r1i]==x+1) neigh^=1<<3;
if(r0i<=r0o && r0[r0i]<x+1) r0i++;
if(r0i<=r0o && r0[r0i]==x+1) neigh^=1<<6;
if(MAPBIT(*rule,neigh)^par^1){
/* cell will exist in next tick */
*p++ = x;
next_pop++;
(*row_start)++;
}
if(gen_to) SETMAPBIT(*gen_to, neigh);
neigh = ((neigh<<1)&0666)|(0111*par); /* transl. neigh by 1 cell */
}
if(row_start!=p-1){ /* row is non-empty */
end = p; /* (trim trailing rows) */
if(row_start==next_w+1) /* first non-empty row */
next_w[0]=y-1;
}
}
*end=0; /* null row terminator */
coord *swp; swp=w; w=next_w; next_w=swp; pop=next_pop;
}
int inc_masked_MAP(map *cntr, map *mask){
int i;
/* Go through the number, ignoring bits that are 0 in mask, toggling every
* bit and advancing if it is a 1, else exiting.
*/
for(i=0;i<512;i++){
if(MAPBIT(*mask, i)){ /* only perform counting on mask bits */
TOGGLEMAPBIT(*cntr, i);
if(MAPBIT(*cntr, i)) break; /* if it was a 0, we're done */
}
}
/* The loop has 2 exit conditions: a 0 was encoutered, breaking the loop
* early, or the value overflowed, and the i<512 check failed.
*/
return i==512; /* return carry out */
}
char *MAP_to_string(map *m){
const char *base64_digs =
"ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";
static char s[90];
int ch, bit, bit6;
s[0]='M', s[1]='A', s[2]='P', s[89]='\0';
for(ch=0;ch<86;ch++){
for(bit=bit6=0;bit<6;bit++){
int mapbit=6*ch+bit;
bit6 = bit6*2 + (mapbit<512?MAPBIT(*m, mapbit):0);
}
s[ch+3] = base64_digs[bit6];
}
return s;
}
int from_base64_digit(unsigned char c){
if((unsigned)(c-'0')<10) return c-'0'+52;
if((unsigned)(c-'A')<26) return c-'A';
if((unsigned)(c-'a')<26) return c-'a'+26;
if(c=='+') return 62;
if(c=='/') return 63;
fprintf(stderr, "Invalid base 64 digit: %c", c);
return 0; }
void string_to_MAP(char *s, map *m){
int bit6, i;
for(i=0;i<64;i++) (*m)[i]=0;
s+=2; /* 'MAP' header */
i = 0;
while(*++s){
bit6 = from_base64_digit(*s);
if(bit6&32) TOGGLEMAPBIT(*m, i+0);
if(bit6&16) TOGGLEMAPBIT(*m, i+1);
if(i==510) break; /* last digit only encodes 2 meaningful bits */
if(bit6&8) TOGGLEMAPBIT(*m, i+2);
if(bit6&4) TOGGLEMAPBIT(*m, i+3);
if(bit6&2) TOGGLEMAPBIT(*m, i+4);
if(bit6&1) TOGGLEMAPBIT(*m, i+5);
i+=6;
}
}
int run_rule(map *rule, map *gen_to){
int t;
/* init world with 1 cell */
pop = 1;
w[0]=0; w[1]=2; w[2]=0; w[3]=0;
for(t=2;t<max_period+2;t+=2){ /* simulate 2 steps at a time */
runstep(0, rule, gen_to);
runstep(1, rule, gen_to);
if(pop<=1) break;
}
if(pop==1) return t;
return 0;
}
/* use linked list to keep track of already known speeds */
typedef struct {
int period;
int x,y;
void *next;
} speed;
speed *spds=0;
void characterise(map *rule){
int t = run_rule(rule, 0);
if(pop==1 && t>=min_period){
/* linear search through list to ensure speed is new */
speed **p;
for(p=&spds;*p;p=(speed**)(&((*p)->next)))
if((*p)->period==t && (*p)->x==w[2] && (*p)->y==w[0])
return;
/* speed is new */
(*p) = malloc(sizeof(speed));
(*p)->period = t;
(*p)->x = w[2];
(*p)->y = w[0];
(*p)->next = 0;
printf("p%d (%d,%d): %s\n", t, w[2], w[0],
MAP_to_string(rule));
}
}
void showhelp(){
printf(
"epsrc: Exhaustive Period Search\nblah 2020\n\n"
"A program that finds alternating 2-state cellular automata defined on the\n"
"range-1 moore neighbourhood in which a single state 1 cell surrounded by 0\n"
"returns to being a single cell in some location after some (even) number of\n"
"steps of the simulation. In this help, a rule and its 1-cell object are\n"
"referred to as if they are the same.\n\n"
"Usage:\n"
" epsrc cha (RULE) [OPTIONS]\n"
" Characterise. Display speed of RULE, and MAP mask of all transitions\n"
" that occur in its evolution.\n"
" epsrc enum (RULE1) (RULE2) [OPTIONS]\n"
" Enumerate. Allow bits that are 1 in RULE1 to vary in RULE2. Characterise\n"
" all rules within this sub-rulespace.\n"
" epsrc ham (NUM) (RULE) [OPTIONS]\n"
" Hamming distance. Characterise all alternating rules with hamming\n"
" distance of NUM from RULE.\n"
" (NUM is an integer between 1 and 510 inclusive.)\n"
" epsrc help\n"
" Display this help and exit.\n"
"A MAP rulestring can have its base64 component truncated; left out digits\n"
"are fillied in with 0s ('A's). For example, in \"epsrc ham 2 MAP\", \"MAP\"\n"
"expands to \"MAP\" with 86 As after it.\n\n"
"Items in parentheses are necessary. Items in square brackets are optional.\n\n"
"Options:\n"
" -p (NUM) Maximum period; disregard rules whose populations do not go below\n"
" 2 in NUM generations. NUM should be an even integer above 0.\n"
" -m (NUM) Minimum period. Do not show rules with periods above this integer.\n"
);
}
int main(int argc, char **argv){
if(argc<2 || strcmp(argv[1],"help")==0) {showhelp(); return 1;}
int i;
/* options */
for(i=0;i<argc;i++){ /* get universal parameters */
if(strcmp(argv[i],"-p")==0 && i+1<argc)
sscanf(argv[i+1], "%d", &max_period);
else if(strcmp(argv[i],"-m")==0 && i+1<argc)
sscanf(argv[i+1], "%d", &min_period);
}
/* initialise world buffers */
#define WORLD_COORDS ((2*max_period+1) * (2*max_period+1) + max_period*2 + 2)
coord *wlds = malloc(2 * WORLD_COORDS * sizeof(coord));
w = wlds; next_w = wlds+WORLD_COORDS;
if(strcmp(argv[1],"cha")==0){
/* characterise */
map to_cha, cha_result;
for(i=0;i<64;i++) to_cha[i]=cha_result[i]=0;
string_to_MAP(argv[2], &to_cha);
int t=run_rule(&to_cha, &cha_result);
if(!t && pop>0){
printf("No periodicity detected.\n");
} else if(pop==0){
printf("Dot dies.\n");
} else { /* spaceship/oscillator */
/* b0 and s8 may not be changed */
RESETMAPBIT(cha_result, 0);
RESETMAPBIT(cha_result, 511);
printf("p%d (%d,%d)\n", t, w[2], w[0]);
printf("%s\n", MAP_to_string(&cha_result));
}
} else if(strcmp(argv[1],"enum")==0) {
/* enumerate */
map base_rule, rule_mask;
for(i=0;i<64;i++) base_rule[i]=rule_mask[i]=0;
string_to_MAP(argv[2], &rule_mask);
string_to_MAP(argv[3], &base_rule);
RESETMAPBIT(rule_mask,0); RESETMAPBIT(rule_mask,511);
MAP_ANDNOT(base_rule, rule_mask, base_rule);
do {
characterise(&base_rule);
} while(!inc_masked_MAP(&base_rule, &rule_mask)); /* until overflow */
} else if(strcmp(argv[1],"ham")==0){
/* test changes to a rule with supplied hamming distance */
int ham_dis; sscanf(argv[2], "%d", &ham_dis);
if(ham_dis==0){
fprintf(stderr,"Hamming distance of 0 not accepted. "
"Use the 'cha' command if you want to evaluate a rule.\n");
return 1;
}
if(ham_dis<0 || ham_dis>510){
fprintf(stderr,"Hamming distance out of range (1-510 inclusive)\n");
return 1;
}
/* The principle is as such: Enumeration is done in-place on two binary
* strings, a counter and the base string. Whenever a bit is toggled in
* the counter it is also toggled in the base string. The task is then
* to loop the counter through every possible value of popcount ham_dis.
*
* We go in an order like 1100 1010 0110 1001 0101 0011; the leftmost 1
* bit is moved to the right and if it moves into any other 1s they get
* pushed over as well.
*
* 10 Start with a sequence of ham_dis 1s on the left, all other bits 0.
* 20 Find the leftmost 1 bit.
* 30 Toggle it and, if the next bit is a 1, set the leftmost 0 bit.
* 40 Move to next bit. If it's 1, GOTO 30, else set it and GOTO 20.
*
* Try running 001101 through that algorithm in your head.
*/
map ctr, base_rule;
for(i=0;i<64;i++) ctr[i]=0;
string_to_MAP(argv[3], &base_rule);
/* note that bit positions 0 and 511 must be 1 and 0, respectively;
* otherwise rule is non-alternating and out of the scope of epsrc */
SETMAPBIT(base_rule,0); RESETMAPBIT(base_rule, 511);
for(i=1;i<=ham_dis;i++) SETMAPBIT(ctr,i), TOGGLEMAPBIT(base_rule,i);
do {
characterise(&base_rule);
/* get next rule */
for(i=1;!MAPBIT(ctr,i);i++); /* find leftmost 1 bit */
int begin = i;
while(1){
/* remove 1 */
TOGGLEMAPBIT(ctr,i);
TOGGLEMAPBIT(base_rule,i);
i++;
if(!MAPBIT(ctr,i)) break;
/* place it back at the beginning */
TOGGLEMAPBIT(ctr,(i-begin)+1);
TOGGLEMAPBIT(base_rule,(i-begin)+1);
}
TOGGLEMAPBIT(ctr,i);
TOGGLEMAPBIT(base_rule,i);
} while(!MAPBIT(ctr,511));
} else {
printf("Command '%s' not recognised. Use '%s help' to list commands.\n",
argv[1], argv[0]);
}
free(wlds);
return 0;
}
I can explain this in more detail if anyone wants. If you have problems using the program, please post about them in this thread.