#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <ctype.h>
#include <string.h>
#include "rubiksmagic.h"

static int polyo[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX];
static int nomirror = 0;
static int antimirror = 0;
static int findpolyogroup = 0;
static int ispolyochiral = 0;

int
main (int argc, char *argv[])
{
  sequence seq, cseq;
  int i, numseq;
  int polyosym, maxpos[3];
  char *argv0;
  int revert, rotate, mirror, inout, spin;

  verbose = 0;
  argv0 = argv[0];

  while (argc >= 2 && *argv[1] == '-')
  {
    if (strcmp (argv[1], "-n") == 0)
    {
      assert (argc >= 3);
      seqlen = atoi (argv[2]);
      assert (seqlen <= SEQDIM);
      argc -= 2;
      argv += 2;
      continue;
    }
    if (strcmp (argv[1], "-p") == 0)
    {
      findpolyogroup = verbose;
      assert (argc >= 3);
      readpolyominoid (argv[2], polyo);
      traslatepolyominoid (polyo, maxpos);
      polyosym = optimalpolyominoid (polyo, maxpos);
      if (findpolyogroup)
      {
        if (ispolyochiral) printf ("Chiral\n");
        printf ("Number of symmetries: %d\n", polyosym);
      }
      printpolyominoid (polyo);
      if (polyosym > 1) printf ("s");
      printf ("\n");
      argc -= 2;
      argv += 2;
      return (0);
    }
    if (strcmp (argv[1], "-c") == 0)
    {
      nocanon = 1;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "-v") == 0)
    {
      verbose++;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "-q") == 0)
    {
      quiet++;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "-d") == 0)
    {
      debug++;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "-M") == 0)
    {
      filtermetric++;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "-T") == 0)
    {
      filterlinking++;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "--nomirror") == 0)
    {
      nomirror++;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "--antimirror") == 0)
    {
      antimirror++;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "--verify") == 0)
    {
      verify++;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "--cycle") == 0)
    {
      cycle = 1;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "--warpcode") == 0 || strcmp (argv[1], "-w") == 0)
    {
      outwarpcode++;
      outbox = 1 - outbox;
      outpolyo = 1 - outpolyo;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "--off") == 0)
    {
      off = 1;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "--nocheckorientable") == 0)
    {
      checkorientable = 0;
      argc--;
      argv++;
      continue;
    }
    if (strcmp (argv[1], "-h") == 0 || strcmp (argv[1], "--help") == 0)
    {
      printf ("%s [-n <n>][-c][-d] [sequence]\n", argv0);
      printf ("  -n <n> number of tiles (default 8)\n");
      printf ("  -c     do not canonify\n");
      printf ("  -M     filter by metric invariant\n");
      printf ("  -T     filter by linking number (not implemented)\n");
      printf ("  -v     increase verbosity\n");
      printf ("  -q     be quiet\n");
      printf ("  -d     increase debug level\n");
      printf ("  --nomirror  avoid mirror images in the canonification process\n");
      printf ("  --antimirror same as 'nomirror' but for a mirror image\n");
      printf ("  --verify    verify invariance of invariants\n");
      printf ("  --cycle     print info on all equivalent strings\n");
      printf ("  --off       use OFF syntax for polyominoid\n");
      printf ("  --nocheckorientable     do not check if sequence is orientable\n");
      printf ("  -p xyz-xyz-xyz-xyz-xyz-xyz-xyz-xyz    compute canonical polyominoid\n");
      printf ("  --warpcode|-w   compute Tom Verhoeff warp code\n");
      printf ("              the second sequence is the canonic representative\n");
      printf ("              (to avoid suppressing other info specify option twice)\n");
      argc--;
      argv++;
      exit (0);
    }
    printf ("Invalid option: %s\n", argv[1]);
    exit (100);
  }

  init_data ();

  if (argc >= 2)
  {
    readseq (&seq, argv[1]);
    verbose++;
    if (off)
    {
      printoff (&seq, stdout);
      return (0);
    }
    if (antimirror)
    {
      action (&seq, &cseq, 0, 0, 1, 0, 0);
      nomirror++;
      seq.type = cseq.type;
      for (i = 0; i < seqlen; i++) seq.seq[i] = cseq.seq[i];
    }
    printseqinfo (&seq);
    if (cycle)
    {
      for (revert = 0; revert <=1; revert++)
      {
        for (rotate = 0; rotate < seqlen; rotate++)
        {
          for (mirror = 0; mirror <= 1; mirror++)
          {
            for (inout = 0; inout <= 1; inout++)
            {
              for (spin = 0; spin < 4; spin++)
              {
                action (&seq, &cseq, revert, rotate, mirror, inout, spin);
                printseqinfo (&cseq);
              }
            }
          }
        }
      }
    }
    if (nocanon == 0)
    {
      canonify (&seq);
      printseqinfo (&seq);
    }
  } else {
    firstsequence (&seq);
    numseq = 0;
    do {
      if (cannotbecanon (&seq)) continue;
      if (!checkvalid (&seq, seqlen)) continue;
      if (checkorientable)
      {
        seqcpy (&cseq, &seq);
        canonify (&cseq);
        if (seqcmp (&seq, &cseq) != 0) continue;
      }
      if (filtermetric && computedelta (&seq) != 0) continue;
      if (filterlinking && computelinking (&seq) != 0) continue;
      minribboncrossings = -1;
      if (verify) verifyinvariance (&seq);
      printseqinfo (&seq);
      numseq++;
    } while (nextsequence (&seq));
    fprintf (stderr, "Found %d sequences\n", numseq);
  }

  return (0);
}

/*
 */

void
seqcpy (sequence *seq1, sequence *seq2)
{
  int i;

  seq1->type = seq2->type;
  for (i = 0; i < seqlen; i++) seq1->seq[i] = seq2->seq[i];
}

void
init_data ()
{
#ifdef OLDDIRS
  chardir[ALTO>>2] = 'N';
  chardir[DX>>2] = 'E';
  chardir[BASSO>>2] = 'S';
  chardir[SX>>2] = 'W';
#else
  chardir[ALTO>>2] = 'U';
  chardir[DX>>2] = 'R';
  chardir[BASSO>>2] = 'D';
  chardir[SX>>2] = 'L';
#endif

  charfold[MF] = 'm';
  charfold[VF] = 'v';

  rotright[ALTO] = DX;
  rotright[DX] = BASSO;
  rotright[BASSO] = SX;
  rotright[SX] = ALTO;
  rotleft[ALTO] = SX;
  rotleft[DX] = ALTO;
  rotleft[BASSO] = DX;
  rotleft[SX] = BASSO;
}

/*
 * ciclo su tutte le sequenze chiuse
 */

void
firstsequence (sequence *seq)
{
  int i;
  seq->type = ISSLASH;
  for (i = 0; i < seqlen; i++) seq->seq[i] = DX;
}

int
nextsequence (sequence *seq)
{
  return (nextsequencerec(seq, 0));
}

int
nextsequencerec (sequence *seq, int start)
{
  int startp1;

  startp1 = start + 1;
  if (startp1 < seqlen)
  {
    if (nextsequencerec (seq, startp1)) return (1);
  }

  while (seq->seq[start] < DIROVER)
  {
    seq->seq[start]++;
    if ((seq->seq[start] & FOLDMASK) == FOLDMASK) seq->seq[start]++;
    if (seq->seq[start] >= DIROVER) break;
    if (start <= 1 && cannotbecanon (seq)) continue;
    if (startp1 >= seqlen) break;
    if (checkvalid (seq, startp1)) break;
  }
  if (seq->seq[start] < DIROVER) return (1);
  seq->seq[start] = DX;
  if (start > 0) return (0);
  if (seq->type == ISBSLASH) return (0);
  firstsequence (seq);
  seq->type = ISBSLASH;
  return (1);
}

