\datethis

@*Introduction. This program determines the number of paths that run from the
source vertex to sink vertices at each level of a given periodic directed
acyclic graph. Although it was designed originally to process the output of
{\mc POLYENUM}, it might also be useful for processing the output of
other programs.

The input file consists of 32-bit words, where the top 5 bits are an
``opcode,'' the next bit is a ``tag,'' and the remaining 26 bits are
the ``address.''

Opcode 0 means: Choose a new vertex as the current vertex from which
successive arcs will emanate. The address is the vertex number.

Opcodes 1--30 mean: Create an arc of this multiplicity from the current
vertex to another vertex specified by the tag and the address.
If the tag is~1, the address gives the serial number of the destination vertex;
all vertices are given a serial number, starting with~1. Address 0 is,
however, legal; it denotes the sink vertex on the current level.
(In the program {\mc POLYENUM}, the sink vertex on
level~$m$ represents all polyominoes that have exactly $m$ cells.)

If the tag is~0,
the address has the form |(a<<8)+b|, where $a$ and $b$ are the
``birth'' and ``death'' levels of a new vertex; this new vertex
is assigned the next serial number. (Levels will be explained
momentarily.) For convenience in English, we
call $a$ and $b$ the birthdate and deathdate. Each new vertex
defined on level~$m$ must have birth and death dates satisfying
$m<a\le b\le n$, where $n$ is the maximum level.

Opcode 31 means: Go to a new level. The address is the new level~$m$, which
should be exactly one greater than the previous level. Initially the
program is at level zero, at which time all arcs emanate from a
special anonymous vertex.

The directed acyclic graph is {\it periodic\/} in the sense that
clones of its vertices can appear almost isomorphically on many different
levels. A vertex $v$ with birthdate~$a$ and deathdate~$b$ actually stands
for $b-a+1$ vertices $v_a$, $v_{a+1}$, \dots,~$v_b$. An arc from
vertex $u_k$ to vertex $v_l$ also implies the existence of arcs from
$u_{k+j}$ to $v_{l+j}$, for all $j$ such that $k+j$ is within the
``lifetime'' of~$u$ and $l+j$ is within the lifetime of~$v$.

We assume that each vertex has a ``weight'' $w$ such that an arc from
$u_k$ to $v_l$ exists only if $l=k+w$. This property holds true in the
output of {\mc POLYENUM}, where vertices represent slices through a
polyomino and $w$ is the number of cells in the slice. I could probably
revise the program to remove this assumption, but I'm keeping it because
it provides an extra sanity check on the integrity of the input file.

The anonymous source vertex lives only at level zero. Successors
to this vertex should not be given with the tag-bit option; they should
always be new vertices.

The final level~$n$ is specified on the command line when this program is
run. When the final level is reached, no further instructions are read;
arcs are created from all vertices that are still alive to the
sink vertex on level~$n$. (Arcs to vertex 0 are implicit, not
read from the input, for all vertices born on the final level.)
Otherwise, every vertex must appear exactly once as a
source vertex (with opcode~0); this appearance must, in fact, occur
on level~$a$, the time of its birth.
Moreover, every clone of every vertex must have at least one outgoing arc.

@d new_pred_code 0 /* high 5 when new vertex is the predecessor */
@d new_level_code 31 /* high 5 when $m$ increases */
@d tag_bit (1<<26)

@c
#include <stdio.h>
@<Type definitions@>@;
@<Global variables@>@;
@<Subroutines@>@;

main(int argc, char *argv[])
{
  @<Local variables@>;
  @<Scan the command line@>;
  @<Initialize@>;
  @<Interpret the instructions in the input@>;
  @<Print statistics@>;
}

@ The program checks frequently that everything in the input file
is legitimate. If not, too bad; the run is terminated (although
a debugger can help diagnose nonobvious problems). Extensive checks like this
have helped the author to detect errors in the program as well as errors in
the input.

@<Sub...@>=
void panic(char *mess)
{
  fprintf(stderr,"%s!\n",mess);
  exit(-1);
}

@ Several gigabytes of input might be needed, so the input file name will
be extended by \.{.0}, \.{.1}, \dots, just as in {\mc POLYENUM}.

