# define SCP 100
# include <stdio.h>
extern char heqg,words,fullg,check,outf[],inf1[],inf2[],inf3[],outf0[];
extern short  perm[],sv[],cp[],actgen[],orb[],
       base[],lorbc[],lorbg[],lorbh[],order[],pno[],
       *pptr[],*svptr1[],*svptr2[],*svptr3[],scpf[],scps[],sadpt[],
       mp,mb;
extern int psp,svsp;
short npt,mxp,sth,stcom,np2,nb,npt1,cps,cpf;
char hing,cthere,hsvth;
FILE *fopen(),*ip,*op;

optprog()
{ short i,j,k,nperms,jobt,*ptr,*iptr;
  int quot;
  char c,err,opt[80];

/* This program stores up to three groups G,H and C. G must always have a s.g.
   set and Schreier Vectors svptr1. H and C may have Schreier vectors
   svptr2 and svptr3. C is not always defined.
   heqg=1 means H=G. cthere=1 means C is defined.
   Perms in G have nos. from 0 to sth-2.
   Perms in H have nos. from sth to stcomm-2 (unless H=G, when H is not stored)
   stcomm-2=np2 if C not defined.
   If C defined, perms in C have nos from stcomm to np2.
   First read inf1 = G
*/
  if ((ip=fopen(inf1,"r"))==0)
  { printf("Cannot open %s.\n",inf1); return(-1);}
  fscanf(ip,"%hd%hd%hd%hd",&npt,&nperms,&nb,&jobt);
  if (jobt<=0) { fprintf(stderr,"Wrong input format for G.\n"); return(-1); }
  npt1=npt+1;
  quot=psp/npt1; if (quot>mp) quot=mp; mxp=quot;
  np2=2*nperms-2;
  if (spacecheck()== -1) return(-1);
  if (svsp/npt<3*nb) {fprintf(stderr,"Not enough sv space.\n"); return(-1); }
  for (i=0;i<mxp;i++) pptr[i]=perm+i*(npt+1)-1;
  for (i=1;i<=nb;i++) svptr1[i]=sv+(i-1)*npt-1;
  for (i=1;i<=nb;i++) svptr2[i]=sv+(i+nb-1)*npt-1;
  for (i=1;i<=nb;i++) svptr3[i]=sv+(i+2*nb-1)*npt-1;
  readbaselo(nb,base,lorbg); readpsv(0,nb,nperms,svptr1);
  sth = (heqg) ? 0 : 2*nperms;
  fclose(ip);
/* if fullg read from gpname.inperm to determine minimal generating set for G */
  if (fullg)
  { if ((ip=fopen(inf3,"r"))==0) 
    { printf("Cannot open %s.\n",inf3); return(-1);}
    fscanf(ip,"%hd%hd",&i,&k); fclose(ip);
  }
  else k=nperms;
/* actgen=1 means we have a genuine generator of G */
  for (i=1;i<=k;i++) actgen[2*i-2]=1;
  for (i=k+1;i<=nperms;i++) actgen[2*i-2]=2;
  if (heqg==0)
/* Read inf2 = H */
  { if ((ip=fopen(inf2,"r"))==0)
    { fprintf(stderr,"Cannot open %s.\n",inf2); return(-1);}
    fscanf(ip,"%hd%hd%hd%hd",&i,&nperms,&j,&k); err=0;
    np2+=2*nperms;
    if (spacecheck()== -1) return(-1);
    if (i!=npt || (k>0 && j!=nb)) err=1;
    if (err==0)
    if (k>0)
    { for (i=1;i<=nb;i++)
      { fscanf(ip,"%hd",&j); if (j!=base[i]) {err=1; break;} }
      if (err==0)
      { for (i=1;i<=nb;i++) fscanf(ip,"%hd",lorbh+i);
        readpsv(sth,nb,nperms,svptr2); hsvth=1;
      }
      for (i=1;i<=nperms;i++) actgen[2*i+sth-2]=1;
    } else
    { seeknln(); if (j!=0) seeknln(); if (k!=0) seeknln();
      for (i=1;i<=nperms;i++)
      { j=sth+2*i-2; ptr=pptr[j]; iptr=ptr+npt1; readvec(ptr,0); seeknln();
        invert(ptr,iptr); actgen[j]=1;
        for (k=1;k<=nb;k++) if (ptr[base[k]]!=base[k]) break;
        ptr[npt1]=k;
      }
      hsvth=0;
    }
    if (err) { fprintf(stderr,"Wrong input format for H.\n"); return(-1); }
    fclose(ip);

/* Check whether H is a subgroup of G and express generators of H as words
   in gens of G if required.
*/
    hing=1;
    for (i=1;i<=nperms;i++)
    { *cp=1; cp[1]=sth+2*i-2;
      if (test(svptr1,words)==0) { hing=0; if (words==0) break; }
    }
  } else
  { sth=0; hing=1; hsvth=1; for (i=1;i<=nb;i++) lorbh[i]=lorbg[i];}
  if (hing)
  { if (check) printf("true\n"); else printf ("H is a subgroup of G.\n");}
  else
  { if (check) printf("false\n"); else printf("H is not a subgroup of G.\n");}
  if (check) return(0);
  cthere=0; stcom=sth+2*nperms; np2=stcom-2;

  while (1)
/* Now start processing according to chosen options. */
  { printf("Choose option. ('l' for list options.)\n");
   scanf("%s",opt);
   if (strcmp(opt,"ex")==0) scanf("%hd",&k);
   while ((c=getchar())!='\n');
    if (strcmp(opt,"l")==0)
    { printf("Options:\nnc    Put C=normal closure of H in G,\n");
      printf("comm  Put C=[G,H],\nint   Put C=G ^ H,\njoin  Put C=<G,H>,\n");
      printf("core  Replace H by its core w.r.t. G,\n");
      printf("oph   Output H,\nopc   Output C,\n");
      printf("hsv   Calculate s.g.set, Schreier vectors and order of H,\n");
      printf("rgh   Replace G by H and set H=G,\n");
      printf("rghc  Replace G by H and H by C,\n");
      printf("rhc   Replace H by C,\n");
      printf("ex n  Test if perm. no n of H lies in G. If so express it in gens of G,\n");
      printf("q     Quit.\n");
      printf("WARNING: core,hsv,rghc and rhc render C undefined.\n");
    }
    else if (strcmp(opt,"nc")==0)
    { if (heqg) printf("H = G. Normal closure = G.\n");
      else if (nc(0,0,sth-2)== -1) return(-1);
    }
    else if (strcmp(opt,"comm")==0) {if (comm()== -1) return(-1);}
    else if (strcmp(opt,"int")==0)
    { if (hing) printf("H <= G. Intersection = H.\n");
      else if (hsvth==0) printf("hsv must be called first.\n");
      else if (intsect(svptr1,svptr2,svptr3,stcom)== -1) return(-1);
    }
    else if (strcmp(opt,"core")==0)
    { if (heqg) printf("H = G. Core = G.\n");
      else if (hsvth==0) printf("hsv must be called first.\n");
      else if (core()== -1) return(-1);
    }
    else if (strcmp(opt,"join")==0)
    { if (hing) printf("H <= G. Join = G.\n"); else join(); }
    else if (strcmp(opt,"oph")==0)
    { if (hsvth==0) fprintf(stderr,"H has no sv. Output not done.\n");
      else if (outp(0)== -1) return(-1);
    }
    else if (strcmp(opt,"opc")==0)
    { if (cthere==0) fprintf(stderr,"No subgroup C to output.\n");
      else if (outp(1)== -1) return(-1);
    }
    else if (strcmp(opt,"hsv")==0)
    { if (hsvth) fprintf(stderr,"H has sv already.\n");
      else
      { cthere=0; np2=stcom-2; if (gp(sth,nb,lorbh,svptr2)== -1) return(-1);
        hsvth=1; stcom=np2+2;
      }
    }
    else if (strcmp(opt,"rhc")==0)
    { if (cthere==0) fprintf(stderr,"No subgroup C.\n"); else rhc(); }
    else if (strcmp(opt,"rgh")==0)
    { if (heqg) fprintf(stderr,"H = G already.\n"); else rgh(0); }
    else if (strcmp(opt,"rghc")==0)
    { if (heqg) fprintf(stderr,"H = G already.\n");
      else if (cthere==0) printf("No subgroup C.\n");
      else rgh(1);
    }
    else if (strcmp(opt,"ex")==0)
    { if (heqg) fprintf(stderr,"H=G, so this option is pointless.\n");
      else if (abs(k)>nperms || k==0)
           fprintf(stderr,"Inappropriate perm no.\n");
      else
      { j= (k>0) ? sth+2*(k-1) : sth+2*(-k-1)+1;
        *cp=1; cp[1]=j;
        if (test(svptr1,0)==0) printf("Perm is not in G.\n");
        else
        { printf("Perm is in G. It is:\n");
          for (i= *cp;i>1;i--)
          { j=cp[i]; k= (j%2==0) ? -j/2-1 : (j-1)/2+1; printf("%d",k);
            if (i==2) printf("\n"); else printf(" * ");
          }
        }
      }
    }
    else if (strcmp(opt,"q")==0) break;
    else printf("Invalid option.\n");
  } /* Main option choosing loop  */
  return(0);
}

