/* file fsacomposite.c  24. 11. 94.
 * 29/9/98 large-scale re-organisation.
 * 3/8/98 - put in code for integer hash tables.
 * 13/1/98 - changes for `gen' type of generator replacing char.
 *
 * This file contains the routines necessary for axiom checking in an
 * automatic group.
 *
 * 13/9/95 - make the length 2 generalised multiplier have state set of
 * type setToLabels, where the labels are lists of words (of length 2),
 * rather than words.
 */

#define LABHTSIZE 8192

#include <stdio.h>
#include "defs.h"
#include "fsa.h"
#include "hash.h"
#include "externals.h"

/* Functions defined in this file: */
fsa  *fsa_genmult2();
fsa  *fsa_genmult2_short();
fsa  *fsa_genmult2_int();
int  fsa_makemult();
int  fsa_makemult2();
fsa  *fsa_composite();
fsa  *fsa_composite_short();
fsa  *fsa_composite_int();

/* Functions used in this file and defined elsewhere */
void compressed_transitions_read();
boolean fsa_equal();
int sparse_target();
void fsa_init();
void fsa_set_is_accepting();
void srec_copy();
void srec_clear();
void fsa_clear();
void hash_init();
void short_hash_init();
int  hash_locate();
int  short_hash_locate();
void hash_clear();
void short_hash_clear();
int* hash_rec();
unsigned short* short_hash_rec();
int hash_rec_len();
int short_hash_rec_len();
int stringlen();

fsa *
fsa_genmult2(genmultptr,op_table_type,destroy,genmult2filename,readback)
	fsa *genmultptr;
	storage_type op_table_type;
	boolean destroy;
	char *genmult2filename;
	boolean readback;
/* *genmultptr should be the general multiplier fsa of an automatic group.
 * This function calculates the transition table of the general multiplier for
 * products of two generators. This is output (in unformatted form) to the
 * file with name genmult2filename, followed by a list of states, which
 * specifies which states are accept-states for which length-two words.
 * It can then be minimised for any such word as required.
 * If readback is false, then
 * the fsa returned does not have its transition table stored internally
 * (since these are in the file genmult2filename) - instead, its table
 * component table->filename contains the name of this file.
 * Otherwise, the transition table is read back in as usual.
 * If genmult2filename is the empty string, then the transition table is
 * not stored at all.
 * (This can be done when all relators of the group have length at most 4.)
 */
{
  if (kbm_print_level>=3)
    printf("    #Calling fsa_genmult2.\n");
  if (genmultptr->states->size < MAXUSHORT)
    return
      fsa_genmult2_short(genmultptr,op_table_type,
                           destroy,genmult2filename,readback);
  else
    return fsa_genmult2_int(genmultptr,op_table_type,
                           destroy,genmult2filename,readback);
}

