Quelle  homos.c   Sprache: C

//
// homos.c               di/graph homomorphisms              Julius Jonusas
//                                                           J. D. Mitchell
//
// Copyright (C) 2014-19 - Julius Jonusas and J. D. Mitchell
//
// This file is free software, see the digraphs/LICENSE.

////////////////////////////////////////////////////////////////////////////////
//
// This file is organised as follows:
//
// 1. Macros
// 2. Forward declarations
// 3. Global variables
// 4. Hook functions
// 5. Static helper functions
// 6. The main recursive functions (and helpers)
// 7. The GAP-level function (and helpers)
//
////////////////////////////////////////////////////////////////////////////////

// TODO(later)
// 0. remove final arg from push_conditions
// 1. Try other bit hacks for iterating through set bits

#include "homos.h"

// C headers
#include <setjmp.h>   // for longjmp, setjmp, jmp_buf
#include <stdbool.h>  // for true, false, bool
#include <stddef.h>   // for NULL
#include <stdint.h>   // for uint16_t, uint64_t
#include <stdlib.h>   // for malloc, NULL

#ifdef DIGRAPHS_ENABLE_STATS
#include <cstdio>  // for printf
#endif

// GAP headers
#include "gap-includes.h"

// Digraphs package headers
#include "bitarray.h"         // for BitArray
#include "bliss-includes.h"   // for bliss stuff
#include "conditions.h"       // for Conditions
#include "digraphs-config.h"  // for DIGRAPHS_HAVE___BUILTIN_CTZLL
#include "digraphs-debug.h"   // for DIGRAPHS_ASSERT
#include "globals.h"          // for UNDEFINED...
#include "homos-graphs.h"     // for Digraph, Graph, . . .
#include "perms.h"            // for UNDEFINED, PermColl, Perm
#include "safemalloc.h"       // for safe_mallov
#include "schreier-sims.h"    // for PermColl, . . .

#ifdef DIGRAPHS_ENABLE_STATS
// This include has to come after digraphs-config.h since that's where
// DIGRAPHS_ENABLE_STATS is defined
#include <time.h>  // for time
#endif

////////////////////////////////////////////////////////////////////////////////
// 1. Macros
////////////////////////////////////////////////////////////////////////////////

#define MAX(a, b) (a < b ? b : a)
#define MIN(a, b) (a < b ? a : b)

#ifndef HOMOS_STRUCTURE_SIZE
uint16_t HOMOS_STRUCTURE_SIZE = 0;
#endif
#ifndef UNDEFINED
uint16_t UNDEFINED = 65535;
#endif

// The next line can be used instead of the first line of STORE_MIN to
// randomise which vertex of minimum degree is used next, but I didn't find any
// cases where this made things faster.

// if (x < current_x || (x == current_x && (rand() & 1))) {

#define STORE_MIN(current_x, current_y, x, y) \
  if (x < current_x) {                        \
    current_x = x;                            \
    current_y = y;                            \
  }

// See the comment before STORE_MIN

// if (x < current_x || (x == current_x && (rand() & 1))) {

#define STORE_MIN_BREAK(current_x, current_y, x, y) \
  if (x < current_x) {                              \
    current_x = x;                                  \
    current_y = y;                                  \
    if (current_x == 1) {                           \
      break;                                        \
    }                                               \
  }

// The following macro is bitmap_decode_ctz_callpack from https://git.io/fho4p
#if SYS_IS_64_BIT && defined(DIGRAPHS_HAVE___BUILTIN_CTZLL)
#define FOR_SET_BITS(__bit_array, __nr_bits, __variable)     \
  for (size_t k = 0; k < number_of_blocks(__nr_bits); ++k) { \
    Block block = __bit_array->blocks[k];                    \
    while (block != 0) {                                     \
      uint64_t t = block & -block;                           \
      int      r = __builtin_ctzll(block);                   \
      __variable = k * 64 + r;                               \
      if (__variable >= __nr_bits) {                         \
        break;                                               \
      }

#define END_FOR_SET_BITS \
  block ^= t;            \
  }                      \
  }

#else
#define FOR_SET_BITS(__bit_array, __nr_bits, __variable)       \
  for (__variable = 0; __variable < __nr_bits; ++__variable) { \
    if (get_bit_array(__bit_array, __variable)) {
#define END_FOR_SET_BITS \
  }                      \
  }
#endif

////////////////////////////////////////////////////////////////////////////////
// 2. Forward declarations
////////////////////////////////////////////////////////////////////////////////

// Defined in digraphs.h
Int DigraphNrVertices(Obj);
Obj FuncOutNeighbours(Obj, Obj);

// GAP level things, imported in digraphs.c
extern Obj IsDigraph;
extern Obj IsMultiDigraph;
extern Obj DIGRAPHS_ValidateVertexColouring;
extern Obj Infinity;
extern Obj IsSymmetricDigraph;
extern Obj GeneratorsOfGroup;
extern Obj AutomorphismGroup;
extern Obj IsPermGroup;
extern Obj IsDigraphAutomorphism;
extern Obj LargestMovedPointPerms;
extern Obj InfoWarning;

////////////////////////////////////////////////////////////////////////////////
// 3. Global variables
////////////////////////////////////////////////////////////////////////////////

static Obj GAP_FUNC;  // Variable to hold a GAP level hook function

static Obj (*HOOK)(void*,  // HOOK function applied to every homo found
                   const uint16_t,
                   const uint16_t*);
static void* USER_PARAM;  // a USER_PARAM for the hook

static jmp_buf OUTOFHERE;  // so we can jump out of the deepest

static bool ORDERED;  // true if the vertices of the domain/source digraph
                      // should be considered in a different order than they are
                      // given, false otherwise.

static BitArray** BIT_ARRAY_BUFFER = NULL;  // A buffer
static BitArray*  IMAGE_RESTRICT;           // Values in MAP must be in this
static BitArray** MAP_UNDEFINED = NULL;     // UNDEFINED positions in MAP
static BitArray*  ORB_LOOKUP;               // points in orbit
static BitArray*  VALS;                     // Values in MAP already

static BitArray** REPS;  // orbit reps organised by depth

static Conditions* CONDITIONS;

static Digraph* DIGRAPH1;  // Digraphs to hold incoming GAP digraphs
static Digraph* DIGRAPH2;

static Graph* GRAPH1;  // Graphs to hold incoming GAP symmetric digraphs
static Graph* GRAPH2;

static BlissGraph** BLISS_GRAPH = NULL;

static uint16_t* MAP           = NULL;  // partial image list
static uint16_t* COLORS2       = NULL;  // colors of range (di)graph
static uint16_t* INVERSE_ORDER = NULL;  // external -> internal
static uint16_t* MAP_BUFFER    = NULL;  // For converting from internal ->
                                        // external and back when calling the
                                        // hook functions.
static uint16_t* ORB   = NULL;  // Array for containing nodes in an orbit.
static uint16_t* ORDER = NULL;  // internal -> external

static PermColl**    STAB_GENS = NULL;  // stabiliser generators
static SchreierSims* SCHREIER_SIMS;

#ifdef DIGRAPHS_ENABLE_STATS
struct homo_stats_struct {
  time_t last_print;
  size_t max_depth;
  double mean_depth;
  size_t nr_calls;
  size_t nr_dead_branches;
  time_t start_time;
};

typedef struct homo_stats_struct HomoStats;

static HomoStats* STATS;

static inline void clear_stats(HomoStats* stats) {
  stats->last_print       = time(0);
  stats->max_depth        = 0;
  stats->mean_depth       = 0;
  stats->nr_calls         = 0;
  stats->nr_dead_branches = 0;
  stats->start_time       = time(0);
}

static void print_stats(HomoStats* stats) {
  printf("Running for %0.0fs . . .\n", difftime(time(0), stats->start_time));
  printf("Number of function calls = %*lu\n", 20, stats->nr_calls);
  printf("Mean depth = %*.2f\n", 20, stats->mean_depth);
  printf("Max depth = %*lu\n", 20, stats->max_depth);
  printf("Number of dead branches = %*lu\n\n", 20, stats->nr_dead_branches);
}

static inline void update_stats(HomoStats* stats, uint64_t depth) {
  if (depth > stats->max_depth) {
    stats->max_depth = depth;
  }
  stats->nr_calls++;
  stats->mean_depth += (depth - stats->mean_depth) / stats->nr_calls;
  if (difftime(time(0), stats->last_print) > 0.9) {
    print_stats(stats);
    stats->last_print = time(0);
  }
}
#endif  // DIGRAPHS_ENABLE_STATS

////////////////////////////////////////////////////////////////////////////////
// 4. Hook functions
////////////////////////////////////////////////////////////////////////////////