/*
 * basic and quick check: 
 * if true, the sequence cannot be canonical
 * it only checks the first sequence element
 */

int
cannotbecanon (sequence *seq)
{
  if (seq->type != ISSLASH) return (1);
  if ((seq->seq[0] & DIRMASK) != DX) return (1);
  if ((seq->seq[0] & FOLDMASK) == VF) return (1);
  return (0);
}

/*
 * check validity on the substring from 0 to 'uptohere'
 * validity: no superposed tiles and no hinge selfcrossings
 */

int
checkvalid (sequence *seq, int uptohere)
{
  int i, j, dir, fold;
  int iposizione[3] = {1,1,0};
  int ifreccia[3] = {0,1,0};
  int inormale[3] = {0,0,1};
  int freccia[3], nfreccia[3], normale[3], nnormale[3];
  int hinges[seqlen][3], posizioni[seqlen+1][3];
  int isstraight[seqlen];
  int diff[3];
  int distsq, maxdist;

  assert (uptohere <= seqlen);
  if (verbose) {printseq (seq); printf ("[%d] - ", uptohere);}
  if (verbose) printf ("F:(%d %d %d)\n", iposizione[0], iposizione[1], iposizione[2]);
  veccpy (posizioni[0], iposizione);
  veccpy (freccia, ifreccia);
  veccpy (normale, inormale);

  for (i = 0; i < uptohere; i++)
  {
    dir = seq->seq[i] & DIRMASK;
    fold = seq->seq[i] & FOLDMASK;
    isstraight[i] = 0;
    if (fold == NF) isstraight[i] = 1;
    newposition (dir, fold, posizioni[i], freccia, normale,
                 posizioni[i+1], nfreccia, nnormale, hinges[i]);
    veccpy (freccia, nfreccia);
    veccpy (normale, nnormale);

    if (verbose) printf ("S:(%d %d %d)\n", hinges[i][0], hinges[i][1], hinges[i][2]);
    if (verbose) printf ("F:(%d %d %d)\n", posizioni[i+1][0], posizioni[i+1][1], posizioni[i+1][2]);
    checkerrno = ERR_SUPERPOSEDTILES;
    if (i + 1 < seqlen)
    {
      for (j = 0; j <= i; j++)
      {
        if (veceq(posizioni[j], posizioni[i+1])) return (0);  // due tasselli sovrapposti
      }
    }
    checkerrno = ERR_SELFINTERSECTION;
    if (isstraight[i])
    {
      for (j = 0; j < i-1; j++)
      {
        if (j == 0 && i + 1 == seqlen) continue; // non devono essere contigui...
        if (!isstraight[j]) continue;
        if (veceq(hinges[j], hinges[i])) return (0);  // attraversamento
      }
    }
  }

  if (uptohere >= seqlen) // check if closed
  {
    checkerrno = ERR_NOTCLOSED;
    if (! veceq(posizioni[0], posizioni[seqlen])) return (0);
    checkerrno = ERR_NOTORIENTABLE;
    if (checkorientable && ! veceq(freccia, ifreccia)) return (0);
    if (checkorientable && ! veceq(normale, inormale)) return (0);
  } else {
    checkerrno = ERR_NOTCLOSED;
    vecdiff (diff, posizioni[uptohere], posizioni[0]);
    distsq = diff[0]*diff[0] + diff[1]*diff[1] + diff[2]*diff[2];
    maxdist = 2*(seqlen - uptohere);
    if (distsq > maxdist*maxdist) return (0);   //impossibile chiudere!
  }
  checkerrno = 0;
  return (1);  // for now
}

/*
 * calcola la posizione e l'orientazione della nuova faccia
 * in base alla direzione e al fold
 */

void
newposition (int dir, int fold, int *posizione, int *freccia, int *normale,
             int *nposizione, int *nfreccia, int *nnormale, int *hinge)
{
  int est[3], nest[3];

  vecprod (est, freccia, normale);
  switch (dir)
  {
    case DX:
      vecsum (hinge, posizione, est);
      break;

    case SX:
      vecdiff (hinge, posizione, est);
      break;

    case ALTO:
      vecsum (hinge, posizione, freccia);
      break;

    case BASSO:
      vecdiff (hinge, posizione, freccia);
      break;
  }

  switch (fold)
  {
    case NF:
    veccpy (nest, est);
    veccpy (nfreccia, freccia);
    veccpy (nnormale, normale);
    break;

    case MF:
    case VF:
    switch (dir)
    {
      case DX: //freccia rimane inalterato
        veccpy (nfreccia, freccia);
        veccpy (nnormale, est);
        if (fold == VF) vecneg (nnormale);
        vecprod (nest, nfreccia, nnormale);
        break;

      case SX:
        veccpy (nfreccia, freccia);
        veccpy (nnormale, est);
        if (fold == MF) vecneg (nnormale);
        vecprod (nest, nfreccia, nnormale);
        break;

      case ALTO: //est rimane inalterato
        veccpy (nest, est);
        veccpy (nnormale, freccia);
        if (fold == VF) vecneg (nnormale);
        vecprod (nfreccia, nnormale, nest);
        break;
        
      case BASSO:
        veccpy (nest, est);
        veccpy (nnormale, freccia);
        if (fold == MF) vecneg (nnormale);
        vecprod (nfreccia, nnormale, nest);
        break;
    }
  }

  switch (dir)
  {
    case DX:
      vecsum (nposizione, hinge, nest);
      break;

    case SX:
      vecdiff (nposizione, hinge, nest);
      break;

    case ALTO:
      vecsum (nposizione, hinge, nfreccia);
      break;

    case BASSO:
      vecdiff (nposizione, hinge, nfreccia);
      break;
  }

  if ((dir == ALTO) || (dir == BASSO))
  {
    vecneg (nfreccia);
    vecneg (nest);
  }
}

/*
 * vector basic operations
 */

void
veccpy (int *res, int *v)
{
  int i;
  for (i = 0; i < 3; i++) res[i] = v[i];
}

void
vecneg (int *res)
{
  int i;
  for (i = 0; i < 3; i++) res[i] = -res[i];
}

void
vecsum (int *res, int *v1, int*v2)
{
  int i;
  for (i = 0; i < 3; i++) res[i] = v1[i] + v2[i];
}
 
void
vecdiff (int *res, int *v1, int*v2)
{
  int i;
  for (i = 0; i < 3; i++) res[i] = v1[i] - v2[i];
}
 
void
vecprod (int *res, int *v1, int*v2)
{
  res[0] = v1[1]*v2[2] - v1[2]*v2[1];
  res[1] = v1[2]*v2[0] - v1[0]*v2[2];
  res[2] = v1[0]*v2[1] - v1[1]*v2[0];
}

int
veceq (int *v1, int *v2)
{
  int i;
  for (i = 0; i < 3; i++) if (v1[i] != v2[i]) return (0);
  return (1);
}

void
vecprint (int *v)
{
  printf ("(%d,%d,%d)", v[0], v[1], v[2]);
}

void
canonify (sequence *seq)
{
  sequence wseq;
  sequence optseq;
  int i;
  int revert, rotate, mirror, inout, spin;

  optseq.type = seq->type;
  for (i = 0; i < seqlen; i++) optseq.seq[i] = seq->seq[i];

  for (revert = 0; revert <=1; revert++)
  {
    for (rotate = 0; rotate < seqlen; rotate++)
    {
      for (mirror = 0; mirror <= 1; mirror++)
      {
        for (inout = 0; inout <= 1; inout++)
        {
          for (spin = 0; spin < 4; spin++)
          {
            action (seq, &wseq, revert, rotate, mirror, inout, spin);
            if (seqcmp (&wseq, &optseq) < 0)
            {
              optseq.type = wseq.type;
              for (i = 0; i < seqlen; i++) optseq.seq[i] = wseq.seq[i];
            }
          }
        }
      }
    }
  }
  seq->type = optseq.type;
  for (i = 0; i < seqlen; i++) seq->seq[i] = optseq.seq[i];
  return;
}

