#include <stdio.h>
extern char  inf[],inf1[],outf1[];
extern short  facexp,tails,stage,depth,no,mng,mexp,marg,
       prime,exp,nng,class,dim,onng,*rpf,*rpb,*eexpnt,*enexpnt,**pcb,
       **opcb,**npcb,**npcb2,*nd1,*nd2,**extno,**subno,chpdim,chsdim,
       rel[],expnt[],nexpnt[],prvec[],pinv[],wt[],
       d1[],d2[],*pcptr[],**powptr[],**comptr[],
       *sspc[],*sspf[],sgen[],sex[],spgen[],spex[],
       spugen[],dpth[];
extern int rsp,wsp,ptrsp;
short  fac;
extern FILE *fopen(),*ip,*op;

ingp(inp) char inp;
/* Input a group with pcp. Normally inp=1, and file inf1 is opened inside
   procedure. inp=0 when file is already open - i.e. when act=1, and pcp is
   the output from scrun.
*/
{ short i,j,k,l,m,sum,*orpf,*orpb,**pcp,**pcpj,jump;
  if (inp)
  { ip=fopen(inf1,"r");
    if (ip==0) { fprintf(stderr,"Cannot open %s\n",inf1); return(-1); }
    fscanf(ip,"%hd%hd%hd%hd%hd%hd",&prime,&exp,&facexp,&no,&class,&m);
    if (exp>=mexp)
    { fprintf(stderr,"exp too big. Increase MEXP.\n"); return(-1); }
    if (m!=0) {fprintf(stderr,"Wrong version of nq.\n"); return(-1);}
  }
  nng=0; sum = (stage==2 || stage==4) ? exp+dim : exp;
  if (inp) for (i=1;i<=sum;i++) fscanf(ip,"%hd",wt+i);
  if (stage<2)
  { if (exp>facexp) for (i=1;i<=exp;i++) fscanf(ip,"%hd",dpth+i);
    else for (i=1;i<=exp;i++) dpth[i]=0;
    depth=dpth[exp];
  }
  for (i=1;i<=sum;i++) fscanf(ip,"%hd",d1+i);
  for (i=1;i<=sum;i++) fscanf(ip,"%hd",d2+i);
/* Now we start to read the pcp into rel, setting the commutator and power
   pointers comptr and powptr, as we go.
   Defs and powers [d[j],d[i]] and d[j]^p are only in the file for j<=no, so
   for the remaining ones we zero them.
   pcp is the current place of the pointer in array of pointers pcptr.
   opcb,pcb,npcb and npcb2 are externals used to record the values of pcp
   at various key places, for use in outgp and prnrel.
   Roughly, opcb=value of pcp at end of PCP of P itself. (=pcb when stage=0).
   pcb=value at end of definition of commutators in [d(M),P]
   (d(M)=dual of M).
   npcb=value at end. npcb2 is only used for ch1.
*/
  rpf=rel; rpb=rel+rsp-1; pcp=pcptr;
  for (i=2;i<=no;i++)
  { comptr[i]=pcp-1-tails;  m= (i>facexp) ? facexp : i-1;
    for (j=1;j<=m;j++)
    { fscanf(ip,"%hd%hd",&k,&l);
      if (l==0) *(pcp++)=0;
      else
      { *(pcp++)=rpf+1; *(rpf++)=k; *rpf=l; orpf=rpf; rpf +=(l+1);
        while ((++orpf)<rpf) fscanf(ip,"%hd",orpf);
      }
      if (tails) *(pcp++)=0;
    }
  }
  if (tails || no<facexp) for (i=no+1;i<=exp;i++)
  { comptr[i]=pcp-1-tails; m= (i>facexp) ? facexp : i-1;
    for (j=1;j<=m;j++) {*(pcp++)=0; if (tails) *(pcp++)=0; }
  }
  m= (no>=facexp) ? facexp : no;
  for (i=1;i<=m;i++)
  { powptr[i]=pcp;
    fscanf(ip,"%hd%hd",&k,&l);
    if (l==0) *(pcp++)=0;
    else
    { *(pcp++)=rpf+1; *(rpf++)=k; *rpf=l; orpf=rpf; rpf +=(l+1);
      while ((++orpf)<rpf) fscanf(ip,"%hd",orpf);
    }
    if (tails) *(pcp++)=0;
  }
  for (i=no+1;i<=facexp;i++)
  { powptr[i]=pcp; *(pcp++)=0; if (tails) *(pcp++)=0; }
  if (stage==3) no=exp;
  if (stage==0) pcb=pcp-1;
/* End of input of PCP of P */
  else
  { opcb=pcp-1; nd1=d1+exp+dim; nd2=d2+exp+dim;
    if (stage==2 || stage==4)
    { fscanf(ip,"%hd",&nng);
      for (i=1;i<=nng;i++) fscanf(ip,"%hd",nd1+i);
      for (i=1;i<=nng;i++) fscanf(ip,"%hd",nd2+i);
      if (stage==4) {jump=exp*dim; pcpj=pcp+jump;}
    }
    for (i=exp+1;i<=exp+dim;i++)
    { comptr[i]=pcp-1;
      if (stage==4) powptr[i]=pcpj-1;
      for (j=1;j<=exp;j++) if (stage==3 || stage==1) *(pcp++)=0;
      else
      { fscanf(ip,"%hd%hd",&k,&l);
        if (l==0) *(pcp++)=0;
        else
        { *(pcp++)=rpf+1; *(rpf++)=k; *rpf=l; orpf=rpf; rpf +=(l+1);
          while ((++orpf)<rpf) fscanf(ip,"%hd",orpf);
        }
        if (stage==4)
        { fscanf(ip,"%hd%hd",&k,&l);
          if (l==0) *(pcpj++)=0;
          else
          { orpf=rpb-l; *(pcpj++)=orpf; *(orpf-1)=k; *orpf=l;
            while ((++orpf)<=rpb) fscanf(ip,"%hd",orpf); rpb-=(l+2);
          }
        }
      }
    }
    pcb=pcp-1;
    if (stage!=4 && stage !=2) for (i=exp+1;i<=exp+dim;i++)
    { powptr[i]=pcp-1; for (j=1;j<=facexp;j++) *(pcp++)=0; }
    else if (stage==4) pcp+=jump;
    npcb=pcp-1; npcb2=npcb;
  }
  wsp=0; eexpnt=expnt+exp;
  if (stage) eexpnt+=dim;
  enexpnt=nexpnt+nng;
  if (stage==2)
  { extno=pcptr+ptrsp-1-nng; subno=extno-nng;
    for (i=1;i<=nng;i++) {extno[i]=0; subno[i]=0;}
    fscanf(ip,"%hd",&chpdim);
    for (i=1;i<=chpdim;i++)
    { fscanf(ip,"%hd",&j); extno[j]=rpf;
      fscanf(ip,"%hd",rpf); l= *rpf; rpf++;
      for (j=1;j<=l;j++) {fscanf(ip,"%hd",rpf); rpf++;}
    }
    fscanf(ip,"%hd",&chsdim);
    for (i=1;i<=chsdim;i++)
    { fscanf(ip,"%hd",&j); subno[j]=rpf;
      fscanf(ip,"%hd",rpf); l= *rpf; rpf++;
      for (j=1;j<=l;j++) {fscanf(ip,"%hd",rpf); rpf++;}
    }
  }
  if (pcp-pcptr>ptrsp)
  { fprintf(stderr,"Not enough pointer space. Increase PTRSP.\n"); return(-1); }
  if (rpf-rel>rsp)
  { fprintf(stderr,"Not enough space. Increase RSP.\n"); return(-1); }
  if (inp) fclose(ip);
  return(0);
}

