#include "defs.h"

extern char inf1[], inf2[], inf3[], outf1[], outf2[], outf3[], temp1[],
    temp2[], expg, exph, cr, dcr, triv;
extern int psp, trsp, svsp;
extern int mpt, mb, tree[], perm[], gorb[], lorbg[], lorbh[], base[], fpt[],
    bpt[], coh_index[], *cp[], *trad[], *sv[], *tailad[], **cpad[],
    **svgptr[], **svhptr[];
int trptr;
int npt, **lcp, **ucp, *stop, dummy, *pst, *pend, nb, fnt, lnt, ind, coind,
    oind, cind, bno, pt, **expp, npg;
FILE *ip, *op, *opx;

void advance(void);

/* The data structures of this program differ from other permutation group
   programs. Permutations are not numbered, but are located by their base
   adresses. Consequently, arrays sv and cp are arrays of pointers. The
   current permutation in cp is stored in the middle of cp (since it has to be
   expanded in both directions), between adresses lcp and ucp of cp.
   The next six procedures are the equivlents of the corresponding ones
   in permfns.c, but addsv can go in both directions.
   In this half of the program, coset reps of H in G are computed.
*/

int image(int pt)
{
  int ** p;
  p = lcp;
  while (p <= ucp) {
    pt = (*p)[pt];
    p++;
  }
  return (pt);
}

void addsvb(int pt, int ** sv)
{
  int * p;
  while ((p = sv[pt]) != stop) {
    pt = p[pt];
    lcp--;
    *lcp = p - npt;
  }
}

void addsvf(int pt, int ** sv)
{
  int * p;
  while ((p = sv[pt]) != stop) {
    pt = p[pt];
    ucp++;
    *ucp = p;
  }
}

void invert(int * a, int * b)
{
  int i;
  for (i = 1; i <= npt; i++)
    b[a[i]] = i;
}

void rdperm(int * a, int * b)
{
  int i, *r;
  for (i = 1; i <= npt; i++)
    fscanf(ip, "%d", a + i);
  r = pst + 1;
  fscanf(ip, "%d", r);
  invert(a, b);
}

void rdsv(int ** sv)
{
  int i, j;
  for (i = 1; i <= npt; i++) {
    fscanf(ip, "%d", &j);
    sv[i] = (j == 0) ? 0 : (j == -1) ? stop : cp[j];
  }
}