/*
 * warning: "revert" reverts the circular order of the tiles *and*
 * performs a mirror image
 * "inout" exchange faceup with facedown, also entails a mirror image
 */

void
action (sequence *seq, sequence *dseq, int revert, int rotate, int mirror,
        int inout, int spin)
{
  sequence tseq;
  int i, destad;
  int dir, fold;
  int typechange;

  assert (rotate >= 0 && rotate < seqlen);
  assert (spin >= 0 && spin < 4);

  if (nomirror)
  {
    mirror = revert;
    if (inout) mirror = 1 - mirror;
  }
  tseq.type = seq->type;
  for (i = 0; i < seqlen; i++) tseq.seq[i] = seq->seq[i];

  if (revert) for (i = 0; i < seqlen; i++) tseq.seq[i] = seq->seq[seqlen - i - 1];
  typechange = rotate % 2;
  if (revert) typechange = 1 - typechange;
  if (mirror) typechange = 1 - typechange;
  if (inout) typechange = 1 - typechange;
  if (spin % 2) typechange = 1 - typechange;
  assert (typechange == 0 || typechange == 1);
  dseq->type = tseq.type;
  if (typechange) dseq->type = ISSLASH + ISBSLASH - tseq.type;
  for (i = 0; i < seqlen; i++)
  {
    destad = rotate + i;
    destad = destad % seqlen;
    dseq->seq[destad] = tseq.seq[i];
  }
  for (i = 0; i < seqlen; i++)
  {
    dir = dseq->seq[i] & DIRMASK;
    fold = dseq->seq[i] & FOLDMASK;
    if (mirror) // exchange DX with SX
    {
      if (dir == DX || dir == SX) dir = DX + SX - dir;
    }
    if (inout) // invert faceup
    {
      if (fold == MF || fold == VF) fold = MF + VF - fold;
    }
    if (spin >= 2)
    {
      if (dir == DX || dir == SX) dir = DX + SX - dir;
      if (dir == ALTO || dir == BASSO) dir = ALTO + BASSO - dir;
    }
    if ((spin % 2))
    {
      if ((i % 2) == 0)
      {
        dir = rotright[dir];
      } else {
        dir = rotleft[dir];
      }
    }
    dseq->seq[i] = dir + fold;
  }
}

int
seqcmp (sequence *seq1, sequence *seq2)
{
  int i;

  if (seq1->type != seq2->type) return ((seq1->type < seq2->type)?(-1):1);
  for (i = 0; i < seqlen; i++)
  {
    if (seq1->seq[i] != seq2->seq[i]) return ((seq1->seq[i] < seq2->seq[i])?(-1):1);
  }
  return (0);
}

void
printseqinfo (sequence *seq)
{
  int box[3];
  int border;
  int polyosym;
  //int polyo[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX];
  int delta, lt, lf, lc;
  int warpcode[SEQDIM];

  printseq (seq);
  if (outbox)
  {
    border = computebox (seq, box);
    printf (" box=%dx%dx%d+%d", box[0], box[1], box[2], border);
  }
  if (outpolyo && (polyosym = polyominoid (seq, polyo)))
  {
    printf (" polyominoid=");
    printpolyominoid (polyo);
    if (polyosym > 1) printf ("s");
  }
  printf (" f=%d", countflaps (seq));
  delta = computedelta (seq);
  printf (" delta=%d", delta);
  ribbontouchings = 0;
  if (checkorientable)
  {
    lt = computelt (seq);
    lf = computelf (seq);
    lc = computelc (seq);
    printf (" linking=%d", lt + lf + lc);
    if (cycle || verify) printf (" lt+lf+lc=%d%+d%+d", lt, lf, lc);
  }
  printf (" symcount=%d", countsymmetries (seq));
  printf (" typeinv=%s", istypeinvariant (seq)?"yes":"no");
  if (minribboncrossings >= 0) printf (" mincrossings=%d", minribboncrossings);
  if (outwarpcode)
  {
    printf (" warpcode:");
    computewarpcode (seq, warpcode);
    printwarpcode (warpcode);
    canonifywarpcode (warpcode);
    printf (":");
    printwarpcode (warpcode);
  }
  printf ("\n");
}

void
printseq (sequence *seq)
{
  int i, fold;

  if (seq->type == ISBSLASH) printf ("\\");
  for (i = 0; i < seqlen; i++)
  {
    printf ("%c", chardir[(seq->seq[i] >> 2)]);
    fold = seq->seq[i] & FOLDMASK;
    if (fold) printf ("%c", charfold[fold]);
  }
}

void
readseq (sequence *seq, char *stringa)
{
  int i;
  char ch, *chpt;

  assert (NF == 0);
  chpt = stringa;
  while (isspace(*chpt)) chpt++;

  seq->type = ISSLASH;
  if (*chpt == '\\' || *chpt == '/')
  {
    if (*chpt == '\\') seq->type = ISBSLASH;
    chpt++;
  }

  for (i = 0; *chpt && i < SEQDIM; i++)
  {
    while (isspace(ch = tolower(*chpt++)));
    assert (isalpha(ch));
    switch (ch)
    {
#ifdef OLDDIRS
      case 'e': seq->seq[i] = DX;
        break;

      case 'n': seq->seq[i] = ALTO;
        break;

      case 'w': seq->seq[i] = SX;
        break;

      case 's': seq->seq[i] = BASSO;
        break;
#else
      case 'r': seq->seq[i] = DX;
        break;

      case 'u': seq->seq[i] = ALTO;
        break;

      case 'l': seq->seq[i] = SX;
        break;

      case 'd': seq->seq[i] = BASSO;
        break;
#endif

      default:
        fprintf (stderr, "invalid character %c\n", ch);
        exit (ERR_INVCHAR);
        break;
    }
    ch = tolower (*chpt);
    if (ch == 'm' || ch == 'v')
    {
      chpt++;
      if (ch == 'm')
      {
        seq->seq[i] += MF;
      } else {
        seq->seq[i] += VF;
      }
    }
  }
  if (i != seqlen)
  {
    fprintf (stderr, "Error: wrong number of tiles in magic code\n");
    exit (ERR_WRONGLEN);
  }
  if (!checkvalid(seq, seqlen))
  {
    fprintf (stderr, "Error: sequence is not valid (errcode %d).\n", checkerrno);
    exit (checkerrno);
  }
}

void
computewarpcode (sequence *seq, int code[SEQDIM])
{
  int i, in, out, type, delta;

  type = seq->type;
  for (i = 0; i < seqlen; i++)
  {
    in = seq->seq[i];
    out = seq->seq[(i+1)%seqlen];
    type = ISSLASH + ISBSLASH - type;
    delta = deltatile (in, out, type);
    if (type == ISSLASH) delta = -delta;
    code[i] = delta;
  }
}

