/* See {OctfEnum.h} */

/* Last edited on 2009-02-10 09:24:41 by stolfi */

#define _GNU_SOURCE
#include <stdio.h>
#include <stdint.h>
#include <assert.h>

#include <fget.h>
#include <vec.h>
#include <affirm.h>

#include <Octf.h>
#include <OctfRep.h>
#include <OctfEnum.h>

#define OctfEnum_C_copyright \
  "Copyright © 1998 Universidade Estadual de Campinas (UNICAMP)"

bool_t EnumRing(Place_t p, StepFunc_t stp, VisitFunc_t visit, Place_vec_t *vP)
  { int nP = 0; /* Visited places are {vP->e[0..nP-1]}, if {vP != NULL}. */ 
    Place_t q = p;
    bool_t stop = FALSE;
    do 
      { stop = ( visit == NULL ? FALSE : visit(q));
        if (vP != NULL) { Place_vec_expand(vP, nP); vP->e[nP] = q; nP++; }
        q = stp(q);
      } 
    while ((q != p) && (! stop));
    if (vP != NULL) { Place_vec_trim(vP, nP); }
    return stop;
  }

#define INIT_QUEUE_SIZE 1024

bool_t EnumOrbits
  ( Place_vec_t root, 
    uint nS, 
    StepFunc_t step[], 
    VisitFunc_t visit[], 
    Place_vec_t *vP
  )
  {
    /* [!!! Rewrite to use a {refset_t} to remember the places that have been seen. !!!] */
    /* [!!! Then delete the {pstack} field from the wedge record. !!!] */
    
    Place_vec_t Q = (vP != NULL ? *vP : Place_vec_new(INIT_QUEUE_SIZE));
    int nQ = 0;
      /* The places visited so far are {Q[0..nQ-1]}, in stratified order. */

    auto void InitMarking(void);
      /* Initializes the tables used by {Marked} and {Mark}. */
    
    auto bool_t Marked(Place_t p);
      /* TRUE iff place {p} has been marked (and therefore visited)
         in this call to {EnumOrbits}. */
    
    auto void Mark(Place_t p);
      /* Marks {p} as visited. Assumes that {p} is not marked. */
    
    auto void DoneMarking(void);
      /* Clears any permanent marks and reclaims storage used by
        {Marked} and {Mark}. */
    
    bool_t stop = FALSE; /* Set to TRUE if any {visit[k]} returns TRUE. */

    auto void EnumOrbitsOfOneRoot(Place_t p, int n);
      /* Enumerates all unvisited places that can be reached from
        {p} by zero or more calls to {step[k]} with {k} in {0..n-1}, and
        appends those places to {Q[0..nQ-1]}. */
    
    /* Initialize the mark tables: */
    InitMarking();
    
    /* Enumerate all roots: */
    int i;
    for (i = 0; (i < root.ne) && (! stop); i++)
      { EnumOrbitsOfOneRoot(root.e[i], nS); }

    /* Release the marking table storage: */
    DoneMarking();
    
    /* Return {Q}, or reclaim its space: */
    if (vP == NULL)
      { free(Q.e); }
    else
      { Place_vec_trim(&Q, nQ); (*vP) = Q; }
      
    /* Return the aborted-enum indicator: */
    return stop;

    /* MAIN PROCEDURE */
    
    auto void EnumOrbitsOfOneRoot(Place_t p, int n);
    /* Like {EnumOrbits}, but restricted to one root,
      and excludes places already reached from 
      previous roots. */
    
    void EnumOrbitsOfOneRoot(Place_t p, int n)
      { if (! Marked(p))
          { /* New place: visit it, mark it: */
            if (visit[n] != NULL) { stop = visit[n](p); } 
            Mark(p); 
            /* Save the queue's state: */
            int nX = nQ; /* Places {Q[nX..nQ-1]} still need {step[n-1]}. */
            /* Add {p} to queue: */
            Place_vec_expand(&Q, nQ);
            Q.e[nQ] = p; nQ++;
            if (stop || (n <= 0)) { return; }
            /* Enumerate lower orbits reachable from {p} by {step[n-1]}: */
            while(TRUE)
              { EnumOrbitsOfOneRoot(p, n-1);
                if (stop || (nX >= nQ)) { return; }
                p = Q.e[nX]; nX++;
                p = step[n-1](p);
              }
          }
      }
             
    /* MARKING PROCEDURES */
    
    int nW;
    uint8_vec_t B;
    Wedge_vec_t W; 
      /* The wedges visited so far are {W[0..nW-1]. For each {i} in
        {0..nW-1}, we must have {W[i]->pstack == i}. A place {p} has been
        visited iff {PWedge(p) = W[i]} for some {i} in {0..nW-1}, and
        the bit number {PBits(p)} of {B[i]} is set. */
    
    void InitMarking(void)
      { B = uint8_vec_new(INIT_QUEUE_SIZE); 
        W = Wedge_vec_new(INIT_QUEUE_SIZE);
        nW = 0; /* No wedges have been seen yet. */
      }
    
    bool_t Marked(Place_t p)
      { Wedge_t pw = PWedge(p);
        uint pn = pw->pstack;
        SRBits_t ps = PBits(p);
        return ((pn < nW) && (W.e[pn] == pw) && (((1 << ps) & B.e[pn]) != 0)); 
      }
    
    void Mark(Place_t p)
      { Wedge_t pw = PWedge(p);
        uint pn = pw->pstack;
        uint8_t wB; /* Places on {pw} that have been visited earlier. */
        if ((pn >= nW) || (W.e[pn] != pw))
          { /* First place visitd on wedge {pw}, add it to the wedge set: */
            pn = nW; nW++; 
            Wedge_vec_expand(&W, pn);
            W.e[pn] = pw;
            /* Set the wedge's {pstack} field to its index in {pw}: */
            pw->pstack = pn;
            wB = 0; 
          }
        else
          { /* Wedge {pw} was visited before (at a different place): */
            wB = B.e[pn]; 
          }
        /* Mark place {p} as visited on {pw}: */
        uint8_vec_expand(&B, pn);
        B.e[pn] = wB | (1 << PBits(p)); 
      }
    
    void DoneMarking(void)
      { /* There is no need to clear the wedge {pstack} fields. */
        free(W.e); free(B.e);
      }
  }

