# include <stdio.h>
extern char inf[],inf1[],outf1[];
extern short facexp,prime,exp,*rpf,*rpb,*eexpnt,**pcb,
       dim,onng,expnt[],**comptr[],*vec[],**mat[],cp[],
       *wf,*wc,**extno,**subno,chsdim,chpdim,gap;
extern FILE *fopen(),*ip,*op;

comprels()
/* Used when crel is true to compute values of relators of P in the chosen
   stable extension of M by P.
*/
{ short i,j,k,l,m,n,k1,*p,*q,*r,stabdim,*stabhom,**stabno,*orpf,
        nb,np,*rno,*covrel,**pgen,*ps,*pf,*v1,*v2,**cbmat;
  short sgn;
  stabdim=chpdim-chsdim; stabno=subno-stabdim;
  stabhom=rpf-1; rpf+=onng; orpf=rpf;
  for (i=1;i<=stabdim;i++) {stabno[i]=rpf-1; rpf+=onng; }
  k=0;
  for (i=1;i<=onng;i++) if (subno[i]==0 && extno[i]!=0)
  { k++; p=stabno[k];
    for (j=1;j<=onng;j++) p[j]=0; p[i]=1;
    for (j=1;j<=onng;j++) if ((q=subno[j])!=0)
    { expand(q,rpf,onng); if ((l=rpf[i])!=0) p[j]=prime-l; }
  }
  if (k!=stabdim) {fprintf(stderr,"stabdim error.\n"); return(-1);}
  if (gap==0 && stabdim>1) fprintf(stderr,"Basis of stable homomorphisms:\n");
  if (gap==0 && stabdim>1) for (i=1;i<=stabdim;i++)
  { for (j=1;j<=onng;j++) fprintf(stderr,"%3d",stabno[i][j]);
    fprintf(stderr,"\n");
  }
  for (i=1;i<=onng;i++) stabhom[i]=0;
  if (gap==0 && stabdim>1)
  fprintf(stderr,"Choose required stable hom as a vector in this basis!\n");
  if (gap)
  /* We just input the number of the basis of stable homs that we want */
  { scanf("%hd",&i); for (j=1;j<=onng;j++) stabhom[j]=stabno[i][j]; }
  else for (i=1;i<=stabdim;i++)
  { if (stabdim>1) scanf("%hd",rpf);
    else *rpf=1;
    if (*rpf<0)
    { for (j=1;j<=onng;j++) scanf("%hd",stabhom+j); break;}
    if (*rpf!=0)
    { p=stabhom; r=p+onng; q=stabno[i]+1;
      while (++p<=r) { *p+= *q*(*rpf); *p%=prime; q++; }
    }
  }
  printf("Chosen hom is:\n");
  for (i=1;i<=onng;i++) printf("%3d",stabhom[i]); printf("\n");
  rpf=orpf;
  strcpy(inf1,inf); strcat(inf1,"psgwds"); ip=fopen(inf1,"r");
  if (ip==0)
  { fprintf(stderr,"Cannot open file %s.\n",inf1); return(-1); }
  pgen=pcb; fscanf(ip,"%hd",&np);
  for (i=1;i<=np;i++)
  { pgen[i]=rpf; fscanf(ip,"%hd",rpf); p=rpf; rpf+=(1+ *p);
    while (++p<rpf) fscanf(ip,"%hd",p);
  }
  fclose(ip);
/* Compute matrix to change basis of module back to the original, by
   using the base change matrices output by matcalc.
*/
  strcpy(inf1,inf); strcat(inf1,"cbmats");
  if ((ip=fopen(inf1,"r"))==0)
  { fprintf(stderr,"Cannot open file %s.\n",inf1); return(-1); }
  fscanf(ip,"%hd%hd%hd",&i,&j,&k);
  if (i!=prime || j!=dim || k!=2)
  { fprintf(stderr,"Error in line 1 of %s.\n",inf1); return(-1);}
  readmat(mat[1]); readmat(mat[2]);
  *cp=2; cp[1]=2; cp[2]=1; prod(cp,mat[3]);
  inv(mat[3],mat[2]); trans(mat[2],mat[1]);
  cbmat=mat[1]; v1=mat[2][1]; v2=mat[3][1];
  fclose(ip);
  strcpy(inf1,inf); strcat(inf1,"psg.rel"); ip=fopen(inf1,"r");
  if (ip==0)
  { fprintf(stderr,"Cannot open file %s.\n",inf1); return(-1); }
  strcpy(outf1,inf); strcat(outf1,"psg.er"); op=fopen(outf1,"w");
  fscanf(ip,"%hd",&nb); rno=rpf-1; rpf+=nb;
  for (i=1;i<=nb;i++) fscanf(ip,"%hd",rno+i);
  fprintf(op,"%4d%4d\n",nb,dim);
  for (i=1;i<=nb;i++) fprintf(op,"%4d",rno[i]); fprintf(op,"\n");
  wf=rpf;
  for (i=1;i<=rno[1];i++)
  { fscanf(ip,"%hd",&l); covrel=rpb-l; p=covrel;
    while (++p<=rpb) fscanf(ip,"%hd",p);
    zero(expnt,eexpnt);
    for (j=l;j>=1;j--)
    { wc=wf-2; k=covrel[j]; k1=k/2+1; m= *pgen[k1];
      if (k%2==0) { sgn=1; ps=pgen[k1]+1; pf=ps+m-2; }
      else { sgn= -1; pf=pgen[k1]+1; ps=pf+m-2; }
      while(1)
      { wc+=2; *wc= *ps; *(wc+1)= *(ps+1)*sgn;
        if (ps==pf) break; ps+=(2*sgn);
      }
      collect(wc,wf,1);
    }
    fprintf(op,"%4d  ",l);
    for (j=1;j<=l;j++) fprintf(op,"%4d",covrel[j]); fprintf(op,"\n");
    zero(v1,v1+dim);
    for (n=1;n<=exp;n++) if ((l=expnt[n])!=0)
    { if (n<=facexp) {fprintf(stderr,"relation error. i,n,l=%d,%d,%d\n",i,n,l);
      return(-1);
    }
      for (j=1;j<=dim;j++)
      { p= *(comptr[exp+j]+n);
        if (p!=0)
        { r=p+ *p;
          while (++p<r) if ((k=stabhom[*p])!=0)
          { v1[j]+=(k*l* *(++p)); v1[j]%=prime; }
          else ++p;
        }
      }
    }
    im(v1,v2,cbmat);
    l=0; p=v2;
    while (++p<=v2+dim) if (*p !=0) l+=2;
    fprintf(op,"%4d  ",l); p=v2;
    while (++p<=v2+dim) if (*p!=0) fprintf(op,"%4d%4d",p-v2,*p);
    fprintf(op,"\n");
    if ((i==rno[1]) && (fscanf(ip,"%hd",&j)>0))
    { fprintf(op,"%4d\n",j); rno[1]+=j;}
  } /* for (i=1;i<=rno[1];... */
  return(0);
}