void
canonifywarpcode (int code[SEQDIM])
{
  int i, j, sign, dir, rot, dest;
  int optcode[SEQDIM], tcode[SEQDIM];

  for (i = 0; i < seqlen; i++) optcode[i] = code[i];

  for (sign = 1; sign >= -1; sign -= 2)
  {
    for (dir = 1; dir >= -1; dir -= 2)
    {
      for (rot = 0; rot < seqlen; rot++)
      {
        for (i = 0; i < seqlen; i++)
        {
          dest = (rot + dir*i + seqlen) % seqlen;
          tcode[dest] = sign*code[i];
        }
        for (i = 0; i < seqlen; i++)  // compare with optimal
        {
          if (tcode[i] < optcode[i]) break;
          if (tcode[i] > optcode[i])
          {
            for (j = 0; j < seqlen; j++) optcode[j] = tcode[j];
            break;
          }
        }
      }
    }
  }
  for (i = 0; i < seqlen; i++) code[i] = optcode[i];
}

void
printwarpcode (int code[SEQDIM])
{
  int i;

  for (i = 0; i < seqlen; i++)
  {
    switch (code[i])
    {
      case 0:
        printf ("0");
        break;

      case 1:
        printf ("+");
        break;

      case (-1):
        printf ("-");
        break;

      case 2:
        printf ("#");
        break;

      case (-2):
        printf ("=");
        break;
    }
  }
}

int
countsymmetries (sequence *seq)
{
  sequence wseq;
  int count = 0;
  int revert, rotate, mirror, inout, spin;

  for (revert = 0; revert <=1; revert++)
  {
    for (rotate = 0; rotate < seqlen; rotate++)
    {
      for (mirror = 0; mirror <= 1; mirror++)
      {
        for (inout = 0; inout <= 1; inout++)
        {
          for (spin = 0; spin < 4; spin++)
          {
            action (seq, &wseq, revert, rotate, mirror, inout, spin);
            if (seqcmp (&wseq, seq) == 0) count++;
          }
        }
      }
    }
  }
  return (count);
}

int
countflaps (sequence *seq)
{
  int i, dir, dir2;
  int flapsnum = 0;

  for (i = 0; i < seqlen; i++)
  {
    dir = seq->seq[i] & DIRMASK;
    dir2 = seq->seq[(i+1)%seqlen] & DIRMASK;
    if (dir == DX || dir == SX)
    {
      if (dir + dir2 == DX + SX) flapsnum++;
    }
    if (dir == ALTO || dir == BASSO)
    {
      if (dir == dir2) flapsnum++;
    }
  }
  return (flapsnum);
}

int
computebox (sequence *seq, int box[3])
{
  int iposizione[3] = {1,1,0};
  int ifreccia[3] = {0,1,0};
  int inormale[3] = {0,0,1};
  int posizione[3];
  int nposizione[3];
  int hinge[3], freccia[3], nfreccia[3], normale[3], nnormale[3];
  int minpos[3];
  int maxpos[3];
  int i, j, dir, fold, countboundary, temp;

  veccpy (posizione, iposizione);
  veccpy (freccia, ifreccia);
  veccpy (normale, inormale);
  veccpy (minpos, iposizione);
  veccpy (maxpos, iposizione);
  for (i = 0; i < seqlen; i++)
  {
    dir = seq->seq[i] & DIRMASK;
    fold = seq->seq[i] & FOLDMASK;
    newposition (dir, fold, posizione, freccia, normale,
                 nposizione, nfreccia, nnormale, hinge);
    veccpy (posizione, nposizione);
    veccpy (freccia, nfreccia);
    veccpy (normale, nnormale);
    for (j = 0; j < 3; j++)
    {
      if (posizione[j] < minpos[j]) minpos[j] = posizione[j];
      if (posizione[j] > maxpos[j]) maxpos[j] = posizione[j];
    }
  }

  for (j = 0; j < 3; j++)
  {
    if ((minpos[j]+2*SEQDIM) % 2) minpos[j]--;  // round to even
    if ((maxpos[j]+2*SEQDIM) % 2) maxpos[j]++;
  }
  countboundary = 0;
  veccpy (posizione, iposizione);
  veccpy (freccia, ifreccia);
  veccpy (normale, inormale);

  for (i = 0; i < seqlen; i++)
  {
    dir = seq->seq[i] & DIRMASK;
    fold = seq->seq[i] & FOLDMASK;
    newposition (dir, fold, posizione, freccia, normale,
                 nposizione, nfreccia, nnormale, hinge);
    veccpy (posizione, nposizione);
    veccpy (freccia, nfreccia);
    veccpy (normale, nnormale);
    for (j = 0; j < 3; j++)
    {
      if (posizione[j] == minpos[j] || posizione[j] == maxpos[j]) countboundary++;
    }
  }

  for (j = 0; j < 3; j++)
  {
    box[j] = (maxpos[j] - minpos[j])/2;
  }
  if (box[1] < box[0])
  {
    temp = box[0];
    box[0] = box[1];
    box[1] = temp;
  }
  if (box[2] < box[0])
  {
    temp = box[0];
    box[0] = box[2];
    box[2] = temp;
  }
  if (box[2] < box[1])
  {
    temp = box[1];
    box[1] = box[2];
    box[2] = temp;
  }

  return (countboundary);
}

/*
 * various types of tiles
 */

int
isstraight (int in, int out)
{
  int indir = in & DIRMASK;
  int outdir = out & DIRMASK;

  switch (indir)
  {
    case ALTO:
    case BASSO:
      if (outdir + indir == ALTO + BASSO) return (1);
      return (0);
      break;

    case DX:
    case SX:
      if (outdir == indir) return (1);
      return (0);
      break;
  }
  assert (0);
  return (1);
}

int
isflap (int in, int out)
{
  int indir = in & DIRMASK;
  int outdir = out & DIRMASK;

  switch (indir)
  {
    case ALTO:
    case BASSO:
      if (outdir == indir) return (1);
      return (0);
      break;

    case DX:
    case SX:
      if (outdir + indir == DX + SX) return (1);
      return (0);
      break;
  }
  assert (0);
  return (1);
}

int
iscurving (int in, int out)
{
  int indir = in & DIRMASK;
  int outdir = out & DIRMASK;

  switch (indir)
  {
    case ALTO:
    case BASSO:
      if (outdir == DX || outdir == SX) return (1);
      return (0);
      break;

    case DX:
    case SX:
      if (outdir == ALTO || outdir == BASSO) return (1);
      return (0);
      break;
  }
  assert (0);
  return (1);
}

int
iscurvingleft (int in, int out)
{
  int indir = in & DIRMASK;
  int outdir = out & DIRMASK;

  switch (indir)
  {
    case ALTO:
      if (outdir == DX) return (1);
      return (0);
      break;

    case BASSO:
      if (outdir == SX) return (1);
      return (0);
      break;

    case DX:
      if (outdir == ALTO) return (1);
      return (0);
      break;

    case SX:
      if (outdir == BASSO) return (1);
      return (0);
      break;
  }
  assert (0);
  return (1);
}

int
iscurvingright (int in, int out)
{
  int indir = in & DIRMASK;
  int outdir = out & DIRMASK;

  switch (indir)
  {
    case ALTO:
      if (outdir == SX) return (1);
      return (0);
      break;

    case BASSO:
      if (outdir == DX) return (1);
      return (0);
      break;

    case DX:
      if (outdir == BASSO) return (1);
      return (0);
      break;

    case SX:
      if (outdir == ALTO) return (1);
      return (0);
      break;
  }
  assert (0);
  return (1);
}

int
ishorizontal (int in)
{
  int indir = in & DIRMASK;

  if (indir == DX || indir == SX) return (1);
  return (0);
}

int
isvertical (int in)
{
  int indir = in & DIRMASK;

  if (indir == ALTO || indir == BASSO) return (1);
  return (0);
}

int
isascending (int in, int out)
{
  int infold = in & FOLDMASK;
  int outfold = out & FOLDMASK;

  if (checkorientable) assert (infold != outfold);
  if (infold == VF) return (0);
  if (infold == MF) return (1);
  if (outfold == VF) return (1);
  return (0);
}