outgp()
/* A group is output to outf1 */
{ short i,k,l,**pcp,*b,*c,*e,*f,sum,nexp,**pcpj,jump;
  char jtl,big;
  op=fopen(outf1,"w");
  sum= (stage) ? exp+dim : (tails) ? exp+nng : exp;
  big = sum>=1000;
  if (tails) no=exp;
  nexp= stage ? exp : sum;
  class=1;
  for (i=1;i<=sum;i++) if (wt[i]>class) class=wt[i];
  fprintf(op,"%3d%5d%4d%5d%4d%4d\n",prime,nexp,facexp,no,class,0);
  for (i=1;i<=sum;i++) fprintf(op,"%4d",wt[i]); fprintf(op,"\n");
  if (stage==0)
  { for (i=1;i<=sum;i++) fprintf(op,"%4d",dpth[i]); fprintf(op,"\n");}
  if (big)
  { for (i=1;i<=sum;i++) fprintf(op,"%5d",d1[i]); fprintf(op,"\n");
    for (i=1;i<=sum;i++) fprintf(op,"%5d",d2[i]); fprintf(op,"\n");
  }
  else
  { for (i=1;i<=sum;i++) fprintf(op,"%4d",d1[i]); fprintf(op,"\n");
    for (i=1;i<=sum;i++) fprintf(op,"%4d",d2[i]); fprintf(op,"\n");
  }
  jtl=0; pcp=pcptr;
  while (pcp<=pcb)
  { b= *pcp; c= tails ? *(pcp+1) : 0;
    if (b!=0) { k= *(*pcp-1); l= **pcp; e=b+l; }
    else { k=0; l=0; }
    if (c!=0) { l+= *c; f=c+ *c; }
    if (k==0 && l==0) fprintf(op,"%d %d",0,0); else fprintf(op,"%3d  %3d",k,l);
    if (b!=0)
    { if (big) while (++b<=e) fprintf(op,"%5d",*b);
      else while (++b<=e) fprintf(op,"%4d",*b);
    }
    if (c!=0)
    { if (big)
      while (++c<=f) {fprintf(op,"%5d",*c+exp); fprintf(op,"%4d",*(++c));}
      else while (++c<=f) {fprintf(op,"%4d",*c+exp); fprintf(op,"%4d",*(++c));}
    }
    if (jtl)
    { for (i=1;i<=10-l;i++) fprintf(op,"    ");
      if (l==0) fprintf(op,"     ");
      pcpj++; b= *pcpj;
      if (b!=0) {l= *b;k= *(b-1);} else {k=0;l=0;}
      e=b+l; fprintf(op,"%3d  %3d",k,l);
      while (++b<=e) fprintf(op,"%4d",*b);
    }
    fprintf(op,"\n");
    if (stage && pcp==opcb)
    { fprintf(op,"%4d\n",onng);
      if (big)
      { for (i=1;i<=onng;i++) fprintf(op,"%5d",nd1[i]); fprintf(op,"\n");
        for (i=1;i<=onng;i++) fprintf(op,"%5d",nd2[i]); fprintf(op,"\n");
      }
      else
      { for (i=1;i<=onng;i++) fprintf(op,"%4d",nd1[i]); fprintf(op,"\n");
        for (i=1;i<=onng;i++) fprintf(op,"%4d",nd2[i]); fprintf(op,"\n");
      }
      if (stage==4)  {jump=dim*exp; pcpj=pcp+jump; jtl=1;}
    }
    pcp +=(1+tails);
  }
  if (stage==2)
  { fprintf(op,"%3d\n",chpdim);
    for (i=1;i<=onng;i++) if ((b=extno[i])!=0)
    { fprintf(op,"%3d   %4d",i,*b); c=b+ *b;
      while (++b<=c) fprintf(op,"%4d",*b); fprintf(op,"\n");
    }
    fprintf(op,"%3d\n",chsdim);
    if (chsdim!=0) for (i=1;i<=onng;i++) if ((b=subno[i])!=0)
    { fprintf(op,"%3d   %4d",i,*b); c=b+ *b;
      while (++b<=c) fprintf(op,"%4d",*b); fprintf(op,"\n");
    }
  }
  fclose(op);
  return(0);
}

