#include <stdio.h>

extern char mult,outft[];
extern short rwd[],***scoeff[],*cst,*cend,*ocst,***coeff,
            wd1[],wd2[],wd3[],wd4[];
extern int rsp;
extern short **simcos[],pinv[],gno[],base[],lorb[],*svptr[],cord[],cp[],bno,
       invg[],**mat[],mwl2,mwdl,mlwdl,maxcos,dim,prime,*spv,**spm,ng,rwl,nb;
extern FILE *fopen(),*ip,*op;

gc()
/* Garbage collection of coefficients by writing and rereading */
{ short i,j;  short l,*p,*q;
  op=fopen(outft,"w");
  for (i=0;i<ng;i++) for (j=1;j<=maxcos;j++) if ((p=coeff[i][j])!=0)
  { l= *p; q=p+l; while (p<=q) {fprintf(op," %d",*p); p++; }
    l= *p; q=p+l; while (p<=q) {fprintf(op," %d",*p); p++; }
  }
  fclose(op);
  ip=fopen(outft,"r");
  cst=ocst;
  for (i=0;i<ng;i++) for (j=1;j<=maxcos;j++) if (coeff[i][j]!=0)
  { coeff[i][j]=cst; fscanf(ip,"%hd",cst); p=cst; q=p+ *p;
    while (++p<=q) fscanf(ip,"%hd",p); cst=q+1;
    fscanf(ip,"%hd",cst); p=cst; q=p+ *p;
    while (++p<=q) fscanf(ip,"%hd",p); cst=q+1;
  }
  fclose(ip); unlink("gctemp");
  rsp=cend+1-cst;
  printf("Performed garbage collection. rsp=%d.\n",rsp);
  return(0);
}

comp(a,dp) short **dp,*a;
/* Compress string + vector in "a" into string in *dp
   Vectors always occupy the last dim entries of words.
*/
{ short *a1,ol,nl,nl2,*p,*q,*r,*ae;
  a1=a+mwdl; ae=a1+dim; p=a1; nl2=0;
  while (++p<=ae) if (*p!=0) nl2+=2;
  ol= (*dp==0) ? 0 : **dp+ *(*dp+ **dp+1);
  nl=nl2+ *a;
  if (nl==0) {*dp=0; return(0); }
  if (nl>ol)
  { nl+=2; if (rsp<nl) gc();
    if (rsp<nl)
    { fprintf(stderr,"Out of cspace. Increase CSPACE.\n"); return(-1); }
    *dp=cst; cst+=nl; rsp-=nl;
  }
  p= *dp; q=a; r=q+ *a;
  while (q<=r) {*p= *q; p++; q++; }
  *p=nl2; q=a1;
  while (++q<=ae) if (*q!=0) {*(++p)=q-a1; *(++p)= *q; }
  return(0);
}

compb(a,dp) short **dp,*a;
/* Similar, but rewrites word at end of short space */
{ short *a1,nl,nl2,*p,*q,*r,*ae;
  a1=a+mwdl; ae=a1+dim; p=a1; nl2=0;
  while (++p<=ae) if (*p!=0) nl2+=2;
  nl=nl2+ *a;
  if (nl==0) {*dp=0; return(0);}
  nl+=2;
  if (rsp<nl) gc();
  if (rsp<nl)
  { fprintf(stderr,"Out of cspace(b). Increase CSPACE.\n"); return(-1); }
  cend-=nl; rsp-=nl; *dp=cend+1;
  p= *dp; q=a; r=q+ *a;
  while (q<=r) {*p= *q; p++; q++; }
  *p=nl2; q=a1;
  while (++q<=ae) if (*q!=0) {*(++p)=q-a1; *(++p)= *q; }
  return(0);
}