uint DegreeOfNode(Place_t p)
  { uint deg = 0;
    auto bool_t CountEdge(Place_t p);
    EnumIncidencesOfElem(Place_vec_make_desc(&p,1), 0, 1, &CountEdge, NULL);
    return deg;
    bool_t CountEdge(Place_t p) { deg++; return FALSE; }
  }

uint DegreeOfWall(Place_t p)
  { uint n = 0;
    Place_t s = p;
    do { n++; s = NextF(s); } while (s != p);
    return n;
  }

uint DegreeOfEdge(Place_t p)
  { uint n = 0;
    Place_t s = p;
    do { n++; s = NextE(s); } while (s != p);
    return n;
  }

bool_t EnumPlaces(Place_vec_t root, VisitFunc_t visit, Place_vec_t *vP)
  { 
    /* Pick three functions {stp[0..2]} that generate 
      all even walks. Then enumerate all places reachable
      by {<stp[0..2]>}. */
      
    StepFunc_t stp[3];
    VisitFunc_t vis[4];
    vis[0] = visit; stp[0] = &ONext; /* == Flip2.Flip1. */
    vis[1] = visit; stp[1] = &SymVF; /* == Flip0.Flip2. */
    vis[2] = visit; stp[2] = &NextF; /* == Flip3.Flip2. */
    vis[3] = visit;                   /* == next root. */
    
    return EnumOrbits(root, 3, stp, vis, vP);
  }

/* ENUMERATING ELEMENTS */