zero(p1,p2) short *p1,*p2;
{ while ((++p1)<=p2) *p1=0; }

setnr(p)  short *p;
/* Really a subprocedure of collect. Adjusts nexpnt by fac times string p. */
{ short *p1,*p2;
  p1=p+ *p;
  while (++p<p1)
  { p2=nexpnt+ *p; *p2+= (fac* *(++p)); *p2%=prime; }
}

collect(spc,spf,sgn)  short *spc,*spf,sgn;
/* The basic collection routine. Taken basically from Canberra NQA.
   There are several complications.
   1) The basic word is collected into an exponent vector expnt[1-exp],
      but collection also generates terms in the new generators, 1-nng.
      These are all central, and these terms come from the tails in the PCP.
      They are collected in the vector nexpnt. The procedure setnr is 
      always used for this.
   2) Negative powers can occur in the string to be collected.
   3) Generators i with i>exp are generators of d(M). These do not need to be
      collected themselves; they are only there to generatoe terms in nng.
      They are not accumulated in expnt, but dealt with in blocks, as they
      occur.
   4) Generators i of P with i>facexp form an elementary abelian group. It
      saves time (particularly for larger primes) to use this fact when
      collecting
*/
{ short gen,ex,pgen,pex,pugen,stkp,exsgn,*e,*p1,*p2,*p3,*p4,**dp,
        i,j,fpg,fpex,*p0,*pe,*q;
  char done;
  stkp=0; pgen=0; pex=0; pugen=0; exsgn=1;

recurse:
  gen= *spc; ex= exsgn*(*(spc+1)); stkp++;
  sspc[stkp]=spc; sspf[stkp]=spf; sgen[stkp]=gen; sex[stkp]=ex;
  spex[stkp]=pex; spugen[stkp]=pugen; spgen[stkp]=pgen;

loop:
  if (gen>exp)
  { p4=prvec+exp+dim; pe=prvec+exp; zero(pe,p4); done=0;
    exsgn= (spc<spf) ? -1 : 1;
/* First we collect all gens i>exp in this block as a vector. We use prvec
   for this purpose, since it is not otherwise in use at present.
   (prvec is fundamentally there for use in prnrel() .)
*/
    while(1)
    { prvec[gen]+=ex;
      if (spc==spf) {done=1; break;}
      spc-=(2*exsgn); gen= *spc; ex= exsgn* *(spc+1);
      if (gen<=exp) break;
    }
    fpg=pgen; fpex=pex;
/* Now we start to generate commutator terms. Those gens i<=facexp may
   generate new terms in prvec, so we must proceed in the correct order.
*/
    while (fpg<=facexp)
    { if (fpex>0) for (i=1;i<=fpex;i++)
      { p3=p4;
        while (p3>pe)
        { if ((j= *p3)!=0)
          { dp=comptr[p3-prvec]+fpg; p1= *dp;
            dp=powptr[p3-prvec]+fpg; p2= *dp;
            if (p2!=0 && *p2!=0) { fac=sgn*j; setnr(p2); }
            if (p1!=0 && *p1!=0)
            { p0=p1+1; q=p1+ *p1;
              while (p0<q)
              { prvec[*p0]+= (j*  *(p0+1)); prvec[*p0]%=prime; p0+=2; }
            }
          }
          p3--;
        }
      }
      fpg++; fpex=expnt[fpg];
    }
    while (fpg<=exp)
    { if (fpex>0)
      { p3=p4;
        while (p3>pe)
        { if ((j= *p3)!=0)
          { dp=comptr[p3-prvec]+fpg; p1= *dp;
            if (p1!=0 && *p1!=0) { fac=sgn*fpex*j; setnr(p1); }
          }
          p3--;
        }
      }
      fpg++; fpex=expnt[fpg];
    }
    p1=expnt+exp;
    while (++pe<=p4)
    { p1++; *p1+= *pe; *p1%=prime; if (*p1<0) *p1+=prime; }
    if (done)
    { stkp--;
      if (stkp==0) return(0);
      gen=sgen[stkp]; ex=sex[stkp];
      if (pex== -1) goto powdone;
      spc=sspc[stkp]; spf=sspf[stkp]; pugen=spugen[stkp];
      goto loop;
    }
    sgen[stkp]=gen; sex[stkp]=ex; sspc[stkp]=spc;
  } /* gen>exp */
  else if (gen>facexp)
  { p4=prvec+exp; pe=prvec+facexp; zero(pe,p4); done=0;
    exsgn= (spc<spf) ? -1 : 1;
/* For facexp<gen<=exp, we again treat generators in blocks, and first add them
   all up together in prvec.
*/
    while(1)
    { prvec[gen]+=ex;
      if (spc==spf) {done=1; break;}
      spc-=(2*exsgn); gen= *spc; ex= exsgn* *(spc+1);
      if (gen<=facexp || gen>exp) break;
    }
    fpg=pgen; fpex=pex;
    while (fpg<=facexp)
    { if (fpex>0) for (i=1;i<=fpex;i++)
      { p3= tails ? p4 : prvec+no;
        while (p3>pe)
        { if ((j= *p3)!=0)
          { dp=comptr[p3-prvec]+(1+tails)*fpg; p1= *dp; p2= *(dp+1);
            if (tails && p2!=0) { fac=sgn*j; setnr(p2); }
            if (p1!=0 && *p1!=0)
            { p0=p1+1; q=p1+ *p1;
              while (p0<q)
              { prvec[*p0]+= (j*  *(p0+1)); prvec[*p0]%=prime; p0+=2;}
            }
          }
          p3--;
        }
      }
      fpg++; fpex=expnt[fpg];
    }
/* Having generated the commutators from this block, we add it on to
   the exponent vector expnt.
*/
    p1=expnt+facexp;
    while (++pe<=p4)
    { p1++; *p1+= *pe; *p1%=prime; if (*p1<0) *p1+=prime; }
    if (done)
    { stkp--;
      if (stkp==0) return(0);
      gen=sgen[stkp]; ex=sex[stkp];
      if (pex== -1) goto powdone;
      spc=sspc[stkp]; spf=sspf[stkp]; pugen=spugen[stkp];
      goto loop;
    }
    sgen[stkp]=gen; sex[stkp]=ex; sspc[stkp]=spc;
  } /* gen>facexp */

/* Most of the remainder of the procedure follows the collect in the NQA */
  e=expnt+pgen;
  if (pgen==gen)
  { if (pugen==gen) { *e += ex; ex=1; sex[stkp]=1; } else (*e)++;
    if (*e < 0) if (gen>facexp) *e+=prime;
    else
/* Deal with negative exponent */
    { *e +=prime; dp=powptr[gen]; p1= *dp; p2= *(dp+1);
      if (tails && p2!=0) {fac= -sgn; setnr(p2); }
      if (p1!=0 && *p1!=0)
      { exsgn= -1; pex= -1; spc= p1+1; spf=spc+ *p1-2; goto recurse;}
    }
    else if (*e>=prime)
    if (gen>facexp) *e-=prime;
    else
/* Deal with exponent >= prime */
    { *e -=prime; dp=powptr[gen]; p1= *dp; p2= *(dp+1);
      if (tails && p2!=0) {fac=sgn; setnr(p2); }
      if (p1!=0 && *p1!=0)
      { exsgn=1; pex= -1; pugen=gen; spf= p1+1; spc=spf+ *p1-2;
        goto recurse;
      }
    }
  }
  else if (pex<=0)
  { if (pugen==pgen) pugen++; pgen++; pex= *(e+1); goto loop; }
  else if (pgen>facexp) {pex=0; goto loop; }
  else
  { dp=comptr[gen]+(1+tails)*pgen; p1= *dp; p2= *(dp+1);
    if (p1!=0 && *p1!=0)
    { if (pugen==pgen)
      { spugen[stkp]= -pgen;
        if (ex<0)
        if (gen>facexp) sex[stkp]=ex+prime;
        else
        { sex[stkp]=ex+prime; dp=powptr[gen]; p3= *dp; p4= *(dp+1);
          if (tails && p4!=0) { fac= -sgn; setnr(p4); }
          if (p3!=0 && *p3!=0)
          { exsgn= -1; spc= p3+1; spf=spc+ *p3-2; goto recurse;}
        }
      }
      else pugen=pgen;
    }
    else
    { if (tails && p2!=0)
      { fac = (pugen==pgen) ? ex*pex*sgn : pex*sgn; setnr(p2); }
      pex=0;
      goto loop;
    }
    spf= p1+1; spc=spf+ *p1-2; exsgn=1;
    if (tails && p2!=0) { fac=sgn; setnr(p2);}
    pex--;
    goto recurse;
  }

powdone:
  pgen=spgen[stkp]; spc=sspc[stkp]; spf=sspf[stkp];pex=spex[stkp];
  pugen=abs(spugen[stkp]);
  ex--;
  if (ex==0)
  { if (spc==spf)
    { stkp--;
      if (stkp==0) return(0);
      gen=sgen[stkp]; ex=sex[stkp];
      if (pex== -1) goto powdone;
      spc=sspc[stkp]; spf=sspf[stkp]; pugen=spugen[stkp];
      goto loop;
    }
    else
    { if (spc<spf) {spc+=2; gen= *spc; ex= -*(spc+1); }
      else { spc-=2; gen= *spc; ex= *(spc+1); }
      sgen[stkp]=gen; sspc[stkp]=spc; sex[stkp]=ex; pugen=pgen;
      goto loop;
    }
  }
  else
  { if (pugen>pgen) { pgen=pugen; pex=expnt[pgen]; }
    sex[stkp]=ex;
    goto loop;
  }
}