@d nmax 30 /* at most this many levels */
@d bad(k,v) sscanf(argv[k],"%d",&v)!=1

@<Scan the command line@>=
if (argc!=6 || bad(1,n) || bad(2,verts) || bad(3,clones) || bad(4,arcs)) {
  fprintf(stderr,"Usage: %s n verts clones arcs infilename\n",argv[0]);
  exit(0);
}
if (n<=0 || n>nmax) {
  fprintf(stderr,"Sorry, n must be between 1 and %d; try again",nmax);
  panic("");
}
if (verts<n) verts=n;
if (clones<n) clones=n;
if (arcs<1) arcs=1;
base_name=argv[5];

@ Memory space for vertices and clones is allocated permanently, when
this program begins, so the command-line parameters |verts| and |clones|
must be large enough to accommodate all of them.
For example, when processing the output of {\mc POLYENUM}, |verts| must be
at least as large as the number of patterns, and |clones| must be at
least as large as the number of slices.

But this program discards arc information as soon as it is no longer needed,
so the user can try to conserve memory by asking for only a limited
amount of arc space. (I could of course do dynamic allocation, using |calloc|
whenever more arc memory was needed; but look, this is not a general-purpose
program, it's a program
for me only, and I don't want to take a chance on thrashing after
the program has run a long time.)

@<Glob...@>=
int n; /* the maximum level */
int verts; /* upper bound on number of proto-vertices */
int clones; /* upper bound on total number of vertices */
int arcs; /* upper bound on memory needed to store the proto-arcs */

@*Input. We might as well start small, with the basic routines that
are needed to read instructions from the input file(s).
As soon as $2^{30}$ bytes of data have been read from file \.{foo.0},
we'll turn to file \.{foo.1}, etc.

@d filelength_threshold 0x10000000 /* in tetrabytes */

@<Glob...@>=
FILE* in_file; /* the input file */
unsigned int buf; /* place for binary input */
int words_in; /* the number of tetrabytes input from the current input file */
int file_extension; /* the number of GGbytes input */
char *base_name, filename[100];

@ @<Sub...@>=
void open_it()
{
  sprintf(filename,"%.90s.%d",base_name,file_extension);
  in_file=fopen(filename,"rb");
  if (!in_file) {
    fprintf(stderr,"I can't open file %s",filename);
    panic(" for input");
  }
  words_in=0;
}

@ @<Sub...@>=
void close_it()
{
  if (fclose(in_file)!=0) panic("I couldn't close the input file");
  printf("[%d bytes read from file %s.]\n",4*words_in,filename);
}

@ @<Sub...@>=
int get_it()
{
  if (words_in==filelength_threshold) {
    close_it();
    file_extension++;
    open_it();
  }
  words_in++;
  return fread(&buf,sizeof(unsigned int),1,in_file)==1;
}

@ @<Init...@>=
open_it();

@* Arithmetic. The integer numbers involved here can get very large,
but we assume that one octabyte (64 bits) is enough to hold them. My
computer deals directly with 32-bit words, so I represent an octabyte
as a pair of unsigned ints:

@<Type...@>=
typedef struct {
  unsigned int h,l; /* high-order and low-order halves */
} octa; /* two tetrabytes make one octabyte */

@ @<Sub...@>=
octa oplus(octa y,octa z) /* compute $(y+z)\bmod2^{64}$ */
{
  octa x;
  x.h=y.h+z.h, x.l=y.l+z.l;
  if (x.l<y.l) x.h++;
  return x;
}

@ @<Sub...@>=
void print_octa(octa o)
{
  register unsigned int hi=o.h, lo=o.l, r, t;
  register int j;
  char dig[20];
  if (lo==0 && hi==0) printf("0");
  else {
    for (j=0;hi;j++) { /* 64-bit division by 10 */
      r=((hi%10)<<16)+(lo>>16);
      hi=hi/10;
      t=((r%10)<<16)+(lo&0xffff);
      lo=((r/10)<<16)+(t/10);
      dig[j]=t%10;
    }
    for (;lo;j++) {
      dig[j]=lo%10;
      lo=lo/10;
    }
    for (j--;j>=0;j--) printf("%c",dig[j]+'0');
  }
}