bool_t EnumElemsEven(Place_vec_t root, uint dimObj, VisitFunc_t visit, Place_vec_t *vP)
  { 
    demand(dimObj <= 3, "invalid dimension for target elements");
    /* Set {stp[0..2]} to flip pairs such that {stp[0..1]}
      generates the even subgroup of {< flp[i] : i != dimObj >}
      and {stp[2]} changes {p.elt[dimObj]}. So, each orbit of {stp[0..1]}
      corresponds to one evenly-reachable total orientation of one
      element of dimension {dimObj}, and vice-versa. */
    StepFunc_t stp[3];
    switch(dimObj)
      {
        case 0: /* Nodes: */
          stp[0] = &Flip21; /* == Flip2.Flip1 = ONext. */
          stp[1] = &Flip32; /* == Flip3.Flip2 = NextF. */
          stp[2] = &Flip03; /* == Flip0.Flip3 = Clock. */
          break;

        case 1: /* Edges: */
          stp[0] = &Flip32; /* == Flip3.Flip2 = NextF. */
          stp[1] = &Flip03; /* == Flip0.Flip3 = Clock. */
          stp[2] = &Flip13; /* == Flip1.Flip3 = SymEC. */
          break;

        case 2: /* Walls: */
          stp[0] = &Flip10; /* == Flip1.Flip0 = PrevE. */
          stp[1] = &Flip03; /* == Flip0.Flip3 = Clock. */
          stp[2] = &Flip02; /* == Flip0.Flip2 = SymVF. */
          break;

        case 3: /* Cells: */
          stp[0] = &Flip21; /* == Flip2.Flip1 = ONext. */
          stp[1] = &Flip10; /* == Flip1.Flip0 = PrevE. */
          stp[2] = &Flip03; /* == Flip0.Flip3 = Clock. */
          break;

        default: assert(FALSE);
      }

    auto bool_t MyVisit(Place_t p);
      /* Appends {p} to {vP} and calls {visitNode(p)}. */
    
    /* We should call {visit(p)} only after a {stp[3]} or a new root: */
    VisitFunc_t vis[4];
    vis[0] = NULL;      
    vis[1] = NULL;      
    vis[2] = MyVisit;
    vis[3] = MyVisit;

    int nP = 0; /* Visited places are {vP->e[0..nP-1]}, if {vP != NULL}. */ 
    
    bool_t stop = EnumOrbits(root, 3, stp, vis, NULL);
    
    /* Trim the visit list to the true size: */
    if (vP != NULL) { Place_vec_trim(vP, nP); }
    
    /* Return the abort indicator: */
    return stop;
    
    bool_t MyVisit(Place_t p)
      { if (vP != NULL) { Place_vec_expand(vP, nP); vP->e[nP] = p; nP++; }
        return (visit != NULL ? visit(p) : FALSE);
      }
  }

bool_t EnumElemsAny(Place_vec_t root, uint dimObj, VisitFunc_t visit, Place_vec_t *vP)
  { demand(dimObj <= 3, "invalid dimension for target elements");
    /* Set {stp[0..3]} to a permutation of {flp[0..3]} so that
      {stp[3]==flp[dimObj]}. So, each orbit of {stp[0..1]} corresponds
      to one element of dimension {dimObj}, and vice-versa. */
    StepFunc_t stp[4];
    stp[0] = &Flip0;
    stp[1] = &Flip1;
    stp[2] = &Flip2;
    stp[3] = &Flip3;
    /* Permute {stp[0..3]} so that {flp[dimObj]} is in {stp[3]}; */
    /* Try to preserve some sort of dual symmetry. */
    uint k = dimObj; /* Invariant: {stp[k] = flp[dimObj]}. */
    if (k < 2) 
      { /* Swap {stp[i]} with {stp[3-i]} for all {i}: */
        StepFunc_t t;
        t = stp[0]; stp[0] = stp[3]; stp[3] = t;
        t = stp[1]; stp[1] = stp[2]; stp[2] = t;
        k = 3 - k;
      }
    /* Now the Invariant holds and {k} is {2} or {3}. */
    if (k == 2)
      { /* Swap even and odd elements of {stp[0..3]}: */
        StepFunc_t t;
        t = stp[0]; stp[0] = stp[1]; stp[1] = t;
        t = stp[3]; stp[3] = stp[2]; stp[2] = t;
        k = 3;
      }
    /* Now the Invariant holds and {k} is {3}. */

    auto bool_t MyVisit(Place_t p);
      /* Appends {p} to {vP} and calls {visitNode(p)}. */
    
    /* We should call {visit(p)} only after a {stp[3]} or a new root: */
    VisitFunc_t vis[5];
    vis[0] = NULL;      
    vis[1] = NULL;      
    vis[2] = NULL;      
    vis[3] = MyVisit;
    vis[4] = MyVisit;

    int nP = 0; /* Visited places are {vP->e[0..nP-1]}, if {vP != NULL}. */ 
    
    bool_t stop = EnumOrbits(root, 4, stp, vis, NULL);
    
    /* Trim the visit list to the true size: */
    if (vP != NULL) { Place_vec_trim(vP, nP); }
    
    /* Return the abort indicator: */
    return stop;
    
    bool_t MyVisit(Place_t p)
      { if (vP != NULL) { Place_vec_expand(vP, nP); vP->e[nP] = p; nP++; }
        return (visit != NULL ? visit(p) : FALSE);
      }
  }