gp(spno,sbno,lorb,svptr) short spno,sbno,*lorb,**svptr;
/* Compute s.g. set and group order of a group K (=H or C).
   Similar to program gprun.
   spno=no. of first perm in K; the last no. is the external np2.
   lorb and svptr for placing orbit lengths and Schreier vectors of K.
   sbno=no. of base points.
   Warning: It is assumed that the given base really is a base for K,
            and no new base points are added. This is the main difference
            from gprun. The answer will be false if this assumption is false.
            However, the call from intsect uses a deliberately smaller value
            of sbno, to save time.
            
*/
{ short i,j,l,m,bno,u,v,w,x,y,z,lsv;
  char trivrel,id;    float grporder;
  bno=sbno;
loop:
  *pno=0; z = (spno==0) ? sth-2 : np2;
  for (i=spno;i<=z;i+=2)
  { if (pptr[i][npt1]>=bno && actgen[i]<=bno)
    { (*pno)++; pno[*pno]=i; }
  }
  lorb[bno]=orbitsv(base[bno],svptr[bno],0);
  if (*pno!=0)
  { np2+=2; y=np2+1;
    if (spacecheck()== -1) return(-1);
    for (i=1;i<bno;i++)
    { u=base[i]; 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=1;
          for (l=bno;l<=nb;l++)
          { v=base[l]; u=image(v); if (svptr[l][u]==0) {id=0; break;}
            addsv(u,svptr[l]); pptr[np2][v]=v; pptr[y][v]=v;
          }
          if (id==0)
          { pptr[np2][npt1]=l; actgen[np2]=bno+1;
            for (m=1;m<=npt;m++) {u=image(m); pptr[np2][m]=u; pptr[y][u]=m;}
            printf("New generator fixing first %d base pts is:\n", 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("*");
            }
            bno=l; goto loop;
          }
        }
      }
    }
    np2-=2;
  }
  bno--;

  if (bno==0)
  { if (spno!=0)
    { printf("Present 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);
    }
    return(0);
  }
  else goto loop;
}