fsa *
fsa_genmult2_short(genmultptr,op_table_type,destroy,genmult2filename,readback)
	fsa *genmultptr;
	storage_type op_table_type;
	boolean destroy;
	char *genmult2filename;
	boolean readback;
{
  int **table, ne, ngens, ngens1, ns, *fsarow, e, e1, e2, es1, ef1, dr,
      nt, cstate, im, csa, csb, csima, csimb,
      i, j, j1, j2, g1, g2, g3, len, reclen, nlab, ct;
  unsigned short *ht_ptr, *ht_chptr, *ht_ptrb, *ht_ptre,
                 *cs_ptr, *cs_ptre, *ptr;
  boolean dense_ip, dense_op, got, leftpad, rightpad, keeptable;
  setToLabelsType *label, l1, l2, *newlabel;
  gen **labellist1, **labellist2;
  short_hash_table ht, labelht;
  fsa *genmult2ptr;
  srec *labelset;
  FILE *tablefile, *fopen();

  if (kbm_print_level>=3)
    printf("    #Calling fsa_genmult2_short.\n");
  if (!genmultptr->flags[DFA]){
    fprintf(stderr,"Error: fsa_genmult2 only applies to DFA's.\n");
    return 0;
  }

  if (genmultptr->alphabet->type!=PRODUCT || genmultptr->alphabet->arity!=2) {
    fprintf(stderr,"Error in fsa_genmult2: fsa must be 2-variable.\n");
    return 0;
  }

  if (genmultptr->states->type!=LABELED) {
    fprintf(stderr,"Error in fsa_genmult2: states of fsa must be labeled.\n");
    return 0;
  }

  keeptable = stringlen(genmult2filename)>0;
  if (readback && !keeptable) {
   fprintf(stderr,"Error in fsa_genmult2: readback true but empty filename.\n");
    return 0;
  }

  tmalloc(genmult2ptr,fsa,1);
  fsa_init(genmult2ptr);
  srec_copy(genmult2ptr->alphabet,genmultptr->alphabet);
  genmult2ptr->num_initial = 1;
  tmalloc(genmult2ptr->initial,int,2);
  genmult2ptr->initial[1] = 1;
  genmult2ptr->flags[DFA] = TRUE;
  genmult2ptr->flags[ACCESSIBLE] = TRUE;
  genmult2ptr->flags[BFS] = TRUE;

  genmult2ptr->table->table_type = op_table_type;
  genmult2ptr->table->denserows = 0;
  genmult2ptr->table->printing_format = op_table_type;
  if (!readback) {
    tmalloc(genmult2ptr->table->filename,char,stringlen(genmult2filename)+1);
    strcpy(genmult2ptr->table->filename,genmult2filename);
  }

  ne = genmultptr->alphabet->size;
  ngens = genmultptr->alphabet->base->size;
  ngens1 = ngens+1;

  if (ne != ngens1*ngens1-1) {
   fprintf(stderr,
     "Error: in a 2-variable fsa, alphabet size should = (ngens+1)^2 - 1.\n");
    return 0;
  }

  fsa_set_is_accepting(genmultptr);

  dense_ip = genmultptr->table->table_type==DENSE;
  dr = genmultptr->table->denserows;
  dense_op = op_table_type==DENSE;
  table = genmultptr->table->table_data_ptr;

  short_hash_init(&ht,FALSE,0,0,0);
  ht_ptr = ht.current_ptr;
  ht_ptr[0] = genmultptr->initial[1];
  ht_ptr[1] = genmultptr->initial[1];
  im = short_hash_locate(&ht,2);
/* Each state in the composite transition table will be represented as a subset
 * of the set of ordered pairs of states of *genmultptr.
 * The initial state contains just one such pair, of which both components are
 * the initial state of *genmultptr.
 * The subsets will be stored as variable-length records in the hash-table,
 * as a list of pairs in increasing order.
 * If a state is reached from a transition ($,x) or (x,$) (with $ the padding
 * symbol), then it needs to be marked as such, since we do not allow a $
 * to be followed by any other generator.
 * We do this by adding a 1 or a 2 to the end of the statelist - this is
 * recognisable by the fact that the length of the statelist then becomes odd.
 */
  if (im!=1) {
    fprintf(stderr,"Hash-initialisation problem in fsa_genmult2.\n");
    return 0;
  }
  if (keeptable)
    if ((tablefile=fopen(genmult2filename,"w"))==0){
      fprintf(stderr,"Error: cannot open file %s\n",genmult2filename);
      return 0;
    }
  if (dense_op)
    tmalloc(fsarow,int,ne)
  else
    tmalloc(fsarow,int,2*ne+1)
 
  label = genmultptr->states->setToLabels;
  cstate = 0;
  if (dense_op)
    len = ne; /* The length of the fsarow output. */
  nt = 0; /* Number of transitions in genmult2ptr */

  while (++cstate <= ht.num_recs) {
    if (kbm_print_level>=3) {
      if ((cstate<=1000 && cstate%100==0)||(cstate<=10000 && cstate%1000==0)||
          (cstate<=100000 && cstate%5000==0) || cstate%50000==0)
       printf("    #cstate = %d;  number of states = %d.\n",cstate,ht.num_recs);
    }
    reclen = short_hash_rec_len(&ht,cstate);
    cs_ptr = short_hash_rec(&ht,cstate);
    cs_ptre = short_hash_rec(&ht,cstate) + reclen - 1;
    if (reclen % 2 == 1) {
      if (*cs_ptre==1) {
        leftpad=TRUE;  rightpad=FALSE;
      } else {
        leftpad=FALSE;  rightpad=TRUE;
      }
      cs_ptre--;
    }
    else {
      leftpad=FALSE; rightpad=FALSE;
    }
      
    if (dense_op)
      for (i=0;i<ne;i++) fsarow[i] = 0;
    else
      len = 0;

    e = 0; /* e is the edge number of generator pair (g1,g3) */
    for (g1=1;g1<=ngens1;g1++)  for (g3=1;g3<=ngens1;g3++) {
/* Calculate action of pair (g1,g3) on state cstate  - to get the image, we
 * have to apply ( (g1,g2), (g2,g3) ) to each ordered pair in the subset
 * corresponding to cstate, * and this for each generator g2 of the
 * base-alphabet (including the padding symbol).
 */
       
      e++;
      if (g1==ngens1 && g3==ngens1)
        continue;

      if ((leftpad && g1<=ngens) || (rightpad && g3<=ngens))
        continue;
      ht_ptrb = ht.current_ptr;
      ht_ptre = ht_ptrb-1;

      es1 = (g1-1)*ngens1 + 1;
      ef1 = g1*ngens1;
/* As g2 ranges from 1 to ngens+1 in the pair (g1,g2), for fixed g1, the
 * corresponding edge number in the fsa ranges from es1 to ef1.
 *
 * As g2 ranges from 1 to ngens+1 in the pair (g2,g3), for fixed g3, the
 * corresponding edge number in the fsa ranges from g3 upwards, with
 * increments of size ngens+1.
 */

      ptr = cs_ptr;
      while (ptr <= cs_ptre) {
        csa = *(ptr++); csb = *(ptr++);
        for (g2=1,e1=es1,e2=g3; e1<=ef1; g2++,e1++,e2+=ngens1){
          csima = g1==ngens1 && g2==ngens1 ?
		   (genmultptr->is_accepting[csa] ? csa : 0) :
                   target(dense_ip,table,e1,csa,dr);
          if (csima==0)
            continue;

          csimb = g2==ngens1 && g3==ngens1 ?
                   (genmultptr->is_accepting[csb] ? csb : 0) :
                   target(dense_ip,table,e2,csb,dr);
          if (csimb==0)
            continue;
          if (ht_ptrb>ht_ptre || csima> *(ht_ptre-1) ||
                               (csima== *(ht_ptre-1) &&  csimb > *ht_ptre) ){
/* We have a new state for the image subset to be added to the end */
            *(++ht_ptre) = csima; *(++ht_ptre) = csimb;
          }
          else {
            ht_chptr = ht_ptrb;
            while (*ht_chptr < csima)
              ht_chptr+=2;
            while (*ht_chptr==csima && *(ht_chptr+1)<csimb)
                ht_chptr+=2;
            if (csima < *ht_chptr || csimb < *(ht_chptr+1)) {
/* we have a new state for the image subset to be added in the middle */
              ht_ptr = ht_ptre;
              ht_ptre += 2;
              while (ht_ptr >= ht_chptr) {
                *(ht_ptr+2) = *ht_ptr;
                ht_ptr--;
              }
              *ht_chptr = csima; *(ht_chptr+1) = csimb;
            }
          }
        }
      }
      if (ht_ptre > ht_ptrb) {
        if (g1==ngens1)
          *(++ht_ptre) = 1;
        else if (g3==ngens1)
          *(++ht_ptre) = 2;
      }
      im = short_hash_locate(&ht,ht_ptre-ht_ptrb+1);
      if (im > 0) {
        if (dense_op)
          fsarow[e-1] = im;
        else  {
         fsarow[++len] = e;
         fsarow[++len] = im;
        }
        nt++;
      }
    }
    if (!dense_op)
      fsarow[0] = len++;
    if (keeptable)
      fwrite((void *)fsarow,sizeof(int),(size_t)len,tablefile);
  }

  if (keeptable)
    fclose(tablefile);
  tfree(fsarow);

  ns =
  genmult2ptr->states->size = ht.num_recs;
  genmult2ptr->table->numTransitions = nt;

  if (kbm_print_level>=3) {
    printf("    #Calculated transitions - %d states, %d transitions.\n",ns,nt);
    printf("    #Now calculating state labels.\n");
  }

/* Locate the accept states for  Mult_(a*b)  each generator pair (a,b).
 * This is slightly cumbersome, since a state S
 * is an accept state if either the  subset representing S contains a
 * pair (s1,s2), where s1 is accept state for Mult_a and s2 for Mult_b,
 * or we can get from to such a state by applying ( ($,g), (g,$) ) to one
 * of the pairs in S, for some generator g.
 * It is most convenient to work out this information taking the states
 * S in reverse order.
 * The information on the accept-states will be stored as labels, which
 * (13/9/95 - change labels from words to lists of words)
 * will be lists of words in the generators. Each such list will have the form
 * [a1*b1,a2*b2,...,ar*br], where the (ai,bi) are the generator pairs for
 * which that particular state is an accept state. The labels will be
 * collected initially in a new hash-table.
 * The identity generator will be numbered ngens+1 and given name '_'
 * rather than represented by the empty word. This is so that we can
 * distinguish between multipliers for "_*a" and "a*_" if necessary.
 */

  genmult2ptr->states->type = LABELED;
  tmalloc(genmult2ptr->states->labels,srec,1);
  labelset = genmult2ptr->states->labels;
  labelset->type = LISTOFWORDS;
  labelset->alphabet_size = ngens1;
  for (i=1;i<=ngens;i++) {
    tmalloc(labelset->alphabet[i],char,
                        stringlen(genmultptr->states->labels->alphabet[i])+1);
    strcpy(labelset->alphabet[i],genmultptr->states->labels->alphabet[i]);
  }
  tmalloc(labelset->alphabet[ngens1],char,2);
  labelset->alphabet[ngens1][0]='_';
  labelset->alphabet[ngens1][1]=0;
  tmalloc(genmult2ptr->states->setToLabels,setToLabelsType,ns+1);
  newlabel = genmult2ptr->states->setToLabels;
  tmalloc(genmult2ptr->is_accepting,boolean,ns+1);
  for (i=1;i<=ns;i++)
    genmult2ptr->is_accepting[i] = FALSE;
  short_hash_init(&labelht,FALSE,0,LABHTSIZE,2*LABHTSIZE);

  es1 = ngens*ngens1 + 1;
  ef1 = ngens1*ngens1-1;
  for (cstate=ns;cstate>=1;cstate--) {
/* We work backwards through the states, since we wish to add on additional
 * elements at the end of the list in the hash-table - this destroys
 * later entries, but that doesn't matter, since we have already dealt
 * with them.
 */
    cs_ptr = short_hash_rec(&ht,cstate);
    reclen = short_hash_rec_len(&ht,cstate);
    if (reclen % 2 ==1)
      reclen--;
       /* The last entry only marks the fact that this is a "padded state"*/
    cs_ptre = short_hash_rec(&ht,cstate) + reclen - 1;
/* Apply generators ( ($,g2), (g2,$) ) and see if we get anything new.
 * We won't bother about having them in increasing order this time.
 */
    ptr = cs_ptr;
    while (ptr <= cs_ptre) {
      csa = *(ptr++); csb = *(ptr++);
      for (e1=es1,e2=ngens1; e1<=ef1; e1++,e2+=ngens1){
        csima = target(dense_ip,table,e1,csa,dr);
        if (csima==0)
          continue;
        csimb = target(dense_ip,table,e2,csb,dr);
        if (csimb==0)
          continue;
  /* see if (csima,csimb) is new */
        ht_chptr = cs_ptr;
        got = FALSE;
        while (ht_chptr < cs_ptre) {
          if (csima == ht_chptr[0] && csimb == ht_chptr[1]) {
            got = TRUE;
            break;
          }
          ht_chptr+=2;
        }
        if (!got) {
  /* add (csima,csimb) to the end */
          *(++cs_ptre) = csima; *(++cs_ptre) = csimb;
        }
      }
    }
  /* Now we see which pairs in the subset are of form (s,t), where s is
   * an accept state for a generator a, and t for a generator b.
   * The list of all such pairs (a,b) will form the label of the state, which
   * will be the list of words [a1*b1,a2*b2,..,ar*br], with the (ai,bi) coming
   * in lex-order.
   */

    ht_ptrb = labelht.current_ptr;
    ht_ptre = ht_ptrb-1;
    ptr = cs_ptr;
    while (ptr <= cs_ptre) {
      csa = *(ptr++); csb = *(ptr++);
      if (((l1=label[csa]) != 0) && ((l2=label[csb]) != 0)) {
        labellist1 = genmultptr->states->labels->wordslist[l1];
        labellist2 = genmultptr->states->labels->wordslist[l2];
        j1=0;
        while (labellist1[j1]) {
          g1 = labellist1[j1][0];
          if (g1==0) g1=ngens1;
          j2=0;
          while (labellist2[j2]) {
            genmult2ptr->is_accepting[cstate]=TRUE;
            g2 = labellist2[j2][0];
            if (g2==0) g2=ngens1;
            if (ht_ptrb>ht_ptre || g1> *(ht_ptre-1) ||
                                 (g1== *(ht_ptre-1) &&  g2 > *ht_ptre) ){
    /* We have a new generator pair to be added to the end */
              *(++ht_ptre) = g1; *(++ht_ptre) = g2;
            }
            else {
              ht_chptr = ht_ptrb;
              while (*ht_chptr < g1)
                ht_chptr+=2;
              while (*ht_chptr==g1 && *(ht_chptr+1)<g2)
                  ht_chptr+=2;
              if (g1 < *ht_chptr || g2 < *(ht_chptr+1)) {
    /* we have a new generator pair to be added in the middle */
                ht_ptr = ht_ptre;
                ht_ptre += 2;
                while (ht_ptr >= ht_chptr) {
                  *(ht_ptr+2) = *ht_ptr;
                  ht_ptr--;
                }
                *ht_chptr = g1; *(ht_chptr+1) = g2;
              }
            }
            j2++;
          }
          j1++;
        }
      }
    }
/* That completes the calculation of the label for cstate */
    newlabel[cstate] = short_hash_locate(&labelht,ht_ptre-ht_ptrb+1);
  } 
  short_hash_clear(&ht);

  ct=0;
  for (i=1;i<=ns;i++) if (genmult2ptr->is_accepting[i])
    ct++;
  genmult2ptr->num_accepting = ct;
  if (ct==1 || ct != ns) {
    tmalloc(genmult2ptr->accepting,int,ct+1);
    ct = 0;
    for (i=1;i<=ns;i++) if (genmult2ptr->is_accepting[i])
      genmult2ptr->accepting[++ct] = i;
  }
  tfree(genmult2ptr->is_accepting);
  tfree(genmultptr->is_accepting);

/* Finally copy the records from the label hash-table into the set of labels */
  nlab =
  labelset->size = labelht.num_recs;
  if (kbm_print_level>=3)
    printf("    #There are %d distinct labels.\n",nlab);
  tmalloc(labelset->wordslist,gen **,nlab+1);
  for (i=1;i<=nlab;i++) {
    len = short_hash_rec_len(&labelht,i)/2;
    tmalloc(labelset->wordslist[i],gen *,len+1);
    ht_ptr = short_hash_rec(&labelht,i);
    for (j=0;j<len;j++) {
      tmalloc(labelset->wordslist[i][j],gen,3);
      labelset->wordslist[i][j][0] = ht_ptr[2*j];
      labelset->wordslist[i][j][1] = ht_ptr[2*j+1];
      labelset->wordslist[i][j][2] = 0;
    }
    labelset->wordslist[i][len] = 0;
  }

  short_hash_clear(&labelht);
  if (destroy)
    fsa_clear(genmultptr);

  if (readback) {
    tablefile = fopen(genmult2filename,"r");
    compressed_transitions_read(genmult2ptr,tablefile);
    fclose(tablefile);
    unlink(genmult2filename);
  }

  return genmult2ptr;
}

