standalone/progs.d/tcp.c

#include "defs.h"

extern char rs,ginrel[];
extern int space,mxc;
extern short mpt,rel[],cosno[],gno[],inv[],gch[],*imcos[];
short ng,ngi,nsg,nr,endsg,endr,maxcos,str,len,ad,*fpt,*bpt,
      ccos,maxd,totd,lastd,cind,nfree,stcr,endcr,fcos,bcos,lcl;
char fullsc,clsd,lkah;
FILE *op;

void inperr (char *s, int n)
{ fprintf(stderr,"Input error in %s no %d\n",s,n); }

int digit (int j)
{ if (j>='0' && j<='9') return(1); else return(0); }

int letter (int j)
{ if ((j>='a' && j<='z') || (j>='A' && j<='Z')) return(1); else return(0);}

void seeknln (void) { while (getchar()!='\n'); }

int readrel (int s, int no)
/* Reads relation (s=1) or subgroup generator (s=0) number "no" into array
  rel, after expanding powers, brackets etc.
   Each relation or subgen is preceded by its length.
  The global variables ad and str are used to mark the current position
  in rel.
*/
{ short stbr,endbr,exp,l,m,n; char ch;
  char gotg,br,clbr,emptybr,name[10];
  if (s) strcpy(name,"relation"); else strcpy(name,"subgen");
  gotg=0; br=0; clbr=0; ad=str; len=0;
  ch=getchar();
  while (ch!='\n')
  { if (ch==' ') ch=getchar();
    else if (ch=='(')
    { if (br || clbr) { inperr(name,no); return(-1); }
      emptybr=1; br=1; stbr=ad+1; gotg=0; ch=getchar();
    }
    else if (ch==')')
    { if (br==0 || emptybr) {inperr(name,no); return(-1); }
      br=0; gotg=0; clbr=1; endbr=ad; ch=getchar();
    }
    else if (letter(ch))
    { if (clbr || gch[ch]==0) {inperr(name,no); return(-1);}
      emptybr=0; gotg=1; len++; ad++; rel[ad]=gch[ch]; ch=getchar();
    }
    else if (ch=='-' || digit(ch))
    { if (gotg==0 && clbr==0) { inperr(name,no); return(-1); }
      gotg=0;
      if (ch=='-')
      { ch=getchar();
        if (digit(ch)==0) exp= -1;
        else
        { exp=0; while (digit(ch)) {exp*=10; exp-=(ch-'0'); ch=getchar();}}
      }
      else
      { exp=0; while (digit(ch)) {exp*=10; exp+=(ch-'0'); ch=getchar();}}
      if (exp==0) {inperr(name,no); return(-1);}
      if (clbr)
      { if (exp<0) for (m=stbr,n=endbr;m<=n;m++,n--)
        { if (m==n) rel[m]= -rel[m];
          else { l= -rel[m]; rel[m]= -rel[n]; rel[n]=l; }
        }
        exp=abs(exp); exp--; clbr=0;
        for (n=1;n<=exp;n++) for (m=stbr;m<=endbr;m++)
        { len++; ad++; rel[ad]=rel[m]; }
      }
      else
      { n=rel[ad];
        if (exp<0) { n= -n; rel[ad]=n; exp= -exp; }
        exp--;
        for (m=1;m<=exp;m++) { len++; ad++; rel[ad]=n; }
      }
    }
    else {inperr(name,no); return(-1); }
  }
  if (br || clbr) {inperr(name,no); return(-1); }
  rel[str]=len; str=ad+1;
  return(0);
}