int
isdescending (int in, int out)
{
  int infold = in & FOLDMASK;
  int outfold = out & FOLDMASK;

  if (checkorientable) assert (infold != outfold);
  if (infold == MF) return (0);
  if (infold == VF) return (1);
  if (outfold == MF) return (1);
  return (0);
}

int
deltatile (int in, int out, int type)
{
  if (isstraight (in, out)) return (0);
  if (iscurvingleft (in, out))
  {
    if (type == ISSLASH) return (-1);
    return (1);
  }
  if (iscurvingright (in, out))
  {
    if (type == ISSLASH) return (1);
    return (-1);
  }
  assert (isflap (in, out));

  if (ishorizontal (in))
  {
    if (isascending (in, out)) return (-2);
    return (2);
  }

  assert (isvertical (in));

  if (isascending (in, out)) return (2);
  return (-2);
}

int
computedelta (sequence *seq)
{
  int i, type, in, out;
  int delta = 0;

  type = seq->type;
  for (i = 0; i < seqlen; i++)
  {
    in = seq->seq[i];
    out = seq->seq[(i+1)%seqlen];
    type = ISSLASH + ISBSLASH - type;
    delta += deltatile (in, out, type);
  }
  return (delta);
}

/*
 * check if the sequence is equivalent
 * to its "inverted" version (change of type)
 */

int istypeinvariant (sequence *seq)
{
  sequence canonseq, canoninvseq;

  seqcpy (&canonseq, seq);
  canonify (&canonseq);
  seqcpy (&canoninvseq, seq);
  canoninvseq.type = ISSLASH + ISBSLASH - seq->type;
  canonify (&canoninvseq);

  if (seqcmp (&canonseq, &canoninvseq) == 0) return (1);
  return (0);
}

int
computelinking (sequence *seq)
{
  return (computelt(seq) + computelf(seq) + computelc(seq));
}

/*
 * compute the lt contribution to L
 */

int
computelt (sequence *seq)
{
  int iposizione[3] = {1,1,0};
  int ifreccia[3] = {0,1,0};
  int inormale[3] = {0,0,1};
  int posizione[3];
  int nposizione[3];
  int hinge[3], freccia[3], nfreccia[3], normale[3], nnormale[3];
  int i, type, in, out, dir, fold, val;
  int twicelt = 0;

  veccpy (posizione, iposizione);
  veccpy (freccia, ifreccia);
  veccpy (normale, inormale);
  type = seq->type;
  for (i = 0; i < seqlen; i++)
  {
    in = seq->seq[i];
    out = seq->seq[(i+1)%seqlen];
    type = ISSLASH + ISBSLASH - type;
    dir = in & DIRMASK;
    fold = in & FOLDMASK;
    newposition (dir, fold, posizione, freccia, normale,
                 nposizione, nfreccia, nnormale, hinge);
    veccpy (posizione, nposizione);
    veccpy (freccia, nfreccia);
    veccpy (normale, nnormale);
    if (isstraight (in, out))
    {
      val = 1;
      if (type == ISBSLASH) val = -val;
      if (ishorizontal (in)) val = -val;
      //if (normale[0] + normale[1] + normale[2] < 0) val = -val;
      twicelt += val;
    }
    /*
     * the flap contribution takes into account:
     * 1. the implicit selfcrossing that would be present
     *    when delta is +2 (tile completely wrapped by
     *    the ribbon).
     * 2. performs a "fake" twist in order to make the
     *    computation of the "fold" contribution easier to
     *    compute (no special treating of flaps).
     */
    if (isflap (in, out))
    {
      val = 1;
      if ((posizione[0] & 1) == 0)  // flap with normal in direction x
      {
        if (hinge[1] > posizione[1] || hinge[2] < posizione[2]) val = -val;
      }
      if ((posizione[1] & 1) == 0)  // flap with normal in direction y
      {
        if (hinge[2] > posizione[2] || hinge[0] < posizione[0]) val = -val;
      }
      if ((posizione[2] & 1) == 0)  // flap with normal in direction z
      {
        if (hinge[0] > posizione[0] || hinge[1] < posizione[1]) val = -val;
      }
//printf ("FLAP: %s - %s - %s\n", (type==ISSLASH)?"slash":"bslash",
//                                isvertical(in)?"vert":"hor",
//                                isascending(in,out)?"asc":"desc");
//printf ("      val = %d\n", val);
      twicelt += val;
    }
  }
  assert ((twicelt & 1) == 0);
  return (twicelt/2);
}

/*
 * Lf is the contribution to L coming from "almost" twists, i.e.
 * mountain or valley folds that revert the "seen" face from a
 * selected point of view.
 *
 * The chosen perspective is as follows
 *
 *    |y
 *    |____x_
 *  z/      /|
 *  /______/ |
 *  | |____|_|
 *  |/     |/
 *  /______|
 *
 * flaps are not treated in a special way, since the fake contribution
 * of flaps to lt balances the result of special treatment of flaps
 * also recall that selfcrossings due to flaps are also included in Lt
 *
 * Moreover the ribbon is deformed such that it connects the center
 * of a face with the center of a hinge.
 */

int
computelf (sequence *seq)
{
  int iposizione[3] = {1,1,0};
  int ifreccia[3] = {0,1,0};
  int inormale[3] = {0,0,1};
  int posizione[3];
  int nposizione[3];
  int hinge[3], freccia[3], nfreccia[3], normale[3], nnormale[3];
  int i, in, dir, fold;
  int f1x, f1y, f2x, f2y, hx, hy, ddd;
  int type, twicelf = 0;

  veccpy (posizione, iposizione);
  veccpy (freccia, ifreccia);
  veccpy (normale, inormale);
  type = seq->type;
  for (i = 0; i < seqlen; i++)
  {
    type = ISSLASH + ISBSLASH - type;
    in = seq->seq[i];
    dir = in & DIRMASK;
    fold = in & FOLDMASK;
    newposition (dir, fold, posizione, freccia, normale,
                 nposizione, nfreccia, nnormale, hinge);
    f1x = 2*posizione[0] - posizione[2];  // perspective projection
    f1y = 2*posizione[1] - posizione[2];
    hx = 2*hinge[0] - hinge[2];
    hy = 2*hinge[1] - hinge[2];
    f2x = 2*nposizione[0] - nposizione[2];
    f2y = 2*nposizione[1] - nposizione[2];
    veccpy (posizione, nposizione);
    veccpy (freccia, nfreccia);
    veccpy (normale, nnormale);
    if (fold == NF) continue;  // straight -> no fold twist
    if (hinge[0] & 1)  /* hinge is in the x direction */
    {
      ddd = (f2y - hy)*(hy - f1y);
      assert (ddd != 0);
      if (ddd > 0) continue;
//printf ("x hinge\n");
      if (f2y < hy) twicelf--;
        else twicelf++;
      continue;
    }
    if (hinge[1] & 1)  /* hinge is in the y direction */
    {
      ddd = (f2x - hx)*(hx - f1x);
      assert (ddd != 0);
      if (ddd > 0) continue;
//printf ("y hinge\n");
      if (f2x < hx) twicelf++;
        else twicelf--;
      continue;
    }
    assert (hinge[2] & 1);  /* hinge in the z direction */
    ddd = (f2y - f2x - hy + hx)*(hy - hx - f1y + f1x);
    assert (ddd != 0);
    if (ddd > 0) continue;
//printf ("z hinge\n");
//printf ("f2y - f2x - hy + hx = %d\n", f2y - f2x - hy + hx);
//printf ("hy - hx - f1y + f1x = %d\n", hy - hx - f1y + f1x);
    if (f2y - f2x - hy + hx < 0) twicelf++;
      else twicelf--;
    continue;
  }
  assert ((twicelf & 1) == 0);
  return (twicelf/2);
}