fsa *
fsa_genmult2_int(genmultptr,op_table_type,destroy,genmult2filename,readback)
	fsa *genmultptr;
	storage_type op_table_type;
	boolean destroy;
	char *genmult2filename;
	boolean readback;
{
  int **table, ne, ngens, ngens1, ns, *fsarow, e, e1, e2, es1, ef1, dr,
      nt, cstate, im, csa, csb, csima, csimb,
      i, j, j1, j2, g1, g2, g3, len, reclen, nlab, ct;
  int *ht_ptr, *ht_chptr, *ht_ptrb, *ht_ptre, *cs_ptr, *cs_ptre, *ptr;
  boolean dense_ip, dense_op, got, leftpad, rightpad, keeptable;
  setToLabelsType *label, l1, l2, *newlabel;
  gen **labellist1, **labellist2;
  hash_table ht, labelht;
  fsa *genmult2ptr;
  srec *labelset;
  FILE *tablefile, *fopen();

  if (kbm_print_level>=3)
    printf("    #Calling fsa_genmult2_int.\n");
  if (!genmultptr->flags[DFA]){
    fprintf(stderr,"Error: fsa_genmult2 only applies to DFA's.\n");
    return 0;
  }

  if (genmultptr->alphabet->type!=PRODUCT || genmultptr->alphabet->arity!=2) {
    fprintf(stderr,"Error in fsa_genmult2: fsa must be 2-variable.\n");
    return 0;
  }

  if (genmultptr->states->type!=LABELED) {
    fprintf(stderr,"Error in fsa_genmult2: states of fsa must be labeled.\n");
    return 0;
  }

  keeptable = stringlen(genmult2filename)>0;
  if (readback && !keeptable) {
   fprintf(stderr,"Error in fsa_genmult2: readback true but empty filename.\n");
    return 0;
  }

  tmalloc(genmult2ptr,fsa,1);
  fsa_init(genmult2ptr);
  srec_copy(genmult2ptr->alphabet,genmultptr->alphabet);
  genmult2ptr->num_initial = 1;
  tmalloc(genmult2ptr->initial,int,2);
  genmult2ptr->initial[1] = 1;
  genmult2ptr->flags[DFA] = TRUE;
  genmult2ptr->flags[ACCESSIBLE] = TRUE;
  genmult2ptr->flags[BFS] = TRUE;

  genmult2ptr->table->table_type = op_table_type;
  genmult2ptr->table->denserows = 0;
  genmult2ptr->table->printing_format = op_table_type;
  if (!readback) {
    tmalloc(genmult2ptr->table->filename,char,stringlen(genmult2filename)+1);
    strcpy(genmult2ptr->table->filename,genmult2filename);
  }

  ne = genmultptr->alphabet->size;
  ngens = genmultptr->alphabet->base->size;
  ngens1 = ngens+1;

  if (ne != ngens1*ngens1-1) {
   fprintf(stderr,
     "Error: in a 2-variable fsa, alphabet size should = (ngens+1)^2 - 1.\n");
    return 0;
  }

  fsa_set_is_accepting(genmultptr);

  dense_ip = genmultptr->table->table_type==DENSE;
  dr = genmultptr->table->denserows;
  dense_op = op_table_type==DENSE;
  table = genmultptr->table->table_data_ptr;

  hash_init(&ht,FALSE,0,0,0);
  ht_ptr = ht.current_ptr;
  ht_ptr[0] = genmultptr->initial[1];
  ht_ptr[1] = genmultptr->initial[1];
  im = hash_locate(&ht,2);
/* Each state in the composite transition table will be represented as a subset
 * of the set of ordered pairs of states of *genmultptr.
 * The initial state contains just one such pair, of which both components are
 * the initial state of *genmultptr.
 * The subsets will be stored as variable-length records in the hash-table,
 * as a list of pairs in increasing order.
 * If a state is reached from a transition ($,x) or (x,$) (with $ the padding
 * symbol), then it needs to be marked as such, since we do not allow a $
 * to be followed by any other generator.
 * We do this by adding a 1 or a 2 to the end of the statelist - this is
 * recognisable by the fact that the length of the statelist then becomes odd.
 */
  if (im!=1) {
    fprintf(stderr,"Hash-initialisation problem in fsa_genmult2.\n");
    return 0;
  }
  if (keeptable)
    if ((tablefile=fopen(genmult2filename,"w"))==0){
      fprintf(stderr,"Error: cannot open file %s\n",genmult2filename);
      return 0;
    }
  if (dense_op)
    tmalloc(fsarow,int,ne)
  else
    tmalloc(fsarow,int,2*ne+1)
 
  label = genmultptr->states->setToLabels;
  cstate = 0;
  if (dense_op)
    len = ne; /* The length of the fsarow output. */
  nt = 0; /* Number of transitions in genmult2ptr */

  while (++cstate <= ht.num_recs) {
    if (kbm_print_level>=3) {
      if ((cstate<=1000 && cstate%100==0)||(cstate<=10000 && cstate%1000==0)||
          (cstate<=100000 && cstate%5000==0) || cstate%50000==0)
       printf("    #cstate = %d;  number of states = %d.\n",cstate,ht.num_recs);
    }
    reclen = hash_rec_len(&ht,cstate);
    cs_ptr = hash_rec(&ht,cstate);
    cs_ptre = hash_rec(&ht,cstate) + reclen - 1;
    if (reclen % 2 == 1) {
      if (*cs_ptre==1) {
        leftpad=TRUE;  rightpad=FALSE;
      } else {
        leftpad=FALSE;  rightpad=TRUE;
      }
      cs_ptre--;
    }
    else {
      leftpad=FALSE; rightpad=FALSE;
    }
      
    if (dense_op)
      for (i=0;i<ne;i++) fsarow[i] = 0;
    else
      len = 0;

    e = 0; /* e is the edge number of generator pair (g1,g3) */
    for (g1=1;g1<=ngens1;g1++)  for (g3=1;g3<=ngens1;g3++) {
/* Calculate action of pair (g1,g3) on state cstate  - to get the image, we
 * have to apply ( (g1,g2), (g2,g3) ) to each ordered pair in the subset
 * corresponding to cstate, * and this for each generator g2 of the
 * base-alphabet (including the padding symbol).
 */
       
      e++;
      if (g1==ngens1 && g3==ngens1)
        continue;

      if ((leftpad && g1<=ngens) || (rightpad && g3<=ngens))
        continue;
      ht_ptrb = ht.current_ptr;
      ht_ptre = ht_ptrb-1;

      es1 = (g1-1)*ngens1 + 1;
      ef1 = g1*ngens1;
/* As g2 ranges from 1 to ngens+1 in the pair (g1,g2), for fixed g1, the
 * corresponding edge number in the fsa ranges from es1 to ef1.
 *
 * As g2 ranges from 1 to ngens+1 in the pair (g2,g3), for fixed g3, the
 * corresponding edge number in the fsa ranges from g3 upwards, with
 * increments of size ngens+1.
 */

      ptr = cs_ptr;
      while (ptr <= cs_ptre) {
        csa = *(ptr++); csb = *(ptr++);
        for (g2=1,e1=es1,e2=g3; e1<=ef1; g2++,e1++,e2+=ngens1){
          csima = g1==ngens1 && g2==ngens1 ?
		   (genmultptr->is_accepting[csa] ? csa : 0) :
                   target(dense_ip,table,e1,csa,dr);
          if (csima==0)
            continue;

          csimb = g2==ngens1 && g3==ngens1 ?
                   (genmultptr->is_accepting[csb] ? csb : 0) :
                   target(dense_ip,table,e2,csb,dr);
          if (csimb==0)
            continue;
          if (ht_ptrb>ht_ptre || csima> *(ht_ptre-1) ||
                               (csima== *(ht_ptre-1) &&  csimb > *ht_ptre) ){
/* We have a new state for the image subset to be added to the end */
            *(++ht_ptre) = csima; *(++ht_ptre) = csimb;
          }
          else {
            ht_chptr = ht_ptrb;
            while (*ht_chptr < csima)
              ht_chptr+=2;
            while (*ht_chptr==csima && *(ht_chptr+1)<csimb)
                ht_chptr+=2;
            if (csima < *ht_chptr || csimb < *(ht_chptr+1)) {
/* we have a new state for the image subset to be added in the middle */
              ht_ptr = ht_ptre;
              ht_ptre += 2;
              while (ht_ptr >= ht_chptr) {
                *(ht_ptr+2) = *ht_ptr;
                ht_ptr--;
              }
              *ht_chptr = csima; *(ht_chptr+1) = csimb;
            }
          }
        }
      }
      if (ht_ptre > ht_ptrb) {
        if (g1==ngens1)
          *(++ht_ptre) = 1;
        else if (g3==ngens1)
          *(++ht_ptre) = 2;
      }
      im = hash_locate(&ht,ht_ptre-ht_ptrb+1);
      if (im > 0) {
        if (dense_op)
          fsarow[e-1] = im;
        else  {
         fsarow[++len] = e;
         fsarow[++len] = im;
        }
        nt++;
      }
    }
    if (!dense_op)
      fsarow[0] = len++;
    if (keeptable)
      fwrite((void *)fsarow,sizeof(int),(size_t)len,tablefile);
  }

  if (keeptable)
    fclose(tablefile);
  tfree(fsarow);

  ns =
  genmult2ptr->states->size = ht.num_recs;
  genmult2ptr->table->numTransitions = nt;

  if (kbm_print_level>=3) {
    printf("    #Calculated transitions - %d states, %d transitions.\n",ns,nt);
    printf("    #Now calculating state labels.\n");
  }

/* Locate the accept states for  Mult_(a*b)  each generator pair (a,b).
 * This is slightly cumbersome, since a state S
 * is an accept state if either the  subset representing S contains a
 * pair (s1,s2), where s1 is accept state for Mult_a and s2 for Mult_b,
 * or we can get from to such a state by applying ( ($,g), (g,$) ) to one
 * of the pairs in S, for some generator g.
 * It is most convenient to work out this information taking the states
 * S in reverse order.
 * The information on the accept-states will be stored as labels, which
 * (13/9/95 - change labels from words to lists of words)
 * will be lists of words in the generators. Each such list will have the form
 * [a1*b1,a2*b2,...,ar*br], where the (ai,bi) are the generator pairs for
 * which that particular state is an accept state. The labels will be
 * collected initially in a new hash-table.
 * The identity generator will be numbered ngens+1 and given name '_'
 * rather than represented by the empty word. This is so that we can
 * distinguish between multipliers for "_*a" and "a*_" if necessary.
 */

  genmult2ptr->states->type = LABELED;
  tmalloc(genmult2ptr->states->labels,srec,1);
  labelset = genmult2ptr->states->labels;
  labelset->type = LISTOFWORDS;
  labelset->alphabet_size = ngens1;
  for (i=1;i<=ngens;i++) {
    tmalloc(labelset->alphabet[i],char,
                        stringlen(genmultptr->states->labels->alphabet[i])+1);
    strcpy(labelset->alphabet[i],genmultptr->states->labels->alphabet[i]);
  }
  tmalloc(labelset->alphabet[ngens1],char,2);
  labelset->alphabet[ngens1][0]='_';
  labelset->alphabet[ngens1][1]=0;
  tmalloc(genmult2ptr->states->setToLabels,setToLabelsType,ns+1);
  newlabel = genmult2ptr->states->setToLabels;
  tmalloc(genmult2ptr->is_accepting,boolean,ns+1);
  for (i=1;i<=ns;i++)
    genmult2ptr->is_accepting[i] = FALSE;
  hash_init(&labelht,FALSE,0,LABHTSIZE,2*LABHTSIZE);

  es1 = ngens*ngens1 + 1;
  ef1 = ngens1*ngens1-1;
  for (cstate=ns;cstate>=1;cstate--) {
/* We work backwards through the states, since we wish to add on additional
 * elements at the end of the list in the hash-table - this destroys
 * later entries, but that doesn't matter, since we have already dealt
 * with them.
 */
    cs_ptr = hash_rec(&ht,cstate);
    reclen = hash_rec_len(&ht,cstate);
    if (reclen % 2 ==1)
      reclen--;
       /* The last entry only marks the fact that this is a "padded state"*/
    cs_ptre = hash_rec(&ht,cstate) + reclen - 1;
/* Apply generators ( ($,g2), (g2,$) ) and see if we get anything new.
 * We won't bother about having them in increasing order this time.
 */
    ptr = cs_ptr;
    while (ptr <= cs_ptre) {
      csa = *(ptr++); csb = *(ptr++);
      for (e1=es1,e2=ngens1; e1<=ef1; e1++,e2+=ngens1){
        csima = target(dense_ip,table,e1,csa,dr);
        if (csima==0)
          continue;
        csimb = target(dense_ip,table,e2,csb,dr);
        if (csimb==0)
          continue;
  /* see if (csima,csimb) is new */
        ht_chptr = cs_ptr;
        got = FALSE;
        while (ht_chptr < cs_ptre) {
          if (csima == ht_chptr[0] && csimb == ht_chptr[1]) {
            got = TRUE;
            break;
          }
          ht_chptr+=2;
        }
        if (!got) {
  /* add (csima,csimb) to the end */
          *(++cs_ptre) = csima; *(++cs_ptre) = csimb;
        }
      }
    }
  /* Now we see which pairs in the subset are of form (s,t), where s is
   * an accept state for a generator a, and t for a generator b.
   * The list of all such pairs (a,b) will form the label of the state, which
   * will be the list of words [a1*b1,a2*b2,..,ar*br], with the (ai,bi) coming
   * in lex-order.
   */

    ht_ptrb = labelht.current_ptr;
    ht_ptre = ht_ptrb-1;
    ptr = cs_ptr;
    while (ptr <= cs_ptre) {
      csa = *(ptr++); csb = *(ptr++);
      if (((l1=label[csa]) != 0) && ((l2=label[csb]) != 0)) {
        labellist1 = genmultptr->states->labels->wordslist[l1];
        labellist2 = genmultptr->states->labels->wordslist[l2];
        j1=0;
        while (labellist1[j1]) {
          g1 = labellist1[j1][0];
          if (g1==0) g1=ngens1;
          j2=0;
          while (labellist2[j2]) {
            genmult2ptr->is_accepting[cstate]=TRUE;
            g2 = labellist2[j2][0];
            if (g2==0) g2=ngens1;
            if (ht_ptrb>ht_ptre || g1> *(ht_ptre-1) ||
                                 (g1== *(ht_ptre-1) &&  g2 > *ht_ptre) ){
    /* We have a new generator pair to be added to the end */
              *(++ht_ptre) = g1; *(++ht_ptre) = g2;
            }
            else {
              ht_chptr = ht_ptrb;
              while (*ht_chptr < g1)
                ht_chptr+=2;
              while (*ht_chptr==g1 && *(ht_chptr+1)<g2)
                  ht_chptr+=2;
              if (g1 < *ht_chptr || g2 < *(ht_chptr+1)) {
    /* we have a new generator pair to be added in the middle */
                ht_ptr = ht_ptre;
                ht_ptre += 2;
                while (ht_ptr >= ht_chptr) {
                  *(ht_ptr+2) = *ht_ptr;
                  ht_ptr--;
                }
                *ht_chptr = g1; *(ht_chptr+1) = g2;
              }
            }
            j2++;
          }
          j1++;
        }
      }
    }
/* That completes the calculation of the label for cstate */
    newlabel[cstate] = hash_locate(&labelht,ht_ptre-ht_ptrb+1);
  } 
  hash_clear(&ht);

  ct=0;
  for (i=1;i<=ns;i++) if (genmult2ptr->is_accepting[i])
    ct++;
  genmult2ptr->num_accepting = ct;
  if (ct==1 || ct != ns) {
    tmalloc(genmult2ptr->accepting,int,ct+1);
    ct = 0;
    for (i=1;i<=ns;i++) if (genmult2ptr->is_accepting[i])
      genmult2ptr->accepting[++ct] = i;
  }
  tfree(genmult2ptr->is_accepting);
  tfree(genmultptr->is_accepting);

/* Finally copy the records from the label hash-table into the set of labels */
  nlab =
  labelset->size = labelht.num_recs;
  if (kbm_print_level>=3)
    printf("    #There are %d distinct labels.\n",nlab);
  tmalloc(labelset->wordslist,gen **,nlab+1);
  for (i=1;i<=nlab;i++) {
    len = hash_rec_len(&labelht,i)/2;
    tmalloc(labelset->wordslist[i],gen *,len+1);
    ht_ptr = hash_rec(&labelht,i);
    for (j=0;j<len;j++) {
      tmalloc(labelset->wordslist[i][j],gen,3);
      labelset->wordslist[i][j][0] = ht_ptr[2*j];
      labelset->wordslist[i][j][1] = ht_ptr[2*j+1];
      labelset->wordslist[i][j][2] = 0;
    }
    labelset->wordslist[i][len] = 0;
  }

  hash_clear(&labelht);
  if (destroy)
    fsa_clear(genmultptr);

  if (readback) {
    tablefile = fopen(genmult2filename,"r");
    compressed_transitions_read(genmult2ptr,tablefile);
    fclose(tablefile);
    unlink(genmult2filename);
  }

  return genmult2ptr;
}