int
tcprog (void)
/* This is the main routine */
{ short i,j,k,ch,*p;  int quot; char redspace,fname[80];
  printf("Todd-Coxeter Coset Enumeration Algorithm. (HLT + Lookahead).\n\n");
  printf("    Total space for relations and coset table = %d.\n",space);
  printf("    Generators should be upper or lower case letters.\n");
  printf("    Relations and subgroup generators are input as strings in the");
  printf(" generators.\n");
  printf("    Generators may be followed by a positive or negative exponent");
  printf(" ('-' = '-1').\n");
  printf("    Brackets (but no nested brackets) are allowed, and a bracketed");
  printf(" expression.\n");
  printf("      must be followed by an exponent.\n");
  printf("    Example:  ab-3c6(ab)-2(XYZ)-ba-\n\n\n");
  printf("Input numbers of generators, subgroup generators and relations.\n");
  scanf("%hd%hd%hd",&ng,&nsg,&nr);
  for (i=0;i<=127;i++) gch[i]=0;
  printf("Input names of generators (upper or lower case letters).\n");
/* Read generators as characters, and assign gen no gch[ch] to generator
   "ch", and -gch[ch] to its inverse. This numbering will be changed later.
*/
  for (i=1;i<=ng;i++)
  { ch= ' '; while (ch==' ' || ch=='\n') ch=getchar();
    if (letter(ch)==0)
    { fprintf(stderr,"Generators must be letters.\n"); return(-1);}
    if (gch[ch]!=0)
    { fprintf(stderr,"Repeated generator.\n"); return(-1); }
    gch[ch]=i; gno[i]=0; ginrel[i]=0;
  }
  seeknln();
  if (nsg>0)
  printf("Input subgroup generators (separated by new lines).\n");
  str=0; ad= -1;
  for (i=1;i<=nsg;i++) if (readrel(0,i)== -1) return(-1);
/* endsg marks address in array rel where subgens stop */
  endsg=ad;
  if (nr>0) printf("Input relators (separated by new lines).\n");
  for (i=1;i<=nr;i++)
  { if (readrel(1,i)== -1) return(-1);
/* Relations of form x^2 are not normally stored. We record the fact at this
   stage by putting gno[x]=1.
   ginrel[x]=1 means that generator x occurs in some relator.
*/
    if (len==2 && rel[ad]==rel[ad-1])
    { gno[abs(rel[ad])]=1; str-=3; ad-=3; }
    else
    { k=ad-len+1; for (j=k;j<=ad;j++) ginrel[abs(rel[j])]=1; }
  }
/* For involutions x with ginrel[x]=0, we store x^2 after all, to avoid
   a number of problems.
*/
  for (i=1;i<=ng;i++) if (gno[i]==1 && ginrel[i]==0)
  { rel[ad+1]=2; rel[ad+2]=i; rel[ad+3]=i; ad+=3; gno[i]=0; }
  endr=ad;
/* endr is the address in rel where relations end. */

  ngi= -1;
/* Now we start the renumbering process, and do the renumbering of the
   generators of the words stored in rel. Numbers will now all be >=0.
   inv[x] is the inverse of generator number x.
   If gen. no. x is an involution, then inv[x]=x. Otherwise inv[x]=x+1 (or x-1).
   gno[i] is the new number of the old generator number  i.
*/
  for (i=1;i<=ng;i++) if (gno[i]==1)
  { ngi++; gno[i]=ngi; gno[i+52]=ngi; inv[ngi]=ngi; }
  else
  { ngi+=2; gno[i]=ngi-1; gno[i+52]=ngi; inv[ngi]=ngi-1; inv[ngi-1]=ngi;}
  i= -1;
  while (++i<=endr)
  { j=rel[i]; k=i+j;
    while (++i<=k) { p=rel+i; *p= (*p<0) ? gno[52- *p] : gno[*p]; }
    i=k;
  }

/* Now we arrange the remaining space in the array rel into the coset table.
   imcos[i][j]=image of coset j under generator i (or 0 when undefined).
   Some space is also required for the arrays fpt and bpt, which are used
   as forward and backward pointers in the linked list of active cosets.
   (But they are used somewhat differently in the coincidence routine.)
   maxcos is the maximum number of cosets that we have room for.
*/
compmaxcos:
  quot= (space-endr-2)/(ngi+3); if (quot>mxc) quot=mxc; maxcos=quot;
  fpt=rel+1+endr; bpt=fpt+maxcos;
  redspace=0;
  for (i=0;i<=ngi;i++) imcos[i]=fpt+maxcos*(2+i);
/* Allow user to reduce maxcos if -r was set. */
  if (rs)
  { printf("Maxcos=%d.  If you wish to reduce this, input smaller number:  ",
            maxcos);
    if ((j=getchar())!='\n')
    { i=0; while (j==' ') j=getchar();
      if (digit(j))
      { i=j-'0'; j=getchar(); while (digit(j)) {i=10*i+j-'0'; j=getchar();}
        if (j!='\n') seeknln();
      }
      else if (j!='\n') seeknln();
      if (i>0 && i<maxcos) { maxcos=i; redspace=1;}
      printf("Maxcos=%d.",maxcos);
    }
    printf("\n\n");
  }
  else printf("Maxcos=%d.\n",maxcos);

/* Now we are ready to start the enumeration.
   ccos=current coset being scanned.
   lastd=last coset defined.
   maxd=max no of cosets that were defined at any one time.
   totd=total number of cosets defined.
   cind=current index=current no of cosets defined
   nfree=next available coset no.
   The defined coset nos. form doubly linked list, using fpt and bpt
   The available coset nos are linked forward with fpt, starting at nfree.
   When nfree=0, there are no further nos. free (cind=maxcos), and then
   lookahead is entered (lkah=1)
*/
  ccos=1; lastd=1; maxd=1; totd=1; cind=1; nfree=2;
  for (i=0;i<=ngi;i++) imcos[i][1]=0;
  for (i=0;i<maxcos;i++) fpt[i]=i+1;
  fpt[1]=0; fpt[maxcos]=0;
  lkah=0; bpt[1]=0; endcr= -1;
/* endcr runs thro' addresses in rel, whilst scanning relations.
   First we scan coset 1 under the subgens.
*/
  while (endcr!=endsg)
  { scanrel();
/*
    The external variable fullsc is set 1 or 0 in scanrel.
    fullsc=1 if the relation or subgen is completely scanned (possibly after
    making new definitions). If not, then cind=maxcos, and no more
    definitions could be made. If this occurs at this stage, then we give up
    straightaway! (I have never known this to happen.)
*/
    if (fullsc==0)
    { fprintf(stderr,"Space overflow in subgroup generation phase.\n");
      return(-1);
    }
  }
/* Now we scan each coset under the relations.
   If fullsc=0, then we enter lookahead. clsd=0 means that a given coset
   was not fully scanned under all relations, so it must be rescanned after
   we come out of lookahead. lcl marks the last fully scanned coset before
   entering lookahead, so this is the return point.
*/
  while (ccos!=0)
  { clsd=1; endcr=endsg;
    while (endcr!=endr)
    { scanrel();
      if (fullsc==0)
      { clsd=0;
        if (lkah==0)
        { lcl=bpt[ccos]; printf("Entering lookahead.\n"); lkah=1; }
      }
    }
    if (lkah)
    { i=fpt[ccos];
/* If the coset was fully scanned, then we change the linking to put it on
   the end of the list of scanned cosets.
*/
      if (clsd)
      { j=bpt[ccos];
        if (j!=lcl)
        { fpt[j]=i;
          if (i==0) lastd=j; else bpt[i]=j;
          j=fpt[lcl]; fpt[lcl]=ccos; bpt[ccos]=lcl; fpt[ccos]=j; bpt[j]=ccos;
        }
        lcl=ccos;
      }
      ccos=i;
      if (ccos==0)
/* End of lookahead */
      { if (cind==maxcos)
        { fprintf(stderr,"Not enough space.\n");
          if (redspace) goto compmaxcos;
          return(-1);
        }
        printf("Exiting lookahead. No. of cosets=%d\n",cind);
        ccos=fpt[lcl]; lkah=0;
      }
    }
    else ccos=fpt[ccos];
  }

/* That ends the enumeration. In case ginrel[i] is 0 for any generator, we
   must check that the coset table is full for that gen. Otherwise index is
   certainly infinite.
*/
  for (i=1;i<=ng;i++) if (ginrel[i]==0)
  { k=gno[i]; j=lastd;
    while (j!=0)
    { if (imcos[k][j]==0)
      { fprintf(stderr,"Coset table incomplete at end of scan. Index is infinite.\n");
       return(-1);
      }
      j=bpt[j];
    }
  }
  printf("Algorithm complete. maxdef,totdef=%d,%d.\n\n",maxd,totd);
  if (nsg==0) printf("The order of the group is %d.\n\n",cind);
  else printf("The index of the subgroup is %d.\n\n",cind);
  if (cind<=mpt)
/* Now we compute and store the permutation action of the generators if required
*/
  { printf("Do you wish to store the permutations? (y/n)    ");
    ch=getchar(); seeknln();
    if (ch=='y')
    { printf("Input filename    "); scanf("%s",fname);
      op=fopen(fname,"w"); fprintf(op,"%4d%4d%4d%4d\n",cind,ng,0,0);
      bpt[1]=1; cosno[1]=1; j=1;
      for (i=2;i<=cind;i++) { j=fpt[j]; cosno[i]=j; bpt[j]=i; }
      for (i=1;i<=ng;i++)
      { k=gno[i];
        if (cind>=1000)
        for (j=1;j<=cind;j++) fprintf(op,"%5d",bpt[imcos[k][cosno[j]]]);
        else
        for (j=1;j<=cind;j++) fprintf(op,"%4d",bpt[imcos[k][cosno[j]]]);
        fprintf(op,"\n");
      }
    }
  }
}

