# include <stdio.h>
extern char mult,inf0[],inf1[],inf2[],outf[],outft[],inf3[];
extern short  csp[],*cpsp[],**cdpsp[],***coeff[];
extern int space,cspace;
extern short sp[],**mat[],*psp[],**imcos[],**cpco[],lorb[],
               wdl,ptrsp,cdptrsp,mb;
short *cst,**cpst,***cdpst,**svptr,*cp,*rel;
short *spst,**pspst,**pptr,**cpptr,npt,nb,nph,nph2,npg,npg2,
      rno,orno,index,*invg;
FILE  *fopen(),*ip,*op;

/* This program differs from most other permutation programs in that perms are
   all stored in the single array sp. Schreier vectors are stored in the short
   array csp. This assumes that there are at most 127 perms, which will
   probably always be the case. The first half of the program (crprog1)
   involves relatively little short space, so it is assumed that there is enough
   In the first half, coefficients as words in H generated by the relators
   of G are computed using the output of conrun.
*/
crprog1()
{ short *pc,*qc,ex,neg; int x;
  short i,j,k,l,m,n,cl,rl,*p,ocl,im,pt,pt1,pn,ipt;
  if ((ip=fopen(inf2,"r"))== 0)
  { fprintf(stderr,"Cannot open %s.\n",inf2); return(-1);}
  fscanf(ip,"%hd%hd%hd%hd",&npt,&nph,&nb,&k);
  if (nb>=mb) {fprintf(stderr,"nb too big. Increase MB.\n"); return(-1);}
  if (k<=2) {fprintf(stderr,"inf2 has illegal format.\n"); return(-1); }
  seeknln(); seeknln();
  for (i=1;i<=nb;i++) fscanf(ip,"%hd",lorb+i);
  pptr=psp-1; pspst=psp+nph; svptr=cpsp-1;  cpst=cpsp+nb;
  invg=sp; nph2=2*nph; spst=sp+nph2;
  for (i=1;i<=nph;i++)
  { pptr[i]=spst+(i-1)*npt-1; p=pptr[i];
    for (j=1;j<=npt;j++) {fscanf(ip,"%hd",&k); p[k]=j; }
    invg[2*i-2]=2*i-1; invg[2*i-1]=2*i-2;
    seeknln();
  }
  spst+=(npt*nph);
  for (i=1;i<=nb;i++)
  { svptr[i]=csp+(i-1)*npt-1; pc=svptr[i];
    for (j=1;j<=npt;j++) {fscanf(ip,"%hd",&k); pc[j]= k>0 ? (k+1)/2 : k;}
  }
  cp=csp+nb*npt; rel=cp+wdl;
  fclose(ip);
/* It is assumed that there is enough space so far. */

  strcpy(inf1,inf0); strcat(inf1,".cp");
  if ((ip=fopen(inf1,"r"))==0)
  { fprintf(stderr,"Cannot open %s.\n",inf1); return(-1);}
  fscanf(ip,"%hd%hd",&index,&npg); npg2=2*npg-1;
  x= (nb+2)*npg*index;
  if (spst+x-sp>space)
  { fprintf(stderr,"Not enough space. Increase SPACE.\n"); return(-1); }
  if (pspst+(nb+2)*npg-psp>ptrsp)
  { fprintf(stderr,"Not enough ptrsp. Increase PTRSP.\n"); return(-1); }
  cpptr=pspst; pspst+=2*npg;
  for (i=0;i<=npg2;i++) { cpptr[i]=spst-1; spst+=index; }
  for (i=1;i<=nb;i++)
  { cpco[i]=pspst-1; pspst+=npg;
    for (j=1;j<=npg;j++) {cpco[i][j]=spst-1; spst+=index;}
  }
  for (i=1;i<=npg;i++)
  { j=2*i-2; k=j+1;
    for (l=1;l<=index;l++)
    { for (m=1;m<=nb;m++) fscanf(ip,"%hd",cpco[m][i]+l);
      fscanf(ip,"%hd",&m); cpptr[j][l]=m; cpptr[k][m]=l;
    }
  }
  fclose(ip);

  strcpy(inf1,inf0); strcat(inf1,".rel");
  if ((ip=fopen(inf1,"r"))==0)
  { fprintf(stderr,"Cannot open %s.\n",inf1); return(-1);}
  fscanf(ip,"%hd%hd",&k,&rno); seeknln();
  op=fopen(outft,"w");
/* Now we have read everything in, and the computation can start */
  orno=0;
  for (i=1;i<=rno;i++)
  { fscanf(ip,"%hd",&rl);
    printf("Scanning rel no %d.\n",i);
    for (j=1;j<=rl;j++) { fscanf(ip,"%hd",&k); rel[j]=k; }
    for (j=1;j<=index;j++)
    { pt=j; cl=0;
      for (k=rl;k>=1;k--)
      { l=rel[k]; im=cpptr[l][pt];
        if (l%2==1) {neg=1;pt1=im; m=(l+1)/2;}
        else {neg=0; pt1=pt; m=(l+2)/2;}
        ocl=cl;
        for (n=1;n<=nb;n++)
        { ipt=cpco[n][m][pt1]; pn=svptr[n][ipt];
          while (pn!= -1)
          { cl++; cp[cl]=2*(pn-1); ipt=pptr[pn][ipt]; pn=svptr[n][ipt]; }
        }
        if (neg)
        { pc=cp+ocl+1; qc=cp+cl;
          while (pc<=qc)
          { if (pc==qc) (*pc)++;
            else {ex= *pc; *pc= *qc+1; *qc= ex+1;}
            pc++; qc--;
          }
        }
        pt=im;
      }
      if (pt!=j) 
      { fprintf(stderr,
                "Relation number %d is not satisfied by the permutations.\n",i);
         return(-1);
      }
      fprintf(op,"%4d",cl);
      for (k=cl;k>=1;k--) fprintf(op," %2d",cp[k]);
      fprintf(op,"\n");
    } /* j loop */
    if (i==rno && orno==0)
    { orno=rno; if (fscanf(ip,"%hd",&j)>0) rno+=j;}
      
  } /* i loop */
  printf("All relations scanned.\n");
  fclose(op); fclose(ip);
  return(0);
}

seeknln()
{ while (getc(ip)!='\n'); }