#include <stdio.h>
#define tmalloc(D,T,N) {D = (T *) malloc(sizeof(T)*(N)); \
  if (D==0) { fprintf(stderr,"Out of space.\n"); return(-1);}}
#define tfree(D) {if (D) free( (char *) D); D=0;}

char *malloc();

extern char mult,inf1[],inf2[],inf3[],outf[];
extern short rwd[],
      ***scoeff[],**cdpsp[],*cpsp[],csp[],wd1[],wd2[],wd3[],wd4[];
extern int cspace,psp,space,cptrsp,svsp;
extern short  perm[],sv[],orb[],imsp[],*ptsp[],**simcos[],pinv[],gno[],
       base[],lorb[],pno[],*pptr[],*svptr[],cord[],
       invg[],rno[],cp[],mspace[],*vec[],**mat[],
       mp,cdptrsp,mwl2,mpt,mb,mdim,mpr,ptrsp,mm,msp,mv,mwdl,mlwdl;
int rsp;
short endr,maxcos,*fpt,*bpt,nelim,bno,lo,**imcos,*spst,**pspst,dim,
      prime,ccos,lastd,cind,nfree,stcr,endcr,fcos,bcos,lcl,npt,np2,*rel,
      *spv,**spm,ng,rwl,nb;
short *cst,*cend,***coeff,**cpst,***cdpst,fullsc,clsd,lkah,conch,**cco,*ocst,
      **def;
FILE *fopen(),*ip,*op;