static Obj
homo_hook_gap(void* user_param, uint16_t const nr, uint16_t const* map) {
  UInt2* ptr;
  Obj    t;
  UInt   i;

  // copy map into new trans2
  t   = NEW_TRANS2(nr);
  ptr = ADDR_TRANS2(t);

  for (i = 0; i < nr; i++) {
    ptr[i] = map[i];
  }
  return CALL_2ARGS(GAP_FUNC, user_param, t);
}

static Obj
homo_hook_collect(void* user_param, uint16_t const nr, uint16_t const* map) {
  UInt2* ptr;
  Obj    t;
  UInt   i;

  // copy map into new trans2
  t   = NEW_TRANS2(nr);
  ptr = ADDR_TRANS2(t);

  for (i = 0; i < nr; i++) {
    ptr[i] = map[i];
  }

  ASS_LIST(user_param, LEN_LIST(user_param) + 1, t);
  return False;
}

////////////////////////////////////////////////////////////////////////////////
// 5. Static helper functions
////////////////////////////////////////////////////////////////////////////////

// print_array is not used, but can be useful for debugging.

// static void print_array(uint16_t const* const array, uint16_t const len) {
//   if (array == NULL) {
//     printf("NULL");
//     return;
//   }
//   printf("<array {");
//   for (uint16_t i = 0; i < len; i++) {
//     printf(" %d", array[i]);
//   }
//   printf(" }>");
// }

bool homos_data_initialized = false;  // did we call this method before?
Obj  FuncDIGRAPHS_FREE_HOMOS_DATA(Obj self) {
  if (homos_data_initialized) {
    free_digraph(DIGRAPH1);
    free_digraph(DIGRAPH2);
    free_graph(GRAPH1);
    free_graph(GRAPH2);
    free_bit_array(IMAGE_RESTRICT);
    free_bit_array(ORB_LOOKUP);
    free(MAP);
    free(COLORS2);
    free(INVERSE_ORDER);
    free(MAP_BUFFER);
    free(ORB);
    free(ORDER);

    for (uint16_t i = 0; i < HOMOS_STRUCTURE_SIZE * 3; i++) {
      bliss_digraphs_release(BLISS_GRAPH[i]);
    }

    for (uint16_t i = 0; i < HOMOS_STRUCTURE_SIZE; i++) {
      free_bit_array(REPS[i]);
      free_bit_array(BIT_ARRAY_BUFFER[i]);
      free_bit_array(MAP_UNDEFINED[i]);
      free_perm_coll(STAB_GENS[i]);
    }

    free(BLISS_GRAPH);
    free(REPS);
    free(BIT_ARRAY_BUFFER);
    free(MAP_UNDEFINED);
    free(STAB_GENS);
    free_bit_array(VALS);
    free_conditions(CONDITIONS);
    free_schreier_sims(SCHREIER_SIMS);
    homos_data_initialized = false;
  }

  return 0L;
}

static void get_automorphism_group_from_gap(Obj digraph_obj, PermColl* out) {
  DIGRAPHS_ASSERT(CALL_1ARGS(IsDigraph, digraph_obj) == True);
  if (CALL_1ARGS(IsMultiDigraph, digraph_obj) == True) {
    ErrorQuit("expected a digraph without multiple edges!", 0L, 0L);
  }
  Obj o = CALL_1ARGS(AutomorphismGroup, digraph_obj);
  o     = CALL_1ARGS(GeneratorsOfGroup, o);
  DIGRAPHS_ASSERT(IS_LIST(o));
  clear_perm_coll(out);
  out->degree = PERM_DEGREE;
  DIGRAPHS_ASSERT(out->capacity >= LEN_LIST(o));
  Int k = 0;  // perm coll index
  for (Int i = 1; i <= LEN_LIST(o); ++i) {
    DIGRAPHS_ASSERT(ISB_LIST(o, i));
    DIGRAPHS_ASSERT(IS_PERM2(ELM_LIST(o, i)) || IS_PERM4(ELM_LIST(o, i)));
    Obj p = ELM_LIST(o, i);
    DIGRAPHS_ASSERT(LargestMovedPointPerm(p) <= PERM_DEGREE);
    if (LargestMovedPointPerm(p) == 0) {
      // the index k must be separate from the index i because of this continue
      continue;
    }
    Perm q = out->perms[k++];
    out->size++;
    DIGRAPHS_ASSERT(out->size == k);
    size_t dep;
    if (IS_PERM2(p)) {
      UInt2 const* ptp2 = CONST_ADDR_PERM2(p);
      dep               = MIN((uint16_t) DEG_PERM2(p), PERM_DEGREE);
      for (uint16_t j = 0; j < dep; ++j) {
        q[j] = (uint16_t) ptp2[j];
      }
    } else {
      DIGRAPHS_ASSERT(IS_PERM4(p));
      UInt4 const* ptp4 = CONST_ADDR_PERM4(p);
      dep               = MIN((uint16_t) DEG_PERM4(p), PERM_DEGREE);
      for (uint16_t j = 0; j < dep; ++j) {
        q[j] = (uint16_t) ptp4[j];
      }
    }
    for (uint16_t j = dep; j < PERM_DEGREE; ++j) {
      q[j] = j;
    }
  }
}

static void init_digraph_from_digraph_obj(Digraph* const digraph,
                                          Obj            digraph_obj,
                                          bool const     reorder) {
  DIGRAPHS_ASSERT(digraph != NULL);
  DIGRAPHS_ASSERT(CALL_1ARGS(IsDigraph, digraph_obj) == True);
  UInt const nr  = DigraphNrVertices(digraph_obj);
  Obj        out = FuncOutNeighbours(0L, digraph_obj);
  DIGRAPHS_ASSERT(nr < MAXVERTS);
  DIGRAPHS_ASSERT(IS_PLIST(out));
  clear_digraph(digraph, nr);

  if (!reorder) {
    for (uint16_t i = 1; i <= nr; i++) {
      Obj nbs = ELM_PLIST(out, i);
      DIGRAPHS_ASSERT(IS_LIST(nbs));
      for (uint16_t j = 1; j <= LEN_LIST(nbs); j++) {
        DIGRAPHS_ASSERT(IS_INTOBJ(ELM_LIST(nbs, j)));
        add_edge_digraph(digraph, i - 1, INT_INTOBJ(ELM_LIST(nbs, j)) - 1);
      }
    }
  } else {
    DIGRAPHS_ASSERT(ORDERED);
    for (uint16_t i = 1; i <= nr; i++) {
      Obj nbs = ELM_PLIST(out, ORDER[i - 1] + 1);
      DIGRAPHS_ASSERT(IS_LIST(nbs));
      for (uint16_t j = 1; j <= LEN_LIST(nbs); j++) {
        DIGRAPHS_ASSERT(IS_INTOBJ(ELM_LIST(nbs, j)));
        add_edge_digraph(
            digraph, i - 1, INVERSE_ORDER[INT_INTOBJ(ELM_LIST(nbs, j)) - 1]);
      }
    }
  }
}

static void init_graph_from_digraph_obj(Graph* const graph,
                                        Obj          digraph_obj,
                                        bool const   reorder) {
  DIGRAPHS_ASSERT(graph != NULL);
  DIGRAPHS_ASSERT(CALL_1ARGS(IsDigraph, digraph_obj) == True);
  DIGRAPHS_ASSERT(CALL_1ARGS(IsSymmetricDigraph, digraph_obj) == True);
  UInt const nr  = DigraphNrVertices(digraph_obj);
  Obj        out = FuncOutNeighbours(0L, digraph_obj);
  DIGRAPHS_ASSERT(nr < MAXVERTS);
  DIGRAPHS_ASSERT(IS_PLIST(out));
  clear_graph(graph, nr);

  if (!reorder) {
    for (uint16_t i = 1; i <= nr; i++) {
      Obj nbs = ELM_PLIST(out, i);
      DIGRAPHS_ASSERT(IS_LIST(nbs));
      for (uint16_t j = 1; j <= LEN_LIST(nbs); j++) {
        DIGRAPHS_ASSERT(IS_INTOBJ(ELM_LIST(nbs, j)));
        add_edge_graph(graph, i - 1, INT_INTOBJ(ELM_LIST(nbs, j)) - 1);
      }
    }
  } else {
    DIGRAPHS_ASSERT(ORDERED);
    for (uint16_t i = 1; i <= nr; i++) {  // Nodes in the new graph
      Obj nbs = ELM_PLIST(out, ORDER[i - 1] + 1);
      DIGRAPHS_ASSERT(IS_LIST(nbs));
      for (uint16_t j = 1; j <= LEN_LIST(nbs); j++) {
        DIGRAPHS_ASSERT(IS_INTOBJ(ELM_LIST(nbs, j)));
        add_edge_graph(
            graph, i - 1, INVERSE_ORDER[INT_INTOBJ(ELM_LIST(nbs, j)) - 1]);
      }
    }
  }
}