int
fsa_makemult(genmultptr,g)
        fsa *genmultptr;
        int g;
/* This procedure takes the fsa *genmultptr produced by fsa_triples, 
 * and builds a particular  multiplier Mult_g1.
 * This merely involves setting the accept states of *genmultptr
 * in accordance with the labels of its states.
 * g can be 0, for the identity generator, which (in shortlex case) inevitably
 * produces the diagonal of the word-acceptor.
 * This procedure alters its arguments and does not return anything.
 */
{
  int ngens, ns, i, j, ct;
  gen ** labellist;
  setToLabelsType *label;

  if (kbm_print_level>=3)
    printf("    #Calling fsa_makemult with generator number %d.\n",g);
  if (!genmultptr->flags[DFA]){
    fprintf(stderr,"Error: fsa_makemult only applies to DFA's.\n");
    return -1;
  }

  if (genmultptr->alphabet->type!=PRODUCT || genmultptr->alphabet->arity!=2) {
    fprintf(stderr,"Error in fsa_makemult: fsa must be 2-variable.\n");
    return -1;
  }

  if (genmultptr->states->type!=LABELED) {
    fprintf(stderr,"Error in fsa_makemult: states of fsa must be labeled.\n");
    return -1;
  }

  ns = genmultptr->states->size;
  ngens = genmultptr->alphabet->base->size;
  if (g<0 || g>ngens+1) {
    fprintf(stderr,"#Error in fsa_makemult: Generator is out of range.\n");
    return -1;
  }
  
  tfree(genmultptr->accepting);
  tfree(genmultptr->is_accepting);
  label = genmultptr->states->setToLabels;
  ct = 0;
  for (i=1;i<=ns;i++) if (label[i]>0) {
    labellist = genmultptr->states->labels->wordslist[label[i]];
    j=0;
    while (labellist[j]) {
      if (labellist[j][0]==g)
        ct++;
      j++;
    }
  }
  genmultptr->num_accepting = ct;
  if (ct<ns || ns==1) { 
    tfree(genmultptr->accepting);
    tmalloc(genmultptr->accepting,int,ct+1);
    ct = 0;
    for (i=1;i<=ns;i++)if (label[i]>0) {
      labellist = genmultptr->states->labels->wordslist[label[i]];
      j=0;
      while (labellist[j]) {
        if (labellist[j][0]==g)
          genmultptr->accepting[++ct] = i;
        j++;
      }
    }
  }
/* The state labelling is no longer relevant, so we clear it */
  genmultptr->states->type=SIMPLE;
  srec_clear(genmultptr->states->labels);
  tfree(genmultptr->states->labels);
  tfree(genmultptr->states->setToLabels);
  return 0;
}