int
scanrel (void)
/* Scan ccos under relation or subgen starting at rel[endcr] */
{ short i,j,k,l,m; char comp;
/* Put endcr to point to next relation */
  fullsc=1; stcr=endcr+2; endcr+=(1+rel[stcr-1]);
  fcos=ccos; bcos=ccos; comp=1;
  for (i=stcr;i<=endcr;i++)
/* Forward scan. If incomplete, try backward scan
   If complete, check for coincidence.
*/
  { k=imcos[rel[i]][fcos];
    if (k==0) { comp=0; break;}
    fcos=k;
  }
  if (comp) { if (fcos!=bcos) coinc(fcos,bcos); return(0); }
  stcr=i;
  for (i=endcr;i>=stcr;i--)
/* Backward scan. If incomplete, make new definitions if possible */
  { l=rel[i]; m=inv[l]; k=imcos[m][bcos];
    if (k==0)
    { if (i==stcr)
      { k=imcos[l][fcos];
        if (k==0) { imcos[l][fcos]=bcos; imcos[m][bcos]=fcos; }
        else if (k!=bcos) coinc(k,bcos);
        return(0);
      }
      if (lkah || nfree==0)
      { fullsc=0; return(0); /* No new definition possible */}
      for (j=0;j<=ngi;j++) imcos[j][nfree]=0;
      totd++; cind++;
      if (maxd<cind) maxd++;
      imcos[m][bcos]=nfree; imcos[l][nfree]=bcos; bcos=nfree;
      bpt[nfree]=lastd; fpt[lastd]=nfree; lastd=nfree;
      nfree=fpt[nfree]; fpt[lastd]=0;
    }
    else bcos=k;
  }
  if (fcos!=bcos) coinc(fcos,bcos);
}

