Quelle digraphs.c Sprache: C

/*******************************************************************************
**
*A  digraphs.c                  GAP package Digraphs          Julius Jonusas
**                                                            James Mitchell
**                                                            Wilf A. Wilson
**                                                            Michael Young
**
**  Copyright (C) 2014-17 - Julius Jonusas, James Mitchell, Wilf A. Wilson,
**  Michael Young
**
**  This file is free software, see the digraphs/LICENSE.
**
*******************************************************************************/

#include "digraphs.h"

#include <stdbool.h>  // for false, true, bool
#include <stdint.h>   // for uint64_t
#include <stdlib.h>   // for NULL, free
#include <string.h>   // for memcpy

#include "bliss-includes.h"   // for bliss stuff
#include "cliques.h"          // for FuncDIGRAPHS_FREE_CLIQUES_DATA
#include "digraphs-config.h"  // for DIGRAPHS_WITH_INCLUDED_BLISS
#include "digraphs-debug.h"   // for DIGRAPHS_ASSERT
#include "homos.h"            // for FuncHomomorphismDigraphsFinder
#include "planar.h"           // for FUNC_IS_PLANAR, . . .
#include "safemalloc.h"       // for safe_malloc

#undef PACKAGE
#undef PACKAGE_BUGREPORT
#undef PACKAGE_NAME
#undef PACKAGE_STRING
#undef PACKAGE_TARNAME
#undef PACKAGE_URL
#undef PACKAGE_VERSION

// GAP level things, imported in InitKernel below
Obj IsDigraph;
Obj IsMultiDigraph;
Obj IsDigraphEdge;
Obj DIGRAPHS_ValidateVertexColouring;
Obj Infinity;
Obj IsSymmetricDigraph;
Obj GeneratorsOfGroup;
Obj AutomorphismGroup;
Obj IsAttributeStoringRepObj;
Obj IsPermGroup;
Obj IsDigraphAutomorphism;
Obj LargestMovedPointPerms;
Obj SmallestMovedPointPerm;
Obj IsClique;
Obj IsTrivial;
Obj Orbit;
Obj Stabilizer;
Obj IsSubset;
Obj OnTuples;
Obj Group;
Obj ClosureGroup;
Obj InfoWarning;

static inline bool IsAttributeStoringRep(Obj o) {
  return (CALL_1ARGS(IsAttributeStoringRepObj, o) == True ? true : false);
}

/*************************************************************************/

static Obj FuncDIGRAPH_NR_VERTICES(Obj self, Obj D) {
  return INTOBJ_INT(DigraphNrVertices(D));
}

Int DigraphNrVertices(Obj D) {
  return LEN_LIST(FuncOutNeighbours(0L, D));
}

static Int RNamOutNeighbours = 0;

Obj FuncOutNeighbours(Obj self, Obj D) {
  if (!RNamOutNeighbours) {
    RNamOutNeighbours = RNamName("OutNeighbours");
  }
  if (CALL_1ARGS(IsDigraph, D) != True) {
    ErrorQuit("expected a digraph, not a %s", (Int) TNAM_OBJ(D), 0L);
  } else if (IsbPRec(D, RNamOutNeighbours)) {
    return ElmPRec(D, RNamOutNeighbours);
  } else {
    ErrorQuit(
        "the `OutNeighbours` component is not set for this digraph,", 0L, 0L);
  }
}

static Obj FuncDIGRAPH_OUT_NEIGHBOURS_FROM_SOURCE_RANGE(Obj self,
                                                        Obj N,
                                                        Obj src,
                                                        Obj ran) {
  DIGRAPHS_ASSERT(LEN_LIST(src) == LEN_LIST(ran));
  Int n = INT_INTOBJ(N);
  if (n == 0) {
    return NEW_PLIST(T_PLIST_EMPTY, 0);
  }
  Obj ret = NEW_PLIST(T_PLIST_TAB, n);
  SET_LEN_PLIST(ret, n);

  for (Int i = 1; i <= n; i++) {
    Obj next = NEW_PLIST(T_PLIST_EMPTY, 0);
    SET_LEN_PLIST(next, 0);
    SET_ELM_PLIST(ret, i, next);
    CHANGED_BAG(ret);
  }

  for (Int i = 1; i <= LEN_LIST(src); i++) {
    Obj list = ELM_PLIST(ret, INT_INTOBJ(ELM_LIST(src, i)));
    ASS_LIST(list, LEN_PLIST(list) + 1, ELM_LIST(ran, i));
    CHANGED_BAG(ret);
  }
  return ret;
}

Int DigraphNrEdges(Obj D) {
  Int nr = 0;
  if (IsbPRec(D, RNamName("DigraphNrEdges"))) {
    return INT_INTOBJ(ElmPRec(D, RNamName("DigraphNrEdges")));
  } else if (IsbPRec(D, RNamName("DigraphSource"))) {
    nr = LEN_LIST(ElmPRec(D, RNamName("DigraphSource")));
  } else {
    Int n    = DigraphNrVertices(D);
    Obj list = FuncOutNeighbours(0L, D);
    for (Int i = 1; i <= n; i++) {
      nr += LEN_LIST(ELM_PLIST(list, i));
    }
  }
  if (IsAttributeStoringRep(D)) {
    AssPRec(D, RNamName("DigraphNrEdges"), INTOBJ_INT(nr));
  }
  return nr;
}

static Obj FuncDIGRAPH_NREDGES(Obj self, Obj D) {
  return INTOBJ_INT(DigraphNrEdges(D));
}

Int DigraphNrAdjacencies(Obj D) {
  Int nr = 0;
  if (IsbPRec(D, RNamName("DigraphNrAdjacencies"))) {
    return INT_INTOBJ(ElmPRec(D, RNamName("DigraphNrAdjacencies")));
  } else {
    Obj const out = FuncOutNeighbours(0L, D);
    for (Int v = 1; v <= LEN_LIST(out); ++v) {
      Obj const out_v = ELM_LIST(out, v);
      for (Int w = 1; w <= LEN_LIST(out_v); ++w) {
        Int u = INT_INTOBJ(ELM_LIST(out_v, w));
        if (v <= u
            || CALL_3ARGS(IsDigraphEdge, D, INTOBJ_INT(u), INTOBJ_INT(v))
                   == False) {
          ++nr;
        }
      }
    }
  }
  if (IsAttributeStoringRep(D)) {
    AssPRec(D, RNamName("DigraphNrAdjacencies"), INTOBJ_INT(nr));
  }
  return nr;
}

static Obj FuncDIGRAPH_NRADJACENCIES(Obj self, Obj D) {
  return INTOBJ_INT(DigraphNrAdjacencies(D));
}

Int DigraphNrAdjacenciesWithoutLoops(Obj D) {
  Int nr = 0;
  if (IsbPRec(D, RNamName("DigraphNrAdjacenciesWithoutLoops"))) {
    return INT_INTOBJ(ElmPRec(D, RNamName("DigraphNrAdjacenciesWithoutLoops")));
  } else {
    Obj const out = FuncOutNeighbours(0L, D);
    for (Int v = 1; v <= LEN_LIST(out); ++v) {
      Obj const out_v = ELM_LIST(out, v);
      for (Int w = 1; w <= LEN_LIST(out_v); ++w) {
        Int u = INT_INTOBJ(ELM_LIST(out_v, w));
        if (v < u
            || CALL_3ARGS(IsDigraphEdge, D, INTOBJ_INT(u), INTOBJ_INT(v))
                   == False) {
          ++nr;
        }
      }
    }
  }
  if (IsAttributeStoringRep(D)) {
    AssPRec(D, RNamName("DigraphNrAdjacenciesWithoutLoops"), INTOBJ_INT(nr));
  }
  return nr;
}

static Obj FuncDIGRAPH_NRADJACENCIESWITHOUTLOOPS(Obj self, Obj D) {
  return INTOBJ_INT(DigraphNrAdjacenciesWithoutLoops(D));
}

/****************************************************************************
**
*F  FuncGABOW_SCC
**
** `digraph' should be a list whose entries and the lists of out-neighbours
** of the vertices. So [[2,3],[1],[2]] represents the graph whose edges are
** 1->2, 1->3, 2->1 and 3->2.
**
** returns a newly constructed record with two components 'comps' and 'id' the
** elements of 'comps' are lists representing the strongly connected components
** of the directed graph, and in the component 'id' the following holds:
** id[i]=PositionProperty(comps, x-> i in x);
** i.e. 'id[i]' is the index in 'comps' of the component containing 'i'.
** Neither the components, nor their elements are in any particular order.
**
** The algorithm is that of Gabow, based on the implementation in Sedgwick:
**   https://algs4.cs.princeton.edu/42directed/GabowSCC.java.html
** (made non-recursive to avoid problems with stack limits) and
** the implementation of STRONGLY_CONNECTED_COMPONENTS_DIGRAPH in listfunc.c.
*/