// Find orbit representatives of the group generated by STAB_GENS[rep_depth]
// and store them in REPS[rep_depth], only values in IMAGE_RESTRICT will be
// chosen as orbit representatives.
static bool compute_stabs_and_orbit_reps(uint16_t const nr_nodes_1,
                                         uint16_t const nr_nodes_2,
                                         uint16_t const rep_depth,
                                         uint16_t const depth,
                                         uint16_t const pt,
                                         bool const     first_call) {
  DIGRAPHS_ASSERT(rep_depth <= depth + 1);
  if (depth == nr_nodes_1 - 1 && !first_call) {
    // first_call is required in the case that nr_nodes_1 is 1, since without
    // first_call this function returns false (in the next line) and REPS is
    // not initialised.
    return false;  // This doesn't really say anything about the stabiliser
  } else if (rep_depth > 0) {
    point_stabilizer(
        SCHREIER_SIMS, STAB_GENS[rep_depth - 1], STAB_GENS[rep_depth], pt);
    if (STAB_GENS[rep_depth]->size == 0) {
      // the stabiliser of pt in STAB_GENS[rep_depth - 1] is trivial
      copy_bit_array(REPS[rep_depth], IMAGE_RESTRICT, nr_nodes_2);
      complement_bit_arrays(REPS[rep_depth], VALS, nr_nodes_2);
      // REPS[rep_depth] is a set of orbit representatives of the stabiliser of
      // the existing values in MAP, which belong to VALS, and since the
      // stabiliser is trivial, we take valid choice as the set of
      // representatives.
      return true;  // the stabiliser is trivial
    }
  }
  init_bit_array(REPS[rep_depth], false, nr_nodes_2);
  copy_bit_array(ORB_LOOKUP, VALS, nr_nodes_2);
  uint16_t fst = 0;
  while (fst < PERM_DEGREE
         && (get_bit_array(ORB_LOOKUP, fst)
             || !get_bit_array(IMAGE_RESTRICT, fst))) {
    fst++;
  }

  while (fst < PERM_DEGREE) {
    ORB[0]     = fst;
    uint16_t n = 1;  // length of ORB

    set_bit_array(REPS[rep_depth], fst, true);
    set_bit_array(ORB_LOOKUP, fst, true);

    for (uint16_t i = 0; i < n; ++i) {
      for (uint16_t j = 0; j < STAB_GENS[rep_depth]->size; ++j) {
        Perm           gen = STAB_GENS[rep_depth]->perms[j];
        uint16_t const img = gen[ORB[i]];
        if (!get_bit_array(ORB_LOOKUP, img)) {
          ORB[n++] = img;
          set_bit_array(ORB_LOOKUP, img, true);
        }
      }
    }
    while (fst < PERM_DEGREE
           && (get_bit_array(ORB_LOOKUP, fst)
               || !get_bit_array(IMAGE_RESTRICT, fst))) {
      fst++;
    }
  }
  return false;  // the stabiliser is not trivial
}

// Rewrite the global variable MAP so that it contains a homomorphism of the
// original GAP level (di)graph, and not the possibly distinct (but isomorphic)
// copy in the homomorphism search.  This should be called before any calls to
// the hook functions (i.e. after an entire homomorphism is found).
static void external_order_map_digraph(Digraph* digraph) {
  if (!ORDERED) {
    return;
  }
  for (uint16_t i = 0; i < digraph->nr_vertices; ++i) {
    MAP_BUFFER[ORDER[i]] = MAP[i];
  }
  for (uint16_t i = 0; i < digraph->nr_vertices; ++i) {
    MAP[i] = MAP_BUFFER[i];
  }
}

static void external_order_map_graph(Graph* graph) {
  if (!ORDERED) {
    return;
  }
  for (uint16_t i = 0; i < graph->nr_vertices; ++i) {
    MAP_BUFFER[ORDER[i]] = MAP[i];
  }
  for (uint16_t i = 0; i < graph->nr_vertices; ++i) {
    MAP[i] = MAP_BUFFER[i];
  }
}

// Rewrite the global variable MAP so that it contains a homomorphism of the
// internal (di)graph, and not the possibly distinct (but isomorphic) GAP level
// (di)graph.  This should be called after any calls to the hook functions
// (i.e. after an entire homomorphism is found).
static void internal_order_map_digraph(Digraph const* const digraph) {
  if (!ORDERED) {
    return;
  }
  for (uint16_t i = 0; i < digraph->nr_vertices; ++i) {
    MAP_BUFFER[INVERSE_ORDER[i]] = MAP[i];
  }
  for (uint16_t i = 0; i < digraph->nr_vertices; ++i) {
    MAP[i] = MAP_BUFFER[i];
  }
}

static void internal_order_map_graph(Graph const* const graph) {
  if (!ORDERED) {
    return;
  }
  for (uint16_t i = 0; i < graph->nr_vertices; ++i) {
    MAP_BUFFER[INVERSE_ORDER[i]] = MAP[i];
  }
  for (uint16_t i = 0; i < graph->nr_vertices; ++i) {
    MAP[i] = MAP_BUFFER[i];
  }
}