/* ENUMERATING THE STAR AND BOUNDARY OF AN ELEMENT */

bool_t EnumPlacesOfElem(Place_vec_t root, uint dimPiv, VisitFunc_t visit, Place_vec_t *vP)
  { 
    demand(dimPiv <= 3, "invalid dimension for pivot elements");

    /* We set {stp[0..1]} to generators of all the 
      even walks that do not change {p.elt[dim]}.
      Then we enumerate all places reachable through
      {<stp[0..1]>}. */
    StepFunc_t stp[2];
    switch(dimPiv)
      {
        case 0: /* Node: */
          stp[0] = &Flip32; /* == Flip3.Flip2. */
          stp[1] = &Flip21; /* == Flip2.Flip1. */
          break;

        case 1: /* Edge */
          stp[0] = &Flip32; /* == Flip3.Flip2. */
          stp[1] = &Flip30; /* == Flip0.Flip3. */
          break;

        case 2: /* Wall */
          stp[0] = &Flip01; /* == Flip0.Flip1. */
          stp[1] = &Flip03; /* == Flip0.Flip3. */
          break;

        case 3: /* Cell */
          stp[0] = &Flip01; /* == Flip0.Flip1. */
          stp[1] = &Flip12; /* == Flip1.Flip2. */
          break;

        default: assert(FALSE);
      }
    
    auto bool_t MyVisit(Place_t p);
      /* Appends {p} to {vP} and calls {visitNode(p)}. */
    
    /* We should visit {p} only after a {stp[2]} or a new root: */
    VisitFunc_t vis[3];
    vis[0] = MyVisit;      
    vis[1] = MyVisit; 
    vis[2] = MyVisit; 
    
    int nP = 0; /* Visited places are {vP->e[0..nP-1]}, if {vP != NULL}. */ 
    
    /* OK, now enumerate the orbits: */
    bool_t stop = EnumOrbits(root, 2, stp, vis, NULL);
    
    /* Trim the visit list to the true size: */
    if (vP != NULL) { Place_vec_trim(vP, nP); }
    
    /* Return the abort indicator: */
    return stop;
    
    bool_t MyVisit(Place_t p)
      { if (vP != NULL) { Place_vec_expand(vP, nP); vP->e[nP] = p; nP++; }
        return (visit != NULL ? visit(p) : FALSE);
      }
  }

bool_t EnumElemsOfElemEven(Place_vec_t root, uint dimPiv, uint dimObj, VisitFunc_t visit, Place_vec_t *vP)
  { 
    demand(dimPiv <= 3, "invalid dimension for pivot elements");
    demand(dimObj <= 3, "invalid dimension for target elements");
    demand(dimPiv != dimObj, "pivot and target dimensions must be distinct");

    /* Set {stp[0..1]} to generators of all the even walks that do 
      not change {p.elt[dimPiv]}, so that only {stp[2]} changes
      {p.elt[dimObj]}. */
    StepFunc_t stp[3];
    switch(dimPiv)
      {
        case 0: /* Pivot is node: */
          switch(dimObj)
            { case 1: /* edges */ stp[0] = &Flip23; stp[1] = &Flip31; break;
              case 2: /* walls */ stp[0] = &Flip31; stp[1] = &Flip12; break;
              case 3: /* cells */ stp[0] = &Flip12; stp[1] = &Flip23;break;
              default: assert(FALSE);
            }
          break;

        case 1: /* Pivot is edge: */
          switch(dimObj)
            { case 0: /* nodes */ stp[0] = &Flip32; stp[1] = &Flip30; break;
              case 2: /* walls */ stp[0] = &Flip03; stp[1] = &Flip02; break;
              case 3: /* cells */ stp[0] = &Flip20; stp[1] = &Flip23; break;
              default: assert(FALSE);
            }
          break;

        case 2: /* Pivot is wall: */
          switch(dimObj)
            { case 0: /* nodes */ stp[0] = &Flip31; stp[1] = &Flip01; break;
              case 1: /* edges */ stp[0] = &Flip03; stp[1] = &Flip13; break;
              case 3: /* cells */ stp[0] = &Flip10; stp[1] = &Flip30; break;
              default: assert(FALSE);
            }
          break;

        case 3: /* Pivot is cell: */
          switch(dimObj)
            { case 0: /* nodes */ stp[0] = &Flip12; stp[1] = &Flip01; break;
              case 1: /* edges */ stp[0] = &Flip20; stp[1] = &Flip12; break;
              case 2: /* walls */ stp[0] = &Flip01; stp[1] = &Flip20; break;
              default: assert(FALSE);
            }
          break;

        default: assert(FALSE);
      }

    /* We should visit {p} only after a {stp[2]} or a new root: */

    auto bool_t MyVisit(Place_t p);
      /* Appends {p} to {vP} and calls {visitNode(p)}. */
    
    VisitFunc_t vis[3];
    vis[0] = NULL;
    vis[1] = &MyVisit;
    vis[2] = &MyVisit;

    int nP = 0; /* Visited places are {vP->e[0..nP-1]}, if {vP != NULL}. */ 

    bool_t stop = EnumOrbits(root, 2, stp, vis, vP);
    
    /* Trim the visit list to the true size: */
    if (vP != NULL) { Place_vec_trim(vP, nP); }
    
    /* Return the abort indicator: */
    return stop;
    
    bool_t MyVisit(Place_t p)
      { if (vP != NULL) { Place_vec_expand(vP, nP); vP->e[nP] = p; nP++; }
        return (visit != NULL ? visit(p) : FALSE);
      }
  }