int cnprg1(void)
{
  int  i, j, k, ad, sad, nph, *svpt, *p1, *p2, **csvg, **csvh, nexp;
  char pthere, w, ok, allok;
  long fac;
  stop = &dummy;
  if ((ip = fopen(inf1, "r")) == 0) {
    fprintf(stderr, "Cannot open %s.\n", inf1);
    return (-1);
  }
  fscanf(ip, "%d%d%d%d", &npt, &npg, &nb, &k);
  if (npt > mpt) {
    fprintf(stderr, "npt too big. Increase mpt.\n");
    return (-1);
  }
  if (nb >= mb) {
    fprintf(stderr, "nb too big. Increase MB.\n");
    return (-1);
  }
  if (2 * nb * npt > svsp) {
    fprintf(stderr, "svsp too small. Increase SVSP.\n");
    return (-1);
  }
  if (k <= 0) {
    fprintf(stderr, "Wrong input format for G.\n");
    return (-1);
  }
  for (i = 1; i <= nb; i++)
    fscanf(ip, "%d", base + i);
  for (i = 1; i <= nb; i++)
    fscanf(ip, "%d", lorbg + i);
  pst = perm - 1;
  pend = perm + psp - 1;
#ifdef __clang__
  /* workaround for a bug in Apple LLVM version 7.0.2 (clang-700.1.81),
    and possibly other versions. See also
    https://github.com/gap-packages/cohomolo/issues/5
  */
  if (psp < 0) /* never should happen */
    fprintf(stderr, "%d", psp);
#endif
  if (pst + 2 * npg * npt > pend) {
    fprintf(stderr, "Out of space.Increase PSP.\n");
    return (-1);
  }
  for (i = 1; i <= npg; i++) {
    p1 = pst;
    p2 = pst + npt;
    pst += (2 * npt);
    rdperm(p1, p2);
    cp[2 * i - 1] = p2;
  }
  for (i = 1; i <= nb; i++) {
    svgptr[i] = sv + npt * (i - 1) - 1;
    rdsv(svgptr[i]);
  }
  fclose(ip);
  /* G is read. Now read the subgroup H */
  if (triv) {
    nph = 0;
    for (i = 1; i <= nb; i++) {
      lorbh[i] = 1;
      svhptr[i] = sv + npt * (nb + i - 1) - 1;
      for (j = 1; j <= npt; j++)
        svhptr[i][j] = 0;
      svhptr[i][base[i]] = stop;
    }
  }
  else {
    if ((ip = fopen(inf2, "r")) == 0) {
      fprintf(stderr, "Cannot open %s.\n", inf2);
      return (-1);
    }
    fscanf(ip, "%d%d%d%d", &i, &nph, &j, &k);
    w = 0;
    if (k <= 0) {
      fprintf(stderr, "Wrong input format for H.\n");
      return (-1);
    }
    if (i != npt || j != nb)
      w = 1;
    else
      for (i = 1; i <= nb; i++) {
        fscanf(ip, "%d", lorbh + i);
        if (lorbh[i] != base[i])
          w = 1;
      }
    if (w) {
      fprintf(stderr, "npt or base for G and H do not agree.\n");
      return (-1);
    }
    for (i = 1; i <= nb; i++)
      fscanf(ip, "%d", lorbh + i);
    if (pst + 2 * npt * nph + 2 > pend) {
      fprintf(stderr, "Out of space.Increase PSP.\n");
      return (-1);
    }
    for (i = 1; i <= nph; i++) {
      p2 = pend - npt;
      p1 = p2 - npt;
      pend -= (2 * npt);
      rdperm(p1, p2);
      cp[2 * i - 1] = p2;
    }
    for (i = 1; i <= nb; i++) {
      svhptr[i] = sv + npt * (nb + i - 1) - 1;
      rdsv(svhptr[i]);
    }
    fclose(ip);
  }
  for (i = 1; i <= nb; i++)
    fpt[i] = 0;
  p1 = pst;
  p2 = p1 + npt;
  pst += (2 * npt);
  lcp = cp;
  expp = cp + 2 * npt - 1;
  cind = 1;
  fac = 1;
  coh_index[nb + 1] = 1;
  lnt = 0;

  /* Now we compute the indices of H in G in the stabilizer chain into the
     array coh_index. The nontrivial indices are linked by fpt and bpt.
  */
  for (i = nb; i >= 1; i--) {
    j = lorbg[i];
    fac = fac * j;
    j = lorbh[i];
    fac = fac / j;
    ind = fac;
    coh_index[i] = ind;
    if (ind > coh_index[i + 1]) {
      if (lnt == 0)
        lnt = i;
      else {
        fpt[i] = fnt;
        bpt[fnt] = i;
      }
      fnt = i;
    }
  }
  bpt[fnt] = -1;
  fpt[lnt] = -1;
  printf("Indices:  ");
  for (i = 1; i <= nb; i++)
    printf("%5d", coh_index[i]);
  printf("\n");
  if (ind == 1) {
    fprintf(stderr, "G=H.\n");
    return (-1);
  }
  nexp = 0;
  ind = 1;
  oind = 1;
  trptr = 0;

  /* Now we start the backtrack search for the left coset reps of H in G. The
     necessary information about these is stored in the array tree. This is
     done so as to make it as easy as possible to compute which coset an
     arbitrary element of G lies in. If expg is true, then those perms in G
     which arise are expanded and stored. If dcr is true, then a small subset
     of the coset reps will need to be output as perms, later. For this
     purpose, the base images of all coset reps are output to a temporary
     file. ind is the current coh_index coh_index[bno], and cind is the no of
     coset reps found so far.
  */
  for (i = fnt; i != -1; i = fpt[i]) {
    trad[i] = tree + trptr;
    tree[trptr] = base[i];
    trptr += 2;
  }
  if (dcr)
    opx = fopen(temp2, "w");
  if (cr) {
    op = fopen(outf3, "w");
    fprintf(op, "%4d %d%4d%4d\n", npt, coh_index[1] - 1, 0, 0);
  }
  sad = trptr;
  tree[trptr] = 1;
  trptr++;
  for (bno = lnt; bno != -1; bno = bpt[bno]) {
    printf("bno=%d.\n", bno);
    *gorb = 1;
    gorb[1] = base[bno];
    ind = coh_index[bno];
    csvg = svgptr[bno];
    csvh = svhptr[bno];
    allok = (lorbh[bno] == 1);
    /* if allok, then all coset reps in this link of the stabilizer chain will
       be coset reps of H in G.
    */
    if (expg) {
      for (i = 1; i <= npt; i++)
        expp[i] = 0;
      expp[base[bno]] = stop;
    }
    for (pt = 1; pt <= npt; pt++) {
      svpt = csvg[pt];
      if (svpt != 0 && svpt != stop) {
        ucp = lcp - 1;
        addsvf(pt, csvg);
        pthere = 0;
        if (expg) {
          for (i = 1; i <= npt; i++)
            p1[i] = image(i);
          invert(p1, p2);
          ucp = lcp;
          *ucp = p1;
        }
        j = sad;
        for (i = fpt[bno]; i != -1; i = fpt[i]) {
          tailad[i] = tree + j;
          cpad[i] = ucp + 1;
          j += 2;
        }
        /* In the search, for each basic coset rep g in the stabilizer chain,
           we have to try out gh as a possible coset rep of H in G, for each h
           in the list of coset reps of H[bno-1] in G[bno-1]. coind and oind
           record the search through this list.
        */
        coind = 1;
        while (coind <= oind) {
          if (coind > 1)
            advance();
          ok = 1;
          if (allok == 0)
            for (i = 1; i <= *gorb && ok; i++)
              if (csvh[image(gorb[i])] != 0)
                ok = 0;
          /* That was the membership test. ok true means we have found a new
           * rep */
          if (ok) {
            if (pthere) {
              ad = fpt[bno];
              while (*tailad[ad] == *trad[ad])
                ad = fpt[ad];
              tree[trptr] = *tailad[ad];
            }
            else {
              pthere = 1;
              ad = bno;
              tree[trptr] = pt;
              if (expg) {
                expp[pt] = p1;
                p1 = pst;
                p2 = p1 + npt;
                pst += (2 * npt);
                nexp++;
                if (pst > pend) {
                  fprintf(stderr, "Out of perm space. Increase PSP.\n");
                  return (-1);
                }
              }
            }
            *(trad[ad] + 1) = trptr;
            trad[ad] = tree + trptr;
            trptr += 2;
            for (i = fpt[ad]; i != -1; i = fpt[i]) {
              *(trad[i] + 1) = -1;
              tree[trptr] = *tailad[i];
              trad[i] = tree + trptr;
              trptr += 2;
            }
            cind++;
            tree[trptr] = cind;
            trptr++;
            if (trptr > trsp) {
              fprintf(stderr, "Out of tree space. Increase TRSP.\n");
              return (-1);
            }
            if (dcr)
              for (i = fnt; i != -1; i = fpt[i])
                fprintf(opx, "%5d", *trad[i]);
            /* if cr, then we output each coset rep */
            if (cr) {
              if (npt >= 1000) {
                for (i = 1; i <= npt; i++)
                  fprintf(op, "%5d", image(i));
                fprintf(op, "\n");
              }
              else {
                for (i = 1; i <= npt; i++)
                  fprintf(op, "%4d", image(i));
                fprintf(op, "\n");
              }
            }
            if (cind == ind)
              goto nextbno;
          } /* if (ok) */
          coind++;
        } /* while (coind<=oind) */
        (*gorb)++;
        gorb[*gorb] = pt;
      } /* if (svpt!=0 ... */
    }   /* for (pt=1;pt<=npt;... */
    fprintf(stderr, "Premature end of search.\n");
    return (-1);
  nextbno:
    if (expg)
      for (i = 1; i <= npt; i++)
        svgptr[bno][i] = expp[i];
    sad -= 2;
    oind = ind;
  } /* main bno loop */
  if (dcr)
    fclose(opx);
  if (cr)
    fclose(op);
  for (i = fnt; i != -1; i = fpt[i])
    *(trad[i] + 1) = -1;
  printf("Tree space used = %d.\n", trptr);
  if (expg)
    printf("Number of perms expanded = %d.\n", nexp);
  return (0);
}

void advance(void)
/* This advances the element h in the search through elements gh */
{
  int ad, k, *p, *q;
  ad = lnt;
  while ((k = *(tailad[ad] + 1)) == -1)
    ad = bpt[ad];
  q = tree + k;
  ucp = cpad[ad] - 1;
  while (ad != -1) {
    tailad[ad] = q;
    cpad[ad] = ucp + 1;
    if (expg) {
      p = svgptr[ad][*q];
      if (p != stop) {
        ucp++;
        *ucp = p;
      }
    }
    else
      addsvf(*q, svgptr[ad]);
    q += 2;
    ad = fpt[ad];
  }
}