static void set_automorphisms(Obj aut_grp_obj, PermColl* out) {
  DIGRAPHS_ASSERT(out != NULL);
  clear_perm_coll(out);
  out->degree = PERM_DEGREE;
  Obj gens    = CALL_1ARGS(GeneratorsOfGroup, aut_grp_obj);
  DIGRAPHS_ASSERT(IS_LIST(gens));
  DIGRAPHS_ASSERT(LEN_LIST(gens) > 0);
  for (Int i = 1; i <= LEN_LIST(gens); ++i) {
    Obj gen_obj = ELM_LIST(gens, i);
    if (LargestMovedPointPerm(gen_obj) > 0) {
      Perm const p = new_perm_from_gap(gen_obj, PERM_DEGREE);
      add_perm_coll(out, p);
      free(p);
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
// 6. The main recursive functions (and helpers)
////////////////////////////////////////////////////////////////////////////////

// Helper for the main recursive homomorphism function.
static ALWAYS_INLINE uint16_t
graph_homo_update_conditions(uint16_t const depth,
                             uint16_t const last_defined,
                             uint16_t const vertex) {
  push_conditions(
      CONDITIONS, depth, vertex, GRAPH2->neighbours[MAP[last_defined]]);
  store_size_conditions(CONDITIONS, vertex);
  return size_conditions(CONDITIONS, vertex);
}

// The main recursive function for homomorphisms of graphs.
//
// The arguments are:
//
// 1. depth             The current depth of the search (i.e. number of
//                      positions in MAP that are assigned).
// 2. pos               The last position in MAP that was filled
// 3. rep_depth         The depth of the stabilisers of points in MAP already.
//                      This is different than depth because we stop increasing
//                      the rep_depth when we first encounter a trivial
//                      stabiliser.
// 4. has_trivial_stab  true if the stabiliser of the values already in MAP is
//                      trivial, false otherwise.
// 5. rank              The current number of distinct values in MAP.
// 6. max_results       The maximum number of results to find.
// 7. hint              The desired number of distinct points in a (full)
//                      homomorphism.
// 8. count             The number of homomorphisms found so far.
static void find_graph_homos(uint16_t        depth,
                             uint16_t        pos,
                             uint16_t        rep_depth,
                             bool            has_trivial_stab,
                             uint16_t        rank,
                             uint64_t const  max_results,
                             uint64_t const  hint,
                             uint64_t* const count) {
#ifdef DIGRAPHS_ENABLE_STATS
  update_stats(STATS, depth);
#endif
  if (depth == GRAPH1->nr_vertices) {
    // Every position in MAP is assigned . . .
    if (hint != UNDEFINED && rank != hint) {
#ifdef DIGRAPHS_ENABLE_STATS
      STATS->nr_dead_branches++;
#endif
      return;
    }
    external_order_map_graph(GRAPH1);
    Obj ret =
        HOOK(USER_PARAM, MAX(GRAPH1->nr_vertices, GRAPH2->nr_vertices), MAP);
    internal_order_map_graph(GRAPH1);
    (*count)++;
    if (*count >= max_results || ret == True) {
      longjmp(OUTOFHERE, 1);
    }
    return;
  }

  uint16_t next = 0;          // the next position to fill
  uint16_t min  = UNDEFINED;  // the minimum number of candidates for MAP[next]
  uint16_t i;

  BitArray* possible = BIT_ARRAY_BUFFER[depth];

  if (depth > 0) {  // this is not the first call of the function
    copy_bit_array(
        MAP_UNDEFINED[depth], MAP_UNDEFINED[depth - 1], GRAPH1->nr_vertices);
    copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
    intersect_bit_arrays(
        possible, GRAPH1->neighbours[pos], GRAPH1->nr_vertices);
    FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
      size_t const n = graph_homo_update_conditions(depth, pos, i);
      if (n == 0) {
#ifdef DIGRAPHS_ENABLE_STATS
        STATS->nr_dead_branches++;
#endif
        pop_conditions(CONDITIONS, depth);
        return;
      }
      STORE_MIN(min, next, n, i);
    }
    END_FOR_SET_BITS
    if (min > 1) {
      copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
      complement_bit_arrays(
          possible, GRAPH1->neighbours[pos], GRAPH1->nr_vertices);
      FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
        size_t const n = size_conditions(CONDITIONS, i);
        STORE_MIN_BREAK(min, next, n, i);
      }
      END_FOR_SET_BITS
    }
  } else {
    for (i = 0; i < GRAPH1->nr_vertices; ++i) {
      size_t const n = size_conditions(CONDITIONS, i);
      STORE_MIN_BREAK(min, next, n, i);
    }
  }
  DIGRAPHS_ASSERT(get_bit_array(MAP_UNDEFINED[depth], next));
  DIGRAPHS_ASSERT(next < GRAPH1->nr_vertices);

  if (rank < hint) {
    copy_bit_array(
        possible, get_conditions(CONDITIONS, next), GRAPH2->nr_vertices);
    complement_bit_arrays(possible, VALS, GRAPH2->nr_vertices);
    intersect_bit_arrays(possible, REPS[rep_depth], GRAPH2->nr_vertices);
    FOR_SET_BITS(possible, GRAPH2->nr_vertices, i) {
      MAP[next] = i;
      set_bit_array(VALS, i, true);
      set_bit_array(MAP_UNDEFINED[depth], next, false);
      if (!has_trivial_stab) {
        find_graph_homos(depth + 1,
                         next,
                         rep_depth + 1,
                         compute_stabs_and_orbit_reps(GRAPH1->nr_vertices,
                                                      GRAPH2->nr_vertices,
                                                      rep_depth + 1,
                                                      depth,
                                                      i,
                                                      false),
                         rank + 1,
                         max_results,
                         hint,
                         count);
      } else {
        find_graph_homos(depth + 1,
                         next,
                         rep_depth,
                         true,
                         rank + 1,
                         max_results,
                         hint,
                         count);
      }
      MAP[next] = UNDEFINED;
      set_bit_array(VALS, i, false);
      set_bit_array(MAP_UNDEFINED[depth], next, true);
    }
    END_FOR_SET_BITS
  }
  copy_bit_array(
      possible, get_conditions(CONDITIONS, next), GRAPH2->nr_vertices);
  intersect_bit_arrays(possible, VALS, GRAPH2->nr_vertices);
  FOR_SET_BITS(possible, GRAPH2->nr_vertices, i) {
    MAP[next] = i;
    set_bit_array(MAP_UNDEFINED[depth], next, false);
    find_graph_homos(depth + 1,
                     next,
                     rep_depth,
                     has_trivial_stab,
                     rank,
                     max_results,
                     hint,
                     count);
    MAP[next] = UNDEFINED;
    set_bit_array(MAP_UNDEFINED[depth], next, true);
  }
  END_FOR_SET_BITS
  pop_conditions(CONDITIONS, depth);
}

// Helper for the main recursive monomorphism function.
static ALWAYS_INLINE uint16_t
graph_mono_update_conditions(uint16_t const depth,
                             uint16_t const last_defined,
                             uint16_t const vertex) {
  push_conditions(
      CONDITIONS, depth, vertex, GRAPH2->neighbours[MAP[last_defined]]);
  store_size_conditions(CONDITIONS, vertex);
  return size_conditions(CONDITIONS, vertex);
}

// The main recursive function for monomorphisms of graphs.
//
// The arguments are:
//
// 1. depth             The current depth of the search (i.e. number of
//                      positions in MAP that are assigned).
// 2. pos               The last position in MAP that was filled
// 3. rep_depth         The depth of the stabilisers of points in MAP already.
//                      This is different than depth because we stop increasing
//                      the rep_depth when we first encounter a trivial
//                      stabiliser.
// 4. has_trivial_stab  true if the stabiliser of the values already in MAP is
//                      trivial, false otherwise.
// 5. max_results       The maximum number of results to find.
// 6. count             The number of embeddings found so far.
static void find_graph_monos(uint16_t        depth,
                             uint16_t        pos,
                             uint16_t        rep_depth,
                             bool            has_trivial_stab,
                             uint64_t const  max_results,
                             uint64_t* const count) {
#ifdef DIGRAPHS_ENABLE_STATS
  update_stats(STATS, depth);
#endif
  if (depth == GRAPH1->nr_vertices) {
    // we've assigned every position in <MAP>
    external_order_map_graph(GRAPH1);
    Obj ret =
        HOOK(USER_PARAM, MAX(GRAPH1->nr_vertices, GRAPH2->nr_vertices), MAP);
    internal_order_map_graph(GRAPH1);
    (*count)++;
    if (*count >= max_results || ret == True) {
      longjmp(OUTOFHERE, 1);
    }
    return;
  }

  uint16_t next = 0;          // the next position to fill
  uint16_t min  = UNDEFINED;  // the minimum number of candidates for MAP[next]
  uint16_t i;

  BitArray* possible = BIT_ARRAY_BUFFER[depth];

  if (depth > 0) {  // this is not the first call of the function
    copy_bit_array(
        MAP_UNDEFINED[depth], MAP_UNDEFINED[depth - 1], GRAPH1->nr_vertices);
    copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
    intersect_bit_arrays(
        possible, GRAPH1->neighbours[pos], GRAPH1->nr_vertices);
    FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
      size_t const n = graph_mono_update_conditions(depth, pos, i);
      if (n == 0) {
#ifdef DIGRAPHS_ENABLE_STATS
        STATS->nr_dead_branches++;
#endif
        pop_conditions(CONDITIONS, depth);
        return;
      }
      STORE_MIN(min, next, n, i);
    }
    END_FOR_SET_BITS
    if (min > 1) {
      copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
      complement_bit_arrays(
          possible, GRAPH1->neighbours[pos], GRAPH1->nr_vertices);
      FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
        size_t const n = size_conditions(CONDITIONS, i);
        STORE_MIN_BREAK(min, next, n, i);
      }
      END_FOR_SET_BITS
    }
  } else {
    for (i = 0; i < GRAPH1->nr_vertices; ++i) {
      size_t const n = size_conditions(CONDITIONS, i);
      STORE_MIN_BREAK(min, next, n, i);
    }
  }
  copy_bit_array(
      possible, get_conditions(CONDITIONS, next), GRAPH2->nr_vertices);
  intersect_bit_arrays(possible, REPS[rep_depth], GRAPH2->nr_vertices);
  complement_bit_arrays(possible, VALS, GRAPH2->nr_vertices);
  FOR_SET_BITS(possible, GRAPH2->nr_vertices, i) {
    MAP[next] = i;
    set_bit_array(VALS, i, true);
    set_bit_array(MAP_UNDEFINED[depth], next, false);
    if (!has_trivial_stab) {
      find_graph_monos(depth + 1,
                       next,
                       rep_depth + 1,
                       compute_stabs_and_orbit_reps(GRAPH1->nr_vertices,
                                                    GRAPH2->nr_vertices,
                                                    rep_depth + 1,
                                                    depth,
                                                    i,
                                                    false),
                       max_results,
                       count);
    } else {
      find_graph_monos(depth + 1, next, rep_depth, true, max_results, count);
    }
    MAP[next] = UNDEFINED;
    set_bit_array(VALS, i, false);
    set_bit_array(MAP_UNDEFINED[depth], next, true);
  }
  END_FOR_SET_BITS
  pop_conditions(CONDITIONS, depth);
}

// Helper for the main recursive embedding function.
static ALWAYS_INLINE uint16_t graph_embed_update_conditions(
    uint16_t const depth,
    uint16_t const last_defined,
    uint16_t const vertex,
    void (*oper)(BitArray* const, BitArray const* const, uint16_t const)) {
  push_conditions(CONDITIONS, depth, vertex, NULL);
  oper(get_conditions(CONDITIONS, vertex),
       GRAPH2->neighbours[MAP[last_defined]],
       GRAPH2->nr_vertices);
  store_size_conditions(CONDITIONS, vertex);
  return size_conditions(CONDITIONS, vertex);
}

