\datethis
@s class ident

@*Intro. Simple playing with sandpiles on the $3\times$3 torus.

@d space (1<<16)
@c
#include <stdio.h>
int link[8][space], cyc[8][space], class[space], next[space];
int rec[space], hit[space];
int classes,recs;
int x[9];
int v[9][9]={
  {4,-1,-1,-1,0,0,-1,0,0},
  {-1,4,-1,0,-1,0,0,-1,0},
  {-1,-1,4,0,0,-1,0,0,-1},@|
  {-1,0,0,4,-1,-1,-1,0,0},
  {0,-1,0,-1,4,-1,0,-1,0},
  {0,0,-1,-1,-1,4,0,0,-1},@|
  {-1,0,0,-1,0,0,4,-1,-1},
  {0,-1,0,0,-1,0,-1,4,-1},
  {0,0,-1,0,0,-1,-1,-1,4}};
@<Sub...@>@;
main()
{
  register int j,k,t,y;
  @<Explore adding one grain at a time@>;
  @<Explore the union-find approach@>;
  @<Continue the grain-at-a-time approach@>;
  @<Find the minimal recurrent states@>;
}

@ @<Sub...@>=
void reduce()
{
  register int g,j,k;
  for (g=0,j=1;g<8;j=(j&7)+1) {
    if (x[j]<4) g++;
    else {
      for (k=1;k<=8;k++) x[k]-=v[j][k];
      g=0;
    }
  }
}

@ In the first experiment, I add a grain of sand and look at the
resulting functional digraph.

Beware: There's a serious bug in this logic, which makes the
results of experiment~1 pretty much meaningless.


@<Explore adding one grain at a time@>=
for (y=0;y<8;y++) {
  for (k=0;k<space;k++) {
    for (t=k,j=8;j;j--) x[j]=t&3, t>>=2;
    x[y+1]++;
    reduce();
    for (t=x[1],j=2;j<=8;j++) t=(t<<2)+x[j];
    link[y][k]=t;
  }
  for (k=0;k<space;k++) if (!cyc[y][k]) {
    for (j=k;!cyc[y][link[y][j]];j=link[y][j])
      cyc[y][j]=1;
    cyc[y][j]=1, j=link[y][j];
    for (;cyc[y][j]==1;j=link[y][j])
      cyc[y][j]=2;
  }
}
for (recs=k=0;k<space;k++) {
  for (y=0;y<8;y++) if (cyc[y][k]!=2) goto nope;
  recs++, rec[k]=1;
nope:@;
}
printf("Adding one grain at a time leads to %d `recurrent' states.\n",recs);
@<Find an example of a false `recurrent state'@>;

@ @<Find an example of a false `recurrent state'@>=
for (k=0;k<space;k++) if (rec[k]) {
  for (y=0;y<8;y++) if (!rec[link[y][k]]) {
    printf("(But e.g. %04x leads to %04x under %d.)\n",k,link[y][k],y);
    goto found;
  }
}
found:@;

@ In the second experiment, I make states equivalent.

@<Sub...@>=
void yunion(int j, int k)
{
  register int i,ii;
  if (class[j]!=class[k]) {
    for (i=j;class[i]!=class[k];ii=i,i=next[i]) class[i]=class[k];
    next[ii]=next[k], next[k]=j, classes--;
  }
}

@ @<Explore the union-find approach@>=
for (k=0;k<space;k++) class[k]=next[k]=k;
classes=space;
for (k=0;k<space;k++) {
  for (t=k,j=8;j;j--) x[j]=t&3, t>>=2;
  x[1]++, x[2]++, x[3]++, x[6]++;
  reduce();
  for (t=x[1],j=2;j<=8;j++) t=(t<<2)+x[j];
  yunion(k,t);
}
printf("Identifying equivalent states leads to %d classes.\n",classes);
@<Print out the zero class@>;

@ @<Print out the zero class@>=
printf("The zero class contains the following members:\n 0000");
for (k=next[0];k;k=next[k]) printf(" %04x",k);
printf("\n");

@ @<Continue the grain-at-a-time approach@>=
loop:@+for (k=0;k<space;k++) hit[k]=0;
for (k=0;k<space;k++) if (rec[k]) {
  for (y=0;y<8;y++) hit[link[y][k]]|=1<<y;
}
for (k=t=0;k<space;k++) if (rec[k]) {
  if (hit[k]!=0xff) rec[k]=0,t=1,recs--;
  else {
    for (y=0;y<8;y++) if (!rec[link[y][k]]) {
      rec[k]=0;t=1;recs--;@+break;
    }
  }
}
if (t) goto loop;
printf("And now I'm down to %d `recurrent' states.\n",recs);
for (k=0;k<space;k++) if (k==class[k]) {
  for (t=0,j=k;;j=next[j]) {
    if (rec[j]) {
      if (t) printf(" and %04x is equiv to %04x!\n",j,t-space);
      else t=space+j;
    }
    if (next[j]==k) break;
  }
}    

@ @<Find the minimal recurrent states@>=
printf("The minimal ones are:\n");
for (t=k=0;k<space;k++) if (rec[k]) {
  for (y=0;y<16;y+=2) {
    if ((k&(3<<y)) && rec[k-(1<<y)]) goto notmin;
  }
  t++;@+printf(" %04x",k);
notmin:@;
}
printf("\n (%d altogether).\n",t);

@*Index.