int
computelc (sequence *seq)
{
  int lc = 0;
  int i, j, ii, jj, rl;
  int a[3], b[3], c[3], d[3], ap[3], cp[3];
  static int xrib[3*SEQDIM];
  static int yrib[3*SEQDIM];
  static int zrib[3*SEQDIM];

  ribboncrossings = 0;
  rl = ribbonpath (seq, xrib, yrib, zrib, 3*SEQDIM);
  if (verbose >= 2) printribbonpath (xrib, yrib, zrib, rl);
  checkribbonpath (xrib, yrib, zrib, rl);
  for (i = 0; i <= rl; i++)
  {
    xrib[i] *= 2;
    yrib[i] *= 2;
    zrib[i] *= 2;
  }
  for (i = 0; i < rl - 1; i++)
  {
    a[0] = xrib[i];
    a[1] = yrib[i];
    a[2] = zrib[i];
    b[0] = xrib[i+1];
    b[1] = yrib[i+1];
    b[2] = zrib[i+1];
    for (j = i + 1; j < rl; j++)
    {
      c[0] = xrib[j];
      c[1] = yrib[j];
      c[2] = zrib[j];
      d[0] = xrib[j+1];
      d[1] = yrib[j+1];
      d[2] = zrib[j+1];
      if (a[0] == c[0] && a[1] == c[1] && a[2] == c[2])
      {
        if ((a[2] & 1)) continue; // vertical touching does not generate crossings
        ribbontouchings++;
        ii = i - 1;
        if (ii < 0) ii += rl;
        jj = j - 1;
        if (jj < 0) jj += rl;
        ap[0] = xrib[ii];
        ap[1] = yrib[ii];
        ap[2] = zrib[ii];
        cp[0] = xrib[jj];
        cp[1] = yrib[jj];
        cp[2] = zrib[jj];
        if (ap[2] + b[2] > 2*c[2])
          c[2]--;
         else
          c[2]++; // lift (or lower) the middle point, this does not create selfcrossings
        lc += ribboncrossing (ap, a, cp, c);
        lc += ribboncrossing (a, b, c, d);
        continue;
      }
      if (j == i+1) continue;
      if (i == 0 && j == rl - 1) continue;
      lc += ribboncrossing (a, b, c, d);
    }
  }

  return (lc);
}

/*
 * the ribbon is deformed to connect the center of faces to the
 * center of hinges
 */

int
ribbonpath (sequence *seq, int *x, int *y, int *z, int bufdim)
{
  int iposizione[3] = {1,1,0};
  int ifreccia[3] = {0,1,0};
  int inormale[3] = {0,0,1};
  int posizione[3];
  int nposizione[3];
  int hinge[3], freccia[3], nfreccia[3], normale[3], nnormale[3];
  int dir, fold, i, in, out, type;
  int ribind = 0;

  type = seq->type;

  veccpy (posizione, iposizione);
  veccpy (freccia, ifreccia);
  veccpy (normale, inormale);
  for (i = 0; i < seqlen; i++)
  {
    in = seq->seq[i];
    out = seq->seq[(i+1)%seqlen];
    type = ISSLASH + ISBSLASH - type;
    dir = in & DIRMASK;
    fold = in & FOLDMASK;
    newposition (dir, fold, posizione, freccia, normale,
                 nposizione, nfreccia, nnormale, hinge);
    veccpy (posizione, nposizione);
    veccpy (freccia, nfreccia);
    veccpy (normale, nnormale);
    if (isflap (in, out)) continue;
    x[ribind] = hinge[0];
    y[ribind] = hinge[1];
    z[ribind] = hinge[2];
    ribind++;
    assert (ribind < bufdim);
    x[ribind] = posizione[0];
    y[ribind] = posizione[1];
    z[ribind] = posizione[2];
    ribind++;
  }
  assert (ribind < bufdim);
  x[ribind] = x[0];
  y[ribind] = y[0];
  z[ribind] = z[0];
  return (ribind);
}

void
printribbonpath (int *x, int *y, int *z, int rlen)
{
  int i;

  for (i = 0; i <= rlen; i++)
  {
    printf ("(%d,%d,%d)-", x[i], y[i], z[i]);
  }
  printf ("\n");
}

void
checkribbonpath (int *x, int *y, int *z, int rlen)
{
  int i, dx, dy, dz;

  assert (2*seqlen >= rlen);
  assert (x[0] == x[rlen]);
  assert (y[0] == y[rlen]);
  assert (z[0] == z[rlen]);

  for (i = 0; i < rlen; i++)
  {
    dx = abs (x[i+1] - x[i]);
    dy = abs (y[i+1] - y[i]);
    dz = abs (z[i+1] - z[i]);
    assert (dx + dy + dz == 1);
  }
}

#define F1 7
#define F2 2
#define F3 1

int
ribboncrossing (int a[3], int b[3], int c[3], int d[3])
{
  int abx, aby, dcx, dcy, acx, acy;
  int result, detsign, det, dett, dets;
  int maxabx, minabx, maxaby, minaby;
  int maxcdx, mincdx, maxcdy, mincdy;
  int qab, qcd;

  if (veceq (a, c)) return (0); // il caso di nastro che si tocca e' gestito a monte
  if (veceq (a, d)) return (0);
  if (veceq (b, c)) return (0);
  if (veceq (b, d)) return (0);

  /* compute perspective */
  abx = F1*b[0] - F2*b[2] - F1*a[0] + F2*a[2];
  aby = F1*b[1] - F3*b[2] - F1*a[1] + F3*a[2];
  dcx = F1*c[0] - F2*c[2] - F1*d[0] + F2*d[2];
  dcy = F1*c[1] - F3*c[2] - F1*d[1] + F3*d[2];
  acx = F1*c[0] - F2*c[2] - F1*a[0] + F2*a[2];
  acy = F1*c[1] - F3*c[2] - F1*a[1] + F3*a[2];

  minabx = maxabx = abx; minaby = maxaby = aby;
  if (abx > 0) minabx = 0;
  if (aby > 0) minaby = 0;
  if (abx < 0) maxabx = 0;
  if (aby < 0) maxaby = 0;

  mincdx = maxcdx = acx; mincdy = maxcdy = acy;
  if (dcx > 0) mincdx -= dcx;
  if (dcy > 0) mincdy -= dcy;
  if (dcx < 0) maxcdx -= dcx;
  if (dcy < 0) maxcdy -= dcy;

  if (minabx > maxcdx) return (0);
  if (minaby > maxcdy) return (0);
  if (maxabx < mincdx) return (0);
  if (maxaby < mincdy) return (0);

  det = abx*dcy - aby*dcx;
  dets  = abx*acy - aby*acx;
  dett  = acx*dcy - acy*dcx;

  if (det == 0)
  {
    assert (dets != 0 || dett != 0);
    return (0);
  }
  detsign = 1;
  if (det < 0)
  {
    det = -det;
    detsign = -1;
    dets = -dets;
    dett = -dett;
  }

  if (dets < 0) return (0);
  if (dett < 0) return (0);
  if (dets > det) return (0);
  if (dett > det) return (0);

  assert (dets > 0 && dets < det);
  assert (dett > 0 && dett < det);

  /* we have a clean crossing */
  ribboncrossings++;
  qab = det*a[2] + dett*(b[2]-a[2]);
  qcd = det*c[2] + dets*(d[2]-c[2]);
  assert (qab != qcd);
  result = detsign;
  if (qab > qcd) result = -detsign;

  return (result);
}