// The main recursive function for embeddings of graphs.
//
// The arguments are:
//
// 1. depth             The current depth of the search (i.e. number of
//                      positions in MAP that are assigned).
// 2. pos               The last position in MAP that was filled
// 3. rep_depth         The depth of the stabilisers of points in MAP already.
//                      This is different than depth because we stop increasing
//                      the rep_depth when we first encounter a trivial
//                      stabiliser.
// 4. has_trivial_stab  true if the stabiliser of the values already in MAP is
//                      trivial, false otherwise.
// 5. max_results       The maximum number of results to find.
// 6. count             The number of embeddings found so far.
static void find_graph_embeddings(uint16_t        depth,
                                  uint16_t        pos,
                                  uint16_t        rep_depth,
                                  bool            has_trivial_stab,
                                  uint64_t const  max_results,
                                  uint64_t* const count) {
#ifdef DIGRAPHS_ENABLE_STATS
  update_stats(STATS, depth);
#endif
  if (depth == GRAPH1->nr_vertices) {
    // we've assigned every position in <MAP>
    external_order_map_graph(GRAPH1);
    Obj ret =
        HOOK(USER_PARAM, MAX(GRAPH1->nr_vertices, GRAPH2->nr_vertices), MAP);
    internal_order_map_graph(GRAPH1);
    (*count)++;
    if (*count >= max_results || ret == True) {
      longjmp(OUTOFHERE, 1);
    }
    return;
  }

  uint16_t next = 0;          // the next position to fill
  uint16_t min  = UNDEFINED;  // the minimum number of candidates for MAP[next]
  uint16_t i;

  BitArray* possible = BIT_ARRAY_BUFFER[depth];

  if (depth > 0) {  // this is not the first call of the function
    copy_bit_array(
        MAP_UNDEFINED[depth], MAP_UNDEFINED[depth - 1], GRAPH1->nr_vertices);
    copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
    intersect_bit_arrays(
        possible, GRAPH1->neighbours[pos], GRAPH1->nr_vertices);
    FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
      size_t const n =
          graph_embed_update_conditions(depth, pos, i, &intersect_bit_arrays);
      if (n == 0) {
#ifdef DIGRAPHS_ENABLE_STATS
        STATS->nr_dead_branches++;
#endif
        pop_conditions(CONDITIONS, depth);
        return;
      }
      STORE_MIN(min, next, n, i);
    }
    END_FOR_SET_BITS
    copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
    complement_bit_arrays(
        possible, GRAPH1->neighbours[pos], GRAPH1->nr_vertices);
    FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
      size_t const n =
          graph_embed_update_conditions(depth, pos, i, &complement_bit_arrays);
      if (n == 0) {
#ifdef DIGRAPHS_ENABLE_STATS
        STATS->nr_dead_branches++;
#endif
        pop_conditions(CONDITIONS, depth);
        return;
      }
      STORE_MIN(min, next, n, i);
    }
    END_FOR_SET_BITS
  } else {
    for (i = 0; i < GRAPH1->nr_vertices; ++i) {
      size_t const n = size_conditions(CONDITIONS, i);
      STORE_MIN_BREAK(min, next, n, i);
    }
  }

  copy_bit_array(
      possible, get_conditions(CONDITIONS, next), GRAPH2->nr_vertices);
  intersect_bit_arrays(possible, REPS[rep_depth], GRAPH2->nr_vertices);
  complement_bit_arrays(possible, VALS, GRAPH2->nr_vertices);

  FOR_SET_BITS(possible, GRAPH2->nr_vertices, i) {
    MAP[next] = i;
    set_bit_array(VALS, i, true);
    set_bit_array(MAP_UNDEFINED[depth], next, false);
    if (!has_trivial_stab) {
      find_graph_embeddings(depth + 1,
                            next,
                            rep_depth + 1,
                            compute_stabs_and_orbit_reps(GRAPH1->nr_vertices,
                                                         GRAPH2->nr_vertices,
                                                         rep_depth + 1,
                                                         depth,
                                                         MAP[next],
                                                         false),
                            max_results,
                            count);
    } else {
      find_graph_embeddings(
          depth + 1, next, rep_depth, true, max_results, count);
    }
    MAP[next] = UNDEFINED;
    set_bit_array(VALS, i, false);
    set_bit_array(MAP_UNDEFINED[depth], next, true);
  }
  END_FOR_SET_BITS
  pop_conditions(CONDITIONS, depth);
}

////////////////////////////////////////////////////////////////////////////////
// The next function should be called before find_graph_homos.  This function
// simulates the first steps of the recursion using the information in
// partial_map_obj (if any) so that the values specified in partial_map_obj are
// installed in MAP first.
////////////////////////////////////////////////////////////////////////////////

static void init_partial_map_and_find_graph_homos(Obj partial_map_obj,
                                                  uint64_t const  max_results,
                                                  uint64_t const  hint,
                                                  uint64_t* const count,
                                                  Obj injective_obj) {
  uint16_t depth                = 0;
  uint16_t rep_depth            = 0;
  uint16_t last_defined         = UNDEFINED;
  bool     last_stab_is_trivial = (STAB_GENS[0]->size == 0 ? true : false);
  uint16_t rank                 = 0;
  uint16_t next                 = UNDEFINED;

  if (partial_map_obj != Fail) {
    for (next = 0; next < LEN_LIST(partial_map_obj); ++next) {
      if (ISB_LIST(partial_map_obj, next + 1)) {
        if (depth > 0) {
          DIGRAPHS_ASSERT(last_defined != UNDEFINED);
          copy_bit_array(MAP_UNDEFINED[depth],
                         MAP_UNDEFINED[depth - 1],
                         GRAPH1->nr_vertices);
          BitArray* possible = BIT_ARRAY_BUFFER[0];
          copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
          intersect_bit_arrays(
              possible, GRAPH1->neighbours[last_defined], GRAPH1->nr_vertices);
          if (INT_INTOBJ(injective_obj) == 0) {
            uint16_t i;  // variable for the next for-loop
            FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
              uint16_t n = graph_homo_update_conditions(depth, last_defined, i);
              if (n == 0) {
                return;
              }
            }
            END_FOR_SET_BITS
          } else if (INT_INTOBJ(injective_obj) == 1) {
            uint16_t i;  // variable for the next for-loop
            FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
              uint16_t n = graph_mono_update_conditions(depth, last_defined, i);
              if (n == 0) {
                return;
              }
            }
            END_FOR_SET_BITS
          } else {
            DIGRAPHS_ASSERT(INT_INTOBJ(injective_obj) == 2);
            uint16_t i;  // variable for the next for-loop
            FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
              uint16_t n = graph_embed_update_conditions(
                  depth, last_defined, i, &intersect_bit_arrays);
              if (n == 0) {
                return;
              }
              END_FOR_SET_BITS
            }
            copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
            complement_bit_arrays(possible,
                                  GRAPH1->neighbours[last_defined],
                                  GRAPH1->nr_vertices);
            FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
              uint16_t n = graph_embed_update_conditions(
                  depth, last_defined, i, &complement_bit_arrays);
              if (n == 0) {
                return;
              }
            }
            END_FOR_SET_BITS
          }
        }
        uint16_t val = INT_INTOBJ(ELM_LIST(partial_map_obj, next + 1)) - 1;
        next         = (ORDERED ? INVERSE_ORDER[next] : next);
        MAP[next]    = val;
        if (!get_bit_array(VALS, MAP[next])) {
          rank++;
          if (rank > hint) {
            return;
          }
        }
        set_bit_array(VALS, MAP[next], true);
        set_bit_array(MAP_UNDEFINED[depth], next, false);
        if (!last_stab_is_trivial) {
          last_stab_is_trivial =
              compute_stabs_and_orbit_reps(GRAPH1->nr_vertices,
                                           GRAPH2->nr_vertices,
                                           rep_depth + 1,
                                           depth,
                                           MAP[next],
                                           false);
          rep_depth++;
        }
        depth++;
        last_defined = next;
        next         = (ORDERED ? ORDER[next] : next);
      }
    }
  }
  if (INT_INTOBJ(injective_obj) == 0) {
    find_graph_homos(depth,
                     last_defined,
                     rep_depth,
                     last_stab_is_trivial,
                     rank,
                     max_results,
                     hint,
                     count);
  } else if (INT_INTOBJ(injective_obj) == 1) {
    find_graph_monos(depth,
                     last_defined,
                     rep_depth,
                     last_stab_is_trivial,
                     max_results,
                     count);
  } else if (INT_INTOBJ(injective_obj) == 2) {
    find_graph_embeddings(depth,
                          last_defined,
                          rep_depth,
                          last_stab_is_trivial,
                          max_results,
                          count);
  }
}