static Obj FuncGABOW_SCC(Obj self, Obj digraph) {
  UInt end1, end2, count, level, w, v, n, nr, idw, *frames, *stack2;
  Obj  id, stack1, out, comp, comps, adj;

  PLAIN_LIST(digraph);
  n = LEN_PLIST(digraph);
  if (n == 0) {
    out = NEW_PREC(2);
    AssPRec(out, RNamName("id"), NEW_PLIST_IMM(T_PLIST_EMPTY, 0));
    AssPRec(out, RNamName("comps"), NEW_PLIST_IMM(T_PLIST_EMPTY, 0));
    return out;
  }

  end1   = 0;
  stack1 = NEW_PLIST(T_PLIST_CYC, n);
  // stack1 is a plist so we can use memcopy below
  SET_LEN_PLIST(stack1, n);

  id = NEW_PLIST_IMM(T_PLIST_CYC, n);
  SET_LEN_PLIST(id, n);

  // init id
  for (v = 1; v <= n; v++) {
    SET_ELM_PLIST(id, v, INTOBJ_INT(0));
  }

  count = n;

  comps = NEW_PLIST_IMM(T_PLIST_TAB, n);

  stack2 = safe_malloc((4 * n + 2) * sizeof(UInt));
  frames = stack2 + n + 1;
  end2   = 0;

  for (v = 1; v <= n; v++) {
    if (INT_INTOBJ(ELM_PLIST(id, v)) == 0) {
      adj = ELM_PLIST(digraph, v);
      PLAIN_LIST(adj);
      level     = 1;
      frames[0] = v;  // vertex
      frames[1] = 1;  // index
      frames[2] = (UInt) adj;
      SET_ELM_PLIST(stack1, ++end1, INTOBJ_INT(v));
      stack2[++end2] = end1;
      SET_ELM_PLIST(id, v, INTOBJ_INT(end1));

      while (1) {
        if (frames[1] > (UInt) LEN_PLIST((Obj) frames[2])) {
          if (stack2[end2] == (UInt) INT_INTOBJ(ELM_PLIST(id, frames[0]))) {
            end2--;
            count++;
            nr = 0;
            do {
              nr++;
              w = INT_INTOBJ(ELM_PLIST(stack1, end1--));
              SET_ELM_PLIST(id, w, INTOBJ_INT(count));
            } while (w != frames[0]);

            comp = NEW_PLIST_IMM(T_PLIST_CYC, nr);
            SET_LEN_PLIST(comp, nr);

            memcpy(ADDR_OBJ(comp) + 1,
                   CONST_ADDR_OBJ(stack1) + (end1 + 1),
                   nr * sizeof(Obj));

            nr = LEN_PLIST(comps) + 1;
            SET_ELM_PLIST(comps, nr, comp);
            SET_LEN_PLIST(comps, nr);
            CHANGED_BAG(comps);
          }
          level--;
          if (level == 0) {
            break;
          }
          frames -= 3;
        } else {
          w   = INT_INTOBJ(ELM_PLIST((Obj) frames[2], frames[1]++));
          idw = INT_INTOBJ(ELM_PLIST(id, w));

          if (idw == 0) {
            adj = ELM_PLIST(digraph, w);
            PLAIN_LIST(adj);
            level++;
            frames += 3;
            frames[0] = w;  // vertex
            frames[1] = 1;  // index
            frames[2] = (UInt) adj;
            SET_ELM_PLIST(stack1, ++end1, INTOBJ_INT(w));
            stack2[++end2] = end1;
            SET_ELM_PLIST(id, w, INTOBJ_INT(end1));
          } else {
            while (stack2[end2] > idw) {
              end2--;  // pop from stack2
            }
          }
        }
      }
    }
  }

  for (v = 1; v <= n; v++) {
    SET_ELM_PLIST(id, v, INTOBJ_INT(INT_INTOBJ(ELM_PLIST(id, v)) - n));
  }

  out = NEW_PREC(2);
  SHRINK_PLIST(comps, LEN_PLIST(comps));
  AssPRec(out, RNamName("id"), id);
  AssPRec(out, RNamName("comps"), comps);
  free(stack2);
  return out;
}

static UInt UF_FIND(UInt* id, UInt i) {
  while (i != id[i])
    i = id[i];
  return i;
}

static void UF_COMBINE_CLASSES(UInt* id, UInt i, UInt j) {
  i = UF_FIND(id, i);
  j = UF_FIND(id, j);
  if (i < j)
    id[j] = i;
  else if (j < i)
    id[i] = j;
  // if i = j then there is nothing to combine
}

static Obj FuncDIGRAPH_CONNECTED_COMPONENTS(Obj self, Obj digraph) {
  UInt n, *id, *nid, i, e, len, f, nrcomps;
  Obj  adj, adji, gid, gcomps, comp, out;

  out = NEW_PREC(2);
  n   = DigraphNrVertices(digraph);
  if (n == 0) {
    gid    = NEW_PLIST_IMM(T_PLIST_EMPTY, 0);
    gcomps = NEW_PLIST_IMM(T_PLIST_EMPTY, 0);
  } else {
    id = safe_malloc(n * sizeof(UInt));
    for (i = 0; i < n; i++) {
      id[i] = i;
    }

    adj = FuncOutNeighbours(self, digraph);
    for (i = 0; i < n; i++) {
      adji = ELM_PLIST(adj, i + 1);
      PLAIN_LIST(adji);
      len = LEN_PLIST(adji);
      for (e = 1; e <= len; e++) {
        UF_COMBINE_CLASSES(id, i, INT_INTOBJ(ELM_PLIST(adji, e)) - 1);
      }
    }

    // "Normalise" id, giving it sensible labels
    nid     = safe_malloc(n * sizeof(UInt));
    nrcomps = 0;
    for (i = 0; i < n; i++) {
      f      = UF_FIND(id, i);
      nid[i] = (f == i) ? ++nrcomps : nid[f];
    }
    free(id);

    // Make GAP object from nid
    gid    = NEW_PLIST_IMM(T_PLIST_CYC, n);
    gcomps = NEW_PLIST_IMM(T_PLIST_CYC, nrcomps);
    SET_LEN_PLIST(gid, n);
    SET_LEN_PLIST(gcomps, nrcomps);
    for (i = 1; i <= nrcomps; i++) {
      SET_ELM_PLIST(gcomps, i, NEW_PLIST_IMM(T_PLIST_CYC, 0));
      CHANGED_BAG(gcomps);
    }
    for (i = 1; i <= n; i++) {
      SET_ELM_PLIST(gid, i, INTOBJ_INT(nid[i - 1]));
      comp = ELM_PLIST(gcomps, nid[i - 1]);
      len  = LEN_PLIST(comp);
      AssPlist(comp, len + 1, INTOBJ_INT(i));
    }
    free(nid);
  }
  AssPRec(out, RNamName("id"), gid);
  AssPRec(out, RNamName("comps"), gcomps);
  return out;
}

static Obj FuncIS_ACYCLIC_DIGRAPH(Obj self, Obj adj) {
  UInt  nr, i, j, k, level;
  Obj   nbs;
  UInt *stack, *ptr;

  nr = LEN_PLIST(adj);

  // init the buf
  ptr   = safe_calloc(nr + 1, sizeof(UInt));
  stack = safe_malloc((2 * nr + 2) * sizeof(UInt));

  for (i = 1; i <= nr; i++) {
    nbs = ELM_PLIST(adj, i);
    if (LEN_LIST(nbs) == 0) {
      ptr[i] = 1;
    } else if (ptr[i] == 0) {
      level    = 1;
      stack[0] = i;
      stack[1] = 1;
      while (1) {
        j = stack[0];
        k = stack[1];
        if (ptr[j] == 2) {
          free(ptr);
          stack -= (2 * level) - 2;  // put the stack back to the start
          free(stack);
          return False;  // We have just travelled around a cycle
        }
        // Check whether:
        // 1. We've previously finished with this vertex, OR
        // 2. Whether we've now investigated all descendant branches
        nbs = ELM_PLIST(adj, j);
        if (ptr[j] == 1 || k > (UInt) LEN_LIST(nbs)) {
          ptr[j] = 1;
          level--;
          if (level == 0) {
            break;
          }
          // Backtrack and choose next available branch
          stack -= 2;
          ptr[stack[0]] = 0;
          stack[1]++;
        } else {  // Otherwise move onto the next available branch
          ptr[j] = 2;
          level++;
          nbs = ELM_PLIST(adj, j);
          stack += 2;
          stack[0] = INT_INTOBJ(CONST_ADDR_OBJ(nbs)[k]);
          stack[1] = 1;
        }
      }
    }
  }
  free(ptr);
  free(stack);
  return True;
}

static Obj FuncDIGRAPH_LONGEST_DIST_VERTEX(Obj self, Obj adj, Obj start) {
  UInt  nr, i, j, k, level, x, prev;
  Obj   nbs;
  UInt *stack, *ptr, *depth;

  nr = LEN_PLIST(adj);
  i  = INT_INTOBJ(start);

  if (i > nr || i < 1) {
    ErrorQuit("DIGRAPH_LONGEST_DIST_VERTEX: the 2nd "
              "arg must be a vertex of the first,",
              0L,
              0L);
  }

  nbs = ELM_PLIST(adj, i);
  if (LEN_LIST(nbs) == 0) {
    return INTOBJ_INT(0);
  }

  ptr   = safe_calloc(nr + 1, sizeof(UInt));
  depth = safe_calloc(nr + 1, sizeof(UInt));
  stack = safe_malloc((2 * nr + 2) * sizeof(UInt));

  level    = 1;
  stack[0] = i;
  stack[1] = 1;
  prev     = 0;
  while (1) {
    j = stack[0];
    k = stack[1];
    if (ptr[j] == 2) {  // we have identified a cycle
      stack -= (2 * level) - 2;
      free(stack);
      free(ptr);
      free(depth);
      return INTOBJ_INT(-2);  // We have just travelled around a cycle
    }

    if (prev > depth[j]) {
      depth[j] = prev;
    }
    // Check whether:
    // 1. We've previously finished with this vertex, OR
    // 2. Whether we've now investigated all descendant branches
    nbs = ELM_PLIST(adj, j);
    if (ptr[j] == 1 || k > (UInt) LEN_LIST(nbs)) {
      ptr[j] = 1;
      level--;
      prev = depth[j];
      if (level == 0) {
        // finished the search
        break;
      }
      // Backtrack and choose next available branch
      stack -= 2;
      ptr[stack[0]] = 0;
      prev++;
      stack[1]++;
    } else {  // Otherwise move onto the next available branch
      ptr[j] = 2;
      level++;
      nbs = ELM_PLIST(adj, j);
      stack += 2;
      stack[0] = INT_INTOBJ(CONST_ADDR_OBJ(nbs)[k]);
      stack[1] = 1;
      prev     = 0;
    }
  }

  x = depth[INT_INTOBJ(start)];
  free(ptr);
  free(depth);
  free(stack);
  return INTOBJ_INT(x);
}

// Forward decl
static Obj FuncDIGRAPH_IN_OUT_NBS(Obj, Obj);

