#include <stdio.h>
extern char  inf1[],outf[], outfm[], gap;
extern short  intexp,mexp,mng,wksp,
       prime,exp,nng,class,*rpf,*rpb,*eexpnt,*enexpnt,**pcb,mnng,mord,
       rel[],expnt[],nexpnt[],cord[],wt[],
       d1[],d2[],*pcptr[],**powptr[],**comptr[],
       *sspc[],*sspf[],sgen[],sex[],spgen[],spex[],
       spugen[],*tlintg[];
extern int ptrsp,rsp;
short  fac;
FILE *fopen(),*ip,*op;

ingp()
/* The PCP for P (output from pcrun) is read in, and all of the associated
   pointers are set. If -a is true, then the input comes from a preceding run
   of nqrun, and tails are also read in.
*/
{ short i,j,k,l,m,x,y,no,*orpf,*orpb,**pcp;  char tails;
  if ((ip=fopen(inf1,"r"))==0)
  { fprintf(stderr,"Cannot open %s.\n",inf1); return(-1); }
  fscanf(ip,"%hd%hd%hd%hd%hd%hd",&prime,&exp,&nng,&no,&class,&m);
  if (exp>mexp) {fprintf(stderr,"exp too big. Increase MEXP.\n"); return(-1);}
  if (nng>mng)  {fprintf(stderr,"nng too big. Increase MNG.\n"); return(-1);}
  if (m!=1) {fprintf(stderr,"Wrong version of nq.\n"); return(-1); }
  tails = exp==no; if (tails==0) nng=0;
  for (i=1;i<=exp;i++) fscanf(ip,"%hd",wt+i);
  for (i=1;i<=exp+nng;i++) fscanf(ip,"%hd",d1+i);
  for (i=1;i<=exp+nng;i++) fscanf(ip,"%hd",d2+i);
  for (i=1;i<=nng;i++) fscanf(ip,"%hd",cord+i);
  rpf=rel; rpb=rel+rsp-1; pcp=pcptr;
  for (i=2;i<=no;i++)
  { comptr[i]=pcp-2;
    for (j=1;j<i;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)
      { fscanf(ip,"%hd%hd",&k,&l);
        if (l==0) *(pcp++)=0;
        else
        { orpb=rpb-nng; *(pcp++)=orpb; *orpb=k;
          while (rpb>orpb) *(rpb--)=0;
          l/=2; for (k=1;k<=l;k++) {fscanf(ip,"%hd%hd",&x,&y); rpb[x]=y;}
          rpb--;
        }
      }
      else *(pcp++)=0;
    }
  }
  for (i=no+1;i<=exp;i++)
  { comptr[i]=pcp-2; for (j=1;j<i;j++) {*(pcp++)=0; *(pcp++)=0; } }
  for (i=1;i<=no;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)
    { fscanf(ip,"%hd%hd",&k,&l);
      if (l==0) *(pcp++)=0;
      else
      { orpb=rpb-nng; *(pcp++)=orpb; *orpb=k;
        while (rpb>orpb) *(rpb--)=0;
        l/=2;
        for (k=1;k<=l;k++) {fscanf(ip,"%hd%hd",&x,&y); rpb[x]=y;}
        rpb--;
      }
    }
    else *(pcp++)=0;
  }
  for (i=no+1;i<=exp;i++)
  { powptr[i]=pcp; *(pcp++)=0; *(pcp++)=0; }
  pcb=pcp-1;
  if (pcb-pcptr>=ptrsp)
  { fprintf(stderr,"Out of ptr space. Increase PTRSP.\n"); return(-1); }
  if (rpb-rpf<wksp)
  { fprintf(stderr,"Out of space. Increase RSP.\n"); return(-1); }
  fclose(ip); return(0);
}

outgp()
/* The PCP is output, together with tails */
{ short i,k,l,**pcp,*b,*e,*c;
  op=fopen(outf,"w");
  fprintf(op,"%4d%4d%4d%4d%4d%4d\n",prime,exp,nng,exp,class,1);
  for (i=1;i<=exp;i++) fprintf(op,"%4d",wt[i]); fprintf(op,"\n");
  for (i=1;i<=exp+nng;i++) fprintf(op,"%4d",d1[i]); fprintf(op,"\n");
  for (i=1;i<=exp+nng;i++) fprintf(op,"%4d",d2[i]); fprintf(op,"\n");
  if (nng>0)
  { for (i=1;i<=nng;i++) fprintf(op,"%4d",cord[i]); fprintf(op,"\n");}
  pcp=pcptr;
  while (pcp<pcb)
  { if (*pcp==0) {fprintf(op,"%3d  %3d",0,0); l=0; }
    else
    { l= **pcp; fprintf(op,"%3d  %3d",*(*pcp-1),l);
      b= *pcp; e=b+l;
      while (++b<=e) fprintf(op,"%4d",*b);
    }
    if (l<10) {l=10-l; for (k=1;k<=l;k++) fprintf(op,"    "); }
    pcp++;
    if (*pcp==0) fprintf(op,"%3d  %3d",0,0);
    else
    { fprintf(op,"%3d  ",**pcp); b= *pcp; e=b+nng; c=cord+1; l=0;
      while (++b<=e)
      { while (*b<0) (*b)+= *c; if (*b!=0) l+=2; c++; }
      fprintf(op,"%3d",l); b= *pcp;
      while (++b<=e) if (*b!=0) fprintf(op,"%4d%4d",b-(*pcp),*b);
    }
    fprintf(op,"\n"); pcp++;
  }
  fclose(op);
}

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