// Helper for the main recursive homomorphism of digraphs function.
static ALWAYS_INLINE uint16_t
digraph_homo_update_conditions(uint16_t const depth,
                               uint16_t const last_defined,
                               uint16_t const vertex) {
  if (is_adjacent_digraph(DIGRAPH1, last_defined, vertex)) {
    push_conditions(
        CONDITIONS, depth, vertex, DIGRAPH2->out_neighbours[MAP[last_defined]]);
    if (is_adjacent_digraph(DIGRAPH1, vertex, last_defined)) {
      intersect_bit_arrays(get_conditions(CONDITIONS, vertex),
                           DIGRAPH2->in_neighbours[MAP[last_defined]],
                           DIGRAPH2->nr_vertices);
    }
    store_size_conditions(CONDITIONS, vertex);
  } else if (is_adjacent_digraph(DIGRAPH1, vertex, last_defined)) {
    push_conditions(
        CONDITIONS, depth, vertex, DIGRAPH2->in_neighbours[MAP[last_defined]]);
    store_size_conditions(CONDITIONS, vertex);
  }
  return size_conditions(CONDITIONS, vertex);
}

// The main recursive function for homomorphisms of digraphs.
//
// The arguments are:
//
// 1. depth             The current depth of the search (i.e. number of
//                      positions in MAP that are assigned).
// 2. pos               The last position in MAP that was filled
// 3. rep_depth         The depth of the stabilisers of points in MAP already.
//                      This is different than depth because we stop increasing
//                      the rep_depth when we first encounter a trivial
//                      stabiliser.
// 4. has_trivial_stab  true if the stabiliser of the values already in MAP is
//                      trivial, false otherwise.
// 5. rank              The current number of distinct values in MAP.
// 6. max_results       The maximum number of results to find.
// 7. hint              The desired number of distinct points in a (full)
//                      homomorphism.
// 8. count             The number of homomorphisms found so far.
static void find_digraph_homos(uint16_t        depth,
                               uint16_t        pos,
                               uint16_t        rep_depth,
                               bool            has_trivial_stab,
                               uint16_t        rank,
                               uint64_t const  max_results,
                               uint64_t const  hint,
                               uint64_t* const count) {
#ifdef DIGRAPHS_ENABLE_STATS
  update_stats(STATS, depth);
#endif
  if (depth == DIGRAPH1->nr_vertices) {
    // we've assigned every position in <MAP>
    if (hint != UNDEFINED && rank != hint) {
#ifdef DIGRAPHS_ENABLE_STATS
      STATS->nr_dead_branches++;
#endif
      return;
    }
    external_order_map_digraph(DIGRAPH1);
    Obj ret = HOOK(
        USER_PARAM, MAX(DIGRAPH1->nr_vertices, DIGRAPH2->nr_vertices), MAP);
    internal_order_map_digraph(DIGRAPH1);
    (*count)++;
    if (*count >= max_results || ret == True) {
      longjmp(OUTOFHERE, 1);
    }
    return;
  }

  uint16_t next = 0;          // the next position to fill
  uint16_t min  = UNDEFINED;  // the minimum number of candidates for MAP[next]
  uint16_t i;

  if (depth > 0) {  // this is not the first call of the function
    copy_bit_array(
        MAP_UNDEFINED[depth], MAP_UNDEFINED[depth - 1], DIGRAPH1->nr_vertices);
    FOR_SET_BITS(MAP_UNDEFINED[depth], DIGRAPH1->nr_vertices, i) {
      uint16_t const n = digraph_homo_update_conditions(depth, pos, i);
      if (n == 0) {
#ifdef DIGRAPHS_ENABLE_STATS
        STATS->nr_dead_branches++;
#endif
        pop_conditions(CONDITIONS, depth);
        return;
      }
      STORE_MIN(min, next, n, i);
    }
    END_FOR_SET_BITS
  } else {  // depth == 0
    for (i = 0; i < DIGRAPH1->nr_vertices; ++i) {
      size_t const n = size_conditions(CONDITIONS, i);
      STORE_MIN_BREAK(min, next, n, i);
    }
  }

  BitArray* possible = BIT_ARRAY_BUFFER[depth];

  if (rank < hint) {
    copy_bit_array(
        possible, get_conditions(CONDITIONS, next), DIGRAPH2->nr_vertices);
    complement_bit_arrays(possible, VALS, DIGRAPH2->nr_vertices);
    intersect_bit_arrays(possible, REPS[rep_depth], DIGRAPH2->nr_vertices);
    FOR_SET_BITS(possible, DIGRAPH2->nr_vertices, i) {
      MAP[next] = i;
      set_bit_array(VALS, i, true);
      set_bit_array(MAP_UNDEFINED[depth], next, false);
      if (!has_trivial_stab) {
        find_digraph_homos(depth + 1,
                           next,
                           rep_depth + 1,
                           compute_stabs_and_orbit_reps(DIGRAPH1->nr_vertices,
                                                        DIGRAPH2->nr_vertices,
                                                        rep_depth + 1,
                                                        depth,
                                                        i,
                                                        false),
                           rank + 1,
                           max_results,
                           hint,
                           count);
      } else {
        find_digraph_homos(depth + 1,
                           next,
                           rep_depth,
                           true,
                           rank + 1,
                           max_results,
                           hint,
                           count);
      }
      MAP[next] = UNDEFINED;
      set_bit_array(VALS, i, false);
      set_bit_array(MAP_UNDEFINED[depth], next, true);
    }
    END_FOR_SET_BITS
  }
  copy_bit_array(
      possible, get_conditions(CONDITIONS, next), DIGRAPH2->nr_vertices);
  intersect_bit_arrays(possible, VALS, DIGRAPH2->nr_vertices);
  FOR_SET_BITS(possible, DIGRAPH2->nr_vertices, i) {
    MAP[next] = i;
    set_bit_array(MAP_UNDEFINED[depth], next, false);
    find_digraph_homos(depth + 1,
                       next,
                       rep_depth,
                       has_trivial_stab,
                       rank,
                       max_results,
                       hint,
                       count);
    MAP[next] = UNDEFINED;
    set_bit_array(MAP_UNDEFINED[depth], next, true);
  }
  END_FOR_SET_BITS
  pop_conditions(CONDITIONS, depth);
}

// Helper for the main recursive monomorphism of digraphs function.
static ALWAYS_INLINE uint16_t
digraph_mono_update_conditions(uint16_t const depth,
                               uint16_t const last_defined,
                               uint16_t const vertex) {
  push_conditions(CONDITIONS, depth, vertex, NULL);
  if (is_adjacent_digraph(DIGRAPH1, last_defined, vertex)) {
    intersect_bit_arrays(get_conditions(CONDITIONS, vertex),
                         DIGRAPH2->out_neighbours[MAP[last_defined]],
                         DIGRAPH2->nr_vertices);
  }
  if (is_adjacent_digraph(DIGRAPH1, vertex, last_defined)) {
    intersect_bit_arrays(get_conditions(CONDITIONS, vertex),
                         DIGRAPH2->in_neighbours[MAP[last_defined]],
                         DIGRAPH2->nr_vertices);
  }
  store_size_conditions(CONDITIONS, vertex);

  return size_conditions(CONDITIONS, vertex);
}