bool_t EnumElemsOfElemAny(Place_vec_t root, uint dimPiv, uint dimObj, VisitFunc_t visit, Place_vec_t *vP)
  { 
    demand(dimPiv <= 3, "invalid dimension for pivot elements");
    demand(dimObj <= 3, "invalid dimension for target elements");
    demand(dimPiv != dimObj, "pivot and target dimensions must be distinct");

    /* Set {stp[0..2]} to a permutation of {{ stp[i] : i != dimPiv }},
      with {flp[dimObj]} in {stp[2]}: */
    StepFunc_t stp[4]; /* We will use {stp[3]} while sorting but not in the enumeration. */
    /* Permute {stp[0..3]} so that {flp[dimPiv]} is in {stp[3]}; */
    /* Try to preserve some sort of dual symmetry. */
    uint kPiv = dimPiv; /* Invariant 1: {stp[kPiv] = flp[dimPiv]}. */
    uint kObj = dimObj; /* Invariant 2: {stp[kObj] = flp[dimObj]}. */
    if (kPiv < 2) 
      { /* Swap {stp[i]} with {stp[3-i]} for all {i}: */
        StepFunc_t t;
        t = stp[0]; stp[0] = stp[3]; stp[3] = t;
        t = stp[1]; stp[1] = stp[2]; stp[2] = t;
        kPiv = 3 - kPiv; kObj = 3 - kObj;
      }
    /* Now the Invariants 1,2 hold and {kPiv} is {2} or {3}. */
    if (kPiv == 2)
      { /* Swap even and odd elements of {stp[0..3]}: */
        StepFunc_t t;
        t = stp[0]; stp[0] = stp[1]; stp[1] = t;
        t = stp[3]; stp[3] = stp[2]; stp[2] = t;
        kPiv = 3; kObj = (kObj ^ 1); 
      }
    /* Now the Invariants 1,2 hold and {kPiv} is {3}. */
    if (kObj == 0)
      { /* Swap {stp[0]} with {stp[2]} */
        StepFunc_t t;
        t = stp[0]; stp[0] = stp[2]; stp[2] = t;
        kObj = 2; 
      }
    else if (kObj == 1)
      { /* Swap {stp[1]} with {stp[2]} */
        StepFunc_t t;
        t = stp[1]; stp[1] = stp[2]; stp[2] = t;
        kObj = 2; 
      }
    /* Now the Invariants 1,2 hold and {kPiv} is {3} and {kObj} is 2. */

    /* We should visit {p} only after a {stp[2]} or a new root: */

    auto bool_t MyVisit(Place_t p);
      /* Appends {p} to {vP} and calls {visitNode(p)}. */
    
    VisitFunc_t vis[4];
    vis[0] = NULL;
    vis[1] = NULL;
    vis[2] = &MyVisit;
    vis[3] = &MyVisit;

    int nP = 0; /* Visited places are {vP->e[0..nP-1]}, if {vP != NULL}. */ 

    bool_t stop = EnumOrbits(root, 2, stp, vis, vP);
    
    /* Trim the visit list to the true size: */
    if (vP != NULL) { Place_vec_trim(vP, nP); }
    
    /* Return the abort indicator: */
    return stop;
    
    bool_t MyVisit(Place_t p)
      { if (vP != NULL) { Place_vec_expand(vP, nP); vP->e[nP] = p; nP++; }
        return (visit != NULL ? visit(p) : FALSE);
      }
  }