setpinv()
{ short i,j;
  for (i=0;i<prime;i++) pinv[i]=0;
  for (i=1;i<prime;i++) if (pinv[i]==0) for (j=1;j<prime;j++)
  if (i*j % prime ==1) { pinv[i]=j; pinv[j]=i; break; }
}

intgen(i,j) short i,j;
/* A new generator is introduced for commutator [i,j] or (if i=j) i^p */
{ short **dp,*p,*nrpb,sum;
  dp= (stage) ? comptr[i]+j : (i==j) ? powptr[i] : comptr[i]+2*j;
  if (*dp != 0 && *(*dp-1) != 0)
  { printf("Reln [%d,%d] is already a defn.\n",i,j); return(1);}
  if (stage>=2) dp=powptr[i]+j;
  if (stage==0) dp++;
  if (*dp!=0)
  { fprintf(stderr,"Reln [%d,%d] already has a tail.\n",i,j); return(1); }
  nng++; enexpnt++;
  if (nng>=mng)
  { fprintf(stderr,"Too many new gens. Increase MNG.\n"); return(-1); }
  if (stage)
  { nd1[nng]=i; nd2[nng]=j; rpb-=4; *(rpb+1)=nng; *(rpb+2)=2;
    *(rpb+3)=nng; *(rpb+4)=1; *dp=rpb+2;
  }
  else
  { sum=exp+nng; d1[sum]=i; d2[sum]=j;dpth[nng+exp]=depth;
    wt[sum]= (i==j) ? wt[i]+1  :  wt[i]+wt[j];
    rpb-=3; *(rpb+1)=2; *(rpb+2)=nng; *(rpb+3)=1; *dp=rpb+1; dp--;
    if (*dp==0) {rpf++; *dp= rpf; *rpf=0; rpf++;}
    *(*dp-1)=sum;
  }
  if (rpb+1-rpf<marg)
{  printf("Running out of space in intgen. Remains,marg,wsp=%d,%d,%d.\n",
                    rpb+1-rpf,marg,wsp);
    if (wsp>marg) bgc(); else return(-1);
  }
  return(0);
}