// The main recursive function for monomorphisms of digraphs.
//
// The arguments are:
//
// 1. depth             The current depth of the search (i.e. number of
//                      positions in MAP that are assigned).
// 2. pos               The last position in MAP that was filled
// 3. rep_depth         The depth of the stabilisers of points in MAP already.
//                      This is different than depth because we stop increasing
//                      the rep_depth when we first encounter a trivial
//                      stabiliser.
// 4. has_trivial_stab  true if the stabiliser of the values already in MAP is
//                      trivial, false otherwise.
// 5. max_results       The maximum number of results to find.
// 6. count             The number of embeddings found so far.
static void find_digraph_monos(uint16_t        depth,
                               uint16_t        pos,
                               uint16_t        rep_depth,
                               bool            has_trivial_stab,
                               uint64_t const  max_results,
                               uint64_t* const count) {
#ifdef DIGRAPHS_ENABLE_STATS
  update_stats(STATS, depth);
#endif
  if (depth == DIGRAPH1->nr_vertices) {
    // we've assigned every position in <MAP>
    external_order_map_digraph(DIGRAPH1);
    Obj ret = HOOK(
        USER_PARAM, MAX(DIGRAPH1->nr_vertices, DIGRAPH2->nr_vertices), MAP);
    internal_order_map_digraph(DIGRAPH1);
    (*count)++;
    if (*count >= max_results || ret == True) {
      longjmp(OUTOFHERE, 1);
    }
    return;
  }

  uint16_t next = 0;          // the next position to fill
  uint16_t min  = UNDEFINED;  // the minimum number of candidates for MAP[next]
  uint16_t i;

  if (depth > 0) {  // this is not the first call of the function
    copy_bit_array(
        MAP_UNDEFINED[depth], MAP_UNDEFINED[depth - 1], DIGRAPH1->nr_vertices);
    FOR_SET_BITS(MAP_UNDEFINED[depth], DIGRAPH1->nr_vertices, i) {
      size_t const n = digraph_mono_update_conditions(depth, pos, i);
      if (n == 0) {
#ifdef DIGRAPHS_ENABLE_STATS
        STATS->nr_dead_branches++;
#endif
        pop_conditions(CONDITIONS, depth);
        return;
      }
      STORE_MIN(min, next, n, i);
    }
    END_FOR_SET_BITS
  } else {
    for (i = 0; i < DIGRAPH1->nr_vertices; i++) {
      size_t const n = size_conditions(CONDITIONS, i);
      STORE_MIN_BREAK(min, next, n, i);
    }
  }

  BitArray* possible = BIT_ARRAY_BUFFER[depth];
  copy_bit_array(
      possible, get_conditions(CONDITIONS, next), DIGRAPH2->nr_vertices);
  intersect_bit_arrays(possible, REPS[rep_depth], DIGRAPH2->nr_vertices);
  complement_bit_arrays(possible, VALS, DIGRAPH2->nr_vertices);
  FOR_SET_BITS(possible, DIGRAPH2->nr_vertices, i) {
    MAP[next] = i;
    set_bit_array(VALS, i, true);
    set_bit_array(MAP_UNDEFINED[depth], next, false);
    if (!has_trivial_stab) {
      find_digraph_monos(depth + 1,
                         next,
                         rep_depth + 1,
                         compute_stabs_and_orbit_reps(DIGRAPH1->nr_vertices,
                                                      DIGRAPH2->nr_vertices,
                                                      rep_depth + 1,
                                                      depth,
                                                      i,
                                                      false),
                         max_results,
                         count);
    } else {
      find_digraph_monos(depth + 1, next, rep_depth, true, max_results, count);
    }
    MAP[next] = UNDEFINED;
    set_bit_array(VALS, i, false);
    set_bit_array(MAP_UNDEFINED[depth], next, true);
  }
  END_FOR_SET_BITS
  pop_conditions(CONDITIONS, depth);
}

// Helper for the main recursive embedding digraphs function.
static ALWAYS_INLINE uint16_t
digraph_embed_update_conditions(uint16_t const depth,
                                uint16_t const last_def,
                                uint16_t const vertex) {
  push_conditions(CONDITIONS, depth, vertex, NULL);
  if (is_adjacent_digraph(DIGRAPH1, last_def, vertex)) {
    intersect_bit_arrays(get_conditions(CONDITIONS, vertex),
                         DIGRAPH2->out_neighbours[MAP[last_def]],
                         DIGRAPH2->nr_vertices);
  } else {
    complement_bit_arrays(get_conditions(CONDITIONS, vertex),
                          DIGRAPH2->out_neighbours[MAP[last_def]],
                          DIGRAPH2->nr_vertices);
  }
  if (is_adjacent_digraph(DIGRAPH1, vertex, last_def)) {
    intersect_bit_arrays(get_conditions(CONDITIONS, vertex),
                         DIGRAPH2->in_neighbours[MAP[last_def]],
                         DIGRAPH2->nr_vertices);
  } else {
    complement_bit_arrays(get_conditions(CONDITIONS, vertex),
                          DIGRAPH2->in_neighbours[MAP[last_def]],
                          DIGRAPH2->nr_vertices);
  }
  store_size_conditions(CONDITIONS, vertex);
  return size_conditions(CONDITIONS, vertex);
}

// The main recursive function for embeddings of digraphs.
//
// The arguments are:
//
// 1. depth             The current depth of the search (i.e. number of
//                      positions in MAP that are assigned).
// 2. pos               The last position in MAP that was filled
// 3. rep_depth         The depth of the stabilisers of points in MAP already.
//                      This is different than depth because we stop increasing
//                      the rep_depth when we first encounter a trivial
//                      stabiliser.
// 4. has_trivial_stab  true if the stabiliser of the values already in MAP is
//                      trivial, false otherwise.
// 5. max_results       The maximum number of results to find.
// 6. count             The number of embeddings found so far.
static void find_digraph_embeddings(uint16_t        depth,
                                    uint16_t        pos,
                                    uint16_t        rep_depth,
                                    bool            has_trivial_stab,
                                    uint64_t const  max_results,
                                    uint64_t* const count) {
#ifdef DIGRAPHS_ENABLE_STATS
  update_stats(STATS, depth);
#endif
  if (depth == DIGRAPH1->nr_vertices) {
    // we've assigned every position in <MAP>
    external_order_map_digraph(DIGRAPH1);
    Obj ret = HOOK(
        USER_PARAM, MAX(DIGRAPH1->nr_vertices, DIGRAPH2->nr_vertices), MAP);
    internal_order_map_digraph(DIGRAPH1);
    (*count)++;
    if (*count >= max_results || ret == True) {
      longjmp(OUTOFHERE, 1);
    }
    return;
  }

  uint16_t next = 0;          // the next position to fill
  uint16_t min  = UNDEFINED;  // the minimum number of candidates for MAP[next]
  uint16_t i;

  if (depth > 0) {  // this is not the first call of the function
    copy_bit_array(
        MAP_UNDEFINED[depth], MAP_UNDEFINED[depth - 1], DIGRAPH1->nr_vertices);
    FOR_SET_BITS(MAP_UNDEFINED[depth], DIGRAPH1->nr_vertices, i) {
      size_t const n = digraph_embed_update_conditions(depth, pos, i);
      if (n == 0) {
#ifdef DIGRAPHS_ENABLE_STATS
        STATS->nr_dead_branches++;
#endif
        pop_conditions(CONDITIONS, depth);
        return;
      }
      STORE_MIN(min, next, n, i);
    }
    END_FOR_SET_BITS
  } else {
    for (i = 0; i < DIGRAPH1->nr_vertices; i++) {
      size_t const n = size_conditions(CONDITIONS, i);
      STORE_MIN_BREAK(min, next, n, i);
    }
  }

  BitArray* possible = get_conditions(CONDITIONS, next);
  copy_bit_array(
      possible, get_conditions(CONDITIONS, next), DIGRAPH2->nr_vertices);
  intersect_bit_arrays(possible, REPS[rep_depth], DIGRAPH2->nr_vertices);
  complement_bit_arrays(possible, VALS, DIGRAPH2->nr_vertices);

  FOR_SET_BITS(possible, DIGRAPH2->nr_vertices, i) {
    MAP[next] = i;
    set_bit_array(VALS, i, true);
    set_bit_array(MAP_UNDEFINED[depth], next, false);
    if (!has_trivial_stab) {
      find_digraph_embeddings(
          depth + 1,
          next,
          rep_depth + 1,
          compute_stabs_and_orbit_reps(DIGRAPH1->nr_vertices,
                                       DIGRAPH2->nr_vertices,
                                       rep_depth + 1,
                                       depth,
                                       i,
                                       false),
          max_results,
          count);
    } else {
      find_digraph_embeddings(
          depth + 1, next, rep_depth, true, max_results, count);
    }
    MAP[next] = UNDEFINED;
    set_bit_array(VALS, i, false);
    set_bit_array(MAP_UNDEFINED[depth], next, true);
  }
  END_FOR_SET_BITS
  pop_conditions(CONDITIONS, depth);
}

////////////////////////////////////////////////////////////////////////////////
// The next function should be called before find_digraph_homos.  This function
// simulates the first steps of the recursion using the information in
// partial_map_obj (if any) so that the values specified in partial_map_obj are
// installed in MAP first.
////////////////////////////////////////////////////////////////////////////////

