#include <stdio.h>
extern char inf1[],inf2[],inf3[],outf1[],outf2[],outf3[],temp1[],temp2[],
     expg,exph,cr,dcr,triv;
extern int psp,trsp,svsp;
extern int  mpt,mb,
       tree[],perm[],gorb[],lorbg[],lorbh[],base[],
       fpt[],bpt[],index[],*cp[],*trad[],*sv[],
       *tailad[],**cpad[],**svgptr[],**svhptr[];
int trptr;
int npt,**lcp,**ucp,*stop,dummy,*pst,*pend,nb,fnt,lnt,
      ind,coind,oind,cind,bno,pt,**expp,npg;
FILE *fopen(),*ip,*op,*opx;

/* The data structures of this program differ from other permutation group
   programs. Permutations are not numbered, but are located by their base
   adresses. Consequently, arrays sv and cp are arrays of pointers. The
   current permutation in cp is stored in the middle of cp (since it has to be
   expanded in both directions), between adresses lcp and ucp of cp.
   The next six procedures are the equivlents of the corresponding ones
   in permfns.c, but addsv can go in both directions.
   In this half of the program, coset reps of H in G are computed.
*/

image(pt) int pt;
{ int **p;
  p=lcp;
  while (p<=ucp) {pt=(*p)[pt]; p++;}
  return(pt);
}

addsvb(pt,sv) int pt,**sv;
{ int *p;
  while ((p=sv[pt])!=stop) {pt=p[pt]; lcp--; *lcp =p-npt; }
}

addsvf(pt,sv) int pt,**sv;
{ int *p;
  while ((p=sv[pt])!=stop) {pt=p[pt]; ucp++; *ucp=p; }
}

invert(a,b) int *a,*b;
{int i; for (i=1;i<=npt;i++) b[a[i]]=i; }

rdperm(a,b) int *a,*b;
{ int i,*r;
  for (i=1;i<=npt;i++) fscanf(ip,"%d",a+i);
  r=pst+1; fscanf(ip,"%d",r); invert(a,b);
}

rdsv(sv) int **sv;
{ int i,j;
  for (i=1;i<=npt;i++)
  { fscanf(ip,"%d",&j); sv[i]= (j==0) ? 0 : (j== -1) ? stop : cp[j];}
}