test(svp,exp) short **svp; char exp;
/* Test whether the current perm (in cp) lies in the group G, using the
   Schreier Vector svp. If exp is true, print out word for element in
   generators of G.
*/
{ short i,k,l;
  for (l=1;l<=nb;l++)
  { k=image(base[l]);
    if (svp[l][k]==0)
    { if (exp) printf("Permutation is not in group.\n"); return(0);}
    addsv(k,svp[l]);
  }
  if (exp) for (i=cp[0];i>1;i--)
  { k=cp[i]; l= (k%2==0) ? -(k+2)/2 : (k+1)/2;
    printf("%2d",l); if (i>2) printf(" *"); else printf("\n");
  }
  return(1);
}

nc(th,sno,eno)  short sno,eno,th;
/* Compute normal closure of a group K under perms from sno to eno.
   If th=0, first copy H to C, put th=stcomm and np2 the last one.
   In any case, K is then generated by perms from th to np2.
*/
{ short i,j,l,m,n,*p,*q,*r;
  if (th==0)
/* Copy H and its Schreier vector into C if necessary. */
  { p=pptr[sth]; q=pptr[stcom]; r=q;
    while (++p<=r) *(++q)= *p;
    np2=2*stcom-sth-2;
    if (spacecheck()== -1) return(-1);
    p=actgen+sth-1; q=actgen+stcom-1; r=q; while (++p<=r) *(++q)= *p;
    if (hsvth)
    { p=svptr2[1]; q=svptr3[1]; r=q; l=stcom-sth;
      while (++p<=r) *(++q) = (*p<=0) ? *p : *p+l;
      p=lorbh; q=lorbc; r=p+nb;
      while (++p<=r) *(++q)= *p;
    }
    else if (gp(stcom,nb,lorbc,svptr3)== -1) return(-1);
    th=stcom;
  }
  for (i=th;i<=np2;i+=2) if (actgen[i]==1)
  for (j=sno;j<=eno;j+=2) if (actgen[j]==1)
  { *cp=3; cp[1]=j+1; cp[2]=i; cp[3]=j;
    if (test(svptr3,0)==0)
/* Add new gen to C */
    { np2+=2; p=pptr[np2]; q=pptr[np2+1]; *cp=3;
      if (spacecheck()== -1) return(-1);
      for (n=1;n<=npt;n++) { m=image(n); p[n]=m; q[m]=n; }
      actgen[np2]=1;
      for (n=1;n<=nb;n++) if (p[base[n]]!=base[n]) break;
      p[npt1]=n; if (gp(stcom,n,lorbc,svptr3)== -1) return(-1);
    }
  }
  cthere=1;
  return(0);
}