// Takes in-neighbours (Obj adj) of a topologically sorted non-multi digraph
// Returns the out-neighbours of its transitive reduction.
// If (Obj loops) == False, loops are removed (transitive reflexive reduction)
static Obj FuncDIGRAPH_TRANS_REDUCTION(Obj self, Obj D) {
  DIGRAPHS_ASSERT(CALL_1ARGS(IsDigraph, D) == True);
  if (!IS_MUTABLE_OBJ(D)) {
    ErrorQuit("the argument (a digraph) must be mutable", 0L, 0L);
  }

  UInt const n = DigraphNrVertices(D);

  // Special case for n = 0
  if (n == 0) {
    return D;
  }

  // Create the GAP out-neighbours structure of the result
  Obj ot_list = NEW_PLIST(T_PLIST_TAB, n);
  SET_LEN_PLIST(ot_list, n);
  for (UInt i = 1; i <= n; i++) {
    SET_ELM_PLIST(ot_list, i, NEW_PLIST(T_PLIST_EMPTY, 0));
    CHANGED_BAG(ot_list);
  }

  Obj const in_list = FuncDIGRAPH_IN_OUT_NBS(self, FuncOutNeighbours(self, D));

  // Create data structures needed for computation
  UInt* ptr   = safe_calloc(n + 1, sizeof(UInt));
  bool* mat   = safe_calloc(n * n, sizeof(bool));
  UInt* stack = safe_malloc((2 * n + 2) * sizeof(UInt));

  // Start a depth-first search from each source of the digraph
  for (UInt i = 1; i <= n; i++) {
    if (ptr[i] == 0) {
      // Remember which vertex was the source
      UInt source = i;
      // not sure if this next line is necessary
      bool backtracking = false;
      UInt level        = 1;

      stack[0] = i;
      stack[1] = 1;
      while (1) {
        UInt j = stack[0];
        UInt k = stack[1];

        // We have found a loop on vertex j
        if (ptr[j] == 2) {
          if (stack[-2] != j) {
            ErrorQuit(
                "the argument (a digraph) must be acyclic except for loops,",
                0L,
                0L);
          }
          backtracking = true;
          level--;
          stack -= 2;
          stack[1]++;
          ptr[j] = 0;
          // Store the loop
          Obj list = ELM_PLIST(ot_list, j);
          ASS_LIST(list, LEN_PLIST(list) + 1, INTOBJ_INT(j));
          CHANGED_BAG(ot_list);
          continue;
        }

        // Calculate if we need to add an edge from j -> (previous vertex)
        if (!backtracking && j != source && !mat[(stack[-2] - 1) * n + j - 1]) {
          Obj list = ELM_PLIST(ot_list, j);
          ASS_LIST(list, LEN_PLIST(list) + 1, INTOBJ_INT(stack[-2]));
          CHANGED_BAG(ot_list);
        }

        Obj nbs = ELM_PLIST(in_list, j);

        // Do we need to backtrack?
        if (ptr[j] == 1 || k > (UInt) LEN_LIST(nbs)) {
          if (level == 1)
            break;

          backtracking = true;
          level--;
          stack -= 2;
          ptr[stack[0]] = 0;
          stack[1]++;
          ptr[j] = 1;

          // w is the vertex we are backtracking to (-1)
          UInt w = stack[0] - 1;
          // Record which vertices we have discovered 'above' w
          for (UInt m = 0; m < n; m++) {
            mat[w * n + m] = mat[w * n + m] + mat[(j - 1) * n + m];
          }
          mat[w * n + (j - 1)] = true;
        } else {
          backtracking = false;
          level++;
          stack += 2;
          stack[0] = INT_INTOBJ(CONST_ADDR_OBJ(nbs)[k]);
          stack[1] = 1;
          ptr[j]   = 2;
        }
      }
    }
  }
  free(mat);
  free(ptr);
  free(stack);
  AssPRec(D, RNamName("OutNeighbours"), ot_list);
  return D;
}

// TODO(*) use generic DFS, when we have one.

static Obj FuncDIGRAPH_PATH(Obj self, Obj adj, Obj u, Obj v) {
  UInt  nr, i, j, k, level, target;
  Obj   nbs, out, path, edge;
  UInt *stack, *ptr;

  i   = INT_INTOBJ(u);
  nbs = ELM_PLIST(adj, i);
  if (LEN_LIST(nbs) == 0) {
    return Fail;
  }
  target = INT_INTOBJ(v);
  nr     = LEN_PLIST(adj);

  // init the buf
  ptr   = safe_calloc(nr + 1, sizeof(UInt));
  stack = safe_malloc((2 * nr + 2) * sizeof(UInt));

  level    = 1;
  stack[0] = i;
  stack[1] = 1;
  while (1) {
    j = stack[0];
    k = stack[1];
    // Check whether:
    // 1. We've previously visited with this vertex, OR
    // 2. Whether we've now investigated all descendant branches
    nbs = ELM_PLIST(adj, j);
    if (ptr[j] != 0 || k > (UInt) LEN_LIST(nbs)) {
      ptr[j] = 1;
      if (--level == 0) {
        break;
      }
      // Backtrack and choose next available branch
      stack -= 2;
      ptr[stack[0]] = 0;
      stack[1]++;
    } else {  // Otherwise move onto the next available branch
      ptr[j] = 2;
      level++;
      nbs = ELM_PLIST(adj, j);
      stack += 2;
      stack[0] = INT_INTOBJ(CONST_ADDR_OBJ(nbs)[k]);
      if (stack[0] == target) {
        // Create output lists
        path = NEW_PLIST_IMM(T_PLIST_CYC, level);
        SET_LEN_PLIST(path, level);
        SET_ELM_PLIST(path, level--, INTOBJ_INT(stack[0]));
        edge = NEW_PLIST_IMM(T_PLIST_CYC, level);
        SET_LEN_PLIST(edge, level);
        out = NEW_PLIST_IMM(T_PLIST_CYC, 2);
        SET_LEN_PLIST(out, 2);

        // Fill output lists
        while (level > 0) {
          stack -= 2;
          SET_ELM_PLIST(edge, level, INTOBJ_INT(stack[1]));
          SET_ELM_PLIST(path, level--, INTOBJ_INT(stack[0]));
        }
        SET_ELM_PLIST(out, 1, path);
        SET_ELM_PLIST(out, 2, edge);
        free(ptr);
        free(stack);
        return out;
      }
      stack[1] = 1;
    }
  }
  free(ptr);
  free(stack);
  return Fail;
}

Obj FuncIS_ANTISYMMETRIC_DIGRAPH(Obj self, Obj adj) {
  Int   nr, i, j, k, l, level, last1, last2;
  Obj   nbs;
  UInt *stack, *ptr;

  nr = LEN_PLIST(adj);
  if (nr <= 1) {
    return True;
  }

  // init the buf (is this correct length?)
  ptr   = safe_calloc(nr + 1, sizeof(UInt));
  stack = safe_malloc((4 * nr + 4) * sizeof(UInt));

  for (i = 1; i <= nr; i++) {
    nbs = ELM_PLIST(adj, i);
    if (LEN_LIST(nbs) == 0) {
      ptr[i] = 1;
    } else if (ptr[i] == 0) {
      level    = 1;
      stack[0] = i;
      stack[1] = 1;
      stack[2] = 0;
      stack[3] = 0;
      while (1) {
        j     = stack[0];
        k     = stack[1];
        last1 = stack[2];
        last2 = stack[3];
        if (j == last2 && j != last1) {
          free(ptr);
          stack -= (4 * level) - 4;  // put the stack back to the start
          free(stack);
          return False;  // Found a non-loop cycle of length 2
        }
        // Check whether:
        // 1. We've previously finished with this vertex, OR
        // 2. Whether we've now investigated all descendant branches
        nbs = ELM_PLIST(adj, j);
        if (ptr[j] == 2) {
          PLAIN_LIST(nbs);
          for (l = 1; l <= LEN_PLIST(nbs); l++) {
            if (last1 != j && INT_INTOBJ(CONST_ADDR_OBJ(nbs)[l]) == last1) {
              free(ptr);
              stack -= (4 * level) - 4;  // put the stack back to the start
              free(stack);
              return False;
            }
          }
        }
        if (k > LEN_LIST(nbs)) {
          ptr[j] = 1;
        }
        if (ptr[j] >= 1) {
          level--;
          if (level == 0) {
            break;
          }
          // Backtrack and choose next available branch
          stack -= 4;
          ptr[stack[0]] = 0;
          stack[1]++;
        } else {  // Otherwise move onto the next available branch
          ptr[j] = 2;
          level++;
          nbs = ELM_PLIST(adj, j);
          stack += 4;
          stack[0] = INT_INTOBJ(ELM_LIST(nbs, k));
          stack[1] = 1;
          stack[2] = j;  // I am wasting memory here, duplicating info
          stack[3] = last1;
        }
      }
    }
  }
  free(ptr);
  free(stack);
  return True;
}

static Obj FuncIS_STRONGLY_CONNECTED_DIGRAPH(Obj self, Obj digraph) {
  UInt n, nextid, *bag, *ptr1, *ptr2, *fptr, *id, w;

  n = LEN_PLIST(digraph);
  if (n == 0) {
    return True;
  }

  nextid = 1;
  bag    = safe_malloc(n * 4 * sizeof(UInt));
  ptr1   = bag;
  ptr2   = bag + n;
  fptr   = bag + n * 2;
  id     = safe_calloc(n + 1, sizeof(UInt));

  // first vertex v=1
  PLAIN_LIST(ELM_PLIST(digraph, 1));
  fptr[0] = 1;  // vertex
  fptr[1] = 1;  // index
  *ptr1   = 1;
  *ptr2   = nextid;
  id[1]   = nextid;

  while (1) {  // we always return before level = 0
    if (fptr[1] > (UInt) LEN_PLIST(ELM_PLIST(digraph, fptr[0]))) {
      if (*ptr2 == id[fptr[0]]) {
        do {
          n--;
        } while (*(ptr1--) != fptr[0]);
        free(bag);
        free(id);
        return n == 0 ? True : False;
      }
      fptr -= 2;
    } else {
      w = INT_INTOBJ(ELM_PLIST(ELM_PLIST(digraph, fptr[0]), fptr[1]++));
      if (id[w] == 0) {
        PLAIN_LIST(ELM_PLIST(digraph, w));
        fptr += 2;
        fptr[0] = w;  // vertex
        fptr[1] = 1;  // index
        nextid++;
        *(++ptr1) = w;
        *(++ptr2) = nextid;
        id[w]     = nextid;
      } else {
        while ((*ptr2) > id[w]) {
          ptr2--;
        }
      }
    }
  }
  // this should never happen, just to keep the compiler happy
  return Fail;
}