@*Data structures.
Basic information about a proto-vertex $v$ is kept in
the vertex node |vtable[v]|. If this vertex has clones $v_a$,
$v_{a+1}$, \dots,~$v_b$, the number of paths from the source
to each clone is accumulated in $b-a+1$ consecutive entries
of |ctable|. The sink vertices for each level are considered to
be clones of |vtable[0]|.

@s arc_struct int

@<Type...@>=
typedef struct vertex_struct {
  char birth; /* the birthdate $a$ */
  char death; /* the deathdate $b$ */
  char weight; /* the level difference from all predecessors */  
  char cleared; /* has the last clone been processed? */
  octa *base; /* the counter for the first clone, $v_a$ */
  struct vertex_struct *next; /* the vertex to process next */
  struct arc_struct *arcs; /* the first arc from this vertex */
} vert; /* 20 bytes per vertex, plus 8 bytes per clone */
@#
typedef struct arc_struct {
  unsigned int tip; /* multiplicity and destination */
  struct arc_struct *next; /* the next arc from the same vertex */
} arc; /* 8 bytes per arc */

@ @<Glob...@>=
vert *vtable; /* master table for all proto-vertices */
vert *bad_vert; /* address at the end of |vtable| */
vert *vtable_ptr; /* the first unused space in |vtable| */
octa *ctable; /* counters for all vertex clones */
octa *bad_count; /* address at the end of |ctable| */
octa *ctable_ptr; /* the first unused space in |ctable| */
arc *atable; /* current pool of arc information */
arc *bad_arc; /* address at the end of |atable| */
arc *atable_ptr; /* the first unused space in |atable| */

@ @<Init...@>=
verts++; /* make room for the proto-sink vertex */
vtable=(vert*)calloc(verts,sizeof(vert));
if (!vtable) panic("I can't allocate the vertex table");
bad_vert=vtable+verts;
clones+=n; /* make room for the sink clones */
ctable=(octa*)calloc(clones,sizeof(octa));
if (!ctable) panic("I can't allocate enough space for counters");
bad_count=ctable+clones;
atable=(arc*)calloc(arcs,sizeof(arc));
if (!atable) panic("I can't allocate that much space for arcs");
bad_arc=atable+arcs;
vtable[0].base=ctable;
vtable[0].birth=1;
vtable[0].death=n;
vtable_ptr=vtable+1;
ctable_ptr=ctable+n; /* the |n| sink vertices occur at the start of |ctable| */
atable_ptr=atable;

@ When a vertex is first mentioned, we place it in a |pending| list of all
vertices having the same birthdate; the |next| field will then point
to the previous vertex sharing that date. After level |m| has been
processed, the list of all vertices born on that level is scanned to
make sure that each one had outgoing arcs. Subsequently the
|next| fields are used to link together all proto-vertices that are
still alive.

@ @<Glob...@>=
vert *pending[nmax+1]; /* lists of vertices to be defined on future levels */
vert *alive; /* list of all vertices still kicking */

@ When an |arc| is no longer needed, it is placed on the |avail| list,
from whence it can be reused. Here is a subroutine that returns a fresh |arc|,
using a familiar allocation ritual:

@<Sub...@>=
arc *new_arc()
{
  register arc*a;
  if (avail) a=avail, avail=a->next, a->next=NULL;
  else if (atable_ptr==bad_arc) panic("Arc memory overflow");
  else a=atable_ptr++;
  return a;
}

@ @<Glob...@>=
arc *avail; /* top of the stack of recycled |arc| nodes */

@* The basic execution cycle. This program is structured like a simple
form of machine simulator, reading instructions one by one.
Luckily, the basic control loop is quite simple,
because we aren't simulating a stored-program computer.

(On the other hand, we do store arc information that is used repeatedly, so
the complexity shows up elsewhere.)

The program will |goto done| when level |n| is reached.

@<Interpret the instructions in the input@>=
while (1) {
  @<Read and unpack the next instruction@>;
  switch(opcode) {
 case new_pred_code: case new_level_code:
  if (u) @<Process all arcs from |u| at level |m|@>;
  if (opcode==new_level_code) @<Begin a new level@>@;
  else @<Define a new predecessor vertex@>;
  break;
 default: @<Create a new proto-arc from the current proto-vertex@>;
  }
  loc++;
}
done: @<Make the final total@>;