expand(a,b) short *a,*b;
/* expand string b to string + vector a */
{ short *a1,*ae,*p,*q,*r;
  a1=a+mwdl; ae=a1+dim; p=a1;
  while (++p<=ae) *p=0;
  if (b==0) { *a=0; return(0);}
  *a= *b; p=a; q=a+ *a; r=b;
  while (++p<=q) *p= *(++r);
  r++; q=r+ *r;
  while (++r<q) {a1[*r]= *(r+1); r++; }
  return(0);
}

action(a,b) short *a,*b;
/* Replace vector a by its image under action by matrices of gens in word b
*/
{ short z,*be,c,*p,*ae; short *spve,*s,*v;
  z=1; p=a; ae=a+dim;
  while (++p<=ae) if (*p!=0) {z=0; break;}
  if (z) return(0);
  spve=spv+dim; be=b+ *b;
  while (++b<=be)
  { s=spv;
    while (++s<=spve) *s=0;
    p=a;
    while (++p<=ae) if ((c= *p)!=0)
    { v=mat[*b][p-a]; s=spv;
      while (++s<=spve) {*s+=(c* *(++v)); *s%=prime; }
    }
    p=a; s=spv;
    while (++p<=ae) *p= *(++s);
  }
  return(0);
}

concat(a,b) short *a,*b;
/* Concatenate strings+vectors "a" and "b" */
{ short *a1,*bb,*be,c,*p,*q,*r,canct,min;
  if (b==0) return(0);
  a1=a+mwdl; bb=b+ *b+1; be=bb+ *bb;
  if (mult)
  { while (++bb<be) {c= *bb; bb++; a1[c]+= *bb; a1[c]%=cord[c]; }}
  else
  { action(a1,b);
    while (++bb<be) {c= *bb; bb++; a1[c]+= *bb; a1[c]%=prime;}
  }
  min= (*a<= *b) ? *a : *b;
  canct=0; p=a+ *a; q=b+1;
  while (canct<min && *p==invg[*q]) {canct++; p--; q++;}
  r=b+ *b;
  while (q<=r) *(++p)= *(q++);
  *a+= (*b-2*canct);
  if (*a>mwl2) if (reduce(a)== -1) return(-1);
  return(0);
}

concatl(b) short *b;
/* Similar, but concatenate b to the long word rwd */
{ short *a1,*bb,*be,c,*p,*q,*r,canct,min;
  if (b==0) return(0);
  a1=rwd+mlwdl; bb=b+ *b+1; be=bb+ *bb;
  if (mult)
  { while (++bb<be) {c= *bb; bb++; a1[c]+= *bb; a1[c]%=cord[c]; }}
  else
  { action(a1,b);
    while (++bb<be) {c= *bb; bb++; a1[c]+= *bb; a1[c]%=prime;}
  }
  min= (rwl<= *b) ? rwl : *b;
  canct=0; p=rwd+ rwl; q=b+1;
  while (canct<min && *p==invg[*q]) {canct++; p--; q++;}
  r=b+ *b;
  rwl+= (*b-2*canct);
  if (rwl>mlwdl)
  { fprintf(stderr,"Word too long. Increase MLWDL.\n"); return(-1);}
  while (q<=r) *(++p)= *(q++);
  return(0);
}

invwd(a,b) short *a,*b;
/* Invert the word+vector a into b */
{ short *a1,*b1,*be,*p,*q;
  a1=a+mwdl; p=a1; b1=b+mwdl; be=b1+dim; q=b1;
  while (++q<=be) *q=0;
  q=b1;
  if (mult)  {while (++q<=be) if (*(++p)!=0) *q= cord[p-a1]- *p; }
  else while (++q<=be) if (*(++p)!=0) *q= prime - *p;
  *b= *a; p=a+ *a; q=b+1;
  while (p>a) {*q=invg[*p]; p--; q++; }
  if (mult==0) action(b1,b);
  return(0);
}

