# include <stdio.h>
extern char wrd,nt,isbase,inf[],outf1[],outf2[],fixed[];
extern int  perm[],sv[],cp[],actgen[],orb[],
       base[],lorb[],order[],pno[], *pptr[],*svptr[],
       mp,mb,mnpt;
extern int psp,svsp;
int npt;
FILE *fopen(),*ip,*op;

/* The data structures in this program are similar to most permutation group
   programs. Permutations are numbered 0,1,2,3,... (where 2x+1 is usually the
   inverse of 2x) and stored in the array perm. Perm no x is pointed to by
   pptr[x]. npt=no. of points. Usually pptr[i][npt+1] gives the number of the
   group in the stabilizer chain in which the perm lies. So this no is i if the
   perm fixes the first i-1 base points. nb=no of base points.
   The base and the lengths of the orbits in the stab chain are stored in base
   and lorb.
   Schreier vectors are stored in sv, and pointed to by svptr[i], i=1,..,nb.
   pno is a list of *pno (=pno[0]) perm nos.
   cp is a similar list of length *cp, but it is used to represent the product
   of the perms  cp[1]cp[2]..cp[*cp]. This product can be evaluated by the
   procedure image in permfns.c
*/
gpprog()
{ int i,j,k,l,m,n,nperms,nb,jobt,np2,seek,cord,ocord,given,ordknown,trivrel,
        lsv,u,v,w,x,y,z,mxp,mnb,bpt,bno,id;
  int quot;
  float grporder;

  if ((ip=fopen(inf,"r"))== 0)
  { fprintf(stderr,"Cannot open file %s\n",inf); return(-1); }
  fscanf(ip,"%d%d%d%d",&npt,&nperms,&nb,&jobt);
  if (npt>mnpt) { fprintf(stderr,"npt too big. Increase NPT.\n"); return(-1); }
  if (jobt>0) { fprintf(stderr,"Wrong input format.\n"); return(-1); }
  quot=psp/(npt+1); if (quot>mp) quot=mp; mxp=quot;
  quot=svsp/npt; if (quot>mb) quot=mb; mnb=quot; if (mnb>=npt) mnb=npt-1;
  if (nb>mnb)
  { fprintf(stderr,"nb to big. Increase SVSP (or MB).\n"); return(-1); }
/* pptr[i] is the i th permutation, svptr[i] is the i the Schreier vector. */
  for (i=0;i<mxp;i++) pptr[i]=perm+i*(npt+1)-1;
  for (i=1;i<=mnb;i++) svptr[i]=sv+(i-1)*npt-1;

/* Next we read any base points and orbit lengths */
  if (nb!=0)
  { for (i=1;i<=npt;i++) orb[i]=0;
    for (i=1;i<=nb;i++)
    { fscanf(ip,"%d",base+i); if (orb[base[i]]!=0)
      { fprintf(stderr,"Repeated base element.\n"); return(-1); }
      orb[base[i]]=1;
    }
  }
  if (jobt<0)
  { jobt= -jobt; if (jobt>nb) {fprintf(stderr,"jobt too big.\n"); return(-1); }
    seek=jobt; given=jobt; ordknown=1;
    for (i=1;i<=jobt;i++)  fscanf(ip,"%d",order+i);
  }
  else ordknown=0;

/* Now we read the generating permutations */
  np2=2*nperms-2;
  for (i=0;i<=np2;i+=2)
  { k=i/2 +1; j=readperm(pptr[i]);
    if (j==2)
    { fprintf(stderr,"Generator no %d is not a permutation.\n",k); return(-1); }
    if (j==1)
    { fprintf(stderr,"Generator no %d is the identity.\n",k);
      if (nt) return(-1);
      nperms-=1; np2-=2; i-=2;
    }
    else
    { invert(pptr[i],pptr[i+1]); actgen[i]=1; x=1;
      for (z=base[x];x<=nb && pptr[i][z]==z;z=base[x]) x++;
      pptr[i][npt+1]=x;
      if (x>nb)
      { if (isbase)
        { fprintf(stderr,"Given base is not a base!\n"); return(-1);}
        nb++;  for (z=1;pptr[i][z]==z;z++); base[nb]=z;
        printf("New base element no %d is %d.\n",nb,z);
      }
      if (x==1) printf("Generator no %d moves first base point.\n",k);
      else printf("Generator no %d fixes first %d base point(s).\n",k,x-1);
    }
    if (nperms==0) { fprintf(stderr,"Trivial group!\n"); return(-1); }
  }
  fclose(ip);

  if (wrd) {op=fopen(outf2,"w"); fprintf(op,"%4d\n",nperms); }
  bno=nb; for (i=0;i<=mnb;i++) lorb[i]=0;

/* Now the main algorithm begins */
loop:
  *pno=0; if (ordknown) ocord= (bno==nb) ? 1 : order[bno+1];
/* We make a list of the perm nos that fix the first bno-1 base pts */
  for (i=0;i<=np2;i+=2)
  { if (pptr[i][npt+1]>=bno && actgen[i]<=bno)
    { (*pno)++; pno[*pno]=i; }
  }
  lorb[bno]=orbitsv(base[bno],svptr[bno],0);
  if (ordknown)
  { cord= (bno>=given) ? ocord*lorb[bno] : lorb[bno];
    if (bno==seek && cord==order[bno])
    { seek--; printf("Order is now as given for bno = %d.\n",bno);
      bno--; if (bno==0) goto foundorder; goto loop;
    }
  }
/* Now we start generating the Schreier generators that fix the firat bno
   base points, test them for membership, and add them as new generators
   if necessary.
*/
  if (*pno!=0)
  { nperms++; np2+=2; y=np2+1;
    if (y>=mxp)
    {fprintf(stderr,"Out of space. Increase PSP (or MP).\n"); return(-1); }
    for (i=1;i<=npt;i++) fixed[i]=0;
    for (i=1;i<bno;i++)
    { u=base[i]; fixed[u]=1; pptr[np2][u]=u; pptr[y][u]=u; }
    for (i=1;i<=lorb[bno];i++)
    { *cp=0; addsv(orb[i],svptr[bno]);
      for (w=1,x= *cp;w<=x;w++,x--)
      { if (w==x) cp[w]-=1; else {u=cp[w]; cp[w]=cp[x]-1; cp[x]=u-1;}}
      lsv= *cp;
      for (j=1;j<=*pno;j++)
      { *cp=lsv; u=pno[j]+1;
        trivrel = (*cp>0) ? cp[*cp]==(u-1) : 0;
        if (trivrel==0)
        { (*cp)++; cp[*cp]=u; id=0;
          for (l=bno;l<=nb;l++) fixed[base[l]]=0;
          for (l=bno;l<=nb;l++)
          { v=base[l]; u=image(v);
            if (svptr[l][u]==0) goto newgen;
            addsv(u,svptr[l]); pptr[np2][v]=v; pptr[y][v]=v; fixed[v]=1;
          }
          l=nb+1; id=1;
newgen:   if (isbase==0)
          for (k=1;k<=npt;k++)
          { if (fixed[k]==0)
            { u=image(k); pptr[np2][k]=u; pptr[y][u]=k;
              if (id && k!=u)
              { id=0; nb++;
  if (nb>mnb) { fprintf(stderr,"nb to big. Increase SVSP (or MB).\n"); return(-1); }
                base[nb]=k;
                printf("New base point no %d is %d.\n",nb,k);
              }
            }
          }
          if (id==0)
          { pptr[np2][npt+1]=l; actgen[np2]=bno+1;
            if (isbase) for (k=1;k<=npt;k++)
            { u=image(k); pptr[np2][k]=u; pptr[y][u]=k;}
            printf("New generator no %d, fixing first %d base pts is:\n",
                    nperms,l-1);
            for (m=1;m<=*cp;m++)
            { z=cp[m]; y=z/2; x= (z==2*y) ? y+1 : -(y+1);
              printf("%3d",x); if (m== *cp) printf("\n"); else printf("*");
            }
            if (wrd)
            {for (m=0;m<=*cp;m++) fprintf(op,"%4d",cp[m]);fprintf(op,"\n");}
            bno=l; goto loop;
          }
        }
      }
    }
    nperms--; np2-=2;
  }
  if (ordknown && bno>given) order[bno]=cord;
  bno--;

foundorder:
  if (bno==0)
  { printf("The order of the group is:\n");
    for (i=1;i<=nb;i++)
    { printf("%3d",lorb[i]);
      if (i==nb) printf(" =\n"); else printf("*");
    }
    grporder=1; for (i=1;i<=nb;i++) grporder *= lorb[i];
    printf("%g\n",grporder);
  }
  else goto loop;
  if (wrd) { fprintf(op,"%d\n",-1); fclose(op); }

/* Now we output the generating perms and Schreier vectors */
  op=fopen(outf1,"w");
  fprintf(op,"%4d%4d%4d%4d\n",npt,nperms,nb,1);
  printbaselo(nb,base,lorb);  *pno=0;
  for (i=1;i<=nperms;i++) {(*pno)++; pno[*pno]=2*(i-1);}
  printpsv(nb,pno,svptr);
  fclose(op); return(0);
}