@ @d address_mask ((1<<27)-1)

@<Read and unpack the next instruction@>=
if (!get_it()) panic("Premature end of input");
opcode=buf>>27;
address=buf&address_mask;
if (verbose) print_inst(0);

@ @<Glob...@>=
int verbose=0; /* set nonzero when debugging */
int loc; /* the number of instructions processed */
unsigned int buf; /* the current instruction from |in_file| */
int opcode; /* its operation code */
int address; /* its address field, including the tag */
int mm; /* global copy of local variable |m| */

@ Sometimes we want to see an instruction in symbolic form.
Each proto-vertex is identified by its decimal serial number,
starting with~1 (not starting with~0 as in the binary file).
Of course this serial number is not the same as the hexadecimal number of the
corresponding pattern in {\mc POLYENUM}; we haven't been told that number.
But the verbose outputs of {\mc DAGENUM} and {\mc POLYENUM} are otherwise
quite similar.

@<Sub...@>=
void print_inst(int err)
{
  register int cur_verts=vtable_ptr-vtable,
     taddress=address^tag_bit;
  register FILE* f=err? stderr: stdout;
  switch (opcode) {
case new_level_code: fprintf(f,"New level %d\n",address);@+break;
case new_pred_code: fprintf(f,"%d:%d\n",address,
  address<cur_verts? vtable[address].birth: -1);@+break;
case 1: fprintf(f,"->"); goto common_ending;
default: fprintf(f,"-%d>",opcode);
common_ending: if (address<tag_bit)
  fprintf(f,"%d:%d..%d\n",cur_verts,address>>8,address&0xff);
else fprintf(f,"%d:%d\n",taddress,
   taddress<cur_verts? vtable[taddress].weight+mm: -1);
  }
}

@ @<Sub...@>=
void ipanic(char *mess)
{
  if (!verbose) fprintf(stderr,"Problem at instruction %d:\n",loc);  
  print_inst(1); /* show context before dying */
  panic(mess);
}

@*Creating vertices and arcs. As this program runs, |m| will be the
current level, |u| will be the current (proto-)vertex from which
arcs are being generated, and |v| will be the current (proto-)vertex
to which a new arc is being directed.

@<Local...@>=
register int m=0; /* the cloning level */
register vert *u=NULL,*v=NULL; /* current vertices of interest */

@ @<Create a new proto-arc from the current proto-vertex@>=
if (address>=tag_bit) {
  v=vtable+(address^tag_bit);
  if (v>=vtable_ptr) ipanic("Vertex number out of range");
  if (m+v->weight>v->death) ipanic("Vertex is already dead");
}@+else {
  if (vtable_ptr==bad_vert) ipanic("Vertex memory overflow");
  v=vtable_ptr++;
  v->birth=address>>8, v->death=address&0xff;
  if (v->birth<=m) ipanic("Birthdate too early");
  if (v->death>n) ipanic("Deathdate too late");
  if (v->death<v->birth) ipanic("Death before birth");
  v->weight=v->birth-m;
  v->base=ctable_ptr;
  ctable_ptr+=v->death-v->birth+1;
  if (ctable_ptr>bad_count) ipanic("Counter memory overflow");
  v->next=pending[v->birth], pending[v->birth]=v;
}
@<Create a new arc from |u| to |v|@>;

@ The value |u=NULL| represents
the source vertex, which is legal and relevant only on level zero.
The value |v=vtable| represents the sink vertex on the current level.

Incidentally, multiple arcs are allowed, in case we need to deal
with graphs that have multiplicities exceeding~30.

@<Create a new arc from |u| to |v|@>=
if (u) {
  buf=(opcode<<27)+(v-vtable);
  s=new_arc();
  s->tip=buf;
  s->next=u->arcs, u->arcs=s;
} else {
  if (m) ipanic("No predecessor specified");
  if (v==vtable) ipanic("There is no sink on level zero");
  if ((*v->base).l) ipanic("Initial arc not to a new vertex");
  (*v->base).l=opcode;
} /* the code now falls through to the end of the master opcode switch */

@ @<Local...@>=
register arc *r,*s; /* current arcs of interest */