bool_t EnumIncidencesOfElem(Place_vec_t root, uint dimPiv, uint dimObj, VisitFunc_t visit, Place_vec_t *vP)
  { 
    demand(dimPiv <= 3, "invalid dimension for pivot elements");
    demand(dimObj <= 3, "invalid dimension for target elements");
    demand(dimPiv != dimObj, "pivot and target dimensions must be distinct");

    /* Get two {stp}s that geneerate all even walks that
      do not change {p.elt[dimPiv]}. The function {stp[0]}
      should also preserve {p.elt[dimObj]}. */
    StepFunc_t stp[2];
    switch(dimPiv)
      {
        case 0: /* Pivot is node: */
          switch(dimObj)
            { case 1: /* edges */ stp[0] = &Flip23; stp[1] = &Flip13; break;
              case 2: /* walls */ stp[0] = &Flip13; stp[1] = &Flip23; break;
              case 3: /* cells */ stp[0] = &Flip12; stp[1] = &Flip32; break;
              default: assert(FALSE);
            }
          break;

        case 1: /* Pivot is edge: */
          switch(dimObj)
            { case 0: /* nodes */ stp[0] = &Flip23; stp[1] = &Flip03; break;
              case 2: /* walls */ stp[0] = &Flip30; stp[1] = &Flip32; break;
              case 3: /* cells */ stp[0] = &Flip02; stp[1] = &Flip32; break;
              default: assert(FALSE);
            }
          break;

        case 2: /* Pivot is wall: */
          switch(dimObj)
            { case 0: /* nodes */ stp[0] = &Flip13; stp[1] = &Flip10; break;
              case 1: /* edges */ stp[0] = &Flip30; stp[1] = &Flip10; break;
              case 3: /* cells */ stp[0] = &Flip01; stp[1] = &Flip03; break;
              default: assert(FALSE);
            }
          break;

        case 3: /* Pivot is cell: */
          switch(dimObj)
            { case 0: /* nodes */ stp[0] = &Flip12; stp[1] = &Flip10; break;
              case 1: /* edges */ stp[0] = &Flip02; stp[1] = &Flip01; break;
              case 2: /* walls */ stp[0] = &Flip01; stp[1] = &Flip02; break;
              default: assert(FALSE);
            }
          break;

        default: assert(FALSE);
      }

    /* We should visit {p} only after a {stp[2]} or a new root: */

    auto bool_t MyVisit(Place_t p);
      /* Appends {p} to {vP} and calls {visitNode(p)}. */

    VisitFunc_t vis[3];
    vis[0] = NULL;      
    vis[1] = MyVisit; 
    vis[2] = MyVisit; 

    int nP = 0; /* Visited places are {vP->e[0..nP-1]}, if {vP != NULL}. */ 

    /* OK, now enumerate the orbits: */
    bool_t stop = EnumOrbits(root, 2, stp, vis, NULL);

    /* Trim the visit list to the true size: */
    if (vP != NULL) { Place_vec_trim(vP, nP); }
    
    /* Return the abort indicator: */
    return stop;
    
    bool_t MyVisit(Place_t p)
      { if (vP != NULL) { Place_vec_expand(vP, nP); vP->e[nP] = p; nP++; }
        return (visit != NULL ? visit(p) : FALSE);
      }
  }

#define OctfEnum_C_author \
  "Various Modula-3 enumeration procedures for 3D maps were written 1998" \
  "by Luis Arturo Perez Lozada, UNICAMP, inspired on the " \
  "similar quad-edge tools that had been implemented in Modula-3 by " \
  "J.Stolfi and R.M.Rosi (ca. 1994).\n" \
  "\n" \
  "The 3D map enumeration procedures were entirely rewritten in\n" \
  "jan/2007 by J.Stolfi, as special cases of a general\n" \
  "{EnumOrbits} procedure. The latter was conceived in 1993,\n" \
  "for the textbook /Fundamentos de Geometria Computacional/\n" \
  "written by P.J.Rezende and J.Stolfi for the Escola de Programação."