comm()
/* Compute commutator subgroup C = [G,H] */
{ short i,j,k,*p,*q,*r,m,n,gnext,hnext;
  p=svptr3[1]; r=p+nb*npt;
  while (++p<=r) *p=0; np2=stcom-2;
  if (spacecheck()== -1) return(-1);
  for (i=1;i<=nb;i++) {svptr3[i][base[i]]= -1; lorbc[i]=1; }
  k= heqg ? stcom : sth;
/* First define C to be generated by commutators of generators */
  for (i=0;i<k;i+=2) if (actgen[i]==1)
  for (j=sth;j<stcom;j+=2) if (actgen[j]==1)
  { *cp=4; cp[1]=i+1; cp[2]=j+1; cp[3]=i; cp[4]=j;
    if (test(svptr3,0)==0)
    { np2+=2; p=pptr[np2]; q=pptr[np2+1]; *cp=4;
      if (spacecheck()== -1) return(-1);
      for (n=1;n<=npt;n++) { m=image(n); p[n]=m; q[m]=n; }
      actgen[np2]=1;
      for (n=1;n<=nb;n++) if (p[base[n]]!=base[n]) break;
      p[npt1]=n; if (gp(stcom,n,lorbc,svptr3)== -1) return(-1);
    }
  }
  if (np2<stcom) {printf("[G,H] is trivial.\n"); cthere=0; return(0); }
  gnext=stcom; hnext=stcom;
/* Now take normal closure under G and H */
  while(1)
  { if (gnext>np2) break;
    n= heqg ? stcom-2 : sth-2;
    if (nc(gnext,0,n)== -1) return(-1); gnext=np2+2;
    if (hing || hnext>np2) break;
    if (nc(hnext,sth,stcom-2)== -1) return(-1);
  }
  return(0);
}

outp(c) char c;
/* Out put a group. H if c=0, C if c=1. */
{ short nperms,i,s,e,**sv;  char sn[12];
  printf("Input name of subgroup to be output:    ");
  scanf("%s",sn); strcpy(outf,outf0); strcat(outf,sn);
  op=fopen(outf,"w");
  nperms = c ? (np2+2-stcom)/2 : (stcom-sth)/2;
  fprintf(op,"%4d%4d%4d%4d\n",npt,nperms,nb,1);
  if (c) printbaselo(nb,base,lorbc); else printbaselo(nb,base,lorbh);
  s= c ? stcom : sth;  e= c ? np2 : stcom-2;
  *pno=0; for (i=s;i<=e;i+=2) pno[++(*pno)]=i;
  sv= c ? svptr3 : heqg ? svptr1 : svptr2;
  printpsv(nb,pno,sv); fclose(op);
  return(0);
}