static Obj FuncDIGRAPH_TOPO_SORT(Obj self, Obj adj) {
  UInt  nr, i, j, k, count;
  UInt  level;
  Obj   nbs, out;
  UInt *stack, *ptr;

  nr = LEN_PLIST(adj);

  if (nr == 0) {
    return NEW_PLIST_IMM(T_PLIST_EMPTY, 0);
  }
  out = NEW_PLIST_IMM(T_PLIST_CYC, nr);
  SET_LEN_PLIST(out, nr);
  if (nr == 1) {
    SET_ELM_PLIST(out, 1, INTOBJ_INT(1));
    return out;
  }

  // init the buf
  ptr   = safe_calloc(nr + 1, sizeof(UInt));
  stack = safe_malloc((2 * nr + 2) * sizeof(UInt));
  count = 0;

  for (i = 1; i <= nr; i++) {
    nbs = ELM_PLIST(adj, i);
    if (LEN_LIST(nbs) == 0) {
      if (ptr[i] == 0) {
        count++;
        SET_ELM_PLIST(out, count, INTOBJ_INT(i));
      }
      ptr[i] = 1;
    } else if (ptr[i] == 0) {
      level    = 1;
      stack[0] = i;
      stack[1] = 1;
      while (1) {
        j = stack[0];
        k = stack[1];
        if (ptr[j] == 2) {
          // SetIsAcyclicDigraph(graph, false);
          stack -= 2;
          level--;
          if (stack[0] != j) {  // Cycle is not just a loop
            free(ptr);
            stack -= (2 * level) - 2;
            free(stack);
            return Fail;
          }
          stack[1]++;
          ptr[j] = 0;
          continue;
        }
        nbs = ELM_PLIST(adj, j);
        if (ptr[j] == 1 || k > (UInt) LEN_LIST(nbs)) {
          if (ptr[j] == 0) {
            // ADD J TO THE END OF OUR LIST
            count++;
            SET_ELM_PLIST(out, count, INTOBJ_INT(j));
          }
          ptr[j] = 1;
          level--;
          if (level == 0) {
            break;
          }
          // Backtrack and choose next available branch
          stack -= 2;
          ptr[stack[0]] = 0;
          stack[1]++;
        } else {  // Otherwise move onto the next available branch
          ptr[j] = 2;
          level++;
          nbs = ELM_PLIST(adj, j);
          stack += 2;
          stack[0] = INT_INTOBJ(ELM_LIST(nbs, k));
          stack[1] = 1;
        }
      }
    }
  }
  free(ptr);
  free(stack);
  return out;
}

// WW. This function performs a depth first search on the digraph defined by
// <adj> and returns the adjacency list <out> of a spanning forest. This is a
// forest rather than a tree since <adj> is not required to be connected. Each
// time a new vertex <j> is discovered (from the previous vertex <i>), the edges
// i <-> j are added to <out>.

// TODO(*) use generic DFS, when we have one.

// Assumes that <adj> is the list of adjacencies of a symmetric digraph
// Multiple edges and loops are allowed
static Obj FuncDIGRAPH_SYMMETRIC_SPANNING_FOREST(Obj self, Obj adj) {
  UInt  nr, i, j, k, next, len, level;
  Obj   nbs, out, out_j, new;
  UInt *stack, *ptr;

  nr = LEN_PLIST(adj);

  if (nr == 0) {
    return NEW_PLIST_IMM(T_PLIST_EMPTY, 0);
  }

  // init the adjacencies of the spanning forest
  out = NEW_PLIST(T_PLIST_TAB, nr);
  SET_LEN_PLIST(out, nr);
  for (i = 1; i <= nr; i++) {
    SET_ELM_PLIST(out, i, NEW_PLIST(T_PLIST_EMPTY, 0));
    CHANGED_BAG(out);
  }

  // init the buffer
  ptr   = safe_calloc(nr + 1, sizeof(UInt));
  stack = safe_malloc((2 * nr + 2) * sizeof(UInt));

  for (i = 1; i <= nr; i++) {
    // perform DFS only on still-undiscovered non-trivial connected components
    if (!ptr[i] && LEN_LIST(ELM_PLIST(adj, i)) > 0) {
      level    = 1;
      stack[0] = i;
      stack[1] = 1;

      while (1) {
        j   = stack[0];
        k   = stack[1];
        nbs = ELM_PLIST(adj, j);

        // idea: is <nbs[k]> a new vertex? add edges <j> <-> <nbs[k]> if so

        // if we're finished with <j>, or <nbs[k]> doesn't exist, then backtrack
        if (ptr[j] || k > (UInt) LEN_LIST(nbs)) {
          ptr[j] = 1;  // vertex <j> and its descendents are now fully explored
          if (--level == 0) {
            break;  // we've explored the whole connected component
          }
          stack -= 2;
          ptr[stack[0]] = 0;
          stack[1]++;
        } else {
          ptr[j] = 1;
          next   = INT_INTOBJ(CONST_ADDR_OBJ(nbs)[k]);
          level++;
          stack += 2;
          stack[0] = next;
          stack[1] = 1;
          // if <next> is a brand new vertex, add it to the spanning tree
          if (ptr[next] == 0) {
            // create the edge <j> -> <next>
            out_j = ELM_PLIST(out, j);
            len   = LEN_PLIST(out_j);
            ASS_LIST(out_j, len + 1, INTOBJ_INT(next));
            // create the edge <next> -> <j>
            // since <next> is new, <out[next]> is still a T_PLIST_EMPTY
            new = ELM_PLIST(out, next);
            ASS_LIST(new, 1, INTOBJ_INT(j));
          }
        }
      }
    }
  }
  free(ptr);
  free(stack);
  return out;
}

static Obj FuncDIGRAPH_SOURCE_RANGE(Obj self, Obj D) {
  Obj src, ran, adj, adji;
  Int i, j, k, m, n, len;

  m   = DigraphNrEdges(D);
  n   = DigraphNrVertices(D);
  adj = FuncOutNeighbours(self, D);

  if (m == 0) {
    src = NEW_PLIST_IMM(T_PLIST_EMPTY, m);
    ran = NEW_PLIST_IMM(T_PLIST_EMPTY, m);
  } else {
    src = NEW_PLIST_IMM(T_PLIST_CYC, m);
    ran = NEW_PLIST_IMM(T_PLIST_CYC, m);
    k   = 0;
    for (i = 1; i <= n; i++) {
      adji = ELM_PLIST(adj, i);
      len  = LEN_LIST(adji);
      for (j = 1; j <= len; j++) {
        k++;
        SET_ELM_PLIST(src, k, INTOBJ_INT(i));
        SET_ELM_PLIST(ran, k, ELM_LIST(adji, j));
      }
    }
  }

  SET_LEN_PLIST(src, m);
  SET_LEN_PLIST(ran, m);
  if (IsAttributeStoringRep(D)) {
    AssPRec(D, RNamName("DigraphSource"), src);
    AssPRec(D, RNamName("DigraphRange"), ran);
    return D;
  } else {
    Obj tmp = NEW_PREC(2);
    SET_LEN_PREC(tmp, 2);
    AssPRec(tmp, RNamName("DigraphSource"), src);
    AssPRec(tmp, RNamName("DigraphRange"), ran);
    return tmp;
  }
}

// Assume we are passed a GAP Int nrvertices
// Two GAP lists of PosInts (all <= nrvertices) of equal length

// Function to change Out-Neighbours to In-Neighbours, and vice versa
static Obj FuncDIGRAPH_IN_OUT_NBS(Obj self, Obj adj) {
  Obj  inn, innk, adji;
  UInt n, i, j, k, len, len2;

  n = LEN_PLIST(adj);
  if (n == 0) {
    return NEW_PLIST_IMM(T_PLIST_EMPTY, 0);
  }
  inn = NEW_PLIST(T_PLIST_TAB, n);
  SET_LEN_PLIST(inn, n);

  // fill adj with empty plists
  for (i = 1; i <= n; i++) {
    SET_ELM_PLIST(inn, i, NEW_PLIST(T_PLIST_EMPTY, 0));
    CHANGED_BAG(inn);
  }

  for (i = 1; i <= n; i++) {
    adji = ELM_PLIST(adj, i);
    PLAIN_LIST(adji);
    len = LEN_PLIST(adji);
    for (j = 1; j <= len; j++) {
      k    = INT_INTOBJ(ELM_PLIST(adji, j));
      innk = ELM_PLIST(inn, k);
      len2 = LEN_PLIST(innk);
      ASS_LIST(innk, len2 + 1, INTOBJ_INT(i));
    }
  }
  return inn;
}

Obj FuncADJACENCY_MATRIX(Obj self, Obj digraph) {
  Int n, i, j, val, len, outj;
  Obj adj, mat, adji, next;

  n = DigraphNrVertices(digraph);
  if (n == 0) {
    return NEW_PLIST_IMM(T_PLIST_EMPTY, 0);
  }

  adj = FuncOutNeighbours(self, digraph);
  mat = NEW_PLIST_IMM(T_PLIST_TAB, n);
  SET_LEN_PLIST(mat, n);

  for (i = 1; i <= n; i++) {
    // Set up the i^th row of the adjacency matrix
    next = NEW_PLIST_IMM(T_PLIST_CYC, n);
    SET_LEN_PLIST(next, n);
    for (j = 1; j <= n; j++) {
      SET_ELM_PLIST(next, j, INTOBJ_INT(0));
    }
    // Fill in the correct values of the matrix
    adji = ELM_PLIST(adj, i);
    len  = LEN_LIST(adji);
    for (j = 1; j <= len; j++) {
      outj = INT_INTOBJ(ELM_LIST(adji, j));
      val  = INT_INTOBJ(ELM_PLIST(next, outj)) + 1;
      SET_ELM_PLIST(next, outj, INTOBJ_INT(val));
    }
    SET_ELM_PLIST(mat, i, next);
    CHANGED_BAG(mat);
  }
  SET_FILT_LIST(mat, FN_IS_RECT);
  return mat;
}

static Obj FuncIS_MULTI_DIGRAPH(Obj self, Obj digraph) {
  Obj  adj, adji;
  UInt n, i, k, j, *seen;

  adj  = FuncOutNeighbours(self, digraph);
  n    = DigraphNrVertices(digraph);
  seen = safe_calloc(n + 1, sizeof(UInt));

  for (i = 1; i <= n; i++) {
    adji = ELM_PLIST(adj, i);
    if ((UInt) LEN_LIST(adji) > n) {
      free(seen);
      return True;
    }
    PLAIN_LIST(adji);
    for (j = 1; j <= (UInt) LEN_PLIST(adji); j++) {
      k = INT_INTOBJ(ELM_PLIST(adji, j));
      if (seen[k] != i) {
        seen[k] = i;
      } else {
        free(seen);
        return True;
      }
    }
  }
  free(seen);
  return False;
}