/*
 * computation or polyominoid
 */

int
polyominoid (sequence *seq, int polyo[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX])
{
  int iposizione[3] = {1,1,0};
  int ifreccia[3] = {0,1,0};
  int inormale[3] = {0,0,1};
  int posizione[3];
  int nposizione[3];
  int hinge[3], freccia[3], nfreccia[3], normale[3], nnormale[3];
  int minpos[3];
  int maxpos[3];
  int dir, fold, i, j, x, y, z, numsym;

  veccpy (posizione, iposizione);
  veccpy (freccia, ifreccia);
  veccpy (normale, inormale);
  veccpy (minpos, iposizione);
  veccpy (maxpos, iposizione);
  for (i = 0; i < seqlen; i++)
  {
    dir = seq->seq[i] & DIRMASK;
    fold = seq->seq[i] & FOLDMASK;
    newposition (dir, fold, posizione, freccia, normale,
                 nposizione, nfreccia, nnormale, hinge);
    veccpy (posizione, nposizione);
    veccpy (freccia, nfreccia);
    veccpy (normale, nnormale);
    for (j = 0; j < 3; j++)
    {
      if (posizione[j] < minpos[j]) minpos[j] = posizione[j];
      if (posizione[j] > maxpos[j]) maxpos[j] = posizione[j];
    }
  }

  for (i = 0; i < 3; i++)
  {
    if ((minpos[i] & 1)) minpos[i]--;
    if ((maxpos[i] & 1)) maxpos[i]++;
    if (maxpos[i] - minpos[i] >= MAXDIMBOX) return (0);
  }

  polyominoidzero (polyo);
  veccpy (posizione, iposizione);
  veccpy (freccia, ifreccia);
  veccpy (normale, inormale);
  for (i = 0; i < seqlen; i++)
  {
    dir = seq->seq[i] & DIRMASK;
    fold = seq->seq[i] & FOLDMASK;
    newposition (dir, fold, posizione, freccia, normale,
                 nposizione, nfreccia, nnormale, hinge);
    veccpy (posizione, nposizione);
    veccpy (freccia, nfreccia);
    veccpy (normale, nnormale);
    x = posizione[0] - minpos[0];
    y = posizione[1] - minpos[1];
    z = posizione[2] - minpos[2];
    assert (polyo[x][y][z] == 0);
    polyo[x][y][z]++;
  }

  for (i = 0; i < 3; i++) maxpos[i] = maxpos[i] - minpos[i] + 1;

  numsym = optimalpolyominoid (polyo, maxpos);

  return (numsym);
}

int
optimalpolyominoid (int polyo[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX], int boxdim[3])
{
  int polyomod[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX];
  int polyoopt[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX];
  int countsym = 0;
  int xneg, yneg, zneg, perm;

  if (findpolyogroup) ispolyochiral = 1;

  polyocpy (polyoopt, polyo);

  for (xneg = 0; xneg < 2; xneg++)
  {
    for (yneg = 0; yneg < 2; yneg++)
    {
      for (zneg = 0; zneg < 2; zneg++)
      {
        for (perm = 0; perm < 6; perm++)
        {
          polyominoidsym (polyomod, polyo, perm, xneg, yneg, zneg, boxdim);
          if (polyominoidcmp (polyomod, polyo) == 0)
          {
            if (findpolyogroup)
            { /* the following works provided the parity of the permutation is consistent
               * with the parity of the variable "perm"
               */
              if ((xneg + yneg + zneg + perm) % 2 == 1) ispolyochiral = 0;
              if (debug >= 2) printf ("%d: %d %d %d %d\n", countsym, xneg, yneg, zneg, perm);
            }
            countsym++;
          }
          if (polyominoidcmp (polyomod, polyoopt) < 0)
          {
            polyocpy (polyoopt, polyomod);
          }
        }
      }
    }
  }

  /* note: boxdim is invalid on new polyominoid */
  polyocpy (polyo, polyoopt);
  // if (findpolyogroup && ischiral) printf ("chiral\n");
  return (countsym);
}

#define PERM_ID 0
#define PERM_BC 1
#define PERM_AC 3
#define PERM_AB 5
#define PERM_ABC 2
#define PERM_ACB 4

void
polyominoidsym (int polyomod[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX], int polyo[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX],
                int perm, int xneg, int yneg, int zneg, int boxdim[3])
{
  int i, j, k, ii, jj, kk, val;

  polyominoidzero (polyomod);

  for (i = 0; i < boxdim[0]; i++)
  {
    for (j = 0; j < boxdim[1]; j++)
    {
      for (k = 0; k < boxdim[2]; k++)
      {
        ii = i; jj = j; kk = k;
        if (xneg) ii = boxdim[0] - i - 1;
        if (yneg) jj = boxdim[1] - j - 1;
        if (zneg) kk = boxdim[2] - k - 1;
        val = polyo[ii][jj][kk];
        switch (perm)
        {
          case PERM_ID:
            polyomod[i][j][k] = val;
            break;

          case PERM_AB:
            polyomod[j][i][k] = val;
            break;

          case PERM_BC:
            polyomod[i][k][j] = val;
            break;

          case PERM_AC:
            polyomod[k][j][i] = val;
            break;

          case PERM_ABC:
            polyomod[j][k][i] = val;
            break;

          case PERM_ACB:
            polyomod[k][i][j] = val;
            break;

          default:
            assert (0);
        }
      }
    }
  }
}

void
polyominoidzero (int polyo[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX])
{
  int i, j, k;

  for (i = 0; i < MAXDIMBOX; i++)
  {
    for (j = 0; j < MAXDIMBOX; j++)
    {
      for (k = 0; k < MAXDIMBOX; k++)
      {
        polyo[i][j][k] = 0;
      }
    }
  }
}

void
polyocpy (int polyores[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX], int polyo[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX])
{
  int i, j, k;

  for (i = 0; i < MAXDIMBOX; i++)
  {
    for (j = 0; j < MAXDIMBOX; j++)
    {
      for (k = 0; k < MAXDIMBOX; k++)
      {
        polyores[i][j][k] = polyo[i][j][k];
      }
    }
  }
}

int
polyominoidcmp (int polyo1[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX], int polyo2[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX])
{
  int i, j, k;

  for (i = 0; i < MAXDIMBOX; i++)
  {
    for (j = 0; j < MAXDIMBOX; j++)
    {
      for (k = 0; k < MAXDIMBOX; k++)
      {
        if (polyo1[i][j][k] > polyo2[i][j][k]) return (-1);
        if (polyo1[i][j][k] < polyo2[i][j][k]) return (1);
      }
    }
  }
  return (0);
}

void
printcharint (int i)
{
  assert (i >= 0);
  if (i <= 9) {printf ("%c", i + '0'); return;}
  i -= 10;
  if (i <= 26) {printf ("%c", i + 'a'); return;}
  printf ("?");
}

void
printpolyominoid (int polyo[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX])
{
  int i, j, k;
  int first = 1;

  for (i = 0; i < MAXDIMBOX; i++)
  {
    for (j = 0; j < MAXDIMBOX; j++)
    {
      for (k = 0; k < MAXDIMBOX; k++)
      {
        if (polyo[i][j][k])
        {
          if (first == 0) printf ("-");
          first = 0;
          printcharint (i);
          printcharint (j);
          printcharint (k);
        }
      }
    }
  }
}