@ @<Define a new predecessor vertex@>=
{
  u=vtable+address;
  if (u>=vtable_ptr) ipanic("Predecessor vertex not yet defined");
  if (u->birth!=m) ipanic("This vertex has the wrong birthdate");
}

@ @<Begin a new level@>=
{
  @<Finish everything that needs to be done on level |m|@>;
  if (m) {
    print_octa(ctable[m-1]);@+printf(" paths to sink[%d]\n",m);
  }
  if (address!=++m) ipanic("Level number is out of sync");
  mm=m;
  printf("Entering level %d (%d)\n",m,atable_ptr-atable);
  if (m==n) goto done;
  u=NULL;
}

@* Accumulating path counts. Now we're ready for the main
action, which is to add $c$ times the counter for $u_m$ to the
counter for $v_{m+w}$ for all arcs that run from $u$ to $v$.

The multiplicity $c$ is usually 1 and almost never very large.
Therefore we don't actually multiply, we simply use addition to
maintain a table of
multiples of |u|'s counter, building that table on the fly as needed.

@<Glob...@>=
octa mcount[31]; /* |mcount[c]=c*mcount[1]| */
int mcount_ptr; /* the number of valid entries of |mcount| */
int recycling; /* is |u| is about to die? */

@ The algorithm would be quite simple if we didn't take the trouble to
recycle arcs that are no longer needed. But hey, that makes things
interesting.

@<Process all arcs from |u| at level |m|@>=
{
  mcount[1]=*(u->base+m-u->birth);
  mcount_ptr=1;
  recycling=(u->death==m);
  if (!u->arcs) ipanic("Vertex has no successor arcs");
  for (r=NULL,s=u->arcs;s;)
    @<Process arc |s| and move to the next, with |r| following |s|@>;
  if (recycling)
    r->next=avail, avail=r, u->cleared=1,  u->arcs=NULL;
}

@ @<Process arc |s|...@>=
{@+register octa *o; register int k;
  k=s->tip>>27; /* the multiplicity */
  while (k>mcount_ptr)
    mcount[mcount_ptr+1]=oplus(mcount[1],mcount[mcount_ptr]),mcount_ptr++;
  v=vtable+(s->tip&address_mask);
  o=v->base+m+v->weight-v->birth;
  *o=oplus(mcount[k],*o);
  if (!recycling && m+v->weight==v->death) { /* recycle this arc only */
    if (!r) u->arcs=s->next, s->next=avail, avail=s, s=u->arcs;
    else r->next=s->next, s->next=avail, avail=s, s=r->next;
  }@+else r=s,s=s->next;
}

@ @<Finish everything that needs to be done on level |m|@>=
@<Check all vertices born at level |m|@>;
@<Process arcs from older vertices still alive@>;
@<Delink vertices that die at level |m|@>;

@ @<Check all vertices born at level |m|@>=
if (pending[m]) {
  for (u=pending[m]; u; v=u,u=u->next)
    if (!u->cleared && !u->arcs)
      ipanic("Missing arcs on previous level");
  v->next=alive;
}

@ @<Process arcs from older vertices still alive@>=
for (u=alive; u; u=u->next)
  @<Process all arcs from |u|...@>;

@ @<Delink vertices that die at level |m|@>=
if (pending[m]) alive=pending[m];
for (v=NULL,u=alive;u;u=u->next)
  if (u->cleared) {
    if (!v) alive=u->next;
    else v->next=u->next;
  }@+else if (!u->arcs) ipanic("An old vertex became a sink");
  else v=u;

@ @<Make the final total@>=
for (u=pending[m];u;u=u->next)
  ctable[m-1]=oplus(*u->base,ctable[m-1]);
for (u=alive;u;u=u->next) {
  s=u->arcs;
  if (!s || s->tip!=(1<<27) || s->next)
     ipanic("Vertex alive at level n should point to sink[n]");
  ctable[m-1]=oplus(*(u->base+m-u->birth),ctable[m-1]);
}
print_octa(ctable[m-1]);
printf(" paths to sink[%d]\n",m);

@ @<Print statistics@>=
printf("All done! (%d arc slots needed)\n",atable_ptr-atable);
close_it();

@*Index.