/***************** GENERAL FLOYD_WARSHALL ALGORITHM***********************
* This function accepts 5 arguments:
*   1. A digraph.
*   2. A special function which takes 5 arguments:
*       - The matrix dist
*       - 3 integers i, j, k
*       - An integer n (the number of vertices of digraph)
*      and modifies the matrix dist according to the values of i, j, k.
*   3. Int val1 s.t initially dist[i][j] = val1 if [ i, j ] isn't an edge.
*   4. Int val2 s.t initially dist[i][j] = val2 if [ i, j ] is an edge.
*   5. bool copy:
*      - If true, FLOYD_WARSHALL stores the initialised dist, and
*        compares it with dist after it has gone through the 3 for-loops,
*        and returns true iff it is unchanged.
*      - If false, proceeds as usual Floyd-Warshall algorithm and returns
*        a GAP object matrix as the result.
*   6. bool diameter: // TODO wouldn't it be better to just take a
*                        "post-processing" function. JDM
*      - If true, FLOYD_WARSHALL goes through dist after the 3 for-loops,
*        returns -1 if dist contains the value -1, else it returns the
*        maximum value of dist
*      - If false, proceeds as usual
*   7. bool shortest:
*      - If true, for each vertex i, dist[i][i] is initially set to 0
*      - If false, this step is skipped
*/
static Obj FLOYD_WARSHALL(Obj digraph,
                          void (*func)(Int** dist, Int i, Int j, Int k, Int n),
                          Int  val1,
                          Int  val2,
                          bool copy,
                          bool diameter,
                          bool shortest) {
  Int  n, i, j, k, *dist;
  Int* adj = NULL;
  Obj  next, out, outi, val;

  n = DigraphNrVertices(digraph);

  // Special case for 0-vertex graph
  if (n == 0) {
    if (diameter) {
      return Fail;
    }
    if (copy) {
      return True;
    }
    return NEW_PLIST_IMM(T_PLIST_EMPTY, 0);
  }

  // Initialise the n x n matrix with val1 and val2
  dist = safe_malloc(n * n * sizeof(Int));
  for (i = 0; i < n * n; i++) {
    dist[i] = val1;
  }
  out = FuncOutNeighbours(0L, digraph);
  for (i = 1; i <= n; i++) {
    outi = ELM_PLIST(out, i);
    PLAIN_LIST(outi);
    for (j = 1; j <= LEN_PLIST(outi); j++) {
      k       = (i - 1) * n + INT_INTOBJ(ELM_PLIST(outi, j)) - 1;
      dist[k] = val2;
    }
  }

  if (shortest) {
    // This is the special case for DIGRAPH_SHORTEST_DIST
    for (i = 0; i < n; i++) {
      dist[i * n + i] = 0;
    }
  }

  if (copy) {
    // This is the special case for IS_TRANSITIVE_DIGRAPH
    adj = safe_malloc(n * n * sizeof(Int));
    for (i = 0; i < n * n; i++) {
      adj[i] = dist[i];
    }
  }

  for (k = 0; k < n; k++) {
    for (i = 0; i < n; i++) {
      for (j = 0; j < n; j++) {
        func(&dist, i, j, k, n);
      }
    }
  }

  // the following is a horrible hack
  if (diameter) {
    // This is the special case for DIGRAPH_DIAMETER
    Int maximum = -1;
    for (i = 0; i < n; i++) {
      for (j = 0; j < n; j++) {
        if (dist[i * n + j] > maximum) {
          maximum = dist[i * n + j];
        } else if (dist[i * n + j] == -1) {
          free(dist);
          free(adj);
          return Fail;
        }
      }
    }
    free(dist);
    return INTOBJ_INT(maximum);
  }

  if (copy) {
    // This is the special case for IS_TRANSITIVE_DIGRAPH
    for (i = 0; i < n * n; i++) {
      if (adj[i] != dist[i]) {
        free(dist);
        free(adj);
        return False;
      }
    }
    free(dist);
    free(adj);
    return True;
  }

  // Create GAP matrix to return
  out = NEW_PLIST(T_PLIST_TAB, n);
  SET_LEN_PLIST(out, n);

  for (i = 1; i <= n; i++) {
    next = NEW_PLIST(T_PLIST_CYC, n);
    SET_LEN_PLIST(next, n);
    for (j = 1; j <= n; j++) {
      val = INTOBJ_INT(dist[(i - 1) * n + (j - 1)]);
      if (val == INTOBJ_INT(-1)) {
        val = Fail;
      }
      SET_ELM_PLIST(next, j, val);
    }
    SET_ELM_PLIST(out, i, next);
    CHANGED_BAG(out);
  }
  SET_FILT_LIST(out, FN_IS_RECT);

  free(dist);
  return out;
}

static void FW_FUNC_SHORTEST_DIST(Int** dist, Int i, Int j, Int k, Int n) {
  if ((*dist)[i * n + k] != -1 && (*dist)[k * n + j] != -1) {
    if ((*dist)[i * n + j] == -1
        || (*dist)[i * n + j] > (*dist)[i * n + k] + (*dist)[k * n + j]) {
      (*dist)[i * n + j] = (*dist)[i * n + k] + (*dist)[k * n + j];
    }
  }
}

static Obj FuncDIGRAPH_SHORTEST_DIST(Obj self, Obj digraph) {
  return FLOYD_WARSHALL(
      digraph, FW_FUNC_SHORTEST_DIST, -1, 1, false, false, true);
}

static Obj FuncDIGRAPH_DIAMETER(Obj self, Obj digraph) {
  return FLOYD_WARSHALL(
      digraph, FW_FUNC_SHORTEST_DIST, -1, 1, false, true, true);
}

static void FW_FUNC_TRANS_CLOSURE(Int** dist, Int i, Int j, Int k, Int n) {
  if ((*dist)[i * n + k] != 0 && (*dist)[k * n + j] != 0) {
    (*dist)[i * n + j] = 1;
  }
}

static Obj FuncIS_TRANSITIVE_DIGRAPH(Obj self, Obj digraph) {
  return FLOYD_WARSHALL(
      digraph, FW_FUNC_TRANS_CLOSURE, 0, 1, true, false, false);
}

static Obj FuncDIGRAPH_TRANS_CLOSURE(Obj self, Obj digraph) {
  return FLOYD_WARSHALL(
      digraph, FW_FUNC_TRANS_CLOSURE, 0, 1, false, false, false);
}

static void
FW_FUNC_REFLEX_TRANS_CLOSURE(Int** dist, Int i, Int j, Int k, Int n) {
  if ((i == j) || ((*dist)[i * n + k] != 0 && (*dist)[k * n + j] != 0)) {
    (*dist)[i * n + j] = 1;
  }
}

static Obj FuncDIGRAPH_REFLEX_TRANS_CLOSURE(Obj self, Obj digraph) {
  return FLOYD_WARSHALL(
      digraph, FW_FUNC_REFLEX_TRANS_CLOSURE, 0, 1, false, false, false);
}

static Obj FuncRANDOM_DIGRAPH(Obj self, Obj nn, Obj pp) {
  UInt    n, i, j;
  Double  p, q;
  Int     len;
  Obj     adj, adji;

  n   = INT_INTOBJ(nn);
  p   = VAL_MACFLOAT(pp);
  adj = NEW_PLIST(T_PLIST_TAB, n);
  SET_LEN_PLIST(adj, n);

  for (i = 1; i <= n; i++) {
    SET_ELM_PLIST(adj, i, NEW_PLIST(T_PLIST_EMPTY, 0));
    CHANGED_BAG(adj);
  }

  for (i = 1; i <= n; i++) {
    for (j = 1; j <= n; j++) {
      q = (double) rand() / RAND_MAX;
      if (q < p) {
        adji = ELM_PLIST(adj, i);
        len  = LEN_PLIST(adji);
        ASS_LIST(adji, len + 1, INTOBJ_INT(j));
      }
    }
  }
  return adj;
}

static Obj FuncRANDOM_MULTI_DIGRAPH(Obj self, Obj nn, Obj mm) {
  UInt n, m, i, j, k;
  Int  len;
  Obj  adj, adjj;

  n   = INT_INTOBJ(nn);
  m   = INT_INTOBJ(mm);
  adj = NEW_PLIST(T_PLIST_TAB, n);
  SET_LEN_PLIST(adj, n);

  for (i = 1; i <= n; i++) {
    SET_ELM_PLIST(adj, i, NEW_PLIST(T_PLIST_EMPTY, 0));
    CHANGED_BAG(adj);
  }

  for (i = 1; i <= m; i++) {
    j    = (rand() % n) + 1;
    k    = (rand() % n) + 1;
    adjj = ELM_PLIST(adj, j);
    len  = LEN_PLIST(adjj);
    ASS_LIST(adjj, len + 1, INTOBJ_INT(k));
  }
  return adj;
}

static bool EqJumbledPlists(Obj l, Obj r, Int nr, Int* buf) {
  bool eq;
  Int  j, jj;

  // Check first whether the lists are identical
  eq = true;
  for (j = 1; j <= nr; j++) {
    jj = INT_INTOBJ(ELM_PLIST(l, j));
    if (jj != INT_INTOBJ(ELM_PLIST(r, j))) {
      eq = false;
      break;
    }
  }

  // Otherwise check that they have equal content
  if (!eq) {
    for (j = 1; j <= nr; j++) {
      jj = INT_INTOBJ(ELM_PLIST(l, j)) - 1;
      buf[jj]++;
      jj = INT_INTOBJ(ELM_PLIST(r, j)) - 1;
      buf[jj]--;
    }

    for (j = 1; j <= nr; j++) {
      jj = INT_INTOBJ(ELM_PLIST(l, j)) - 1;
      if (buf[jj] != 0) {
        return false;
      }
    }
  }
  return true;
}

//-----------------------------------------------------------------------------
// MurmurHash3 was written by Austin Appleby, and is placed in the public
// domain. The author hereby disclaims copyright to this source code.

/* Minor modifications to get it to compile in C rather than C++ and
integrate with GAP  SL*/

/* Digraphs takes parts of this source from GAP. */

#define BIG_CONSTANT(x) (x##LLU)

static inline UInt fmix(UInt h) {
#ifndef SYS_IS_64_BIT
  h ^= h >> 16;
  h *= 0x85ebca6b;
  h ^= h >> 13;
  h *= 0xc2b2ae35;
  h ^= h >> 16;
#else
  h ^= h >> 33;
  h *= BIG_CONSTANT(0xff51afd7ed558ccd);
  h ^= h >> 33;
  h *= BIG_CONSTANT(0xc4ceb9fe1a85ec53);
  h ^= h >> 33;
#endif

  return h;
}