int
fsa_makemult2(genmult2ptr,g1,g2)
        fsa *genmult2ptr;
        int g1, g2;
/* This procedure takes the fsa *genmult2ptr produced by fsa_genmult2, 
 * and builds a particular length-2 multiplier Mult_g1g2.
 * This merely involves locating the accept states.
 * This procedure alters its arguments and does not return anything.
 */
{
  int  ngens, ns, nlabs, i, j, ct;
  gen ***labellist;
  boolean *accept;
  setToLabelsType *labelnumber;

  if (kbm_print_level>=3)
    printf("    #Calling fsa_makemult2 with generators %d and %d.\n",g1,g2);
  if (!genmult2ptr->flags[DFA]){
    fprintf(stderr,"Error: fsa_makemult2 only applies to DFA's.\n");
    return -1;
  }

  if (genmult2ptr->alphabet->type!=PRODUCT || genmult2ptr->alphabet->arity!=2) {
    fprintf(stderr,"Error in fsa_makemult2: fsa must be 2-variable.\n");
    return -1;
  }

  if (genmult2ptr->states->type!=LABELED) {
    fprintf(stderr,"Error in fsa_makemult2: states of fsa must be labeled.\n");
    return -1;
  }

  if (genmult2ptr->states->labels->type!=LISTOFWORDS) {
    fprintf(stderr,"Error in fsa_makemult2: labels must be lists of words.\n");
    return -1;
  }

  ns = genmult2ptr->states->size;
  nlabs = genmult2ptr->states->labels->size;
  ngens = genmult2ptr->states->labels->alphabet_size;
  if (g1<=0 || g1>ngens || g2<=0 || g2>ngens) {
    fprintf(stderr,"#Error in fsa_makemult2: A generator is out of range.\n");
    return -1;
  }
  
  tmalloc(accept,boolean,nlabs+1);
  labellist = genmult2ptr->states->labels->wordslist;
/* Now we see which labels are accept-labels for the pair (g1,g2) - the
 * label is an accept-label if the list of words which is its name
 * contains g1*g2.
 */
  accept[0] = FALSE;
  for (i=1;i<=nlabs;i++) {
    accept[i] = FALSE;
    j=0;
    while (labellist[i][j]) {
      if (labellist[i][j][0]==g1 && labellist[i][j][1]==g2) {
        accept[i] = TRUE;
        break;
      }
      j++;
    }
  }

/* Now we can see which states are accept-states. A state is an accept-state
 * iff its label is an accept-label.
 * First we count the number.
 */
  ct = 0;
  labelnumber = genmult2ptr->states->setToLabels;
  for (i=1;i<=ns;i++) if (accept[labelnumber[i]])
    ct++;

  genmult2ptr->num_accepting = ct;
  if (ct<ns || ns==1) { 
    tfree(genmult2ptr->accepting);
    tmalloc(genmult2ptr->accepting,int,ct+1);
    ct = 0;
    for (i=1;i<=ns;i++) if (accept[labelnumber[i]])
      genmult2ptr->accepting[++ct] = i;
  }

  tfree(accept);

/* The state labelling is no longer relevant, so we clear it */
  genmult2ptr->states->type=SIMPLE;
  srec_clear(genmult2ptr->states->labels);
  tfree(genmult2ptr->states->labels);
  tfree(genmult2ptr->states->setToLabels);
  return 0;
}

fsa *
fsa_composite(mult1ptr,mult2ptr,op_table_type,destroy,compfilename,readback)
	fsa *mult1ptr, *mult2ptr;
	storage_type op_table_type;
	boolean destroy;
	char *compfilename;
        boolean readback;
/* *mult1ptr and *mult2ptr should be multiplier fsa's of an automatic group.
 * This function calculates the composite of these two multipliers.
 * So if *mult1ptr is the mutiplier for the word w1 and *mult2ptr for w2,
 * then the returned fsa is the multiplier for w1*w2.
 * In fact, *mult1ptr and *mult2ptr can be any 2-variable automata over the
 * same alphabet, and the composite is returned.
 */
{
  if (kbm_print_level>=3)
    printf("    #Calling fsa_composite.\n");
  if (mult1ptr->states->size < MAXUSHORT && mult2ptr->states->size < MAXUSHORT)
    return fsa_composite_short(mult1ptr,mult2ptr,op_table_type,
                                                destroy,compfilename,readback);
  else
    return fsa_composite_int(mult1ptr,mult2ptr,op_table_type,
                                                destroy,compfilename,readback);
}