void
readpolyominoid (char *s, int polyo[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX])
{
  while (isspace(*s)) s++;
  int i, x, y, z;

  assert (MAXDIMBOX >= 9);
  polyominoidzero (polyo);
  for (i = 0; i < seqlen; i++)
  {
    if (i > 0) assert (*s++ == '-');
    assert (*s);
    x = *s++ - '0';
    assert (x >= 0 && x < 9);
    assert (*s);
    y = *s++ - '0';
    assert (y >= 0 && y < 9);
    assert (*s);
    z = *s++ - '0';
    assert (z >= 0 && z < 9);
    assert (((x + y + z) & 1) == 0);
    polyo[x][y][z] = 1;
  }
  if (*s == 's') s++;
  assert (*s == 0 || isspace(*s));
}

void
traslatepolyominoid (int polyo[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX], int maxpos[3])
{
  int i, x, xx, y, yy, z, zz;
  int minpos[3];

  minpos[0] = minpos[1] = minpos[2] = MAXDIMBOX;
  maxpos[0] = maxpos[1] = maxpos[2] = 0;
  for (x = 0; x < MAXDIMBOX; x++)
  {
    for (y = 0; y < MAXDIMBOX; y++)
    {
      for (z = 0; z < MAXDIMBOX; z++)
      {
        if (polyo[x][y][z])
        {
          if (x > maxpos[0]) maxpos[0] = x;
          if (x < minpos[0]) minpos[0] = x;
          if (y > maxpos[1]) maxpos[1] = y;
          if (y < minpos[1]) minpos[1] = y;
          if (z > maxpos[2]) maxpos[2] = z;
          if (z < minpos[2]) minpos[2] = z;
        }
      }
    }
  }

  for (i = 0; i < 3; i++)
  {
    if ((minpos[i] & 1)) minpos[i]--;
    if ((maxpos[i] & 1)) maxpos[i]++;
  }

  for (i = 0; i < 3; i++)
  {
    maxpos[i] = maxpos[i] - minpos[i] + 1;
  }

  if (minpos[0] + minpos[1] + minpos[2] == 0) return;

  for (x = 0; x < MAXDIMBOX; x++)
  {
    for (y = 0; y < MAXDIMBOX; y++)
    {
      for (z = 0; z < MAXDIMBOX; z++)
      {
        xx = x+minpos[0];
        yy = y+minpos[1];
        zz = z+minpos[2];
        polyo[x][y][z] = 0;
        if (xx >= MAXDIMBOX || yy >= MAXDIMBOX || zz >= MAXDIMBOX) continue;
        polyo[x][y][z] = polyo[xx][yy][zz];
      }
    }
  }

  return;
}

#define COLORREDR 255
#define COLORREDG 0
#define COLORREDB 0
#define COLORCYANR 0
#define COLORCYANG 255
#define COLORCYANB 255
#define COLORSILVERR 192
#define COLORSILVERG 192
#define COLORSILVERB 192
#define COLORYELLOWR 255
#define COLORYELLOWG 255
#define COLORYELLOWB 0
#define COLORGREENR 0
#define COLORGREENG 255
#define COLORGREENB 0
#define COLORBLUER 0
#define COLORBLUEG 0
#define COLORBLUEB 255
#define COLORMAGENTAR 255
#define COLORMAGENTAG 0
#define COLORMAGENTAB 255
#define COLORGRAYR 128
#define COLORGRAYG 128
#define COLORGRAYB 128
#define COLORFLATR 255
#define COLORFLATG 128
#define COLORFLATB 128

void
printoff (sequence *seq, FILE *out)
{
  int polyo[MAXDIMBOX][MAXDIMBOX][MAXDIMBOX];
  int box[3];
  int i, j, k, ii, jj, kk;
  int count = 0;
  int colorr, colorg, colorb;

  polyominoid (seq, polyo);
  computebox (seq, box);

  colorr = COLORFLATR;
  colorg = COLORFLATG;
  colorb = COLORFLATB;
  if (box[0] == 1 && box[1] == 2 && box[2] == 2)
  {
    colorr = COLORREDR;
    colorg = COLORREDG;
    colorb = COLORREDB;
  }
  if (box[0] == 2 && box[1] == 2 && box[2] == 2)
  {
    colorr = COLORCYANR;
    colorg = COLORCYANG;
    colorb = COLORCYANB;
  }
  if (box[0] == 1 && box[1] == 2 && box[2] == 3)
  {
    colorr = COLORSILVERR;
    colorg = COLORSILVERG;
    colorb = COLORSILVERB;
  }
  if (box[0] == 2 && box[1] == 2 && box[2] == 3)
  {
    colorr = COLORYELLOWR;
    colorg = COLORYELLOWG;
    colorb = COLORYELLOWB;
  }
  if (box[0] == 1 && box[1] == 1 && box[2] == 3)
  {
    colorr = COLORGREENR;
    colorg = COLORGREENG;
    colorb = COLORGREENB;
  }
  if (box[0] == 1 && box[1] == 1 && box[2] == 4)
  {
    colorr = COLORBLUER;
    colorg = COLORBLUEG;
    colorb = COLORBLUEB;
  }
  if (box[0] == 1 && box[1] == 1 && box[2] == 2)
  {
    colorr = COLORMAGENTAR;
    colorg = COLORMAGENTAG;
    colorb = COLORMAGENTAB;
  }
  if (box[0] == 1 && box[1] == 3 && box[2] == 3)
  {
    colorr = COLORGRAYR;
    colorg = COLORGRAYG;
    colorb = COLORGRAYB;
  }
  if (box[0] == 0)
  {
    colorr = COLORFLATR;
    colorg = COLORFLATG;
    colorb = COLORFLATB;
  }

  fprintf (out, "OFF\n");
  fprintf (out, "# Polyominoid Rubik's Magic\n");
  fprintf (out, "%d %d 0\n", 4*seqlen, seqlen);
  for (i = 0; i < MAXDIMBOX; i++)
  {
    for (j = 0; j < MAXDIMBOX; j++)
    {
      for (k = 0; k < MAXDIMBOX; k++)
      {
        if (polyo[i][j][k])
        {
          count++;
          for (ii = i-1; ii <= i+1; ii++)
          {
            if ((ii & 1) != 0) continue;
            for (jj = j-1; jj <= j+1; jj++)
            {
              if ((jj & 1) != 0) continue;
              for (kk = k-1; kk <= k+1; kk++)
              {
                if ((kk & 1) != 0) continue;
                fprintf (out, "%d.0 %d.0 %d.0\n", ii/2, jj/2, kk/2);
              }
            }
          }
        }
      }
    }
  }
  assert (count == seqlen);

  for (i = 0; i < seqlen; i++)
  {
    fprintf (out, "4 %d %d %d %d %d %d %d\n", 4*i, 4*i+1, 4*i+3, 4*i+2, colorr, colorg, colorb);
  }
  fprintf (stderr, "You can load file with 'geomview file -c ./geomview.gcl'\n");
}

/*
 * verify invariance of the invariants
 */

int
verifyinvariance (sequence *seq)
{
  int revert, rotate, mirror, inout, spin;
  int delta, linking;
  sequence cseq;

  delta = abs(computedelta (seq));
  linking = abs(computelinking (seq));
  minribboncrossings = ribboncrossings;

  for (revert = 0; revert <=1; revert++)
  {
    for (rotate = 0; rotate < seqlen; rotate++)
    {
      for (mirror = 0; mirror <= 1; mirror++)
      {
        for (inout = 0; inout <= 1; inout++)
        {
          for (spin = 0; spin < 4; spin++)
          {
            action (seq, &cseq, revert, rotate, mirror, inout, spin);
            assert (delta == abs(computedelta (&cseq)));
            assert (linking == abs(computelinking (&cseq)));
            if (ribboncrossings < minribboncrossings) minribboncrossings = ribboncrossings;
          }
        }
      }
    }
  }
  return (1);
}