cnprg1()
{ int i,j,k,ad,sad,nph,*svpt,*p1,*p2,**csvg,**csvh,nexp;
  char pthere,w,ok,allok;  long fac;
  stop = &dummy;
  if ((ip=fopen(inf1,"r"))==0)
  { fprintf(stderr,"Cannot open %s.\n",inf1); return(-1); }
  fscanf(ip,"%d%d%d%d",&npt,&npg,&nb,&k);
  if (npt>mpt) {fprintf(stderr,"npt too big. Increase mpt.\n"); return(-1); }
  if (nb>=mb) {fprintf(stderr,"nb too big. Increase MB.\n"); return(-1); }
  if (2*nb*npt>svsp)
  { fprintf(stderr,"svsp too small. Increase SVSP.\n"); return(-1); }
  if (k<=0) {fprintf(stderr,"Wrong input format for G.\n"); return(-1); }
  for (i=1;i<=nb;i++) fscanf(ip,"%d",base+i);
  for (i=1;i<=nb;i++) fscanf(ip,"%d",lorbg+i);
  pst=perm-1; pend=perm+psp-1;
  if (pst+2*npg*npt>pend)
  { fprintf(stderr,"Out of space.Increase PSP.\n"); return(-1);}
  for (i=1;i<=npg;i++)
  { p1=pst; p2=pst+npt; pst+=(2*npt); rdperm(p1,p2); cp[2*i-1]=p2; }
  for (i=1;i<=nb;i++) {svgptr[i]=sv+npt*(i-1)-1; rdsv(svgptr[i]); }
  fclose(ip);
/* G is read. Now read the subgroup H */
  if (triv)
  { nph=0; for (i=1;i<=nb;i++)
    { lorbh[i]=1; svhptr[i]=sv+npt*(nb+i-1)-1;
      for (j=1;j<=npt;j++) svhptr[i][j]=0; svhptr[i][base[i]]=stop;
    }
  }
  else
  { if ((ip=fopen(inf2,"r"))==0)
    { fprintf(stderr,"Cannot open %s.\n",inf2); return(-1); }
    fscanf(ip,"%d%d%d%d",&i,&nph,&j,&k); w=0;
    if (k<=0) { fprintf(stderr,"Wrong input format for H.\n"); return(-1); }
    if (i!=npt || j!=nb) w=1;
    else for (i=1;i<=nb;i++)
    { fscanf(ip,"%d",lorbh+i); if (lorbh[i]!=base[i]) w=1; }
    if (w)
    { fprintf(stderr,"npt or base for G and H do not agree.\n");return(-1);}
    for (i=1;i<=nb;i++) fscanf(ip,"%d",lorbh+i);
    if (pst+2*npt*nph+2>pend)
    { fprintf(stderr,"Out of space.Increase PSP.\n"); return(-1);}
    for (i=1;i<=nph;i++)
    { p2=pend-npt; p1=p2-npt; pend-=(2*npt); rdperm(p1,p2); cp[2*i-1]=p2;}
    for (i=1;i<=nb;i++) {svhptr[i]=sv+npt*(nb+i-1)-1; rdsv(svhptr[i]); }
    fclose(ip);
  }
  for (i=1;i<=nb;i++) fpt[i]=0;
  p1=pst; p2=p1+npt; pst+=(2*npt); lcp=cp; expp= cp+2*npt-1;
  cind=1; fac=1; index[nb+1]=1; lnt=0;

/* Now we compute the indices of H in G in the stabilizer chain into the array
   index. The nontrivial indices are linked by fpt and bpt.
*/
  for (i=nb;i>=1;i--)
  { j=lorbg[i]; fac=fac*j; j=lorbh[i]; fac=fac/j; ind=fac; index[i]=ind;
    if (ind>index[i+1])
    { if (lnt==0) lnt=i;
      else { fpt[i]=fnt; bpt[fnt]=i; }
      fnt=i;
    }
  }
  bpt[fnt]= -1; fpt[lnt]= -1;
  printf("Indices:  ");
  for (i=1;i<=nb;i++) printf("%5d",index[i]); printf("\n");
  if (ind==1) { fprintf(stderr,"G=H.\n"); return(-1); }
  nexp=0; ind=1; oind=1; trptr=0;
 
/* Now we start the backtrack search for the left coset reps of H in G. The
   necessary information about these is stored in the array tree. This is done
   so as to make it as easy as possible to compute which coset an arbitrary
   element of G lies in. If expg is true, then those perms in G which arise
   are expanded and stored.
   If dcr is true, then a small subset of the coset reps will need to be
   output as perms, later. For this purpose, the base images of all coset reps
   are output to a temporary file.
   ind is the current index index[bno], and cind is the no of coset reps
   found so far.
*/
  for (i=fnt;i!= -1;i=fpt[i])
  { trad[i]=tree+trptr; tree[trptr]=base[i]; trptr+=2; }
  if (dcr) opx=fopen(temp2,"w");
  if (cr)
  { op=fopen(outf3,"w");
    fprintf(op,"%4d %d%4d%4d\n",npt,index[1]-1,0,0);
  }
  sad=trptr; tree[trptr]=1; trptr++;
  for (bno=lnt;bno!= -1;bno=bpt[bno])
  { printf("bno=%d.\n",bno);
    *gorb=1; gorb[1]=base[bno];
    ind=index[bno]; csvg=svgptr[bno]; csvh=svhptr[bno]; allok=(lorbh[bno]==1);
/* if allok, then all coset reps in this link of the stabilizer chain will be
   coset reps of H in G.
*/
    if (expg) {for (i=1;i<=npt;i++) expp[i]=0; expp[base[bno]]=stop; }
    for (pt=1;pt<=npt;pt++)
    { svpt=csvg[pt];
      if (svpt!=0 && svpt!=stop)
      { ucp=lcp-1; addsvf(pt,csvg); pthere=0;
        if (expg)
        { for (i=1;i<=npt;i++) p1[i]=image(i);
          invert(p1,p2); ucp=lcp; *ucp=p1;
        }
        j=sad;
        for (i=fpt[bno];i!= -1;i=fpt[i])
        { tailad[i]=tree+j; cpad[i]=ucp+1; j+=2; }
/* In the search, for each basic coset rep g in the stabilizer chain, we have
   to try out gh as a possible coset rep of H in G, for each h in the list of
   coset reps of H[bno-1] in G[bno-1]. coind and oind record the search
   through this list.
*/
        coind=1;
        while (coind<=oind)
        { if (coind>1) advance(); ok=1;
          if (allok==0)
          for (i=1;i<=*gorb && ok;i++) if (csvh[image(gorb[i])]!=0) ok=0;
/* That was the membership test. ok true means we have found a new rep */
          if (ok)
          { if (pthere)
            { ad=fpt[bno];
              while (*tailad[ad]== *trad[ad]) ad=fpt[ad];
              tree[trptr]= *tailad[ad];
            }
            else
            { pthere=1; ad=bno; tree[trptr]=pt;
              if (expg)
              { expp[pt]=p1; p1=pst; p2=p1+npt; pst+=(2*npt); nexp++;
                if (pst>pend)
                { fprintf(stderr,"Out of perm space. Increase PSP.\n");
                  return(-1);
                }
              }
            }
            *(trad[ad]+1)=trptr; trad[ad]=tree+trptr; trptr+=2;
            for (i=fpt[ad];i!= -1;i=fpt[i])
            { *(trad[i]+1)= -1; tree[trptr]= *tailad[i];
              trad[i]=tree+trptr; trptr+=2;
            }
            cind++; tree[trptr]=cind; trptr++;
            if (trptr>trsp)
            { fprintf(stderr,"Out of tree space. Increase TRSP.\n");
              return(-1);
            }
            if (dcr) for (i=fnt;i!= -1;i=fpt[i]) fprintf(opx,"%5d",*trad[i]);
/* if cr, then we output each coset rep */
            if (cr)
            if (npt>=1000)
            { for (i=1;i<=npt;i++) fprintf(op,"%5d",image(i));fprintf(op,"\n");}
            else
            { for (i=1;i<=npt;i++) fprintf(op,"%4d",image(i));fprintf(op,"\n");}
            if (cind==ind) goto nextbno;
          } /* if (ok) */
          coind++;
        }  /* while (coind<=oind) */
        (*gorb)++; gorb[*gorb]=pt;
      }  /* if (svpt!=0 ... */
    } /* for (pt=1;pt<=npt;... */
    fprintf(stderr,"Premature end of search.\n"); return(-1);
nextbno:
    if (expg) for (i=1;i<=npt;i++) svgptr[bno][i]=expp[i];
    sad-=2; oind=ind;
  } /* main bno loop */
  if (dcr) fclose(opx); if (cr) fclose(op);
  for (i=fnt;i!= -1;i=fpt[i]) *(trad[i]+1)= -1;
  printf("Tree space used = %d.\n",trptr);
  if (expg) printf("Number of perms expanded = %d.\n",nexp);
  return(0);
}

advance()
/* This advances the element h in the search through elements gh */
{ int ad,k,*p,*q;
  ad=lnt;
  while ((k= *(tailad[ad]+1)) == -1) ad=bpt[ad];
  q=tree+k; ucp=cpad[ad]-1;
  while (ad!= -1)
  { tailad[ad]=q; cpad[ad]=ucp+1;
    if (expg) { p=svgptr[ad][*q]; if (p!=stop) {ucp++; *ucp=p; }}
    else addsvf(*q,svgptr[ad]);
    q+=2; ad=fpt[ad];
  }
}