subrel(i,j) short i,j;
/* Value of [i,j] or i^p is computed using an associativity condition. */
{ short **dp,*p,*nrpb,sum,x,y,z;
  dp= (stage) ? comptr[i]+j : (i==j) ? powptr[i] : comptr[i]+2*j;
  if (*dp != 0 && *(*dp-1) != 0)
  { printf("Reln [%d,%d] is already a defn.\n",i,j); return(0);}
  if (stage>=2) dp=powptr[i]+j;
  if (stage==0) dp++;
  if (*dp!=0)
  { fprintf(stderr,"Reln [%d,%d] already has a tail.\n",i,j); return(-1); }
  z=0;
  for (x=nng;x>=1;x--) if ((y=nexpnt[x])!=0)
  { if (y<0) y+=prime; *rpb=y; *(rpb-1)=x; rpb-=2; z+=2; }
  *rpb=z;
  if (stage) {*dp=rpb; rpb-=2; *(rpb+1)=0; } else { *dp=rpb;rpb--;}
  if (rpb+1-rpf<marg)
{ printf("Running out of space in subrel. Remains,marg,wsp=%d,%d,%d.\n",
                    rpb+1-rpf,marg,wsp);
    if (wsp>marg) bgc(); else return(-1);
  }
  return(0);
}