fsa *
fsa_composite_short(mult1ptr,mult2ptr,op_table_type,
                                                  destroy,compfilename,readback)
	fsa *mult1ptr, *mult2ptr;
	storage_type op_table_type;
	boolean destroy;
	char *compfilename;
        boolean readback;
{
  int **table1, **table2, ne, ngens, ngens1, ns, *fsarow,
      e, e1, e2, es1, ef1, dr1, dr2, nt, cstate, im, csa, csb, csima, csimb,
      i, ct, g1, g2, g3, len, reclen;
  unsigned short *ht_ptr, *ht_chptr, *ht_ptrb, *ht_ptre,
                 *cs_ptr, *cs_ptre, *ptr;
  boolean dense_ip1, dense_ip2, dense_op, got, leftpad, rightpad;
  short_hash_table ht;
  fsa *compositeptr;
  FILE *tempfile, *fopen();

  if (kbm_print_level>=3)
    printf("    #Calling fsa_composite_short.\n");
  if (!mult1ptr->flags[DFA] || !mult2ptr->flags[DFA]){
    fprintf(stderr,"Error: fsa_composite only applies to DFA's.\n");
    return 0;
  }

  if (mult1ptr->alphabet->type!=PRODUCT || mult1ptr->alphabet->arity!=2) {
    fprintf(stderr,"Error in fsa_composite: fsa must be 2-variable.\n");
    return 0;
  }
  if (mult2ptr->alphabet->type!=PRODUCT || mult2ptr->alphabet->arity!=2) {
    fprintf(stderr,"Error in fsa_composite: fsa must be 2-variable.\n");
    return 0;
  }

  if (!srec_equal(mult1ptr->alphabet,mult2ptr->alphabet)) {
    fprintf(stderr,"Error in fsa_composite: fsa's must have same alphabet.\n");
    return 0;
  }

  tmalloc(compositeptr,fsa,1);
  fsa_init(compositeptr);
  srec_copy(compositeptr->alphabet,mult1ptr->alphabet);
  compositeptr->states->type = SIMPLE;
  compositeptr->num_initial = 1;
  tmalloc(compositeptr->initial,int,2);
  compositeptr->initial[1] = 1;
  compositeptr->flags[DFA] = TRUE;
  compositeptr->flags[ACCESSIBLE] = TRUE;
  compositeptr->flags[BFS] = TRUE;

  compositeptr->table->table_type = op_table_type;
  compositeptr->table->denserows = 0;
  compositeptr->table->printing_format = op_table_type;
  if (!readback) {
    tmalloc(compositeptr->table->filename,char,stringlen(compfilename)+1);
    strcpy(compositeptr->table->filename,compfilename);
  }

  ne = mult1ptr->alphabet->size;
  ngens = mult1ptr->alphabet->base->size;
  ngens1 = ngens+1;

  if (ne != ngens1*ngens1-1) {
   fprintf(stderr,
     "Error: in a 2-variable fsa, alphabet size should = (ngens+1)^2 - 1.\n");
    return 0;
  }

  fsa_set_is_accepting(mult1ptr);
  fsa_set_is_accepting(mult2ptr);

  dense_ip1 = mult1ptr->table->table_type==DENSE;
  dr1 = mult1ptr->table->denserows;
  dense_ip2 = mult2ptr->table->table_type==DENSE;
  dr2 = mult2ptr->table->denserows;
  dense_op = op_table_type==DENSE;
  table1 = mult1ptr->table->table_data_ptr;
  table2 = mult2ptr->table->table_data_ptr;

  short_hash_init(&ht,FALSE,0,0,0);
  ht_ptr = ht.current_ptr;
  ht_ptr[0] = mult1ptr->initial[1];
  ht_ptr[1] = mult2ptr->initial[1];
  im = short_hash_locate(&ht,2);
/* Each state in the composite transition table will be represented as a subset
 * of the set of ordered pairs in which the first component is a state in
 * *multptr1 and the second a state in *multptr2.
 * The initial state contains just one such pair, of which the components are
 * the initial states of *mult1ptr and *multptr2.
 * The subsets will be stored as variable-length records in the hash-table,
 * as a list of pairs in increasing order.
 * If a state is reached from a transition ($,x) or (x,$) (with $ the padding
 * symbol), then it needs to be marked as such, since we do not allow a $
 * to be followed by any other generator.
 * We do this by adding a 1 or a 2 to the end of the statelist - this is
 * recognisable by the fact that the length of the statelist then becomes odd.
 */
  if (im!=1) {
    fprintf(stderr,"Hash-initialisation problem in fsa_composite.\n");
    return 0;
  }
  if ((tempfile=fopen(compfilename,"w"))==0){
    fprintf(stderr,"Error: cannot open file %s\n",compfilename);
    return 0;
  }
  if (dense_op)
    tmalloc(fsarow,int,ne)
  else
    tmalloc(fsarow,int,2*ne+1)
 
  cstate = 0;
  if (dense_op)
    len = ne; /* The length of the fsarow output. */
  nt = 0; /* Number of transitions in compositeptr */

  while (++cstate <= ht.num_recs) {
    if (kbm_print_level>=3) {
      if ((cstate<=1000 && cstate%100==0)||(cstate<=10000 && cstate%1000==0)||
          (cstate<=100000 && cstate%5000==0) || cstate%50000==0)
       printf("    #cstate = %d;  number of states = %d.\n",cstate,ht.num_recs);
    }
    reclen = short_hash_rec_len(&ht,cstate);
    cs_ptr = short_hash_rec(&ht,cstate);
    cs_ptre = short_hash_rec(&ht,cstate) + reclen - 1;
    /* for (i=0;i<reclen;i++) printf(" %d",cs_ptr[i]); printf("\n"); */
    if (reclen % 2 == 1) {
      if (*cs_ptre==1) {
        leftpad=TRUE;  rightpad=FALSE;
      } else {
        leftpad=FALSE;  rightpad=TRUE;
      }
      cs_ptre--;
    }
    else {
      leftpad=FALSE; rightpad=FALSE;
    }
      
    if (dense_op)
      for (i=0;i<ne;i++) fsarow[i] = 0;
    else
      len = 0;

    e = 0; /* e is the edge number of generator pair (g1,g3) */
    for (g1=1;g1<=ngens1;g1++)  for (g3=1;g3<=ngens1;g3++) {
/* Calculate action of pair (g1,g3) on state cstate  - to get the image, we
 * have to apply ( (g1,g2), (g2,g3) ) to each ordered pair in the subset
 * corresponding to cstate, and this for each generator g2 of the
 * base-alphabet (including the padding symbol).
 */
      e++;
      if (g1==ngens1 && g3==ngens1)
        continue;
      if ((leftpad && g1<=ngens) || (rightpad && g3<=ngens))
        continue;
      ht_ptrb = ht.current_ptr;
      ht_ptre = ht_ptrb-1;
      es1 = (g1-1)*ngens1 + 1;
      ef1 = g1*ngens1;
/* As g2 ranges from 1 to ngens+1 in the pair (g1,g2), for fixed g1, the
 * corresponding edge number in the fsa ranges from es1 to ef1.
 *
 * As g2 ranges from 1 to ngens+1 in the pair (g2,g3), for fixed g3, the
 * corresponding edge number in the fsa ranges from g3 upwards, with
 * increments of size ngens+1.
 */

      ptr = cs_ptr;
      while (ptr <= cs_ptre) {
        csa = *(ptr++); csb = *(ptr++);
        for (g2=1,e1=es1,e2=g3; e1<=ef1; g2++,e1++,e2+=ngens1){
          csima = g1==ngens1 && g2==ngens1 ?
                   (mult1ptr->is_accepting[csa]>0 ? csa : 0) :
                   target(dense_ip1,table1,e1,csa,dr1);
          if (csima==0)
            continue;

          csimb = g2==ngens1 && g3==ngens1 ?
                   (mult2ptr->is_accepting[csb]>0 ? csb : 0) :
                   target(dense_ip2,table2,e2,csb,dr2);
          if (csimb==0)
            continue;

          if (ht_ptrb>ht_ptre || csima> *(ht_ptre-1) ||
                               (csima== *(ht_ptre-1) &&  csimb > *ht_ptre) ){
/* We have a new state for the image subset to be added to the end */
            *(++ht_ptre) = csima; *(++ht_ptre) = csimb;
          }
          else {
            ht_chptr = ht_ptrb;
            while (*ht_chptr < csima)
              ht_chptr+=2;
            while (*ht_chptr==csima && *(ht_chptr+1)<csimb)
                ht_chptr+=2;
            if (csima < *ht_chptr || csimb < *(ht_chptr+1)) {
/* we have a new state for the image subset to be added in the middle */
              ht_ptr = ht_ptre;
              ht_ptre += 2;
              while (ht_ptr >= ht_chptr) {
                *(ht_ptr+2) = *ht_ptr;
                ht_ptr--;
              }
              *ht_chptr = csima; *(ht_chptr+1) = csimb;
            }
          }
        }
      }
      if (ht_ptre > ht_ptrb) {
        if (g1==ngens1)
          *(++ht_ptre) = 1;
        else if (g3==ngens1)
          *(++ht_ptre) = 2;
      }
      im = short_hash_locate(&ht,ht_ptre-ht_ptrb+1);
      if (im > 0) {
        if (dense_op)
          fsarow[e-1] = im;
        else  {
          fsarow[++len] = e;
          fsarow[++len] = im;
        }
        nt++;
      }
    }
    if (!dense_op)
      fsarow[0] = len++;
    fwrite((void *)fsarow,sizeof(int),(size_t)len,tempfile);
  }
  fclose(tempfile);

  ns = compositeptr->states->size = ht.num_recs;
  compositeptr->table->numTransitions = nt;

  if (kbm_print_level>=3)
    printf("    #Calculated transitions - %d states, %d transitions.\n",ns,nt);

/* Locate the accept states. This is slightly cumbersome, since a state
 * is an accept state if either the corresponding subset contains a
 * a pair (s1,s2), where s1 is and accept state of *mult1ptr and s2 an
 * accept state of *mult2ptr, or we can get from some such state to an
 * accept state by applying elements ( ($,x), (x,$ ).
 * We will need to use the array compositeptr->is_accepting.
 */

  tmalloc(compositeptr->is_accepting,boolean,ns+1);
  for (i=1;i<=ns;i++)
    compositeptr->is_accepting[i] = FALSE;
  ct = 0;
  es1 = ngens*ngens1 + 1;
  ef1 = ngens1*ngens1-1;
  for (cstate=ns;cstate>=1;cstate--) {
/* We work backwards through the states, since we wish to add on additional
 * elements at the end of the list in the hash-table - this destroys
 * later entries, but that doesn't matter, since we have already dealt
 * with them.
 */
    cs_ptr = short_hash_rec(&ht,cstate);
    reclen = short_hash_rec_len(&ht,cstate);
    if (reclen % 2 ==1)
      reclen--; /* The last entry marks the fact that this is a "padded state"*/
    cs_ptre = short_hash_rec(&ht,cstate) + reclen - 1;
/* First see if the set itself contains an accept-state */
    ptr = cs_ptr;
    while (ptr <= cs_ptre) {
      csa = *(ptr++); csb = *(ptr++);
      if (mult1ptr->is_accepting[csa] && mult2ptr->is_accepting[csb]) {
        compositeptr->is_accepting[cstate] = TRUE;
        ct++;
        break;
      }
    }
    if (compositeptr->is_accepting[cstate])
      continue;
/* Next apply generators ( ($,g2), (g2,$) ) and see if we get anything new.
 * We won't bother about having them in increasing order this time.
 */
    ptr = cs_ptr;
    while (ptr <= cs_ptre) {
      csa = *(ptr++); csb = *(ptr++);
      for (e1=es1,e2=ngens1; e1<=ef1; e1++,e2+=ngens1){
        csima = target(dense_ip1,table1,e1,csa,dr1);
        if (csima==0)
          continue;
        csimb = target(dense_ip2,table2,e2,csb,dr2);
        if (csimb==0)
          continue;

  /* see if (csima,csimb) is accepting */
        if (mult1ptr->is_accepting[csima] && mult2ptr->is_accepting[csimb]) {
          compositeptr->is_accepting[cstate] = TRUE;
          ct++;
          break;
        }
  /* now see if it is new */
        ht_chptr = cs_ptr;
        got = FALSE;
        while (ht_chptr < cs_ptre) {
          if (csima == ht_chptr[0] && csimb == ht_chptr[1]) {
            got = TRUE;
            break;
          }
          ht_chptr+=2;
        }
        if (!got) {
  /* add (csima,csimb) to the end */
          *(++cs_ptre) = csima; *(++cs_ptre) = csimb;
        }
      }
      if (compositeptr->is_accepting[cstate])
        continue;
    }
  }

  compositeptr->num_accepting = ct;
  if (ct==1 || ct != ns) {
    tmalloc(compositeptr->accepting,int,ct+1);
    ct = 0;
    for (i=1;i<=ns;i++) if (compositeptr->is_accepting[i])
      compositeptr->accepting[++ct] = i;
  }
  tfree(compositeptr->is_accepting);
  tfree(mult1ptr->is_accepting);
  tfree(mult2ptr->is_accepting);
  short_hash_clear(&ht);
  tfree(fsarow);
  if (destroy) {
    fsa_clear(mult1ptr); fsa_clear(mult2ptr);
  }

  if (readback) {
    tempfile = fopen(compfilename,"r");
    compressed_transitions_read(compositeptr,tempfile);
    fclose(tempfile);
    unlink(compfilename);
  }

  return compositeptr;
}