extpprog()
{ short i,j,k,l,m,n,b,l1,nperms,npt1,mxp,maxv,maxm,pos,*p,*q,*r;
  int quot;
  short **dp,**dq,**dr,*pc; float a;
  if ((ip=fopen(inf1,"r"))==0)
  { fprintf(stderr,"Cannot open %s.\n",inf1); return(-1); }
  fscanf(ip,"%hd%hd%hd%hd",&npt,&nperms,&nb,&k);
  if (npt>mpt) {fprintf(stderr,"npt too big. Increase NPT.\n"); return(-1); }
  if (nb>mb) {fprintf(stderr,"nb too big. Increase MB.\n"); return(-1); }
  if (nb*npt>svsp)
  { fprintf(stderr,"Not enough sv space. Increase SVSP.\n"); return(-1); }
  if (k<=2) {fprintf(stderr,"Wrong input format.\n"); return(-1); }
  npt1=npt+1; np2=2*nperms-1;
  quot=psp/npt1; if (quot>mp) quot=mp; mxp=quot;
  if (np2>=mxp)
  { fprintf(stderr,"Too many permutations. Increase MP and/or PSP.\n");
    return(-1);
  }
  for (i=0;i<mxp;i++) pptr[i]=perm+i*npt1-1;
  for (i=1;i<=nb;i++) svptr[i]=sv+(i-1)*npt-1;
  readbaselo(nb,base,lorb); readpsv(0,nb,nperms,svptr);
  fclose(ip);
  for (i=0;i<np2;i+=2) { invg[i]=i+1; invg[i+1]=i; }

/* We need to read in matrices for generators when action is nontrivial */
  if (mult==0)
  { if ((ip=fopen(inf3,"r"))==0)
    { fprintf(stderr,"Cannot open %s.\n",inf3); return(-1); }
    fscanf(ip,"%hd%hd%hd",&prime,&dim,&k);
    if (prime>mpr)
    { fprintf(stderr,"Prime too big. Increase MPR.\n"); return(-1); }
    if (dim>mdim)
    { fprintf(stderr,"Dim too big. Increase MDIM.\n"); return(-1); }
    if (k!=nperms) { fprintf(stderr,"No of mats wrong.\n"); return(-1); }
    setpinv();
    maxv=msp/dim; if (maxv>mv) maxv=mv;
    for (i=0;i<maxv;i++) vec[i]=mspace-1+i*dim;
    maxm=maxv/dim-1; if (maxm>=mm) maxm=mm-1;
    if (np2>=maxm)
    { fprintf(stderr,"Not enough mat space. Increase MSP and/or MV,MM.\n");
      return(-1);
    }
    for (i=0;i<=maxm;i++) mat[i]=vec-1+i*dim;
    spm=mat[maxm]; spv=spm[1];
    for (i=0;i<np2;i+=2) { readmat(mat[i]); inv(mat[i],mat[i+1]);}
    fclose(ip);
  }

  printf("Input a,b,l1.\nmaxcos = a*lorb+b\n");
  printf("Exit lookahead after eliminating abs(l1) cosets.\n");
  printf("l1 negative for consistency checking.\n");
  scanf("%f%hd%hd",&a,&b,&l1);
  if (l1<0) {l1= -l1; conch=1;} else conch=0;

/* Now we read in the relations */
  if ((ip=fopen(inf2,"r"))==0)
  { fprintf(stderr,"Cannot open %s.\n",inf2); return(-1); }
  fscanf(ip,"%hd%hd",&i,&j);
  if (i!=nb || (mult==0 && j!=dim))
  { fprintf(stderr,"Error in %s, line 1.\n",inf2); return(-1); }
  for (i=1;i<=nb;i++) fscanf(ip,"%hd",rno+i);
  if (mult) {dim=j; for (i=1;i<=dim;i++) fscanf(ip,"%hd",cord+i); }
  rel=imsp; spst=rel;
  for (i=1;i<=2*rno[1];i++)
  { fscanf(ip,"%hd",spst); l= *spst;
    for (j=1;j<=l;j++) fscanf(ip,"%hd",spst+j);
    spst+=(l+1);
  }
  fclose(ip);

/* Now we are ready to start the modified Todd-Coxeter procedure. This
   follows the ordinary tcrun with procedures added in the appropriate
   places to compute the coefficients as words in the generators of G and M.
   The Tood-Coxeter procedures are in exb.c and the word operation
   procedures in exc.c
*/
  pos=0; endr= -1; pspst=ptsp; cst=csp; ocst=cst; cend=csp+cspace-1; cpst=cpsp;
  cdpst=cdpsp; rsp=cspace;
  for (bno=nb;bno>=1;bno--)
  { printf("bno=%d.\n",bno);
    l=rno[bno]-pos;
    for (i=1;i<=2*l;i++) {endr++; endr+= rel[endr];}
    pos=rno[bno];
    lo=lorb[bno];
    if (lo==1) continue;
    maxcos=lo*a+b; ng=0;
    tmalloc(def,short *,maxcos+1);
    for (i=0;i<np2;i+=2) if (pptr[i][npt1]>=bno) ng+=2; else break;
    k=(ng+1)*maxcos; l=k+maxcos+1; k-=lo;
    if (cpst+k-cpsp>cptrsp)
    { fprintf(stderr,"Out of cptrsp. Increase CPTRSP.\n"); return(-1); }
    if (spst+l-imsp>space)
    { fprintf(stderr,"Out of im space. Increase SPACE.\n"); return(-1); }
    if (pspst+ng-ptsp>ptrsp)
    { fprintf(stderr,"Out of ptrsp. Increase PTRSP.\n"); return(-1); }
    if (cdpst+ng-cdpsp>cdptrsp)
    { fprintf(stderr,"Out of cdptr Increase CDPTRSP.\n"); return(-1); }
    imcos=pspst; pspst+=ng; gno[bno]=ng;
    for (i=0;i<ng;i++)
    {imcos[i]=spst-1; p=spst-1; spst+=maxcos; while (++p<spst) *p=0; }
    fpt=spst; bpt=spst+maxcos;
    coeff=cdpst; cdpst+=ng;
    for (i=0;i<ng;i++)
    {coeff[i]=cpst-1; dp=cpst-1; cpst+=maxcos; while (++dp<cpst) *dp=0;}
    cco=cpst-lo-1;
    *pno=0;
    for (i=0;i<np2;i+=2) if (pptr[i][npt1]>=bno)
    { (*pno)++; pno[*pno]=i;
      if (pptr[i][npt1]>bno)
      { imcos[i][1]=1; imcos[i+1][1]=1;
        coeff[i][1]=cst; *cst=1; *(++cst)=i; *(++cst)=0; cst++;
        coeff[i+1][1]=cst; *cst=1; *(++cst)=i+1; *(++cst)=0; cst++; rsp-=6;
      }
    }
    if (orbitsv(base[bno],svptr[bno],0)!=lo)
    { fprintf(stderr,"Orbit length wrong!\n"); return(-1); }
    if (mult==0)
    { tmalloc(def[1],short,1); def[1][0]=0;}
    for (l=1;l<=lo;l++)
    { n=orb[l]; cp[n]=l;
      if ((k=svptr[bno][n])!= -1)
      { j=pptr[k][n]; m=cp[j]; imcos[k-1][m]=l; imcos[k][l]=m;
        if (mult==0)
        { tmalloc(def[l],short,def[m][0]+2);
          def[l][0]=def[m][0]+1; def[l][1]=k;
          for (i=1;i<=def[m][0];i++) def[l][i+1]=def[m][i];
        }
      }
    }

    lastd=lo; cind=lo; nfree=lo+1;
    for (i=1;i<=lo;i++) bpt[i]=i-1;
    for (i=0;i<maxcos;i++) fpt[i]=i+1;
    fpt[lo]=0; fpt[maxcos]=0;
    ccos=1; lkah=0;
/* At last we are ready to proceed with the enumeration */
    while (ccos!=0)
    { clsd=1; endcr= -1;
      while (endcr!=endr)
      { if (scanrel()== -1)  return(-1);
        if (fullsc==0)
        { clsd=0;
          if (lkah==0)
          { lcl=bpt[ccos]; lkah=1; nelim=0;
            printf("Entering lookahead.\n");
          }
        }
      }
      ccos=fpt[ccos];
      if (lkah)
      { if (nelim>=l1 || ccos==0)
        { if (cind!=maxcos)
          { printf("Exiting lookahead. cind=%d.\n",cind);
            ccos=fpt[lcl]; lkah=0;
          }
          else { fprintf(stderr,"maxcos too small.\n"); return(-1); }
        }
      }
    } /* while ccos!=0 */
/* End of enumeration. Tidy up and prepare for next bno */
    p=imcos[0]+lo; q=imcos[1]; k=maxcos-lo; dp=coeff[0]+lo; dq=coeff[1];
    for (i=1;i<ng;i++)
    { imcos[i]=p; r=p+lo;
      while (++p<=r) *p= *(++q);
      q+=k;
      coeff[i]=dp; dr=dp+lo;
      while (++dp<=dr) *dp= *(++dq);
      dq+=k;
    }
    spst=p+1; cpst=dp+1; simcos[bno]=imcos; scoeff[bno]=coeff;
    if (bno>1) { maxcos=lo; gc(); ocst=cst; }
    if (mult==0)
      for (i=1;i<=lo;i++) tfree(def[i]);
    tfree(def);
  } /* for (bno=;... */

/* End of program. Now output */
  op=fopen(outf,"w");
  fprintf(op,"%4d%4d\n",nb,dim);
  for (i=1;i<=nb;i++) fprintf(op," %3d",lorb[i]); fprintf(op,"\n");
  if (mult) {for (i=1;i<=dim;i++) fprintf(op,"%4d",cord[i]); fprintf(op,"\n");}
  for (i=1;i<=nb;i++) fprintf(op,"%4d",gno[i]); fprintf(op,"\n");
  for (i=1;i<=nb;i++) if ((l=lorb[i])>1) for (j=0;j<gno[i];j++)
  { p=simcos[i][j];
    for (k=1;k<=l;k++) fprintf(op," %3d",p[k]); fprintf(op,"\n");
    dp=scoeff[i][j];
    for (k=1;k<=l;k++)
    { if ((pc=dp[k])==0) fprintf(op,"%3d %3d ",0,0);
      else
      { fprintf(op,"%3d ",*pc);
        for (m=1;m<= *pc;m++) fprintf(op," %2d",pc[m]);
        pc+= (1+ *pc); fprintf(op," %3d ",*pc);
        for (m=1;m<= *pc;m++) fprintf(op," %2d",pc[m]);
      }
    }
    fprintf(op,"\n");
  }
  fclose(op); return(0);
}