static Obj FuncDIGRAPH_HASH(Obj self, Obj digraph) {
  UInt n, i, h;
  Obj  out, a;
  Int  nr, j;

  h   = 0;
  n   = DigraphNrVertices(digraph);
  out = FuncOutNeighbours(self, digraph);

  for (i = 1; i <= n; i++) {
    a = ELM_PLIST(out, i);
    PLAIN_LIST(a);
    nr = LEN_PLIST(a);
    for (j = 1; j <= nr; j++)
      h += fmix(INT_INTOBJ(ELM_PLIST(a, j)));
    h = fmix(h + 0x9e3779b9);
  }

  return INTOBJ_INT(h);
}

static Obj FuncDIGRAPH_EQUALS(Obj self, Obj digraph1, Obj digraph2) {
  UInt i, n1, n2, m1, m2;
  Obj  out1, out2, a, b;
  Int  nr, *buf;

  // Check NrVertices is equal
  n1 = DigraphNrVertices(digraph1);
  n2 = DigraphNrVertices(digraph2);
  if (n1 != n2) {
    return False;
  }

  // Check NrEdges is equal
  m1 = DigraphNrEdges(digraph1);
  m2 = DigraphNrEdges(digraph2);

  if (m1 != m2) {
    return False;
  }

  out1 = FuncOutNeighbours(self, digraph1);
  out2 = FuncOutNeighbours(self, digraph2);

  buf = safe_calloc(n1, sizeof(Int));

  // Compare OutNeighbours of each vertex in turn
  for (i = 1; i <= n1; i++) {
    a = ELM_PLIST(out1, i);
    b = ELM_PLIST(out2, i);
    PLAIN_LIST(a);
    PLAIN_LIST(b);

    nr = LEN_PLIST(a);
    // Check that the OutDegrees match
    if (nr != LEN_PLIST(b)) {
      free(buf);
      return False;
    }

    if (!EqJumbledPlists(a, b, nr, buf)) {
      free(buf);
      return False;
    }
  }
  free(buf);
  return True;
}

static Int LTJumbledPlists(Obj l, Obj r, Int nr1, Int nr2, Int* buf, Int n) {
  bool eq;
  Int  j, jj, min;

  // Check first whether the lists are identical
  if (nr1 == nr2) {
    eq = true;
    for (j = 1; j <= nr1; j++) {
      jj = INT_INTOBJ(ELM_PLIST(l, j));
      if (jj != INT_INTOBJ(ELM_PLIST(r, j))) {
        eq = false;
        break;
      }
    }
  } else {
    eq = false;
  }

  // Otherwise compare their content
  if (!eq) {
    min = nr1 < nr2 ? nr1 : nr2;

    for (j = 1; j <= min; j++) {
      jj = INT_INTOBJ(ELM_PLIST(l, j)) - 1;
      buf[jj]++;
      jj = INT_INTOBJ(ELM_PLIST(r, j)) - 1;
      buf[jj]--;
    }

    for (j = min + 1; j <= nr1; j++) {
      jj = INT_INTOBJ(ELM_PLIST(l, j)) - 1;
      buf[jj]++;
    }

    for (j = min + 1; j <= nr2; j++) {
      jj = INT_INTOBJ(ELM_PLIST(r, j)) - 1;
      buf[jj]--;
    }

    for (j = 0; j < n; j++) {
      if (buf[j] < 0) {
        // Pr("Found difference: range: %d, number: %d\n", j + 1, buf[j]);
        return 2;
      } else if (buf[j] > 0) {
        // Pr("Found difference: range: %d, number: %d\n", j + 1, buf[j]);
        return 1;
      }
    }
  }
  return 0;
  // Return 0: l = r (as multisets)
  // Return 1: l < r
  // Return 2: r < l
}

static Obj FuncDIGRAPH_LT(Obj self, Obj digraph1, Obj digraph2) {
  UInt i, n1, n2, m1, m2;
  Obj  out1, out2, a, b;
  Int  nr1, nr2, *buf, comp, max;

  // Compare NrVertices
  n1 = DigraphNrVertices(digraph1);
  n2 = DigraphNrVertices(digraph2);
  if (n1 < n2) {
    return True;
  } else if (n2 < n1) {
    return False;
  }

  // Compare NrEdges
  m1 = DigraphNrEdges(digraph1);
  m2 = DigraphNrEdges(digraph2);

  if (m1 < m2) {
    return True;
  } else if (m2 < m1) {
    return False;
  }

  out1 = FuncOutNeighbours(self, digraph1);
  out2 = FuncOutNeighbours(self, digraph2);

  buf = safe_calloc(n1, sizeof(Int));

  // Compare Sorted(out1[i]) and Sorted(out2[i]) for each vertex i
  for (i = 1; i <= n1; i++) {
    a = ELM_PLIST(out1, i);
    b = ELM_PLIST(out2, i);
    PLAIN_LIST(a);
    PLAIN_LIST(b);

    nr1 = LEN_PLIST(a);
    nr2 = LEN_PLIST(b);
    max = nr1 < nr2 ? nr2 : nr1;

    // Check whether both vertices have 0 out-degree
    if (max != 0) {
      if (nr1 == 0) {
        free(buf);
        return False;
      } else if (nr2 == 0) {
        free(buf);
        return True;
      }
      // Both vertices have positive out-degree

      // Compare out1[i] and out2[i]
      comp = LTJumbledPlists(a, b, nr1, nr2, buf, n1);
      if (comp == 1) {
        free(buf);
        return True;
      } else if (comp == 2) {
        free(buf);
        return False;
      }
      // if comp == 0 then the lists are equal, so continue
    }
  }
  free(buf);
  return False;
}

// bliss

static BlissGraph* buildBlissMultiDigraph(Obj digraph) {
  UInt        n, i, j, k, l, nr;
  Obj         adji, adj;
  BlissGraph* graph;

  n     = DigraphNrVertices(digraph);
  graph = bliss_digraphs_new(n);

  adj = FuncOutNeighbours(0L, digraph);
  for (i = 1; i <= n; i++) {
    adji = ELM_PLIST(adj, i);
    nr   = LEN_PLIST(adji);
    for (j = 1; j <= nr; j++) {
      k = bliss_digraphs_add_vertex(graph, 1);
      l = bliss_digraphs_add_vertex(graph, 2);
      bliss_digraphs_add_edge(graph, i - 1, k);
      bliss_digraphs_add_edge(graph, k, l);
      bliss_digraphs_add_edge(graph, l, INT_INTOBJ(ELM_PLIST(adji, j)) - 1);
    }
  }
  return graph;
}

// TODO: document mult (and everything else)
static BlissGraph*
buildBlissDigraph(Obj digraph, Obj vert_colours, Obj edge_colours) {
  uint64_t    colour, mult, num_vc, num_ec, n, i, j, k, nr;
  Obj         adjj, adj;
  BlissGraph* graph;

  n      = DigraphNrVertices(digraph);
  num_vc = 0;
  num_ec = 0;

  // TODO: make a decision about this
  // mult = (orientation_double == True) ? 2 : 1;

  mult = 2;

  if (vert_colours != Fail) {
    DIGRAPHS_ASSERT(n == (uint64_t) LEN_LIST(vert_colours));
    for (i = 1; i <= n; i++) {
      num_vc = MAX(num_vc, (uint64_t) INT_INTOBJ(ELM_LIST(vert_colours, i)));
    }
  }

  adj = FuncOutNeighbours(0L, digraph);

  if (edge_colours != Fail) {
    DIGRAPHS_ASSERT(n == (uint64_t) LEN_LIST(edge_colours));
    for (i = 1; i <= n; i++) {
      Int len = LEN_LIST(ELM_PLIST(edge_colours, i));
      DIGRAPHS_ASSERT(LEN_LIST(ELM_PLIST(adj, i)) == len);
      for (Int l = 1; l <= len; l++) {
        uint64_t x = INT_INTOBJ(ELM_LIST(ELM_LIST(edge_colours, i), l));
        num_ec     = MAX(num_ec, x);
      }
    }
  } else if (DigraphNrEdges(digraph) > 0) {
    num_ec = 1;
  }

  graph = bliss_digraphs_new(0);

  // TODO: make this safe
  uint64_t num_layers = 64 - __builtin_clzll(num_ec);

  // Take care of the case where there are no edges in the digraph
  if (DigraphNrEdges(digraph) == 0) {
    num_layers = 1;
    mult       = 1;
  }

  if (vert_colours == Fail) {
    num_vc = 1;
  }

  // TODO: is duplicating the best idea here?
  for (i = 1; i <= mult * num_layers; i += mult) {
    for (j = 1; j <= n; j++) {
      colour = (vert_colours != Fail)
                   ? (i - 1) * num_vc + INT_INTOBJ(ELM_LIST(vert_colours, j))
                   : i - 1;
      bliss_digraphs_add_vertex(graph, colour);
    }
    if (mult == 2) {
      for (j = 1; j <= n; j++) {
        colour = (vert_colours != Fail)
                     ? i * num_vc + INT_INTOBJ(ELM_LIST(vert_colours, j))
                     : i;
        bliss_digraphs_add_vertex(graph, colour);
      }
    }
  }

  if (mult == 2) {
    for (i = 0; i < n; i++) {
      j = bliss_digraphs_add_vertex(graph, num_vc * num_layers * mult + 2);
      bliss_digraphs_add_edge(graph, j, i);
      bliss_digraphs_add_edge(graph, j, i + n);
      for (k = 0; k < num_layers; k++) {
        bliss_digraphs_add_edge(graph, j, i + k * mult * n);
        bliss_digraphs_add_edge(graph, j, i + (k * mult + 1) * n);
      }
    }
  }

  for (i = 1; i < num_layers; i++) {
    for (j = 1; j <= mult * n; j++) {
      bliss_digraphs_add_edge(
          graph, (i - 1) * mult * n + (j - 1), i * mult * n + (j - 1));
    }
  }

  for (j = 1; j <= n; j++) {
    adjj = ELM_PLIST(adj, j);
    nr   = LEN_PLIST(adjj);
    for (k = 1; k <= nr; k++) {
      uint64_t w = INT_INTOBJ(ELM_PLIST(adjj, k));
      for (i = 0; i < num_layers; i++) {
        colour = edge_colours != Fail
                     ? INT_INTOBJ(ELM_LIST(ELM_LIST(edge_colours, j), k))
                     : 1;
        if ((1 << i) & colour) {
          bliss_digraphs_add_edge(graph,
                                  i * mult * n + (j - 1),
                                  ((i + 1) * mult - 1) * n + (w - 1));
        }
      }
    }
  }
  return graph;
}