intsect(sv1,sv2,sv3,stint) short **sv1,**sv2,**sv3,stint;
/* Compute intersection of groups with Schreier vectors sv1 and sv2.
   This uses a becktrack search thro' the sv1 group.
   The result has S. vector sv3, and s.g. set from stint to np2.
   It is necessary for the current perm to expand in both directions, so
   we start it at position SCP defined above. Chaos would result if this
   should be too small. In general it goes from cps to cpf. in the array cp.
*/
{ short i,m,n,k,adno,adpt,bno,*p,*q,*r;
  np2=stint-2; p=sv3[1]; r=p+nb*npt; while (++p<=r) *p=0;
  if (spacecheck()== -1) return(-1);
  for (i=1;i<=nb;i++) { sv3[i][base[i]]= -1; lorbc[i]=1;}
  scps[1]=SCP; scpf[1]=scps[1]-1;
  for (bno=nb;bno>=1;bno--)
  { adno=bno; cps=scps[1]; cpf=scpf[1]; scps[bno]=cps; scpf[bno]=cpf;
    adpt=0;
    while(1)
/* First advance thro' search */
    { adpt++; while (adpt<=npt && sv1[adno][adpt]==0) adpt++;
      if (adpt>npt)
      { if (adno==bno) break;
        adno--; cps=scps[adno]; cpf=scpf[adno]; adpt=sadpt[adno];
        continue;
      }
      if (adno==bno && (sv2[bno][adpt]==0 || sv3[bno][adpt]!=0)) continue;
      addsvb(adpt,sv1[adno]);
      k=im(base[adno]);
      if (sv2[adno][k]==0) { cps=scps[adno]; continue;}
      addsvf(k,sv2[adno]); sadpt[adno]=adpt;
      if (adno==nb)
/* We have a new generator in the intersection. */
      { np2+=2;
        if (spacecheck()== -1) return(-1);
        p=pptr[np2]; q=pptr[np2+1]; p[npt1]=bno; actgen[np2]=1;
        cps=scps[bno];
        for (n=1;n<=npt;n++) {m=im(n); p[n]=m; q[m]=n; }
        if (gp(stint,bno,lorbc,sv3)== -1) return(-1);
        adno=bno; adpt=sadpt[bno]; cpf=scpf[bno]; continue;
      }
      adno++; scps[adno]=cps; scpf[adno]=cpf; adpt=0;
    }
  }
  if (np2<stint)
  {printf("Intersection is trivial.\n"); cthere=0; return(0); }
   else cthere=1; return(1);
}

im(pt)  short pt;
/* Image under current perm used in intsect */
{ short i;
  for (i=cps;i<=cpf;i++) pt=pptr[cp[i]][pt];
  return(pt);
}

addsvb(pt,sv) short pt,*sv;
/* Similar to addsv in permfns.c, but goes backwards, using current perm 
   in intsect
*/
{ short pn;
  pn=sv[pt]; while (pn!= -1)
  { cps--; cp[cps]=pn-1; pt=pptr[pn][pt]; pn=sv[pt]; }
}

addsvf(pt,sv) short pt,*sv;
/* Same but goes forwards */
{ short pn;
  pn=sv[pt]; while (pn!= -1)
  { cpf++; cp[cpf]=pn; pt=pptr[pn][pt]; pn=sv[pt]; }
}

