#include "defs.h"

extern char  mult, inf0[], inf1[], inf2[], outf[], outft[], inf3[];
extern short csp[], *cpsp[], **cdpsp[], ***coeff[];
extern int   space, cspace;
extern short sp[], **mat[], *psp[], **imcos[], **cpco[], lorb[], wdl, ptrsp,
    cdptrsp, mb;
short *cst, **cpst, ***cdpst, **svptr, *cp, *rel;
short *spst, **pspst, **pptr, **cpptr, npt, nb, nph, nph2, npg, npg2, rno,
    orno, coh_index, *invg;

void seeknln(FILE * ip)
{
  while (getc(ip) != '\n')
    ;
}

/* This program differs from most other permutation programs in that perms are
   all stored in the single array sp. Schreier vectors are stored in the short
   array csp. This assumes that there are at most 127 perms, which will
   probably always be the case. The first half of the program (crprog1)
   involves relatively little short space, so it is assumed that there is
   enough In the first half, coefficients as words in H generated by the
   relators of G are computed using the output of conrun.
*/
int crprog1(void)
{
  short *pc, *qc, ex, neg;
  int    x;
  FILE * ip, *op;
  short  i, j, k, l, m, n, cl, rl, *p, ocl, im, pt, pt1, pn, ipt;
  if ((ip = fopen(inf2, "r")) == 0) {
    fprintf(stderr, "Cannot open %s.\n", inf2);
    return (-1);
  }
  fscanf(ip, "%hd%hd%hd%hd", &npt, &nph, &nb, &k);
  if (nb >= mb) {
    fprintf(stderr, "nb too big. Increase MB.\n");
    return (-1);
  }
  if (k <= 2) {
    fprintf(stderr, "inf2 has illegal format.\n");
    return (-1);
  }
  seeknln(ip);
  seeknln(ip);
  for (i = 1; i <= nb; i++)
    fscanf(ip, "%hd", lorb + i);
  pptr = psp - 1;
  pspst = psp + nph;
  svptr = cpsp - 1;
  cpst = cpsp + nb;
  invg = sp;
  nph2 = 2 * nph;
  spst = sp + nph2;
  for (i = 1; i <= nph; i++) {
    pptr[i] = spst + (i - 1) * npt - 1;
    p = pptr[i];
    for (j = 1; j <= npt; j++) {
      fscanf(ip, "%hd", &k);
      p[k] = j;
    }
    invg[2 * i - 2] = 2 * i - 1;
    invg[2 * i - 1] = 2 * i - 2;
    seeknln(ip);
  }
  spst += (npt * nph);
  for (i = 1; i <= nb; i++) {
    svptr[i] = csp + (i - 1) * npt - 1;
    pc = svptr[i];
    for (j = 1; j <= npt; j++) {
      fscanf(ip, "%hd", &k);
      pc[j] = k > 0 ? (k + 1) / 2 : k;
    }
  }
  cp = csp + nb * npt;
  rel = cp + wdl;
  fclose(ip);
  /* It is assumed that there is enough space so far. */

  strcpy(inf1, inf0);
  strcat(inf1, ".cp");
  if ((ip = fopen(inf1, "r")) == 0) {
    fprintf(stderr, "Cannot open %s.\n", inf1);
    return (-1);
  }
  fscanf(ip, "%hd%hd", &coh_index, &npg);
  npg2 = 2 * npg - 1;
  x = (nb + 2) * npg * coh_index;
  if (spst + x - sp > space) {
    fprintf(stderr, "Not enough space. Increase SPACE.\n");
    return (-1);
  }
  if (pspst + (nb + 2) * npg - psp > ptrsp) {
    fprintf(stderr, "Not enough ptrsp. Increase PTRSP.\n");
    return (-1);
  }
  cpptr = pspst;
  pspst += 2 * npg;
  for (i = 0; i <= npg2; i++) {
    cpptr[i] = spst - 1;
    spst += coh_index;
  }
  for (i = 1; i <= nb; i++) {
    cpco[i] = pspst - 1;
    pspst += npg;
    for (j = 1; j <= npg; j++) {
      cpco[i][j] = spst - 1;
      spst += coh_index;
    }
  }
  for (i = 1; i <= npg; i++) {
    j = 2 * i - 2;
    k = j + 1;
    for (l = 1; l <= coh_index; l++) {
      for (m = 1; m <= nb; m++)
        fscanf(ip, "%hd", cpco[m][i] + l);
      fscanf(ip, "%hd", &m);
      cpptr[j][l] = m;
      cpptr[k][m] = l;
    }
  }
  fclose(ip);

  strcpy(inf1, inf0);
  strcat(inf1, ".rel");
  if ((ip = fopen(inf1, "r")) == 0) {
    fprintf(stderr, "Cannot open %s.\n", inf1);
    return (-1);
  }
  fscanf(ip, "%hd%hd", &k, &rno);
  seeknln(ip);
  op = fopen(outft, "w");
  /* Now we have read everything in, and the computation can start */
  orno = 0;
  for (i = 1; i <= rno; i++) {
    fscanf(ip, "%hd", &rl);
    printf("Scanning rel no %d.\n", i);
    for (j = 1; j <= rl; j++) {
      fscanf(ip, "%hd", &k);
      rel[j] = k;
    }
    for (j = 1; j <= coh_index; j++) {
      pt = j;
      cl = 0;
      for (k = rl; k >= 1; k--) {
        l = rel[k];
        im = cpptr[l][pt];
        if (l % 2 == 1) {
          neg = 1;
          pt1 = im;
          m = (l + 1) / 2;
        }
        else {
          neg = 0;
          pt1 = pt;
          m = (l + 2) / 2;
        }
        ocl = cl;
        for (n = 1; n <= nb; n++) {
          ipt = cpco[n][m][pt1];
          pn = svptr[n][ipt];
          while (pn != -1) {
            cl++;
            cp[cl] = 2 * (pn - 1);
            ipt = pptr[pn][ipt];
            pn = svptr[n][ipt];
          }
        }
        if (neg) {
          pc = cp + ocl + 1;
          qc = cp + cl;
          while (pc <= qc) {
            if (pc == qc)
              (*pc)++;
            else {
              ex = *pc;
              *pc = *qc + 1;
              *qc = ex + 1;
            }
            pc++;
            qc--;
          }
        }
        pt = im;
      }
      if (pt != j) {
        fprintf(stderr,
                "Relation number %d is not satisfied by the permutations.\n",
                i);
        return (-1);
      }
      fprintf(op, "%4d", cl);
      for (k = cl; k >= 1; k--)
        fprintf(op, " %2d", cp[k]);
      fprintf(op, "\n");
    } /* j loop */
    if (i == rno && orno == 0) {
      orno = rno;
      if (fscanf(ip, "%hd", &j) > 0)
        rno += j;
    }

  } /* i loop */
  printf("All relations scanned.\n");
  fclose(op);
  fclose(ip);
  return (0);
}