static void init_partial_map_and_find_digraph_homos(Obj partial_map_obj,
                                                    uint64_t const  max_results,
                                                    uint64_t const  hint,
                                                    uint64_t* const count,
                                                    Obj injective_obj) {
  uint16_t depth                = 0;
  uint16_t rep_depth            = 0;
  uint16_t last_defined         = UNDEFINED;
  bool     last_stab_is_trivial = (STAB_GENS[0]->size == 0 ? true : false);
  uint16_t rank                 = 0;
  uint16_t next                 = UNDEFINED;

  if (partial_map_obj != Fail) {
    for (next = 0; next < LEN_LIST(partial_map_obj); ++next) {
      if (ISB_LIST(partial_map_obj, next + 1)) {
        if (depth > 0) {
          DIGRAPHS_ASSERT(last_defined != UNDEFINED);
          copy_bit_array(MAP_UNDEFINED[depth],
                         MAP_UNDEFINED[depth - 1],
                         DIGRAPH1->nr_vertices);
          uint16_t i;  // variable for the next for-loop
          FOR_SET_BITS(MAP_UNDEFINED[depth], DIGRAPH1->nr_vertices, i) {
            uint16_t n;
            if (INT_INTOBJ(injective_obj) == 0) {
              n = digraph_homo_update_conditions(depth, last_defined, i);
            } else if (INT_INTOBJ(injective_obj) == 1) {
              n = digraph_mono_update_conditions(depth, last_defined, i);
            } else {
              DIGRAPHS_ASSERT(INT_INTOBJ(injective_obj) == 2);
              n = digraph_embed_update_conditions(depth, last_defined, i);
            }
            if (n == 0) {
              return;
            }
          }
          END_FOR_SET_BITS
        }
        uint16_t val = INT_INTOBJ(ELM_LIST(partial_map_obj, next + 1)) - 1;
        next         = (ORDERED ? INVERSE_ORDER[next] : next);
        MAP[next]    = val;
        if (!get_bit_array(VALS, MAP[next])) {
          rank++;
          if (rank > hint) {
            return;
          }
        }
        set_bit_array(VALS, MAP[next], true);
        set_bit_array(MAP_UNDEFINED[depth], next, false);
        if (!last_stab_is_trivial) {
          last_stab_is_trivial =
              compute_stabs_and_orbit_reps(DIGRAPH1->nr_vertices,
                                           DIGRAPH2->nr_vertices,
                                           rep_depth + 1,
                                           depth,
                                           MAP[next],
                                           false);
          rep_depth++;
        }
        depth++;
        last_defined = next;
        next         = (ORDERED ? ORDER[next] : next);
      }
    }
  }
  if (INT_INTOBJ(injective_obj) == 0) {
    find_digraph_homos(depth,
                       last_defined,
                       rep_depth,
                       last_stab_is_trivial,
                       rank,
                       max_results,
                       hint,
                       count);
  } else if (INT_INTOBJ(injective_obj) == 1) {
    find_digraph_monos(depth,
                       last_defined,
                       rep_depth,
                       last_stab_is_trivial,
                       max_results,
                       count);
  } else {
    DIGRAPHS_ASSERT(INT_INTOBJ(injective_obj) == 2);
    find_digraph_embeddings(depth,
                            last_defined,
                            rep_depth,
                            last_stab_is_trivial,
                            max_results,
                            count);
  }
  return;
}

////////////////////////////////////////////////////////////////////////////////
// 7. The GAP-level function (and helpers)
////////////////////////////////////////////////////////////////////////////////

// Initialises the data structures required by the recursive functions for
// finding homomorphisms. If true is returned everything was initialised ok, if
// false is returned, then the arguments already imply that there can be no
// homomorphisms.
static bool init_data_from_args(Obj digraph1_obj,
                                Obj digraph2_obj,
                                Obj hook_obj,
                                Obj user_param_obj,
                                Obj max_results_obj,
                                Obj hint_obj,
                                Obj injective_obj,
                                Obj image_obj,
                                Obj partial_map_obj,
                                Obj colors1_obj,
                                Obj colors2_obj,
                                Obj order_obj,
                                Obj aut_grp_obj) {
  uint16_t calculated_max_verts =
      MAX(DigraphNrVertices(digraph1_obj), DigraphNrVertices(digraph2_obj));
  if (!homos_data_initialized
      || (calculated_max_verts >= HOMOS_STRUCTURE_SIZE)) {
    FuncDIGRAPHS_FREE_HOMOS_DATA(0L);
    homos_data_initialized = true;

    // Rather arbitrary, but we multiply by 1.2 to avoid
    // n = 1,2,3,4,5... causing constant reallocation
    HOMOS_STRUCTURE_SIZE =
        (calculated_max_verts + calculated_max_verts / 5) + 1;
    // The previous line includes "+ 1" because below we do:
    // "BLISS_GRAPH[3 * PERM_DEGREE]" but we only allocate bliss graphs in
    // BLISS_GRAPH up to but not including 3 * HOMOS_STRUCTURE_SIZE. So if
    // PERM_DEGREE = (# nodes in digraph2) = HOMOS_STRUCTURE_SIZE (the
    // first equality always holds, the second if  # nodes in digraph2 >
    // # nodes in digraph1), then this is out of bounds.
#ifdef DIGRAPHS_ENABLE_STATS
    STATS = safe_malloc(sizeof(HomoStats));
#endif
    DIGRAPH1 = new_digraph(HOMOS_STRUCTURE_SIZE);
    DIGRAPH2 = new_digraph(HOMOS_STRUCTURE_SIZE);

    GRAPH1 = new_graph(HOMOS_STRUCTURE_SIZE);
    GRAPH2 = new_graph(HOMOS_STRUCTURE_SIZE);

    IMAGE_RESTRICT = new_bit_array(HOMOS_STRUCTURE_SIZE);
    ORB_LOOKUP     = new_bit_array(HOMOS_STRUCTURE_SIZE);
    REPS           = safe_malloc(HOMOS_STRUCTURE_SIZE * sizeof(BitArray*));
    BIT_ARRAY_BUFFER =
        (BitArray**) safe_calloc(HOMOS_STRUCTURE_SIZE, sizeof(BitArray*));
    MAP_UNDEFINED =
        (BitArray**) safe_calloc(HOMOS_STRUCTURE_SIZE, sizeof(BitArray*));
    BLISS_GRAPH = (BlissGraph**) safe_calloc(3 * HOMOS_STRUCTURE_SIZE,
                                             sizeof(BlissGraph*));
    MAP     = (uint16_t*) safe_calloc(HOMOS_STRUCTURE_SIZE, sizeof(uint16_t));
    COLORS2 = (uint16_t*) safe_calloc(HOMOS_STRUCTURE_SIZE, sizeof(uint16_t));
    INVERSE_ORDER =
        (uint16_t*) safe_calloc(HOMOS_STRUCTURE_SIZE, sizeof(uint16_t));
    MAP_BUFFER =
        (uint16_t*) safe_calloc(HOMOS_STRUCTURE_SIZE, sizeof(uint16_t));
    ORB   = (uint16_t*) safe_calloc(HOMOS_STRUCTURE_SIZE, sizeof(uint16_t));
    ORDER = (uint16_t*) safe_calloc(HOMOS_STRUCTURE_SIZE, sizeof(uint16_t));
    STAB_GENS =
        (PermColl**) safe_calloc(HOMOS_STRUCTURE_SIZE, sizeof(PermColl*));

    for (uint16_t i = 0; i < HOMOS_STRUCTURE_SIZE * 3; i++) {
      BLISS_GRAPH[i] = bliss_digraphs_new(i);
    }

    for (uint16_t i = 0; i < HOMOS_STRUCTURE_SIZE; i++) {
      REPS[i]             = new_bit_array(HOMOS_STRUCTURE_SIZE);
      BIT_ARRAY_BUFFER[i] = new_bit_array(HOMOS_STRUCTURE_SIZE);
      MAP_UNDEFINED[i]    = new_bit_array(HOMOS_STRUCTURE_SIZE);
      STAB_GENS[i] = new_perm_coll(HOMOS_STRUCTURE_SIZE, HOMOS_STRUCTURE_SIZE);
    }
    VALS       = new_bit_array(HOMOS_STRUCTURE_SIZE);
    CONDITIONS = new_conditions(HOMOS_STRUCTURE_SIZE, HOMOS_STRUCTURE_SIZE);

    SCHREIER_SIMS = new_schreier_sims();
  }
#ifdef DIGRAPHS_ENABLE_STATS
  clear_stats(STATS);
#endif

  uint16_t nr1 = DigraphNrVertices(digraph1_obj);
  uint16_t nr2 = DigraphNrVertices(digraph2_obj);

  init_bit_array(MAP_UNDEFINED[0], true, nr1);
  init_bit_array(VALS, false, nr2);

  if (IS_LIST(order_obj)) {
    ORDERED = true;
    for (uint16_t i = 0; i < nr1; i++) {
      ORDER[i]                = INT_INTOBJ(ELM_LIST(order_obj, i + 1)) - 1;
      INVERSE_ORDER[ORDER[i]] = i;
    }
  } else {
    ORDERED = false;
  }

  bool is_undirected;
  if (CALL_1ARGS(IsSymmetricDigraph, digraph1_obj) == True
      && CALL_1ARGS(IsSymmetricDigraph, digraph2_obj) == True) {
    init_graph_from_digraph_obj(GRAPH1, digraph1_obj, ORDERED);
--> --------------------

--> maximum size reached

--> --------------------