core()
/* Replace H by its core under G */
{ short i,j,l,m,n,stint,*p,*q,*r;
  np2=stcom-2;
  if (spacecheck()== -1) return(-1);
restart:
  for (i=0;i<sth;i+=2) if (actgen[i]==1) for (j=sth;j<stcom;j+=2)
  { *cp=3; cp[1]=i+1; cp[2]=j; cp[3]=i;
    if (test(svptr2,0)==0)
/* The conjugate of the generator of H under the generator of G does not
   lie in H, so we must take an intersection.
   First compute g-1 H g and then its intersection with H
*/
    { for (j=sth;j<stcom;j+=2)
      { *cp=3; cp[2]=j; np2+=2; p=pptr[np2]; q=pptr[np2+1];
        for (n=1;n<=npt;n++) { m=image(n); p[n]=m; q[m]=n; }
        actgen[np2]=1;
        for (n=1;n<=nb;n++) if (p[base[n]]!=base[n]) break; p[npt1]=n;
      }
      if (gp(stcom,nb,lorbc,svptr3)== -1) return(-1);
      stint=np2+2;
      if ((n=intsect(svptr2,svptr3,svptr1,stint))==-1) return(-1);
      if (n==0)
/* If core is trivial, we exit after putting heqg=1 for good measure. */
      { printf("Core is trivial. Setting H=G.\n");
        if (gp(0,nb,lorbg,svptr1)== -1) return(-1);
        stcom=sth; np2=sth-2; sth=0;
        hing=1; heqg=1; cthere=0; return(0);
      }
      l=stint-sth; p=svptr1[1]; q=svptr2[1]; r=p+nb*npt;
      while (++p<=r) *(++q)= (*p<=0) ? *p : *p-l;
      for (i=1;i<=nb;i++) lorbh[i]=lorbc[i];
      p=pptr[stint]; q=pptr[sth]; r=pptr[np2]+2*npt1;
      while (++p<=r) *(++q)= *p;
      np2=sth+np2-stint; stcom=np2+2;
      for (i=sth;i<=np2;i++) actgen[i]=1; goto restart;
    }
  }
/* We have messed up svptr1 during this computation, so we recompute it */
  cthere=0; gp(0,nb,lorbg,svptr1);
  *cp=1; hing=1; for (i=sth;i<=np2;i+=2)
  { cp[1]=i; if (test(svptr1,0)==0) { hing=0; break; } }
  if (hing) printf ("H is a subgroup of G.\n");
  else printf("H is not a subgroup of G.\n");
  return(0);
}

rhc()
/* Replace H by C */
{ short i,l,*p,*q,*r;
  cthere=0;
  if (heqg) { sth=stcom; heqg=0; }
  l=stcom-sth; p=svptr3[1]; q=svptr2[1]; r=p+nb*npt;
  while (++p<=r) *(++q)= (*p<=0) ? *p : *p-l;
  for (i=1;i<=nb;i++) lorbh[i]=lorbc[i];
  if (heqg==0)
  { p=pptr[stcom]; q=pptr[sth]; r=pptr[np2]+2*npt1;
    while (++p<=r) *(++q)= *p;
    p=actgen+stcom-1; q=actgen+sth-1; r=actgen+np2+1;
    while (++p<=r) *(++q)= *p;
  }
  np2=sth+np2-stcom; stcom=np2+2; hsvth=1;
  *cp=1; hing=1; for (i=sth;i<=np2;i+=2)
  { cp[1]=i; if (test(svptr1,0)==0) { hing=0; break; } }
  if (hing) printf ("H is a subgroup of G.\n");
  else printf("H is not a subgroup of G.\n");
}

rgh(c) char c;
/* Replace G by H, and is c=1 also H by C */
{ short i,l,*p,*q,*r;
  l=sth; p=svptr2[1]; q=svptr1[1]; r=p+nb*npt;
  while (++p<=r) *(++q)= (*p<=0) ? *p : *p-l;
  for (i=1;i<=nb;i++) lorbg[i]=lorbh[i];
  p=pptr[sth]; q=pptr[0]; r=pptr[stcom];
  while (++p<=r) *(++q)= *p;
  p=actgen+sth-1; q=actgen-1; r=actgen+stcom-1;
  while (++p<=r) *(++q)= *p;
  if (c)
  { sth=stcom-sth; rhc(); }  else
  { stcom-=sth; np2=stcom-2; sth=0; cthere=0; hsvth=1; hing=1; heqg=1;}
}

join()
/* Compute C=<G,H> */
{ short i,*p,*q,*r;
  np2=stcom-2;
  for (i=0;i<stcom;i+=2) if (actgen[i]==1)
  { if (spacecheck()== -1) return(-1);
    np2+=2; p=pptr[i]; q=pptr[np2]; r=p+2*npt1;
    if (spacecheck()== -1) return(-1);
    while (++p<=r) *(++q) = *p; actgen[np2]=1;
  }
  gp(stcom,nb,lorbc,svptr3); cthere=1;
  return(0);
}

spacecheck()
{ if (np2+1>=mxp) 
  { fprintf(stderr,"Out of perm space. Increase PSP (or MP).\n"); return(-1);}
  return(0);
}