setnr(p)  short *p;
/* This is called only by collect, which follows */
{ short *p1,*p2,*c;
  p1=nexpnt; p2=p1+nng; c=cord+1;
  while ((++p1)<=p2)
  { p++;  *p1 += fac*(*p); *p1 %= *c; c++; }
}

collect(spc,spf,sgn)  short *spc,*spf,sgn;
/* The basic collecting algorithm. The word to be collected is stored as a
   gen-pow string. spc points to the last gen, and spf to the first, in this
   string. The collected word is put into expnt, as a vector, and its tail
   is added to nexpnt (or subtracted if sgn= -1). The algorithm is similar to
   that used in the Canberra NQA, but negative exponents are allowed.
*/
{ short gen,ex,pgen,pex,pugen,stkp,exsgn,*e,*p1,*p2,*p3,*p4,**dp;
  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:
  e=expnt+pgen;
  if (pgen==gen)
  { if (pugen==gen) { *e += ex; ex=1; sex[stkp]=1; } else (*e)++;
    if (*e < 0)
    { *e +=prime; dp=powptr[gen]; p1= *dp; p2= *(dp+1);
      if (p2!=0) {fac= -sgn; setnr(p2); }
      if (p1!=0)
      { exsgn= -1; pex= -1; spc= p1+1; spf=spc+ *p1-2; goto recurse;}
    }
    else if (*e>=prime)
    { *e -=prime; dp=powptr[gen]; p1= *dp; p2= *(dp+1);
      if (p2!=0) {fac=sgn; setnr(p2); }
      if (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
  { dp=comptr[gen]+2*pgen; p1= *dp; p2= *(dp+1);
    if (p1==0)
    { if (p2!=0)
      { fac = (pugen==pgen) ? ex*pex*sgn : pex*sgn; setnr(p2); }
       pex=0; goto loop;
    }
    else
    { if (pugen==pgen)
      { spugen[stkp]= -pgen;
        if (ex<0)
        { sex[stkp]=ex+prime; dp=powptr[gen]; p3= *dp; p4= *(dp+1);
          if (p4!=0) { fac= -sgn; setnr(p4); }
          if (p3!=0)
          { exsgn= -1; spc= p3+1; spf=spc+ *p3-2; goto recurse;}
        }
      }
      else pugen=pgen;
    }
    spf= p1+1; spc=spf+ *p1-2; exsgn=1;
    if (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;
  }
}

intgen(i,j) short i,j;
/* A new generator is introduced into the tail of the relation [i,j], or
   i^p if i=j.
*/
{ short **dp,*p,*nrpb,sum;
  if (i==j) dp=powptr[i];  else dp=comptr[i]+2*j;
  if (*dp != 0 && *(*dp-1) != 0)
  { printf("Reln [%d,%d] is already a defn.\n",i,j);  return(0);}
  p= *(dp+1);
  if (p!=0)
  { fprintf(stderr,"Reln [%d,%d] already has a tail.\n",i,j); return(-1); }
  nng++; sum=exp+nng;
  if (sum>mng)
  { fprintf(stderr,"Too many newgens. Increase MNG.\n"); return(-1);}
  d1[sum]=i; d2[sum]=j; cord[nng]=mord; enexpnt++;
  nrpb=rpb-mnng;
  if (nrpb-rpf<wksp)
  { fprintf(stderr,"Out of space. Increase RSP.\n"); return(-1);}
  zero(nrpb,rpb); *(dp+1)=nrpb;
  nrpb[nng]=1; *nrpb=sum; rpb=nrpb-1;
  return(0);
}

subrel(i,j) short i,j;
/* The element currently in nexpnt, is introduced as the tail of the relation
   [i,j]
*/
{ short **dp,*p,*nrpb,sum;
  if (i==j) dp=powptr[i];  else dp=comptr[i]+2*j;
  if (*dp != 0 && *(*dp-1) != 0)
  { printf("Reln [%d,%d] is already a defn.\n",i,j); return(0);}
  p= *(dp+1);
  if (p!=0)
  { fprintf(stderr,"Reln [%d,%d] already has a tail.\n",i,j); return(-1); }
  nrpb=rpb-mnng;
  if (nrpb-rpf<wksp)
  { fprintf(stderr,"Out of space. Increase RSP.\n"); return(-1);}
  zero(nrpb,rpb); *(dp+1)=nrpb; *nrpb=0; rpb=nrpb-1;
  for (i=1;i<=nng;i++) nrpb[i]=nexpnt[i];
  return(0);
}

assoc(g1,g2,g3)  short g1,g2,g3;
/* The difference between the tails of the elements (g1g2)g3 and g1(g2g3) is
   computed into nexpnt. Thisdifference will of course have to be made trivial
   later. If g1=g2, we use (g1^p)g2 and g1(g1^(p-1)g2), and similarly if
   g2=g3.
*/
{ char eq12,eq23,prnt,triv; short i,l,e,*p;
  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<=exp;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(corrtl) char corrtl;
/* The element stored in nexpnt is made trivial. This involves eliminating or
   reducing the order of one or more new generators. corrtl is true when -a is
   set, and images of generators in the Sylow intersection have their own tails,
   which require updating.
   When orders involve larger powers of primes, sometimes more than one run
   is necessary to complete the reductions. Hence label "restart".
*/
{ short i,j,k,len,hp,ct,fac,x,y,a,b,c,hpi,gno,*p,**dp,term[30];
   char elim,sub,rep,triv;
restart:
  hp=exp/2; len=0; rep=0;
  for (i=nng;i>=1;i--)
  { hpi=0; x=nexpnt[i];
    if (x!=0)
    { len++;
      while (x%prime==0) {hpi++; x/=prime; }
      if (hpi<hp || (hpi==hp && cord[i]<b))
      { gno=i; hp=hpi; a=x; b=cord[i]; }
    }
  }
  if (len==0) return(0);
  elim= hp==0; sub= len>1;
  x=1;
  for (i=1;i<=hp;i++) x*=prime;
  if (sub)
  { ct=0; c=cord[gno]/x;
    if (a<0) a+=cord[gno];
    while (a!=1) { ct++; term[ct]=c/a; b=c%a; c=a; a=b; }
    if (ct==0) fac=1;
    else
    { term[ct+1]=1; term[ct+2]=0; p=term+ct+1;
      while (--p>term) *p = -(*p) * *(p+1) + *(p+2);
      fac= term[1];
    }
    p=nexpnt;
    while (++p<=enexpnt) {*p *=fac; *p /=x; *p %=cord[p-nexpnt]; }
    dp=pcptr;
    while (++dp<=pcb)
    { p= *dp; if (p!=0)
      { y=p[gno];
        if (y!=0) for (j=1;j<=nng;j++) if (j==gno) continue;
        else if ((k=nexpnt[j])!=0) {p[j]-=(y*k); p[j]%=cord[j]; }
      }
      dp++;
    }
    if (corrtl) for (i=3;i<=intexp;i++) if ((p=tlintg[i])!=0)
    { y=p[gno];
      if (y!=0) for (j=1;j<=nng;j++) if (j==gno) continue;
      else
      if ((k=nexpnt[j])!=0) {p[j]-=(y*k); p[j]%=cord[j]; }
    }
    a=cord[gno]; p=nexpnt;
    while (++p<=enexpnt)
    { *p *=a; *p %=cord[p-nexpnt]; if (*p!=0) rep=1; }
  }
  if (elim)
  { printf("New gen no %d is eliminated.\n",gno);
    dp=pcptr;
    while (++dp<=pcb)
    { p= *dp;
      if (p!=0)
      { if (*p>(exp+gno)) (*p)--; else if (*p==(exp+gno)) *p=0;
        triv=1;
        for (j=1;j<gno;j++) if (p[j]!=0) triv=0;
        for (j=gno;j<nng;j++)
        { p[j]=p[j+1];
          if (triv && p[j]!=0) triv=0;
        }
        p[nng]=0;
        if (triv) *dp=0;
      }
      dp++;
    }
    for (j=gno;j<nng;j++)
    { cord[j]=cord[j+1]; d1[exp+j]=d1[exp+j+1]; d2[exp+j]=d2[exp+j+1];
      nexpnt[j]=nexpnt[j+1];
    }
    if (corrtl) for (i=3;i<=intexp;i++) if ((p=tlintg[i])!=0)
    for (j=gno;j<nng;j++) p[j]=p[j+1];
    nng--; mnng--; enexpnt--;
    if (nng==0)
    { if (gap)
      { op=fopen(outfm,"w"); fprintf(op,"COHOMOLO.Multiplier:=[];\n");
        fclose(op);
        printf("All new generators eliminated. Multiplier is trivial.\n");
      }
      else fprintf(stderr,
                  "All new generators eliminated. Multiplier is trivial.\n");
      return(-1);
    }
  }
  else
  { printf("New order of new gen no %d is %d.\n",gno,x);
    cord[gno]=x; dp=pcptr;
    while (++dp<=pcb)
    { p= *dp; if (p!=0) p[gno]%=x; dp++; }
    if (corrtl) for (i=3;i<=intexp;i++) if ((p=tlintg[i])!=0) p[gno] %= x;
  }
  if (rep) goto restart;
  return(1);
}