assoc(g1,g2,g3)  short g1,g2,g3;
/* Associativity relation is collected. */
{ char eq12,eq23,prnt,triv; short i,l,e,*p,sum;
  sum= stage ? exp+dim: exp;
  if (g3<0) {prnt=1; g3= -g3;} else prnt=0;
  eq12= g1==g2; eq23= g2==g3 && g1!=g2;
  zero(nexpnt,enexpnt); zero(expnt,eexpnt);
  p=rpf; *(p++)=g1; *(p++)= eq12 ? prime-1 : 1; *p=g2; *(p+1)=1;
  collect(p,rpf,1);
  p=rpf;
  for (i=1;i<=sum;i++)
  { e=expnt[i]; if (e!=0) {*(p++)=i; *(p++)=e;}}
  if (p!=rpf)
  { zero(expnt,eexpnt); expnt[g3]= eq23 ? prime-1 : 1;
    collect(p-2,rpf,1);
  }
  zero(expnt,eexpnt);
  if (eq12 && g2==g3) expnt[g2]=2;
  else
  { p=rpf; *(p++)=g2; *(p++)=1; *p=g3; *(p+1)= eq23 ? prime-1 : 1;
    collect(p,rpf,-1);
  }
  p=rpf; *p=g1; *(p+1)= eq12 ? prime-1 : 1;
  collect(p,rpf,-1);
  if (prnt)
  { p=nexpnt; while (++p<=enexpnt) printf("%3d",*p); printf("\n"); }
  triv=1;
  for (i=1;i<=nng;i++) if (triv && nexpnt[i]!=0) triv=0;
  if (triv) return(0);
  else return(1);
}

