Quelle  nqp3+.c

#include "defs.h"

/* This version contains the routine subact() called at end of each main
   cycle in calcfm.
*/
extern char  act, inf[], inf1[], inf3[], outf1[], safilech;
extern short facexp, prime, exp, depth, *rpf, *rpb, *eexpnt, *enexpnt, **pcb,
    dim, nng, onng, expnt[], nexpnt[], *pcptr[], **comptr[], *vec[], **mat[],
    cp[], dpth[], pinv[], wt[], *wf, *wc, **extno, **subno, chsdim, chpdim,
    marg;
extern int ptrsp, wsp;
FILE *     ip, *ips, *op;

int comprels(void)
/* 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;
  char 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 (stabdim > 1)
    fprintf(stderr, "Basis of stable homomorphisms:\n");
  if (stabdim > 1)
    for (i = 1; i <= stabdim; i++) {
      for (j = 1; j <= onng; j++)
        fprintf(stderr, "%3d", stabno[i][j]);
      printf("\n");
    }
  for (i = 1; i <= onng; i++)
    stabhom[i] = 0;
  if (stabdim > 1)
    fprintf(stderr,
            "Choose required stable hom as a vector in this basis!\n");
  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.\n");
          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;
    }
  }
  return (0);
}

int subact(void)
{
  short substeps, subexp, *orpf, *commno, **commer, **comm, **subg, **pcp,
      subdp, subdp0, ncommer, ncomm, commst, commlim, i, j, k, l, m, n, *ptr,
      *ptrlim, *ptre, *ptr2;
  char subfile[3], infc[80], str[2];
  bgc();
  strcpy(subfile, "xa");
  orpf = rpf;
restart:
  pcp = pcb + 1;
  rpf = orpf;
  if (act) {
    strcpy(infc, inf3);
    str[0] = safilech;
    str[1] = '\0';
    strcat(infc, str);
  }
  else
    strcpy(infc, inf);
  strcat(infc, subfile);
  if ((ips = fopen(infc, "r")) == 0)
    return (0);
  fscanf(ips, "%hd%hd", &subexp, &substeps);
  if (depth <= substeps) {
    fclose(ips);
    subfile[1]++;
    goto restart;
  }
  subg = pcp - 1;
  pcp += facexp;
  if (pcp - pcptr > ptrsp) {
    printf("Not enough pointer space. Increase PTRSP.\n");
    return (-1);
  }
  subdp0 = 0;
  commst = facexp + 1;
  while (subexp == facexp) {
    subdp0++;
    fscanf(ips, "%hd", &subexp);
    if (subexp == 0)
      subexp = facexp;
    while (commst <= exp && dpth[commst] <= subdp0)
      commst++;
  }

  for (i = 1; i <= subexp; i++) {
    subg[i] = rpf;
    fscanf(ips, "%hd", rpf);
    l = *rpf;
    rpf++;
    for (j = 1; j <= l; j++) {
      fscanf(ips, "%hd", rpf);
      rpf++;
    }
  }
  if (rpb + 1 - rpf < marg) {
    printf("Running out of space in subact.\n");
    if (wsp > marg)
      bgc();
    else
      return (-1);
  }
  commer = pcp - 1;
  commlim = commst;
  i = 0;
  while (commlim <= exp && dpth[commlim] <= depth - substeps + subdp0) {
    *(pcp++) = rpf;
    *(rpf++) = 2;
    *(rpf++) = commlim;
    *(rpf++) = 1;
    if (pcp - pcptr > ptrsp) {
      printf("Not enough pointer space. Increase PTRSP.\n");
      return (-1);
    }
    i++;
    commlim++;
  }
  if (rpb + 1 - rpf < marg) {
    printf("Running out of space.\n");
    return (-1);
  }
  ncommer = i;
  for (subdp = subdp0 + 1; subdp <= substeps; subdp++) {
    if (subdp < substeps) {
      while (dpth[commst] <= subdp)
        commst++;
      while (commlim <= exp && dpth[commlim] <= depth - substeps + subdp)
        commlim++;
      commno = rpf - commst;
      rpf += (commlim - commst);
      if (rpb + 1 - rpf < marg) {
        printf("Running out of space.\n");
        return (-1);
      }
      for (i = commst; i < commlim; i++)
        commno[i] = 0;
      comm = pcp - 1;
      pcp += (commlim - commst);
      ncomm = 0;
      if (pcp - pcptr > ptrsp) {
        printf("Not enough pointer space. Increase PTRSP.\n");
        return (-1);
      }
    }

    for (i = 1; i <= ncommer; i++)
      for (j = 1; j <= subexp; j++) {
        zero(expnt, eexpnt);
        zero(nexpnt, enexpnt);
        wf = rpf;
        wc = wf - 2;
        enter(commer[i], -1);
        enter(subg[j], -1);
        enter(commer[i], 1);
        enter(subg[j], 1);
        collect(wc, wf, 1);
        if (subdp == substeps) {
          if (prnrel() == -1)
            return (-1);
        }
        else {
          ptr = expnt + commst - 1;
          ptrlim = expnt + commlim;
          while (++ptr < ptrlim)
            if ((l = *ptr) != 0) {
              n = pinv[l];
              k = ptr - expnt;
              while (ptr < ptrlim) {
                *ptr *= n;
                *(ptr++) %= prime;
              }
              if ((n = commno[k]) != 0) {
                ptr = comm[n];
                l = *ptr;
                ptre = ptr + l;
                while (++ptr < ptre) {
                  ptr2 = expnt + *ptr;
                  *ptr2 -= *(++ptr);
                  if (*ptr2 < 0)
                    *ptr2 += prime;
                }
                ptr = expnt + k;
              }
              else {
                ptr = expnt + k - 1;
                ncomm++;
                commno[k] = ncomm;
                comm[ncomm] = rpf;
                l = 0;
                while (++ptr < ptrlim)
                  if ((n = *ptr) != 0) {
                    l += 2;
                    *(++rpf) = ptr - expnt;
                    *(++rpf) = n;
                  }
                *comm[ncomm] = l;
                rpf++;
                if (rpb + 1 - rpf < marg) {
                  printf("Running out of space in subact.\n");
                  if (wsp > marg)
                    bgc();
                  else
                    return (-1);
                }
                break;
              }
            }
        }
      }

    fscanf(ips, "%hd", &i);
    if (i != 0) {
      subexp = i;
      rpf = orpf;
      if (subexp == facexp) {
        k = 0;
        for (i = 1; i <= facexp; i++)
          if (wt[i] == 1) {
            k++;
            subg[k] = rpf;
            *(rpf++) = 2;
            *(rpf++) = i;
            *(rpf++) = 1;
          }
        subexp = k;
      }
      else
        for (i = 1; i <= subexp; i++) {
          subg[i] = rpf;
          fscanf(ips, "%hd", rpf);
          l = *rpf;
          rpf++;
          for (j = 1; j <= l; j++) {
            fscanf(ips, "%hd", rpf);
            rpf++;
          }
        }
      commer[1] = rpf;
    }

    if (subdp < substeps) {
      rpf = commer[1];
      ptr = comm[1];
      for (i = 1; i <= ncomm; i++) {
        commer[i] = rpf;
        *(rpf++) = l = *(ptr++);
        for (j = 1; j <= l; j++)
          *(rpf++) = *(ptr++);
      }
      ncommer = ncomm;
      pcp = commer + ncommer + 1;
    }
  }
  printf("File %s.   ", infc);
  printf("Reduced order at depth %d is: %d ^ %d\n", depth, prime, exp + nng);
  fflush(stdout);
  fclose(ips);
  subfile[1]++;
  goto restart;
}