ainvb(a,b,c) short *a,*b,*c;
/* Concatenate word+vector "a" with inverse of word+vector "b" and write into c
*/
{ short min,canct,n,*p,*q,*r,*c1,*ce,*a1;
  min= (*a<= *b) ? *a : *b;
  p=a+1; q=b+1; canct=0;
  while (canct<min && *p== *q) {canct++; p++; q++; }
  p= a+ *a+1; q=a+canct+1; r=c; *r= *a+ *b-2*canct;
  while (--p>=q) *(++r)= invg[*p];
  p= b+canct; q=b+ *b;
  while (++p<=q) *(++r)= *p;
  a1=a+mwdl; c1=c+mwdl; ce=c1+dim; p=c1;
  while (++p<=ce) *p=0;
  p=c1; q=a1; r= b+ mwdl;
  if (mult)
  { while (++p<=ce) if (*(++q)!=0) *p=cord[p-c1]- *q;
    p=c1;
    while (++p<=ce) if (*(++r)!=0) {*p+= *r; *p%= cord[p-c1]; }
  }
  else
  { while (++p<=ce) if (*(++q)!=0) *p=prime- *q;
    action(c1,c);
    p=c1;
    while (++p<=ce) if (*(++r)!=0) {*p+= *r; *p%= prime; }
  }
  if (*c>mwl2) if (reduce(c)== -1) return(-1);
  return(0);
}

reduce(a) short *a;
/* Attempt to reduce length of word+vector "a", using coefficients computed
   from higher values of bno
*/
{ short *q,*r,*a1,*rwd1,*rwde; short ol,nl,cos,i,j,*p,*s,*cp1,*cpe;
  ol= *a; *cp=ol; p=cp; q=a; r=a+ol; cp1=cp+mlwdl; cpe=cp1+dim;
  while (++q<=r) *(++p)= *q;
  p=cp1;
  while (++p<=cpe) *p=0;
  for (i=bno+1;i<=nb;i++) addsv(image(base[i]),svptr[i]);
  nl= *cp- ol;
  if (nl<ol)
  { *a=nl; p=cp+ *cp+1; q=a; r=a+nl;
    while (++q<=r) *q= invg[*(--p)];
    rwd1=rwd+mlwdl; rwde=rwd1+dim;
    for (i=bno+1;i<=nb;i++) if (lorb[i]>1)
    { q=rwd1;
      while (++q<=rwde) *q=0;
      cos=1; rwl=0; p=cp; s=p+ *cp;
      while (++p<=s)
      { if (concatl(scoeff[i][*p][cos])== -1) return(-1);
        cos= simcos[i][*p][cos];
      }
      if (cos!=1) {fprintf(stderr,"Reduction error.\n"); return(-1); }
      *cp= rwl; p=cp; s=p+rwl; q=rwd;
      while (++p<=s) *p= *(++q);
      q=rwd1; r=rwde; p=cp1;
      if (mult) while (++q<=r) {*(++p)+= *q; j=cord[q-rwd1]; *p%=j; }
      else while (++q<=r) {*(++p)+= *q; *p%=prime; }
    }
    p=cp1; r=rwd1;
    while (++p<=cpe) *(++r)= *p;
    a1=a+mwdl; r=rwd1; q=a1;
    if (mult) {while (++r<=rwde) {*(++q)+= *r; *q%=cord[q-a1]; }}
    else
    { action(rwd1,a);
      while (++r<=rwde) {*(++q)+= *r; *q%=prime;}
    }
  }
  if (*a>mwl2)
  { fprintf(stderr,"Reduction unsuccessful.\n"); return(-1); }
  return(0);
}

concheck()
/* Consistency check */
{ short *p,*q,ok;
  if (ainvb(wd1,wd2,wd3)== -1) return(-1);
  if (reduce(wd3)== -1) return(-1);
  ok= *wd3==0;
  if (ok)
  { p=wd3+mwdl; q=p+dim;
    while (++p<=q) if (*p!=0) {ok=0; break;}
  }
  if (ok==0)
  { fprintf(stderr,"Inconsistent relation.\n"); return(-1);}
  return(0);
}

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; }
  return(0);
}