fsa *
fsa_composite_int(mult1ptr,mult2ptr,op_table_type,destroy,compfilename,readback)
	fsa *mult1ptr, *mult2ptr;
	storage_type op_table_type;
	boolean destroy;
	char *compfilename;
        boolean readback;
{
  int **table1, **table2, ne, ngens, ngens1, ns, *fsarow,
      e, e1, e2, es1, ef1, dr1, dr2, nt, cstate, im, csa, csb, csima, csimb,
      i, ct, g1, g2, g3, len, reclen;
  int *ht_ptr, *ht_chptr, *ht_ptrb, *ht_ptre,
                 *cs_ptr, *cs_ptre, *ptr;
  boolean dense_ip1, dense_ip2, dense_op, got, leftpad, rightpad;
  hash_table ht;
  fsa *compositeptr;
  FILE *tempfile, *fopen();

  if (kbm_print_level>=3)
    printf("    #Calling fsa_composite_int.\n");
  if (!mult1ptr->flags[DFA] || !mult2ptr->flags[DFA]){
    fprintf(stderr,"Error: fsa_composite only applies to DFA's.\n");
    return 0;
  }

  if (mult1ptr->alphabet->type!=PRODUCT || mult1ptr->alphabet->arity!=2) {
    fprintf(stderr,"Error in fsa_composite: fsa must be 2-variable.\n");
    return 0;
  }
  if (mult2ptr->alphabet->type!=PRODUCT || mult2ptr->alphabet->arity!=2) {
    fprintf(stderr,"Error in fsa_composite: fsa must be 2-variable.\n");
    return 0;
  }

  if (!srec_equal(mult1ptr->alphabet,mult2ptr->alphabet)) {
    fprintf(stderr,"Error in fsa_composite: fsa's must have same alphabet.\n");
    return 0;
  }

  tmalloc(compositeptr,fsa,1);
  fsa_init(compositeptr);
  srec_copy(compositeptr->alphabet,mult1ptr->alphabet);
  compositeptr->states->type = SIMPLE;
  compositeptr->num_initial = 1;
  tmalloc(compositeptr->initial,int,2);
  compositeptr->initial[1] = 1;
  compositeptr->flags[DFA] = TRUE;
  compositeptr->flags[ACCESSIBLE] = TRUE;
  compositeptr->flags[BFS] = TRUE;

  compositeptr->table->table_type = op_table_type;
  compositeptr->table->denserows = 0;
  compositeptr->table->printing_format = op_table_type;
  if (!readback) {
    tmalloc(compositeptr->table->filename,char,stringlen(compfilename)+1);
    strcpy(compositeptr->table->filename,compfilename);
  }

  ne = mult1ptr->alphabet->size;
  ngens = mult1ptr->alphabet->base->size;
  ngens1 = ngens+1;

  if (ne != ngens1*ngens1-1) {
   fprintf(stderr,
     "Error: in a 2-variable fsa, alphabet size should = (ngens+1)^2 - 1.\n");
    return 0;
  }

  fsa_set_is_accepting(mult1ptr);
  fsa_set_is_accepting(mult2ptr);

  dense_ip1 = mult1ptr->table->table_type==DENSE;
  dr1 = mult1ptr->table->denserows;
  dense_ip2 = mult2ptr->table->table_type==DENSE;
  dr2 = mult2ptr->table->denserows;
  dense_op = op_table_type==DENSE;
  table1 = mult1ptr->table->table_data_ptr;
  table2 = mult2ptr->table->table_data_ptr;

  hash_init(&ht,FALSE,0,0,0);
  ht_ptr = ht.current_ptr;
  ht_ptr[0] = mult1ptr->initial[1];
  ht_ptr[1] = mult2ptr->initial[1];
  im = hash_locate(&ht,2);
/* Each state in the composite transition table will be represented as a subset
 * of the set of ordered pairs in which the first component is a state in
 * *multptr1 and the second a state in *multptr2.
 * The initial state contains just one such pair, of which the components are
 * the initial states of *mult1ptr and *multptr2.
 * The subsets will be stored as variable-length records in the hash-table,
 * as a list of pairs in increasing order.
 * If a state is reached from a transition ($,x) or (x,$) (with $ the padding
 * symbol), then it needs to be marked as such, since we do not allow a $
 * to be followed by any other generator.
 * We do this by adding a 1 or a 2 to the end of the statelist - this is
 * recognisable by the fact that the length of the statelist then becomes odd.
 */
  if (im!=1) {
    fprintf(stderr,"Hash-initialisation problem in fsa_composite.\n");
    return 0;
  }
  if ((tempfile=fopen(compfilename,"w"))==0){
    fprintf(stderr,"Error: cannot open file %s\n",compfilename);
    return 0;
  }
  if (dense_op)
    tmalloc(fsarow,int,ne)
  else
    tmalloc(fsarow,int,2*ne+1)
 
  cstate = 0;
  if (dense_op)
    len = ne; /* The length of the fsarow output. */
  nt = 0; /* Number of transitions in compositeptr */

  while (++cstate <= ht.num_recs) {
    if (kbm_print_level>=3) {
      if ((cstate<=1000 && cstate%100==0)||(cstate<=10000 && cstate%1000==0)||
          (cstate<=100000 && cstate%5000==0) || cstate%50000==0)
       printf("    #cstate = %d;  number of states = %d.\n",cstate,ht.num_recs);
    }
    reclen = hash_rec_len(&ht,cstate);
    cs_ptr = hash_rec(&ht,cstate);
    cs_ptre = hash_rec(&ht,cstate) + reclen - 1;
    /* for (i=0;i<reclen;i++) printf(" %d",cs_ptr[i]); printf("\n"); */
    if (reclen % 2 == 1) {
      if (*cs_ptre==1) {
        leftpad=TRUE;  rightpad=FALSE;
      } else {
        leftpad=FALSE;  rightpad=TRUE;
      }
      cs_ptre--;
    }
    else {
      leftpad=FALSE; rightpad=FALSE;
    }
      
    if (dense_op)
      for (i=0;i<ne;i++) fsarow[i] = 0;
    else
      len = 0;

    e = 0; /* e is the edge number of generator pair (g1,g3) */
    for (g1=1;g1<=ngens1;g1++)  for (g3=1;g3<=ngens1;g3++) {
/* Calculate action of pair (g1,g3) on state cstate  - to get the image, we
 * have to apply ( (g1,g2), (g2,g3) ) to each ordered pair in the subset
 * corresponding to cstate, and this for each generator g2 of the
 * base-alphabet (including the padding symbol).
 */
      e++;
      if (g1==ngens1 && g3==ngens1)
        continue;
      if ((leftpad && g1<=ngens) || (rightpad && g3<=ngens))
        continue;
      ht_ptrb = ht.current_ptr;
      ht_ptre = ht_ptrb-1;
      es1 = (g1-1)*ngens1 + 1;
      ef1 = g1*ngens1;
/* As g2 ranges from 1 to ngens+1 in the pair (g1,g2), for fixed g1, the
 * corresponding edge number in the fsa ranges from es1 to ef1.
 *
 * As g2 ranges from 1 to ngens+1 in the pair (g2,g3), for fixed g3, the
 * corresponding edge number in the fsa ranges from g3 upwards, with
 * increments of size ngens+1.
 */

      ptr = cs_ptr;
      while (ptr <= cs_ptre) {
        csa = *(ptr++); csb = *(ptr++);
        for (g2=1,e1=es1,e2=g3; e1<=ef1; g2++,e1++,e2+=ngens1){
          csima = g1==ngens1 && g2==ngens1 ?
                   (mult1ptr->is_accepting[csa]>0 ? csa : 0) :
                   target(dense_ip1,table1,e1,csa,dr1);
          if (csima==0)
            continue;

          csimb = g2==ngens1 && g3==ngens1 ?
                   (mult2ptr->is_accepting[csb]>0 ? csb : 0) :
                   target(dense_ip2,table2,e2,csb,dr2);
          if (csimb==0)
            continue;

          if (ht_ptrb>ht_ptre || csima> *(ht_ptre-1) ||
                               (csima== *(ht_ptre-1) &&  csimb > *ht_ptre) ){
/* We have a new state for the image subset to be added to the end */
            *(++ht_ptre) = csima; *(++ht_ptre) = csimb;
          }
          else {
            ht_chptr = ht_ptrb;
            while (*ht_chptr < csima)
              ht_chptr+=2;
            while (*ht_chptr==csima && *(ht_chptr+1)<csimb)
                ht_chptr+=2;
            if (csima < *ht_chptr || csimb < *(ht_chptr+1)) {
/* we have a new state for the image subset to be added in the middle */
              ht_ptr = ht_ptre;
              ht_ptre += 2;
              while (ht_ptr >= ht_chptr) {
                *(ht_ptr+2) = *ht_ptr;
                ht_ptr--;
              }
              *ht_chptr = csima; *(ht_chptr+1) = csimb;
            }
          }
        }
      }
      if (ht_ptre > ht_ptrb) {
        if (g1==ngens1)
          *(++ht_ptre) = 1;
        else if (g3==ngens1)
          *(++ht_ptre) = 2;
      }
      im = hash_locate(&ht,ht_ptre-ht_ptrb+1);
      if (im > 0) {
        if (dense_op)
          fsarow[e-1] = im;
        else  {
          fsarow[++len] = e;
          fsarow[++len] = im;
        }
        nt++;
      }
    }
    if (!dense_op)
      fsarow[0] = len++;
    fwrite((void *)fsarow,sizeof(int),(size_t)len,tempfile);
  }
  fclose(tempfile);

  ns = compositeptr->states->size = ht.num_recs;
  compositeptr->table->numTransitions = nt;

  if (kbm_print_level>=3)
    printf("    #Calculated transitions - %d states, %d transitions.\n",ns,nt);

/* Locate the accept states. This is slightly cumbersome, since a state
 * is an accept state if either the corresponding subset contains a
 * a pair (s1,s2), where s1 is and accept state of *mult1ptr and s2 an
 * accept state of *mult2ptr, or we can get from some such state to an
 * accept state by applying elements ( ($,x), (x,$ ).
 * We will need to use the array compositeptr->is_accepting.
 */

  tmalloc(compositeptr->is_accepting,boolean,ns+1);
  for (i=1;i<=ns;i++)
    compositeptr->is_accepting[i] = FALSE;
  ct = 0;
  es1 = ngens*ngens1 + 1;
  ef1 = ngens1*ngens1-1;
  for (cstate=ns;cstate>=1;cstate--) {
/* We work backwards through the states, since we wish to add on additional
 * elements at the end of the list in the hash-table - this destroys
 * later entries, but that doesn't matter, since we have already dealt
 * with them.
 */
    cs_ptr = hash_rec(&ht,cstate);
    reclen = hash_rec_len(&ht,cstate);
    if (reclen % 2 ==1)
      reclen--; /* The last entry marks the fact that this is a "padded state"*/
    cs_ptre = hash_rec(&ht,cstate) + reclen - 1;
/* First see if the set itself contains an accept-state */
    ptr = cs_ptr;
    while (ptr <= cs_ptre) {
      csa = *(ptr++); csb = *(ptr++);
      if (mult1ptr->is_accepting[csa] && mult2ptr->is_accepting[csb]) {
        compositeptr->is_accepting[cstate] = TRUE;
        ct++;
        break;
      }
    }
    if (compositeptr->is_accepting[cstate])
      continue;
/* Next apply generators ( ($,g2), (g2,$) ) and see if we get anything new.
 * We won't bother about having them in increasing order this time.
 */
    ptr = cs_ptr;
    while (ptr <= cs_ptre) {
      csa = *(ptr++); csb = *(ptr++);
      for (e1=es1,e2=ngens1; e1<=ef1; e1++,e2+=ngens1){
        csima = target(dense_ip1,table1,e1,csa,dr1);
        if (csima==0)
          continue;
        csimb = target(dense_ip2,table2,e2,csb,dr2);
        if (csimb==0)
          continue;

  /* see if (csima,csimb) is accepting */
        if (mult1ptr->is_accepting[csima] && mult2ptr->is_accepting[csimb]) {
          compositeptr->is_accepting[cstate] = TRUE;
          ct++;
          break;
        }
  /* now see if it is new */
        ht_chptr = cs_ptr;
        got = FALSE;
        while (ht_chptr < cs_ptre) {
          if (csima == ht_chptr[0] && csimb == ht_chptr[1]) {
            got = TRUE;
            break;
          }
          ht_chptr+=2;
        }
        if (!got) {
  /* add (csima,csimb) to the end */
          *(++cs_ptre) = csima; *(++cs_ptre) = csimb;
        }
      }
      if (compositeptr->is_accepting[cstate])
        continue;
    }
  }

  compositeptr->num_accepting = ct;
  if (ct==1 || ct != ns) {
    tmalloc(compositeptr->accepting,int,ct+1);
    ct = 0;
    for (i=1;i<=ns;i++) if (compositeptr->is_accepting[i])
      compositeptr->accepting[++ct] = i;
  }
  tfree(compositeptr->is_accepting);
  tfree(mult1ptr->is_accepting);
  tfree(mult2ptr->is_accepting);
  hash_clear(&ht);
  tfree(fsarow);
  if (destroy) {
    fsa_clear(mult1ptr); fsa_clear(mult2ptr);
  }

  if (readback) {
    tempfile = fopen(compfilename,"r");
    compressed_transitions_read(compositeptr,tempfile);
    fclose(tempfile);
    unlink(compfilename);
  }

  return compositeptr;
}