prnrel()
/* The new relation in nexpnt is processed, and a new generator eleiminated
   or made redundant if necessary.
*/
{ short i,j,k,l,nl,w,x,y,ct,elno,fac,pow,len,gno,*p,*q,*eprvec,**dp,**edp;
  char sub,triv,gth,u;
  u= (stage) ? 1 : 0; pow=0; len=0;
  for (i=nng;i>=1;i--) if ((x=nexpnt[i])!=0)
  { len++;
    if (pow==0 || (stage==0 && wt[i+exp]<w)) {pow=x; gno=i; w=wt[i+exp];}
  }
  if (len==0) return(1);
/* When stage=2 and gno<=onng, we do not eliminate a generator. In this case,
   the relation is in fact a generator of the dual of H^2(P,M), and it is
   stored as such.
   Otherwise we will be eliminating gno.
*/
  if (stage==2 && gno<=onng) return(1+gno);
  sub= len>1;
  if (sub)
  { if (pow<0) pow+=prime;
    fac=pinv[pow];
    p=nexpnt;
    while (++p<=enexpnt) {*p *=fac; *p %=prime; }
    /* printf("New gen no %d is redundant.\n",gno); */
  }
  /* else printf("New gen no %d is eliminated.\n",gno); */
  eprvec=prvec+nng;  elno= (u) ? gno : exp+gno;
  if (stage>=2) { dp=pcb+1; edp= (stage>2) ? npcb2 : npcb;}
  else
  if (stage==1) { dp=opcb+facexp+1; ct=facexp+1; edp=pcb;}
  else { dp=pcptr; edp=pcb; }
/* Now we edit the PCP relations, wasting as little space as possible */
  while (dp<=edp)
  { p= *dp;
    if (p!=0)
    {p--; if (*p>elno) (*p)--; else if (*p==elno) *p=0; p++; }
    if (stage==0) {dp++; p= *dp;}
    if (p!=0)
    { l= *p; q=p+l; gth=0;
      while (++p<q)
      { if (gth) { *(p-2)= *p-1; *(p-1)= *(p+1); }
        else if (*p>gno) (*p)--;
        else if (*p==gno)
        { gth=1;
          if (sub) break;
          (**dp)-=2;
          if (**dp==0) {*dp=0; wsp+=(3+u);}
          else wsp+=(2+u);
        }
        p++;
      }
      if (gth && sub)
      { zero(prvec,eprvec); p= *dp;
        while (++p<q) {prvec[*p]= *(p+1); p++;}
        y=prvec[gno]; nl=l-2;
        for (j=1;j<=nng;j++) if (j!=gno && (k=nexpnt[j])!=0)
        { q=prvec+j; if (*q==0) nl+=2;
          *q-=y*k;
          *q%=prime;
          if (*q<0) *q+=prime;
          if (*q==0) nl-=2;
        }
        if (nl==0) {*dp=0; wsp+=(l+1+u); }
        else
        { if (nl>l)
          { if (rpb+1-rpf<marg)
     { printf("Running out of space in prnrel. Remains,marg,wsp=%d,%d,%d.\n",
                    rpb+1-rpf,marg,wsp);
              if (wsp>marg) bgc();
              else return(-1);
            }
            if (stage)
            { rpb-=(nl+2); *(rpb+1)= *(*dp-1); *dp=rpb+2; wsp+=(l+2); }
            else { rpb-=(nl+1); *dp=rpb+1; wsp+=(l+1);}
          }
          for (j=gno;j<nng;j++) prvec[j]=prvec[j+1];
          p= *dp; *p=nl;
          for (j=1;j<nng;j++) if ((y=prvec[j])>0) {*(++p)=j; *(++p)=y; }
        }
      }
    }
    dp++;
    if (stage==1)
    { if (ct==exp) {dp+=facexp; ct=facexp+1; } else ct++; }
  }
  for (i=gno;i<nng;i++) if (stage)
  { nd1[i]=nd1[i+1]; nd2[i]=nd2[i+1]; }
  else
  { j=exp+i; k=j+1;  d1[j]=d1[k]; d2[j]=d2[k]; wt[j]=wt[k];}
  nng--; enexpnt--;
  if (nng==0)
  { printf("All new generators eliminated.\n"); return(0); }
  return(1);
}