int
coinc (int c1, int c2)
/* Process coincidence c1 = c2.
   fpt and bpt are used differently in this routine.
   For pairs d1,d2 of coincidences in the queue waiting to be processed,
   (d1>d2), bpt[d1] is the next coset in the queue, and fpt[d1]= -d2.
   d1 will be eliminated eventually.
*/
{ short lc,hc,qh,qt,i,j,x,fhc,bhc,lim,him;
  if (c1<c2) { lc=c1; hc=c2; }
  else { lc=c2; hc=c1; }
/* hc will be eliminated. qh,qt are head and tail of coincidence queue */
  qh=0; qt=0;
/* Unlink hc from linked list of live cosets */
  fhc=fpt[hc]; bhc=bpt[hc]; fpt[bhc]=fhc;
  if (fhc==0) lastd=bhc; else bpt[fhc]=bhc;
  if (ccos==hc) { ccos=bhc; endcr=endr; clsd=0; }
  if (lkah && lcl==hc) lcl=bhc;
  while (1)
/* Each loop corresponds to one coincidence pair hc=lc */
  { fpt[hc]=nfree; nfree=hc; cind--;
/* Look at images of hc,lc under each generator, to deduce further possible
   coincidences
*/
    for (i=0;i<=ngi;i++)
    { him=imcos[i][hc];
      if (him!=0)
      { j=inv[i]; lim=imcos[i][lc];
        if (him==hc) him=lc; else imcos[j][him]=0;
/* imcos[j][him] was previously hc, so we remove this entry. It will either be
   put equal to lc at the end of this loop, or later filled in as a
   consequence of another coincidence involving lc.
*/
        if (lim==0) imcos[i][lc]=him;
        else
        { if (lim==hc) { imcos[i][lc]=lc; lim=lc; }
/* him,lim is a new coincident pair. First we reduce it using the queue */
          while (fpt[him]<0) him= -fpt[him];
          while (fpt[lim]<0) lim= -fpt[lim];
          if (him!=lim)
          { if (lim>him) { x=lim; lim=him; him=x; }
/* Eliminate him from live cosets */
            fhc=fpt[him]; bhc=bpt[him]; fpt[bhc]=fhc;
            if (fhc==0) lastd=bhc; else bpt[fhc]=bhc;
            if (ccos==him) { ccos=bhc; endcr=endr; clsd=0; }
            if (lkah && lcl==him) lcl=bhc;
/* Add him,lim to coincidence queue */
            fpt[him]= -lim;
            if (qh==0) qh=him; else bpt[qt]=him;
            qt=him; bpt[qt]=0;
          }
        }
        x=imcos[i][lc];
        if (imcos[j][x]==0) imcos[j][x]=lc;
      }
    }
/* Get next coincident pair from head of queue, if any */
    if (qh==0) break;
    hc=qh; qh=bpt[qh]; lc= -fpt[hc]; bpt[hc]=0;
  }
}