static BlissGraph* buildBlissMultiDigraphWithColours(Obj digraph, Obj colours) {
  UInt        n, i, j, k, l, nr;
  Obj         adji, adj;
  BlissGraph* graph;

  n = DigraphNrVertices(digraph);
  DIGRAPHS_ASSERT(n == (UInt) LEN_LIST(colours));
  graph = bliss_digraphs_new(0);
  adj   = FuncOutNeighbours(0L, digraph);

  for (i = 1; i <= n; i++) {
    bliss_digraphs_add_vertex(graph, INT_INTOBJ(ELM_LIST(colours, i)));
  }
  for (i = 1; i <= n; i++) {
    bliss_digraphs_add_vertex(graph, n + 1);
  }
  for (i = 1; i <= n; i++) {
    bliss_digraphs_add_vertex(graph, n + 2);
  }

  for (i = 1; i <= n; i++) {
    bliss_digraphs_add_edge(graph, i - 1, n + i - 1);
    bliss_digraphs_add_edge(graph, i - 1, 2 * n + i - 1);
    adji = ELM_PLIST(adj, i);
    nr   = LEN_PLIST(adji);
    for (j = 1; j <= nr; j++) {
      k = bliss_digraphs_add_vertex(graph, n + 3);
      l = bliss_digraphs_add_vertex(graph, n + 4);
      bliss_digraphs_add_edge(graph, n + i - 1, k);
      bliss_digraphs_add_edge(graph, k, l);
      bliss_digraphs_add_edge(
          graph, l, 2 * n + INT_INTOBJ(ELM_PLIST(adji, j)) - 1);
    }
  }

  return graph;
}

static void digraph_hook_function(void*               user_param,
                                  unsigned int        N,
                                  const unsigned int* aut) {
  UInt4* ptr;
  Obj    p, gens;
  UInt   i, n;

  n = INT_INTOBJ(ELM_PLIST(user_param, 2));  // the degree
  DIGRAPHS_ASSERT(n <= N);
  p   = NEW_PERM4(n);
  ptr = ADDR_PERM4(p);

  for (i = 0; i < n; i++) {
    ptr[i] = aut[i];
  }

  gens = ELM_PLIST(user_param, 1);
  AssPlist(gens, LEN_PLIST(gens) + 1, p);
}

// Take a list of C integers, and multiply them together into a GAP int
static Obj MultiplyList(int* vals, int length) {
  Obj res = INTOBJ_INT(1);
  for (int i = 0; i < length; ++i) {
    res = ProdInt(res, INTOBJ_INT(vals[i]));
  }
  return res;
}

static Obj FuncDIGRAPH_AUTOMORPHISMS(Obj self,
                                     Obj digraph,
                                     Obj vert_colours,
                                     Obj edge_colours) {
  Obj                 autos, p, n;
  BlissGraph*         graph;
  UInt4*              ptr;
  const unsigned int* canon;
  Int                 i;

  graph = buildBlissDigraph(digraph, vert_colours, edge_colours);

  autos = NEW_PLIST(T_PLIST, 3);
  n     = INTOBJ_INT(DigraphNrVertices(digraph));

  SET_ELM_PLIST(autos, 1, NEW_PLIST(T_PLIST, 0));  // perms of the vertices
  CHANGED_BAG(autos);
  SET_ELM_PLIST(autos, 2, n);
  SET_LEN_PLIST(autos, 2);

  BlissStats stats;

  canon = bliss_digraphs_find_canonical_labeling(
      graph, digraph_hook_function, autos, &stats);

  p   = NEW_PERM4(INT_INTOBJ(n));
  ptr = ADDR_PERM4(p);

  for (i = 0; i < INT_INTOBJ(n); i++) {
    ptr[i] = canon[i];
  }
  SET_ELM_PLIST(autos, 2, p);
  CHANGED_BAG(autos);

  bliss_digraphs_release(graph);
  if (LEN_PLIST(ELM_PLIST(autos, 1)) != 0) {
    SortDensePlist(ELM_PLIST(autos, 1));
    RemoveDupsDensePlist(ELM_PLIST(autos, 1));
  }

#ifdef DIGRAPHS_WITH_INCLUDED_BLISS
  Obj size = MultiplyList(stats.group_size, stats.group_size_len);
  bliss_digraphs_free_blissstats(&stats);

  SET_LEN_PLIST(autos, 3);
  SET_ELM_PLIST(autos, 3, size);
#endif

  return autos;
}

// user_param = [vertex perms, nr vertices, edge perms, nr edges]
static void multidigraph_hook_function(void*               user_param,
                                       unsigned int        N,
                                       const unsigned int* aut) {
  UInt4* ptr;
  Obj    p, gens;
  UInt   i, n, m;
  bool   stab;

  m = INT_INTOBJ(ELM_PLIST(user_param, 2));  // the nr of vertices
  DIGRAPHS_ASSERT(m <= N);

  stab = true;
  for (i = 0; i < m; i++) {
    if (aut[i] != i) {
      stab = false;
    }
  }
  if (stab) {                                    // permutation of the edges
    n   = INT_INTOBJ(ELM_PLIST(user_param, 4));  // the nr of edges
    p   = NEW_PERM4(n);
    ptr = ADDR_PERM4(p);
    DIGRAPHS_ASSERT(2 * (n - 1) + m <= N);
    for (i = 0; i < n; i++) {
      ptr[i] = (aut[2 * i + m] - m) / 2;
    }
    gens = ELM_PLIST(user_param, 3);
  } else {  // permutation of the vertices
    p   = NEW_PERM4(m);
    ptr = ADDR_PERM4(p);
    DIGRAPHS_ASSERT(m <= N);
    for (i = 0; i < m; i++) {
      ptr[i] = aut[i];
    }
    gens = ELM_PLIST(user_param, 1);
  }

  AssPlist(gens, LEN_PLIST(gens) + 1, p);
}

// user_param = [vertex perms, nr vertices, edge perms, nr edges]
static void multidigraph_colours_hook_function(void*               user_param,
                                               unsigned int        N,
                                               const unsigned int* aut) {
  UInt4* ptr;
  Obj    p, gens;
  UInt   i, n, m;
  bool   stab;

  m = INT_INTOBJ(ELM_PLIST(user_param, 2));  // the nr of vertices
  DIGRAPHS_ASSERT(m <= N);

  stab = true;
  for (i = 0; i < m; i++) {
    if (aut[i] != i) {
      stab = false;
    }
  }
  if (stab) {                                    // permutation of the edges
    n   = INT_INTOBJ(ELM_PLIST(user_param, 4));  // the nr of edges
    p   = NEW_PERM4(n);
    ptr = ADDR_PERM4(p);
    DIGRAPHS_ASSERT(2 * (n - 1) + 3 * m <= N);
    for (i = 0; i < n; i++) {
      ptr[i] = (aut[2 * i + 3 * m] - 3 * m) / 2;
    }
    gens = ELM_PLIST(user_param, 3);
  } else {  // permutation of the vertices
    p   = NEW_PERM4(m);
    ptr = ADDR_PERM4(p);
    DIGRAPHS_ASSERT(m < N);
    for (i = 0; i < m; i++) {
      ptr[i] = aut[i];
    }
    gens = ELM_PLIST(user_param, 1);
  }

  AssPlist(gens, LEN_PLIST(gens) + 1, p);
}

static Obj FuncMULTIDIGRAPH_AUTOMORPHISMS(Obj self, Obj digraph, Obj colours) {
  Obj                 autos, p, q, out;
  BlissGraph*         graph;
  UInt4*              ptr;
  const unsigned int* canon;
  Int                 i, m, n;

  if (colours == False) {
    graph = buildBlissMultiDigraph(digraph);
  } else {
    graph = buildBlissMultiDigraphWithColours(digraph, colours);
  }
  autos = NEW_PLIST(T_PLIST, 4);
  SET_ELM_PLIST(autos, 1, NEW_PLIST(T_PLIST, 0));  // perms of the vertices
  CHANGED_BAG(autos);
  SET_ELM_PLIST(autos, 2, INTOBJ_INT(DigraphNrVertices(digraph)));
  SET_ELM_PLIST(autos, 3, NEW_PLIST(T_PLIST, 0));  // perms of the edges
  CHANGED_BAG(autos);
  SET_ELM_PLIST(autos, 4, INTOBJ_INT(DigraphNrEdges(digraph)));

  BlissStats stats;

  if (colours == False) {
    canon = bliss_digraphs_find_canonical_labeling(
        graph, multidigraph_hook_function, autos, &stats);
  } else {
    canon = bliss_digraphs_find_canonical_labeling(
        graph, multidigraph_colours_hook_function, autos, &stats);
  }

  // Get canonical labeling as GAP perms
  m   = DigraphNrVertices(digraph);
  p   = NEW_PERM4(m);  // perm of vertices
  ptr = ADDR_PERM4(p);

  for (i = 0; i < m; i++) {
    ptr[i] = canon[i];
  }

  n   = DigraphNrEdges(digraph);
  q   = NEW_PERM4(n);  // perm of edges
  ptr = ADDR_PERM4(q);

  if (colours == False) {
    for (i = 0; i < n; i++) {
      ptr[i] = canon[2 * i + m] - m;
    }
  } else {
    for (i = 0; i < n; i++) {
      ptr[i] = canon[2 * i + 3 * m] - 3 * m;
    }
  }

  bliss_digraphs_release(graph);

  // put the canonical labeling (as a list of two perms) into autos[2]
  out = NEW_PLIST(T_PLIST, 2);
  SET_ELM_PLIST(out, 1, p);
  SET_ELM_PLIST(out, 2, q);
  SET_LEN_PLIST(out, 2);
  CHANGED_BAG(out);

  SET_ELM_PLIST(autos, 2, out);
  CHANGED_BAG(autos);

  // remove 4th entry of autos (the number of edges) . . .
  SET_LEN_PLIST(autos, 3);

  if (LEN_PLIST(ELM_PLIST(autos, 1)) != 0) {
    SortDensePlist(ELM_PLIST(autos, 1));
    RemoveDupsDensePlist(ELM_PLIST(autos, 1));
  }
  if (LEN_PLIST(ELM_PLIST(autos, 3)) != 0) {
    SortDensePlist(ELM_PLIST(autos, 3));
    RemoveDupsDensePlist(ELM_PLIST(autos, 3));
  }

#ifdef DIGRAPHS_WITH_INCLUDED_BLISS
  Obj size = MultiplyList(stats.group_size, stats.group_size_len);
  bliss_digraphs_free_blissstats(&stats);

  SET_LEN_PLIST(autos, 4);
  SET_ELM_PLIST(autos, 4, size);
#endif

  return autos;
}