bgc()
/* Garbage collection by outputting to a temporary file and reading in again 
   Occurs after wasted space accumulates from eleiemnating generators.
*/
{ short **pcp,*p,*q,l,k,d,**epcp,ct;
  char gcfile[80];
  FILE *fopen(),*wt;
  strcpy(gcfile,inf); strcat(gcfile,"wrd.t");
  printf("Back garbage collection. wsp=%d.\n",wsp);
  fflush(stdout); fflush(stderr);
  wt=fopen(gcfile,"w");
  if (stage>=2) { pcp=pcb+1; epcp= (stage>2) ? npcb2 : npcb; d=1;}
  else
  if (stage==1) { pcp=opcb+facexp+1; ct=facexp+1; epcp=pcb; d=1;}
  else { pcp=pcptr; epcp=pcb; d=0;}
  while (pcp<=epcp)
  { if (stage==0) pcp++; if ((p=  *pcp)!=0)
    { q=p+ *p; p-=(1+d); while (++p<=q) fprintf(wt,"%d ",*p); }
    pcp++;
    if (stage==1)
    { if (ct==exp) { pcp+=facexp; ct=facexp+1; } else  ct++; }
  }
  rpb=rel+rsp-1;
  printf("Written to file.\n"); fflush(stdout);
  fclose(wt); wt=fopen(gcfile,"r");
  if (stage>=2) { pcp=pcb+1; epcp= (stage>2) ? npcb2 : npcb;}
  else
  if (stage==1) { pcp=opcb+facexp+1; ct=facexp+1; epcp=pcb; }
  else { pcp=pcptr; epcp=pcb; }
  while (pcp<=epcp)
  { if (stage==0) pcp++;
    if (*pcp!=0)
    { if (stage) { fscanf(wt,"%hd%hd",&k,&l); rpb-=(l+2); p=rpb+2; *(p-1)=k;}
      else { fscanf(wt,"%hd",&l); rpb-=(l+1); p=rpb+1;}
      *pcp=p; *p=l; q=p+l;
      while (++p<=q) fscanf(wt,"%hd",p);
    }
    pcp++;
    if (stage==1)
    { if (ct==exp) { pcp+=facexp; ct=facexp+1; } else  ct++; }
  }
  wsp=0;
  fclose(wt); unlink(gcfile);
  printf("Reread file.\n"); fflush(stdout);
}