static Obj FuncDIGRAPH_CANONICAL_LABELLING(Obj self, Obj digraph, Obj colours) {
  Obj                 p;
  UInt4*              ptr;
  BlissGraph*         graph;
  Int                 n, i;
  const unsigned int* canon;

  if (colours == Fail) {
    graph = buildBlissDigraph(digraph, NULL, NULL);
  } else {
    graph = buildBlissDigraph(digraph, colours, NULL);
  }

  canon = bliss_digraphs_find_canonical_labeling(graph, 0, 0, 0);

  n   = DigraphNrVertices(digraph);
  p   = NEW_PERM4(n);
  ptr = ADDR_PERM4(p);

  for (i = 0; i < n; i++) {
    ptr[i] = canon[i];
  }
  bliss_digraphs_release(graph);

  return p;
}

static Obj
FuncMULTIDIGRAPH_CANONICAL_LABELLING(Obj self, Obj digraph, Obj colours) {
  Obj                 p, q, out;
  UInt4*              ptr;
  BlissGraph*         graph;
  Int                 m, n, i;
  const unsigned int* canon;

  if (colours == Fail) {
    graph = buildBlissMultiDigraph(digraph);
  } else {
    graph = buildBlissMultiDigraphWithColours(digraph, colours);
  }

  canon = bliss_digraphs_find_canonical_labeling(graph, 0, 0, 0);

  m   = DigraphNrVertices(digraph);
  p   = NEW_PERM4(m);  // perm of vertices
  ptr = ADDR_PERM4(p);

  for (i = 0; i < m; i++) {
    ptr[i] = canon[i];
  }

  n   = DigraphNrEdges(digraph);
  q   = NEW_PERM4(n);  // perm of edges
  ptr = ADDR_PERM4(q);

  if (colours == Fail) {
    for (i = 0; i < n; i++) {
      ptr[i] = canon[2 * i + m] - m;
    }
  } else {
    for (i = 0; i < n; i++) {
      ptr[i] = canon[2 * i + 3 * m] - 3 * m;
    }
  }

  bliss_digraphs_release(graph);

  out = NEW_PLIST(T_PLIST, 2);
  SET_ELM_PLIST(out, 1, p);
  SET_ELM_PLIST(out, 2, q);
  SET_LEN_PLIST(out, 2);
  CHANGED_BAG(out);

  return out;
}

/*F * * * * * * * * * * * * * initialize package * * * * * * * * * * * * * * */

/******************************************************************************
*V  GVarFuncs . . . . . . . . . . . . . . . . . . list of functions to export
*/

static StructGVarFunc GVarFuncs[] = {
    GVAR_FUNC(DIGRAPH_NREDGES, 1, "digraph"),
    GVAR_FUNC(DIGRAPH_NRADJACENCIES, 1, "digraph"),
    GVAR_FUNC(DIGRAPH_NRADJACENCIESWITHOUTLOOPS, 1, "digraph"),
    GVAR_FUNC(GABOW_SCC, 1, "adj"),
    GVAR_FUNC(DIGRAPH_CONNECTED_COMPONENTS, 1, "digraph"),
    GVAR_FUNC(IS_ACYCLIC_DIGRAPH, 1, "adj"),
    GVAR_FUNC(DIGRAPH_LONGEST_DIST_VERTEX, 2, "adj, start"),
    GVAR_FUNC(DIGRAPH_TRANS_REDUCTION, 1, "list"),
    GVAR_FUNC(IS_ANTISYMMETRIC_DIGRAPH, 1, "adj"),
    GVAR_FUNC(IS_STRONGLY_CONNECTED_DIGRAPH, 1, "adj"),
    GVAR_FUNC(DIGRAPH_TOPO_SORT, 1, "adj"),
    GVAR_FUNC(DIGRAPH_SYMMETRIC_SPANNING_FOREST, 1, "adj"),
    GVAR_FUNC(DIGRAPH_SOURCE_RANGE, 1, "digraph"),
    GVAR_FUNC(OutNeighbours, 1, "D"),
    GVAR_FUNC(DIGRAPH_OUT_NEIGHBOURS_FROM_SOURCE_RANGE, 3, "N, source, range"),
    GVAR_FUNC(DIGRAPH_NR_VERTICES, 1, "D"),
    GVAR_FUNC(DIGRAPH_IN_OUT_NBS, 1, "adj"),
    GVAR_FUNC(ADJACENCY_MATRIX, 1, "digraph"),
    GVAR_FUNC(IS_MULTI_DIGRAPH, 1, "digraph"),
    GVAR_FUNC(DIGRAPH_SHORTEST_DIST, 1, "digraph"),
    GVAR_FUNC(DIGRAPH_DIAMETER, 1, "digraph"),
    GVAR_FUNC(IS_TRANSITIVE_DIGRAPH, 1, "digraph"),
    GVAR_FUNC(DIGRAPH_TRANS_CLOSURE, 1, "digraph"),
    GVAR_FUNC(DIGRAPH_REFLEX_TRANS_CLOSURE, 1, "digraph"),
    GVAR_FUNC(RANDOM_DIGRAPH, 2, "nn, limm"),
    GVAR_FUNC(RANDOM_MULTI_DIGRAPH, 2, "nn, mm"),
    GVAR_FUNC(DIGRAPH_EQUALS, 2, "digraph1, digraph2"),
    GVAR_FUNC(DIGRAPH_LT, 2, "digraph1, digraph2"),
    GVAR_FUNC(DIGRAPH_HASH, 1, "digraph"),
    GVAR_FUNC(DIGRAPH_PATH, 3, "digraph, u, v"),
    GVAR_FUNC(DIGRAPH_AUTOMORPHISMS, 3, "digraph, vert_colours, edge_colours"),
    GVAR_FUNC(MULTIDIGRAPH_AUTOMORPHISMS, 2, "digraph, colours"),
    GVAR_FUNC(DIGRAPH_CANONICAL_LABELLING, 2, "digraph, colours"),
    GVAR_FUNC(MULTIDIGRAPH_CANONICAL_LABELLING, 2, "digraph, colours"),
    GVAR_FUNC(HomomorphismDigraphsFinder,
              -1,
              "digraph1, digraph2, hook, user_param, max_results, hint, "
              "injective, image, partial_map, colors1, colors2"),
    GVAR_FUNC(DigraphsCliquesFinder,
              -1,
              "digraph, hook, user_param, limit, "
              "include, exclude, max, size"),
    GVAR_FUNC(IS_PLANAR, 1, "digraph"),
    GVAR_FUNC(PLANAR_EMBEDDING, 1, "digraph"),
    GVAR_FUNC(KURATOWSKI_PLANAR_SUBGRAPH, 1, "digraph"),
    GVAR_FUNC(IS_OUTER_PLANAR, 1, "digraph"),
    GVAR_FUNC(OUTER_PLANAR_EMBEDDING, 1, "digraph"),
    GVAR_FUNC(KURATOWSKI_OUTER_PLANAR_SUBGRAPH, 1, "digraph"),
    GVAR_FUNC(SUBGRAPH_HOMEOMORPHIC_TO_K23, 1, "digraph"),
    GVAR_FUNC(SUBGRAPH_HOMEOMORPHIC_TO_K33, 1, "digraph"),
    GVAR_FUNC(SUBGRAPH_HOMEOMORPHIC_TO_K4, 1, "digraph"),
    GVAR_FUNC(DIGRAPHS_FREE_HOMOS_DATA, 0, ""),
    GVAR_FUNC(DIGRAPHS_FREE_CLIQUES_DATA, 0, ""),

    {0, 0, 0, 0, 0} /* Finish with an empty entry */
};

/******************************************************************************
*F  InitKernel( <module> )  . . . . . . . . initialise kernel data structures
*/
static Int InitKernel(StructInitInfo* module) {
  /* init filters and functions                                          */
  InitHdlrFuncsFromTable(GVarFuncs);
  ImportGVarFromLibrary("IsDigraph", &IsDigraph);
  ImportGVarFromLibrary("IsMultiDigraph", &IsMultiDigraph);
  ImportGVarFromLibrary("IsDigraphEdge", &IsDigraphEdge);
  ImportGVarFromLibrary("DIGRAPHS_ValidateVertexColouring",
                        &DIGRAPHS_ValidateVertexColouring);
  ImportGVarFromLibrary("infinity", &Infinity);
  ImportGVarFromLibrary("IsSymmetricDigraph", &IsSymmetricDigraph);
  ImportGVarFromLibrary("AutomorphismGroup", &AutomorphismGroup);
  ImportGVarFromLibrary("GeneratorsOfGroup", &GeneratorsOfGroup);
  ImportGVarFromLibrary("IsAttributeStoringRep", &IsAttributeStoringRepObj);
  ImportGVarFromLibrary("IsPermGroup", &IsPermGroup);
  ImportGVarFromLibrary("IsDigraphAutomorphism", &IsDigraphAutomorphism);
  ImportGVarFromLibrary("LargestMovedPointPerms", &LargestMovedPointPerms);
  ImportGVarFromLibrary("SmallestMovedPointPerm", &SmallestMovedPointPerm);
  ImportGVarFromLibrary("IsClique", &IsClique);
  ImportGVarFromLibrary("IsTrivial", &IsTrivial);
  ImportGVarFromLibrary("Orbit", &Orbit);
  ImportGVarFromLibrary("Stabilizer", &Stabilizer);
  ImportGVarFromLibrary("IsSubset", &IsSubset);
  ImportGVarFromLibrary("OnTuples", &OnTuples);
  ImportGVarFromLibrary("Group", &Group);
  ImportGVarFromLibrary("ClosureGroup", &ClosureGroup);
  ImportGVarFromLibrary("InfoWarning", &InfoWarning);
  /* return success                                                      */
  return 0;
}

/******************************************************************************
*F  InitLibrary( <module> ) . . . . . . .  initialise library data structures
*/
static Int InitLibrary(StructInitInfo* module) {
  /* init filters and functions */
  InitGVarFuncsFromTable(GVarFuncs);

  /* return success                                                      */
  return 0;
}

/******************************************************************************
*F  InitInfopl()  . . . . . . . . . . . . . . . . . table of init functions
*/
static StructInitInfo module = {
    .type        = MODULE_DYNAMIC,
    .name        = "digraphs",
    .initKernel  = InitKernel,
    .initLibrary = InitLibrary,
};

extern StructInitInfo* Init__Dynamic(void);  // prototype to avoid warnings
StructInitInfo*        Init__Dynamic(void) {
  return &module;
}

quality92%

¤ Dauer der Verarbeitung: 0.78 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.