1 //
2 // homos.c               di/graph homomorphisms              Julius Jonusas
3 //                                                           J. D. Mitchell
4 //
5 // Copyright (C) 2014-19 - Julius Jonusas and J. D. Mitchell
6 //
7 // This file is free software, see the digraphs/LICENSE.
8 
9 ////////////////////////////////////////////////////////////////////////////////
10 //
11 // This file is organised as follows:
12 //
13 // 1. Macros
14 // 2. Forward declarations
15 // 3. Global variables
16 // 4. Hook functions
17 // 5. Static helper functions
18 // 6. The main recursive functions (and helpers)
19 // 7. The GAP-level function (and helpers)
20 //
21 ////////////////////////////////////////////////////////////////////////////////
22 
23 // TODO(later)
24 // 0. remove final arg from push_conditions
25 // 1. Try other bit hacks for iterating through set bits
26 
27 #include "homos.h"
28 
29 // C headers
30 #include <limits.h>   // for CHAR_BIT
31 #include <setjmp.h>   // for longjmp, setjmp, jmp_buf
32 #include <stdbool.h>  // for true, false, bool
33 #include <stddef.h>   // for NULL
34 #include <stdint.h>   // for uint16_t, uint64_t
35 #include <stdlib.h>   // for malloc, NULL
36 #include <time.h>     // for time
37 
38 // GAP headers
39 #include "src/compiled.h"
40 
41 // Digraphs package headers
42 #include "bitarray.h"         // for BitArray
43 #include "conditions.h"       // for Conditions
44 #include "digraphs-config.h"  // for DIGRAPHS_HAVE___BUILTIN_CTZLL
45 #include "digraphs-debug.h"   // for DIGRAPHS_ASSERT
46 #include "homos-graphs.h"     // for Digraph, Graph, . . .
47 #include "perms.h"            // for MAXVERTS, UNDEFINED, PermColl, Perm
48 #include "schreier-sims.h"    // for PermColl, . . .
49 
50 #ifdef DIGRAPHS_WITH_INCLUDED_BLISS
51 #include "bliss-0.73/bliss_C.h"  // for bliss_digraphs_release, . . .
52 #else
53 #include "bliss/bliss_C.h"
54 #define bliss_digraphs_add_edge bliss_add_edge
55 #define bliss_digraphs_new bliss_new
56 #define bliss_digraphs_add_vertex bliss_add_vertex
57 #define bliss_digraphs_find_canonical_labeling bliss_find_canonical_labeling
58 #define bliss_digraphs_release bliss_release
59 #endif
60 
61 ////////////////////////////////////////////////////////////////////////////////
62 // 1. Macros
63 ////////////////////////////////////////////////////////////////////////////////
64 
65 #define MAX(a, b) (a < b ? b : a)
66 #define MIN(a, b) (a < b ? a : b)
67 
68 #define NUMBER_BITS_PER_BLOCK (sizeof(Block) * CHAR_BIT)
69 
70 // The next line can be used instead of the first line of STORE_MIN to
71 // randomise which vertex of minimum degree is used next, but I didn't find any
72 // cases where this made things faster.
73 
74 // if (x < current_x || (x == current_x && (rand() & 1))) {
75 
76 #define STORE_MIN(current_x, current_y, x, y) \
77   if (x < current_x) {                        \
78     current_x = x;                            \
79     current_y = y;                            \
80   }
81 
82 // See the comment before STORE_MIN
83 
84 // if (x < current_x || (x == current_x && (rand() & 1))) {
85 
86 #define STORE_MIN_BREAK(current_x, current_y, x, y) \
87   if (x < current_x) {                              \
88     current_x = x;                                  \
89     current_y = y;                                  \
90     if (current_x == 1) {                           \
91       break;                                        \
92     }                                               \
93   }
94 
95 #if SIZEOF_VOID_P == 4 && DIGRAPHS_HAVE___BUILTIN_CTZLL
96 #undef DIGRAPHS_HAVE___BUILTIN_CTZLL
97 #endif
98 
99 // The following macro is bitmap_decode_ctz_callpack from https://git.io/fho4p
100 #ifdef DIGRAPHS_HAVE___BUILTIN_CTZLL
101 #define FOR_SET_BITS(__bit_array, __nr_bits, __variable)     \
102   for (size_t k = 0; k < NR_BLOCKS_LOOKUP[__nr_bits]; ++k) { \
103     Block block = __bit_array->blocks[k];                    \
104     while (block != 0) {                                     \
105       uint64_t t = block & -block;                           \
106       int      r = __builtin_ctzll(block);                   \
107       __variable = k * 64 + r;                               \
108       if (__variable >= __nr_bits) {                         \
109         break;                                               \
110       }
111 
112 #define END_FOR_SET_BITS \
113   block ^= t;            \
114   }                      \
115   }
116 
117 #else
118 #define FOR_SET_BITS(__bit_array, __nr_bits, __variable)       \
119   for (__variable = 0; __variable < __nr_bits; ++__variable) { \
120     if (get_bit_array(__bit_array, __variable)) {
121 #define END_FOR_SET_BITS \
122   }                      \
123   }
124 #endif
125 
126 ////////////////////////////////////////////////////////////////////////////////
127 // 2. Forward declarations
128 ////////////////////////////////////////////////////////////////////////////////
129 
130 // Defined in digraphs.h
131 Int DigraphNrVertices(Obj);
132 Obj FuncOutNeighbours(Obj, Obj);
133 
134 // GAP level things, imported in digraphs.c
135 extern Obj IsDigraph;
136 extern Obj DIGRAPHS_ValidateVertexColouring;
137 extern Obj Infinity;
138 extern Obj IsSymmetricDigraph;
139 extern Obj GeneratorsOfGroup;
140 extern Obj AutomorphismGroup;
141 extern Obj IsPermGroup;
142 extern Obj IsDigraphAutomorphism;
143 extern Obj LargestMovedPointPerms;
144 
145 ////////////////////////////////////////////////////////////////////////////////
146 // 3. Global variables
147 ////////////////////////////////////////////////////////////////////////////////
148 
149 static Obj GAP_FUNC;  // Variable to hold a GAP level hook function
150 
151 static Obj (*HOOK)(void*,  // HOOK function applied to every homo found
152                    const uint16_t,
153                    const uint16_t*);
154 static void* USER_PARAM;  // a USER_PARAM for the hook
155 
156 // Values in MAP are restricted to those positions in IMAGE_RESTRICT
157 static jmp_buf OUTOFHERE;  // so we can jump out of the deepest
158 
159 static bool ORDERED;  // true if the vertices of the domain/source digraph
160                       // should be considered in a different order than they are
161                       // given, false otherwise.
162 
163 static BitArray* BIT_ARRAY_BUFFER[MAXVERTS];  // A buffer
164 static BitArray* IMAGE_RESTRICT;              // Values in MAP must be in this
165 static BitArray* MAP_UNDEFINED[MAXVERTS];     // Undefined positions in MAP
166 static BitArray* ORB_LOOKUP;                  // points in orbit
167 static BitArray* VALS;                        // Values in MAP already
168 
169 static BitArray** REPS;  // orbit reps organised by depth
170 
171 static Conditions* CONDITIONS;
172 
173 static Digraph* DIGRAPH1;  // Digraphs to hold incoming GAP digraphs
174 static Digraph* DIGRAPH2;
175 
176 static Graph* GRAPH1;  // Graphs to hold incoming GAP symmetric digraphs
177 static Graph* GRAPH2;
178 
179 static BlissGraph* BLISS_GRAPH[3 * MAXVERTS];
180 
181 static uint16_t MAP[MAXVERTS];            // partial image list
182 static uint16_t COLORS2[MAXVERTS];        // colors of range (di)graph
183 static uint16_t INVERSE_ORDER[MAXVERTS];  // external -> internal
184 static uint16_t MAP_BUFFER[MAXVERTS];     // For converting from internal ->
185                                           // external and back when calling the
186                                           // hook functions.
187 static uint16_t ORB[MAXVERTS];    // Array for containing nodes in an orbit.
188 static uint16_t ORDER[MAXVERTS];  // internal -> external
189 
190 static PermColl*     STAB_GENS[MAXVERTS];  // stabiliser generators
191 static SchreierSims* SCHREIER_SIMS;
192 
193 #ifdef DIGRAPHS_ENABLE_STATS
194 struct homo_stats_struct {
195   time_t last_print;
196   size_t max_depth;
197   double mean_depth;
198   size_t nr_calls;
199   size_t nr_dead_branches;
200   time_t start_time;
201 };
202 
203 typedef struct homo_stats_struct HomoStats;
204 
205 static HomoStats* STATS;
206 
clear_stats(HomoStats * stats)207 static inline void clear_stats(HomoStats* stats) {
208   stats->last_print       = time(0);
209   stats->max_depth        = 0;
210   stats->mean_depth       = 0;
211   stats->nr_calls         = 0;
212   stats->nr_dead_branches = 0;
213   stats->start_time       = time(0);
214 }
215 
print_stats(HomoStats * stats)216 static void print_stats(HomoStats* stats) {
217   printf("Running for %0.0fs . . .\n", difftime(time(0), stats->start_time));
218   printf("Number of function calls = %*lu\n", 20, stats->nr_calls);
219   printf("Mean depth               = %*.2f\n", 20, stats->mean_depth);
220   printf("Max depth                = %*lu\n", 20, stats->max_depth);
221   printf("Number of dead branches  = %*lu\n\n", 20, stats->nr_dead_branches);
222 }
223 
update_stats(HomoStats * stats,uint64_t depth)224 static inline void update_stats(HomoStats* stats, uint64_t depth) {
225   if (depth > stats->max_depth) {
226     stats->max_depth = depth;
227   }
228   stats->nr_calls++;
229   stats->mean_depth += (depth - stats->mean_depth) / stats->nr_calls;
230   if (difftime(time(0), stats->last_print) > 0.9) {
231     print_stats(stats);
232     stats->last_print = time(0);
233   }
234 }
235 #endif  // DIGRAPHS_ENABLE_STATS
236 
237 ////////////////////////////////////////////////////////////////////////////////
238 // 4. Hook functions
239 ////////////////////////////////////////////////////////////////////////////////
240 
241 static Obj
homo_hook_gap(void * user_param,uint16_t const nr,uint16_t const * map)242 homo_hook_gap(void* user_param, uint16_t const nr, uint16_t const* map) {
243   UInt2* ptr;
244   Obj    t;
245   UInt   i;
246 
247   // copy map into new trans2
248   t   = NEW_TRANS2(nr);
249   ptr = ADDR_TRANS2(t);
250 
251   for (i = 0; i < nr; i++) {
252     ptr[i] = map[i];
253   }
254   return CALL_2ARGS(GAP_FUNC, user_param, t);
255 }
256 
257 static Obj
homo_hook_collect(void * user_param,uint16_t const nr,uint16_t const * map)258 homo_hook_collect(void* user_param, uint16_t const nr, uint16_t const* map) {
259   UInt2* ptr;
260   Obj    t;
261   UInt   i;
262 
263   if (TNUM_OBJ((Obj) user_param) == T_PLIST_EMPTY) {
264     RetypeBag(user_param, T_PLIST);
265   }
266 
267   // copy map into new trans2
268   t   = NEW_TRANS2(nr);
269   ptr = ADDR_TRANS2(t);
270 
271   for (i = 0; i < nr; i++) {
272     ptr[i] = map[i];
273   }
274 
275   ASS_LIST(user_param, LEN_LIST(user_param) + 1, t);
276   return False;
277 }
278 
279 ////////////////////////////////////////////////////////////////////////////////
280 // 5. Static helper functions
281 ////////////////////////////////////////////////////////////////////////////////
282 
283 // print_array is not used, but can be useful for debugging.
284 
285 // static void print_array(uint16_t const* const array, uint16_t const len) {
286 //   if (array == NULL) {
287 //     printf("NULL");
288 //     return;
289 //   }
290 //   printf("<array {");
291 //   for (uint16_t i = 0; i < len; i++) {
292 //     printf(" %d", array[i]);
293 //   }
294 //   printf(" }>");
295 // }
296 
get_automorphism_group_from_gap(Obj digraph_obj,PermColl * out)297 static void get_automorphism_group_from_gap(Obj digraph_obj, PermColl* out) {
298   DIGRAPHS_ASSERT(CALL_1ARGS(IsDigraph, digraph_obj) == True);
299   Obj o = CALL_1ARGS(AutomorphismGroup, digraph_obj);
300   o     = CALL_1ARGS(GeneratorsOfGroup, o);
301   DIGRAPHS_ASSERT(IS_LIST(o));
302   clear_perm_coll(out);
303   out->degree = PERM_DEGREE;
304   DIGRAPHS_ASSERT(out->capacity >= LEN_LIST(o));
305   for (Int i = 1; i <= LEN_LIST(o); ++i) {
306     DIGRAPHS_ASSERT(ISB_LIST(o, i));
307     DIGRAPHS_ASSERT(IS_PERM2(ELM_LIST(o, i)) || IS_PERM4(ELM_LIST(o, i)));
308     Obj p = ELM_LIST(o, i);
309     DIGRAPHS_ASSERT(LargestMovedPointPerm(p) <= PERM_DEGREE);
310     if (LargestMovedPointPerm(p) == 0) {
311       continue;
312     }
313     Perm q = out->perms[i - 1];
314     out->size++;
315     DIGRAPHS_ASSERT(out->size == i);
316     size_t dep;
317     if (IS_PERM2(p)) {
318       UInt2 const* ptp2 = CONST_ADDR_PERM2(p);
319       dep               = MIN((uint16_t) DEG_PERM2(p), PERM_DEGREE);
320       for (uint16_t j = 0; j < dep; ++j) {
321         q[j] = (uint16_t) ptp2[j];
322       }
323     } else {
324       DIGRAPHS_ASSERT(IS_PERM4(p));
325       UInt4 const* ptp4 = CONST_ADDR_PERM4(p);
326       dep               = MIN((uint16_t) DEG_PERM4(p), PERM_DEGREE);
327       for (uint16_t j = 0; j < dep; ++j) {
328         q[j] = (uint16_t) ptp4[j];
329       }
330     }
331     for (uint16_t j = dep; j < PERM_DEGREE; ++j) {
332       q[j] = j;
333     }
334   }
335 }
336 
init_digraph_from_digraph_obj(Digraph * const digraph,Obj digraph_obj,bool const reorder)337 static void init_digraph_from_digraph_obj(Digraph* const digraph,
338                                           Obj            digraph_obj,
339                                           bool const     reorder) {
340   DIGRAPHS_ASSERT(digraph != NULL);
341   DIGRAPHS_ASSERT(CALL_1ARGS(IsDigraph, digraph_obj) == True);
342   UInt const nr  = DigraphNrVertices(digraph_obj);
343   Obj        out = FuncOutNeighbours(0L, digraph_obj);
344   DIGRAPHS_ASSERT(nr < MAXVERTS);
345   DIGRAPHS_ASSERT(IS_PLIST(out));
346   clear_digraph(digraph, nr);
347 
348   if (!reorder) {
349     for (uint16_t i = 1; i <= nr; i++) {
350       Obj nbs = ELM_PLIST(out, i);
351       DIGRAPHS_ASSERT(IS_LIST(nbs));
352       for (uint16_t j = 1; j <= LEN_LIST(nbs); j++) {
353         DIGRAPHS_ASSERT(IS_INTOBJ(ELM_LIST(nbs, j)));
354         add_edge_digraph(digraph, i - 1, INT_INTOBJ(ELM_LIST(nbs, j)) - 1);
355       }
356     }
357   } else {
358     DIGRAPHS_ASSERT(ORDERED);
359     for (uint16_t i = 1; i <= nr; i++) {
360       Obj nbs = ELM_PLIST(out, ORDER[i - 1] + 1);
361       DIGRAPHS_ASSERT(IS_LIST(nbs));
362       for (uint16_t j = 1; j <= LEN_LIST(nbs); j++) {
363         DIGRAPHS_ASSERT(IS_INTOBJ(ELM_LIST(nbs, j)));
364         add_edge_digraph(
365             digraph, i - 1, INVERSE_ORDER[INT_INTOBJ(ELM_LIST(nbs, j)) - 1]);
366       }
367     }
368   }
369 }
370 
init_graph_from_digraph_obj(Graph * const graph,Obj digraph_obj,bool const reorder)371 static void init_graph_from_digraph_obj(Graph* const graph,
372                                         Obj          digraph_obj,
373                                         bool const   reorder) {
374   DIGRAPHS_ASSERT(graph != NULL);
375   DIGRAPHS_ASSERT(CALL_1ARGS(IsDigraph, digraph_obj) == True);
376   DIGRAPHS_ASSERT(CALL_1ARGS(IsSymmetricDigraph, digraph_obj) == True);
377   UInt const nr  = DigraphNrVertices(digraph_obj);
378   Obj        out = FuncOutNeighbours(0L, digraph_obj);
379   DIGRAPHS_ASSERT(nr < MAXVERTS);
380   DIGRAPHS_ASSERT(IS_PLIST(out));
381   clear_graph(graph, nr);
382 
383   if (!reorder) {
384     for (uint16_t i = 1; i <= nr; i++) {
385       Obj nbs = ELM_PLIST(out, i);
386       DIGRAPHS_ASSERT(IS_LIST(nbs));
387       for (uint16_t j = 1; j <= LEN_LIST(nbs); j++) {
388         DIGRAPHS_ASSERT(IS_INTOBJ(ELM_LIST(nbs, j)));
389         add_edge_graph(graph, i - 1, INT_INTOBJ(ELM_LIST(nbs, j)) - 1);
390       }
391     }
392   } else {
393     DIGRAPHS_ASSERT(ORDERED);
394     for (uint16_t i = 1; i <= nr; i++) {  // Nodes in the new graph
395       Obj nbs = ELM_PLIST(out, ORDER[i - 1] + 1);
396       DIGRAPHS_ASSERT(IS_LIST(nbs));
397       for (uint16_t j = 1; j <= LEN_LIST(nbs); j++) {
398         DIGRAPHS_ASSERT(IS_INTOBJ(ELM_LIST(nbs, j)));
399         add_edge_graph(
400             graph, i - 1, INVERSE_ORDER[INT_INTOBJ(ELM_LIST(nbs, j)) - 1]);
401       }
402     }
403   }
404 }
405 
406 // Find orbit representatives of the group generated by STAB_GENS[rep_depth]
407 // and store them in REPS[rep_depth], only values in IMAGE_RESTRICT will be
408 // chosen as orbit representatives.
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)409 static bool compute_stabs_and_orbit_reps(uint16_t const nr_nodes_1,
410                                          uint16_t const nr_nodes_2,
411                                          uint16_t const rep_depth,
412                                          uint16_t const depth,
413                                          uint16_t const pt,
414                                          bool const     first_call) {
415   DIGRAPHS_ASSERT(rep_depth <= depth + 1);
416   if (depth == nr_nodes_1 - 1 && !first_call) {
417     // first_call is required in the case that nr_nodes_1 is 1, since without
418     // first_call this function returns false (in the next line) and REPS is
419     // not initialised.
420     return false;  // This doesn't really say anything about the stabiliser
421   } else if (rep_depth > 0) {
422     point_stabilizer(
423         SCHREIER_SIMS, STAB_GENS[rep_depth - 1], STAB_GENS[rep_depth], pt);
424     if (STAB_GENS[rep_depth]->size == 0) {
425       // the stabiliser of pt in STAB_GENS[rep_depth - 1] is trivial
426       copy_bit_array(REPS[rep_depth], IMAGE_RESTRICT, nr_nodes_2);
427       complement_bit_arrays(REPS[rep_depth], VALS, nr_nodes_2);
428       // REPS[rep_depth] is a set of orbit representatives of the stabiliser of
429       // the existing values in MAP, which belong to VALS, and since the
430       // stabiliser is trivial, we take valid choice as the set of
431       // representatives.
432       return true;  // the stabiliser is trivial
433     }
434   }
435   init_bit_array(REPS[rep_depth], false, nr_nodes_2);
436   copy_bit_array(ORB_LOOKUP, VALS, nr_nodes_2);
437   uint16_t fst = 0;
438   while (fst < PERM_DEGREE
439          && (get_bit_array(ORB_LOOKUP, fst)
440              || !get_bit_array(IMAGE_RESTRICT, fst))) {
441     fst++;
442   }
443 
444   while (fst < PERM_DEGREE) {
445     ORB[0]     = fst;
446     uint16_t n = 1;  // length of ORB
447 
448     set_bit_array(REPS[rep_depth], fst, true);
449     set_bit_array(ORB_LOOKUP, fst, true);
450 
451     for (uint16_t i = 0; i < n; ++i) {
452       for (uint16_t j = 0; j < STAB_GENS[rep_depth]->size; ++j) {
453         Perm           gen = STAB_GENS[rep_depth]->perms[j];
454         uint16_t const img = gen[ORB[i]];
455         if (!get_bit_array(ORB_LOOKUP, img)) {
456           ORB[n++] = img;
457           set_bit_array(ORB_LOOKUP, img, true);
458         }
459       }
460     }
461     while (fst < PERM_DEGREE
462            && (get_bit_array(ORB_LOOKUP, fst)
463                || !get_bit_array(IMAGE_RESTRICT, fst))) {
464       fst++;
465     }
466   }
467   return false;  // the stabiliser is not trivial
468 }
469 
470 // Rewrite the global variable MAP so that it contains a homomorphism of the
471 // original GAP level (di)graph, and not the possibly distinct (but isomorphic)
472 // copy in the homomorphism search.  This should be called before any calls to
473 // the hook functions (i.e. after an entire homomorphism is found).
external_order_map_digraph(Digraph * digraph)474 static void external_order_map_digraph(Digraph* digraph) {
475   if (!ORDERED) {
476     return;
477   }
478   for (uint16_t i = 0; i < digraph->nr_vertices; ++i) {
479     MAP_BUFFER[ORDER[i]] = MAP[i];
480   }
481   for (uint16_t i = 0; i < digraph->nr_vertices; ++i) {
482     MAP[i] = MAP_BUFFER[i];
483   }
484 }
485 
external_order_map_graph(Graph * graph)486 static void external_order_map_graph(Graph* graph) {
487   if (!ORDERED) {
488     return;
489   }
490   for (uint16_t i = 0; i < graph->nr_vertices; ++i) {
491     MAP_BUFFER[ORDER[i]] = MAP[i];
492   }
493   for (uint16_t i = 0; i < graph->nr_vertices; ++i) {
494     MAP[i] = MAP_BUFFER[i];
495   }
496 }
497 
498 // Rewrite the global variable MAP so that it contains a homomorphism of the
499 // internal (di)graph, and not the possibly distinct (but isomorphic) GAP level
500 // (di)graph.  This should be called after any calls to the hook functions
501 // (i.e. after an entire homomorphism is found).
internal_order_map_digraph(Digraph const * const digraph)502 static void internal_order_map_digraph(Digraph const* const digraph) {
503   if (!ORDERED) {
504     return;
505   }
506   for (uint16_t i = 0; i < digraph->nr_vertices; ++i) {
507     MAP_BUFFER[INVERSE_ORDER[i]] = MAP[i];
508   }
509   for (uint16_t i = 0; i < digraph->nr_vertices; ++i) {
510     MAP[i] = MAP_BUFFER[i];
511   }
512 }
513 
internal_order_map_graph(Graph const * const graph)514 static void internal_order_map_graph(Graph const* const graph) {
515   if (!ORDERED) {
516     return;
517   }
518   for (uint16_t i = 0; i < graph->nr_vertices; ++i) {
519     MAP_BUFFER[INVERSE_ORDER[i]] = MAP[i];
520   }
521   for (uint16_t i = 0; i < graph->nr_vertices; ++i) {
522     MAP[i] = MAP_BUFFER[i];
523   }
524 }
525 
set_automorphisms(Obj aut_grp_obj,PermColl * out)526 static void set_automorphisms(Obj aut_grp_obj, PermColl* out) {
527   DIGRAPHS_ASSERT(out != NULL);
528   clear_perm_coll(out);
529   out->degree = PERM_DEGREE;
530   Obj gens    = CALL_1ARGS(GeneratorsOfGroup, aut_grp_obj);
531   DIGRAPHS_ASSERT(IS_LIST(gens));
532   DIGRAPHS_ASSERT(LEN_LIST(gens) > 0);
533   for (UInt i = 1; i <= LEN_LIST(gens); ++i) {
534     Obj gen_obj = ELM_LIST(gens, i);
535     if (LargestMovedPointPerm(gen_obj) > 0) {
536       add_perm_coll(out, new_perm_from_gap(gen_obj, PERM_DEGREE));
537     }
538   }
539 }
540 
541 ////////////////////////////////////////////////////////////////////////////////
542 // 6. The main recursive functions (and helpers)
543 ////////////////////////////////////////////////////////////////////////////////
544 
545 // Helper for the main recursive homomorphism function.
546 static ALWAYS_INLINE uint16_t
graph_homo_update_conditions(uint16_t const depth,uint16_t const last_defined,uint16_t const vertex)547 graph_homo_update_conditions(uint16_t const depth,
548                              uint16_t const last_defined,
549                              uint16_t const vertex) {
550   push_conditions(
551       CONDITIONS, depth, vertex, GRAPH2->neighbours[MAP[last_defined]]);
552   store_size_conditions(CONDITIONS, vertex);
553   return size_conditions(CONDITIONS, vertex);
554 }
555 
556 // The main recursive function for homomorphisms of graphs.
557 //
558 // The arguments are:
559 //
560 // 1. depth             The current depth of the search (i.e. number of
561 //                      positions in MAP that are assigned).
562 // 2. pos               The last position in MAP that was filled
563 // 3. rep_depth         The depth of the stabilisers of points in MAP already.
564 //                      This is different than depth because we stop increasing
565 //                      the rep_depth when we first encounter a trivial
566 //                      stabiliser.
567 // 4. has_trivial_stab  true if the stabiliser of the values already in MAP is
568 //                      trivial, false otherwise.
569 // 5. rank              The current number of distinct values in MAP.
570 // 6. max_results       The maximum number of results to find.
571 // 7. hint              The desired number of distinct points in a (full)
572 //                      homomorphism.
573 // 8. count             The number of homomorphisms found so far.
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)574 static void find_graph_homos(uint16_t        depth,
575                              uint16_t        pos,
576                              uint16_t        rep_depth,
577                              bool            has_trivial_stab,
578                              uint16_t        rank,
579                              uint64_t const  max_results,
580                              uint64_t const  hint,
581                              uint64_t* const count) {
582 #ifdef DIGRAPHS_ENABLE_STATS
583   update_stats(STATS, depth);
584 #endif
585   if (depth == GRAPH1->nr_vertices) {
586     // Every position in MAP is assigned . . .
587     if (hint != UNDEFINED && rank != hint) {
588 #ifdef DIGRAPHS_ENABLE_STATS
589       STATS->nr_dead_branches++;
590 #endif
591       return;
592     }
593     external_order_map_graph(GRAPH1);
594     Obj ret =
595         HOOK(USER_PARAM, MAX(GRAPH1->nr_vertices, GRAPH2->nr_vertices), MAP);
596     internal_order_map_graph(GRAPH1);
597     (*count)++;
598     if (*count >= max_results || ret == True) {
599       longjmp(OUTOFHERE, 1);
600     }
601     return;
602   }
603 
604   uint16_t next = 0;          // the next position to fill
605   uint16_t min  = UNDEFINED;  // the minimum number of candidates for MAP[next]
606   uint16_t i;
607 
608   BitArray* possible = BIT_ARRAY_BUFFER[depth];
609 
610   if (depth > 0) {  // this is not the first call of the function
611     copy_bit_array(
612         MAP_UNDEFINED[depth], MAP_UNDEFINED[depth - 1], GRAPH1->nr_vertices);
613     copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
614     intersect_bit_arrays(
615         possible, GRAPH1->neighbours[pos], GRAPH1->nr_vertices);
616     FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
617       size_t const n = graph_homo_update_conditions(depth, pos, i);
618       if (n == 0) {
619 #ifdef DIGRAPHS_ENABLE_STATS
620         STATS->nr_dead_branches++;
621 #endif
622         pop_conditions(CONDITIONS, depth);
623         return;
624       }
625       STORE_MIN(min, next, n, i);
626     }
627     END_FOR_SET_BITS
628     if (min > 1) {
629       copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
630       complement_bit_arrays(
631           possible, GRAPH1->neighbours[pos], GRAPH1->nr_vertices);
632       FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
633         size_t const n = size_conditions(CONDITIONS, i);
634         STORE_MIN_BREAK(min, next, n, i);
635       }
636       END_FOR_SET_BITS
637     }
638   } else {
639     for (i = 0; i < GRAPH1->nr_vertices; ++i) {
640       size_t const n = size_conditions(CONDITIONS, i);
641       STORE_MIN_BREAK(min, next, n, i);
642     }
643   }
644   DIGRAPHS_ASSERT(get_bit_array(MAP_UNDEFINED[depth], next));
645 
646   if (rank < hint) {
647     copy_bit_array(
648         possible, get_conditions(CONDITIONS, next), GRAPH2->nr_vertices);
649     complement_bit_arrays(possible, VALS, GRAPH2->nr_vertices);
650     intersect_bit_arrays(possible, REPS[rep_depth], GRAPH2->nr_vertices);
651     FOR_SET_BITS(possible, GRAPH2->nr_vertices, i) {
652       MAP[next] = i;
653       set_bit_array(VALS, i, true);
654       set_bit_array(MAP_UNDEFINED[depth], next, false);
655       if (!has_trivial_stab) {
656         find_graph_homos(depth + 1,
657                          next,
658                          rep_depth + 1,
659                          compute_stabs_and_orbit_reps(GRAPH1->nr_vertices,
660                                                       GRAPH2->nr_vertices,
661                                                       rep_depth + 1,
662                                                       depth,
663                                                       i,
664                                                       false),
665                          rank + 1,
666                          max_results,
667                          hint,
668                          count);
669       } else {
670         find_graph_homos(depth + 1,
671                          next,
672                          rep_depth,
673                          true,
674                          rank + 1,
675                          max_results,
676                          hint,
677                          count);
678       }
679       MAP[next] = UNDEFINED;
680       set_bit_array(VALS, i, false);
681       set_bit_array(MAP_UNDEFINED[depth], next, true);
682     }
683     END_FOR_SET_BITS
684   }
685   copy_bit_array(
686       possible, get_conditions(CONDITIONS, next), GRAPH2->nr_vertices);
687   intersect_bit_arrays(possible, VALS, GRAPH2->nr_vertices);
688   FOR_SET_BITS(possible, GRAPH2->nr_vertices, i) {
689     MAP[next] = i;
690     set_bit_array(MAP_UNDEFINED[depth], next, false);
691     find_graph_homos(depth + 1,
692                      next,
693                      rep_depth,
694                      has_trivial_stab,
695                      rank,
696                      max_results,
697                      hint,
698                      count);
699     MAP[next] = UNDEFINED;
700     set_bit_array(MAP_UNDEFINED[depth], next, true);
701   }
702   END_FOR_SET_BITS
703   pop_conditions(CONDITIONS, depth);
704 }
705 
706 // Helper for the main recursive monomorphism function.
707 static ALWAYS_INLINE uint16_t
graph_mono_update_conditions(uint16_t const depth,uint16_t const last_defined,uint16_t const vertex)708 graph_mono_update_conditions(uint16_t const depth,
709                              uint16_t const last_defined,
710                              uint16_t const vertex) {
711   push_conditions(
712       CONDITIONS, depth, vertex, GRAPH2->neighbours[MAP[last_defined]]);
713   store_size_conditions(CONDITIONS, vertex);
714   return size_conditions(CONDITIONS, vertex);
715 }
716 
717 // The main recursive function for monomorphisms of graphs.
718 //
719 // The arguments are:
720 //
721 // 1. depth             The current depth of the search (i.e. number of
722 //                      positions in MAP that are assigned).
723 // 2. pos               The last position in MAP that was filled
724 // 3. rep_depth         The depth of the stabilisers of points in MAP already.
725 //                      This is different than depth because we stop increasing
726 //                      the rep_depth when we first encounter a trivial
727 //                      stabiliser.
728 // 4. has_trivial_stab  true if the stabiliser of the values already in MAP is
729 //                      trivial, false otherwise.
730 // 5. max_results       The maximum number of results to find.
731 // 6. count             The number of embeddings found so far.
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)732 static void find_graph_monos(uint16_t        depth,
733                              uint16_t        pos,
734                              uint16_t        rep_depth,
735                              bool            has_trivial_stab,
736                              uint64_t const  max_results,
737                              uint64_t* const count) {
738 #ifdef DIGRAPHS_ENABLE_STATS
739   update_stats(STATS, depth);
740 #endif
741   if (depth == GRAPH1->nr_vertices) {
742     // we've assigned every position in <MAP>
743     external_order_map_graph(GRAPH1);
744     Obj ret =
745         HOOK(USER_PARAM, MAX(GRAPH1->nr_vertices, GRAPH2->nr_vertices), MAP);
746     internal_order_map_graph(GRAPH1);
747     (*count)++;
748     if (*count >= max_results || ret == True) {
749       longjmp(OUTOFHERE, 1);
750     }
751     return;
752   }
753 
754   uint16_t next = 0;          // the next position to fill
755   uint16_t min  = UNDEFINED;  // the minimum number of candidates for MAP[next]
756   uint16_t i;
757 
758   BitArray* possible = BIT_ARRAY_BUFFER[depth];
759 
760   if (depth > 0) {  // this is not the first call of the function
761     copy_bit_array(
762         MAP_UNDEFINED[depth], MAP_UNDEFINED[depth - 1], GRAPH1->nr_vertices);
763     copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
764     intersect_bit_arrays(
765         possible, GRAPH1->neighbours[pos], GRAPH1->nr_vertices);
766     FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
767       size_t const n = graph_mono_update_conditions(depth, pos, i);
768       if (n == 0) {
769 #ifdef DIGRAPHS_ENABLE_STATS
770         STATS->nr_dead_branches++;
771 #endif
772         pop_conditions(CONDITIONS, depth);
773         return;
774       }
775       STORE_MIN(min, next, n, i);
776     }
777     END_FOR_SET_BITS
778     if (min > 1) {
779       copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
780       complement_bit_arrays(
781           possible, GRAPH1->neighbours[pos], GRAPH1->nr_vertices);
782       FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
783         size_t const n = size_conditions(CONDITIONS, i);
784         STORE_MIN_BREAK(min, next, n, i);
785       }
786       END_FOR_SET_BITS
787     }
788   } else {
789     for (i = 0; i < GRAPH1->nr_vertices; ++i) {
790       size_t const n = size_conditions(CONDITIONS, i);
791       STORE_MIN_BREAK(min, next, n, i);
792     }
793   }
794   copy_bit_array(
795       possible, get_conditions(CONDITIONS, next), GRAPH2->nr_vertices);
796   intersect_bit_arrays(possible, REPS[rep_depth], GRAPH2->nr_vertices);
797   complement_bit_arrays(possible, VALS, GRAPH2->nr_vertices);
798   FOR_SET_BITS(possible, GRAPH2->nr_vertices, i) {
799     MAP[next] = i;
800     set_bit_array(VALS, i, true);
801     set_bit_array(MAP_UNDEFINED[depth], next, false);
802     if (!has_trivial_stab) {
803       find_graph_monos(depth + 1,
804                        next,
805                        rep_depth + 1,
806                        compute_stabs_and_orbit_reps(GRAPH1->nr_vertices,
807                                                     GRAPH2->nr_vertices,
808                                                     rep_depth + 1,
809                                                     depth,
810                                                     i,
811                                                     false),
812                        max_results,
813                        count);
814     } else {
815       find_graph_monos(depth + 1, next, rep_depth, true, max_results, count);
816     }
817     MAP[next] = UNDEFINED;
818     set_bit_array(VALS, i, false);
819     set_bit_array(MAP_UNDEFINED[depth], next, true);
820   }
821   END_FOR_SET_BITS
822   pop_conditions(CONDITIONS, depth);
823 }
824 
825 // Helper for the main recursive embedding function.
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))826 static ALWAYS_INLINE uint16_t graph_embed_update_conditions(
827     uint16_t const depth,
828     uint16_t const last_defined,
829     uint16_t const vertex,
830     void (*oper)(BitArray* const, BitArray const* const, uint16_t const)) {
831   push_conditions(CONDITIONS, depth, vertex, NULL);
832   oper(get_conditions(CONDITIONS, vertex),
833        GRAPH2->neighbours[MAP[last_defined]],
834        GRAPH2->nr_vertices);
835   store_size_conditions(CONDITIONS, vertex);
836   return size_conditions(CONDITIONS, vertex);
837 }
838 
839 // The main recursive function for embeddings of graphs.
840 //
841 // The arguments are:
842 //
843 // 1. depth             The current depth of the search (i.e. number of
844 //                      positions in MAP that are assigned).
845 // 2. pos               The last position in MAP that was filled
846 // 3. rep_depth         The depth of the stabilisers of points in MAP already.
847 //                      This is different than depth because we stop increasing
848 //                      the rep_depth when we first encounter a trivial
849 //                      stabiliser.
850 // 4. has_trivial_stab  true if the stabiliser of the values already in MAP is
851 //                      trivial, false otherwise.
852 // 5. max_results       The maximum number of results to find.
853 // 6. count             The number of embeddings found so far.
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)854 static void find_graph_embeddings(uint16_t        depth,
855                                   uint16_t        pos,
856                                   uint16_t        rep_depth,
857                                   bool            has_trivial_stab,
858                                   uint64_t const  max_results,
859                                   uint64_t* const count) {
860 #ifdef DIGRAPHS_ENABLE_STATS
861   update_stats(STATS, depth);
862 #endif
863   if (depth == GRAPH1->nr_vertices) {
864     // we've assigned every position in <MAP>
865     external_order_map_graph(GRAPH1);
866     Obj ret =
867         HOOK(USER_PARAM, MAX(GRAPH1->nr_vertices, GRAPH2->nr_vertices), MAP);
868     internal_order_map_graph(GRAPH1);
869     (*count)++;
870     if (*count >= max_results || ret == True) {
871       longjmp(OUTOFHERE, 1);
872     }
873     return;
874   }
875 
876   uint16_t next = 0;          // the next position to fill
877   uint16_t min  = UNDEFINED;  // the minimum number of candidates for MAP[next]
878   uint16_t i;
879 
880   BitArray* possible = BIT_ARRAY_BUFFER[depth];
881 
882   if (depth > 0) {  // this is not the first call of the function
883     copy_bit_array(
884         MAP_UNDEFINED[depth], MAP_UNDEFINED[depth - 1], GRAPH1->nr_vertices);
885     copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
886     intersect_bit_arrays(
887         possible, GRAPH1->neighbours[pos], GRAPH1->nr_vertices);
888     FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
889       size_t const n =
890           graph_embed_update_conditions(depth, pos, i, &intersect_bit_arrays);
891       if (n == 0) {
892 #ifdef DIGRAPHS_ENABLE_STATS
893         STATS->nr_dead_branches++;
894 #endif
895         pop_conditions(CONDITIONS, depth);
896         return;
897       }
898       STORE_MIN(min, next, n, i);
899     }
900     END_FOR_SET_BITS
901     copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
902     complement_bit_arrays(
903         possible, GRAPH1->neighbours[pos], GRAPH1->nr_vertices);
904     FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
905       size_t const n =
906           graph_embed_update_conditions(depth, pos, i, &complement_bit_arrays);
907       if (n == 0) {
908 #ifdef DIGRAPHS_ENABLE_STATS
909         STATS->nr_dead_branches++;
910 #endif
911         pop_conditions(CONDITIONS, depth);
912         return;
913       }
914       STORE_MIN(min, next, n, i);
915     }
916     END_FOR_SET_BITS
917   } else {
918     for (i = 0; i < GRAPH1->nr_vertices; ++i) {
919       size_t const n = size_conditions(CONDITIONS, i);
920       STORE_MIN_BREAK(min, next, n, i);
921     }
922   }
923 
924   copy_bit_array(
925       possible, get_conditions(CONDITIONS, next), GRAPH2->nr_vertices);
926   intersect_bit_arrays(possible, REPS[rep_depth], GRAPH2->nr_vertices);
927   complement_bit_arrays(possible, VALS, GRAPH2->nr_vertices);
928 
929   FOR_SET_BITS(possible, GRAPH2->nr_vertices, i) {
930     MAP[next] = i;
931     set_bit_array(VALS, i, true);
932     set_bit_array(MAP_UNDEFINED[depth], next, false);
933     if (!has_trivial_stab) {
934       find_graph_embeddings(depth + 1,
935                             next,
936                             rep_depth + 1,
937                             compute_stabs_and_orbit_reps(GRAPH1->nr_vertices,
938                                                          GRAPH2->nr_vertices,
939                                                          rep_depth + 1,
940                                                          depth,
941                                                          MAP[next],
942                                                          false),
943                             max_results,
944                             count);
945     } else {
946       find_graph_embeddings(
947           depth + 1, next, rep_depth, true, max_results, count);
948     }
949     MAP[next] = UNDEFINED;
950     set_bit_array(VALS, i, false);
951     set_bit_array(MAP_UNDEFINED[depth], next, true);
952   }
953   END_FOR_SET_BITS
954   pop_conditions(CONDITIONS, depth);
955 }
956 
957 ////////////////////////////////////////////////////////////////////////////////
958 // The next function should be called before find_graph_homos.  This function
959 // simulates the first steps of the recursion using the information in
960 // partial_map_obj (if any) so that the values specified in partial_map_obj are
961 // installed in MAP first.
962 ////////////////////////////////////////////////////////////////////////////////
963 
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)964 static void init_partial_map_and_find_graph_homos(Obj partial_map_obj,
965                                                   uint64_t const  max_results,
966                                                   uint64_t const  hint,
967                                                   uint64_t* const count,
968                                                   Obj injective_obj) {
969   uint16_t depth                = 0;
970   uint16_t rep_depth            = 0;
971   uint16_t last_defined         = UNDEFINED;
972   bool     last_stab_is_trivial = (STAB_GENS[0]->size == 0 ? true : false);
973   uint16_t rank                 = 0;
974   uint16_t next                 = UNDEFINED;
975 
976   if (partial_map_obj != Fail) {
977     for (next = 0; next < LEN_LIST(partial_map_obj); ++next) {
978       if (ISB_LIST(partial_map_obj, next + 1)) {
979         if (depth > 0) {
980           DIGRAPHS_ASSERT(last_defined != UNDEFINED);
981           copy_bit_array(MAP_UNDEFINED[depth],
982                          MAP_UNDEFINED[depth - 1],
983                          GRAPH1->nr_vertices);
984           BitArray* possible = BIT_ARRAY_BUFFER[0];
985           copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
986           intersect_bit_arrays(
987               possible, GRAPH1->neighbours[last_defined], GRAPH1->nr_vertices);
988           if (INT_INTOBJ(injective_obj) == 0) {
989             uint16_t i;  // variable for the next for-loop
990             FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
991               uint16_t n = graph_homo_update_conditions(depth, last_defined, i);
992               if (n == 0) {
993                 return;
994               }
995             }
996             END_FOR_SET_BITS
997           } else if (INT_INTOBJ(injective_obj) == 1) {
998             uint16_t i;  // variable for the next for-loop
999             FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
1000               uint16_t n = graph_mono_update_conditions(depth, last_defined, i);
1001               if (n == 0) {
1002                 return;
1003               }
1004             }
1005             END_FOR_SET_BITS
1006           } else {
1007             DIGRAPHS_ASSERT(INT_INTOBJ(injective_obj) == 2);
1008             uint16_t i;  // variable for the next for-loop
1009             FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
1010               uint16_t n = graph_embed_update_conditions(
1011                   depth, last_defined, i, &intersect_bit_arrays);
1012               if (n == 0) {
1013                 return;
1014               }
1015               END_FOR_SET_BITS
1016             }
1017             copy_bit_array(possible, MAP_UNDEFINED[depth], GRAPH1->nr_vertices);
1018             complement_bit_arrays(possible,
1019                                   GRAPH1->neighbours[last_defined],
1020                                   GRAPH1->nr_vertices);
1021             FOR_SET_BITS(possible, GRAPH1->nr_vertices, i) {
1022               uint16_t n = graph_embed_update_conditions(
1023                   depth, last_defined, i, &complement_bit_arrays);
1024               if (n == 0) {
1025                 return;
1026               }
1027             }
1028             END_FOR_SET_BITS
1029           }
1030         }
1031         uint16_t val = INT_INTOBJ(ELM_LIST(partial_map_obj, next + 1)) - 1;
1032         next         = (ORDERED ? INVERSE_ORDER[next] : next);
1033         MAP[next]    = val;
1034         if (!get_bit_array(VALS, MAP[next])) {
1035           rank++;
1036           if (rank > hint) {
1037             return;
1038           }
1039         }
1040         set_bit_array(VALS, MAP[next], true);
1041         set_bit_array(MAP_UNDEFINED[depth], next, false);
1042         if (!last_stab_is_trivial) {
1043           last_stab_is_trivial =
1044               compute_stabs_and_orbit_reps(GRAPH1->nr_vertices,
1045                                            GRAPH2->nr_vertices,
1046                                            rep_depth + 1,
1047                                            depth,
1048                                            MAP[next],
1049                                            false);
1050           rep_depth++;
1051         }
1052         depth++;
1053         last_defined = next;
1054         next         = (ORDERED ? ORDER[next] : next);
1055       }
1056     }
1057   }
1058   if (INT_INTOBJ(injective_obj) == 0) {
1059     find_graph_homos(depth,
1060                      last_defined,
1061                      rep_depth,
1062                      last_stab_is_trivial,
1063                      rank,
1064                      max_results,
1065                      hint,
1066                      count);
1067   } else if (INT_INTOBJ(injective_obj) == 1) {
1068     find_graph_monos(depth,
1069                      last_defined,
1070                      rep_depth,
1071                      last_stab_is_trivial,
1072                      max_results,
1073                      count);
1074   } else if (INT_INTOBJ(injective_obj) == 2) {
1075     find_graph_embeddings(depth,
1076                           last_defined,
1077                           rep_depth,
1078                           last_stab_is_trivial,
1079                           max_results,
1080                           count);
1081   }
1082 }
1083 
1084 // Helper for the main recursive homomorphism of digraphs function.
1085 static ALWAYS_INLINE uint16_t
digraph_homo_update_conditions(uint16_t const depth,uint16_t const last_defined,uint16_t const vertex)1086 digraph_homo_update_conditions(uint16_t const depth,
1087                                uint16_t const last_defined,
1088                                uint16_t const vertex) {
1089   if (is_adjacent_digraph(DIGRAPH1, last_defined, vertex)) {
1090     push_conditions(
1091         CONDITIONS, depth, vertex, DIGRAPH2->out_neighbours[MAP[last_defined]]);
1092     if (is_adjacent_digraph(DIGRAPH1, vertex, last_defined)) {
1093       intersect_bit_arrays(get_conditions(CONDITIONS, vertex),
1094                            DIGRAPH2->in_neighbours[MAP[last_defined]],
1095                            DIGRAPH2->nr_vertices);
1096     }
1097     store_size_conditions(CONDITIONS, vertex);
1098   } else if (is_adjacent_digraph(DIGRAPH1, vertex, last_defined)) {
1099     push_conditions(
1100         CONDITIONS, depth, vertex, DIGRAPH2->in_neighbours[MAP[last_defined]]);
1101     store_size_conditions(CONDITIONS, vertex);
1102   }
1103   return size_conditions(CONDITIONS, vertex);
1104 }
1105 
1106 // The main recursive function for homomorphisms of digraphs.
1107 //
1108 // The arguments are:
1109 //
1110 // 1. depth             The current depth of the search (i.e. number of
1111 //                      positions in MAP that are assigned).
1112 // 2. pos               The last position in MAP that was filled
1113 // 3. rep_depth         The depth of the stabilisers of points in MAP already.
1114 //                      This is different than depth because we stop increasing
1115 //                      the rep_depth when we first encounter a trivial
1116 //                      stabiliser.
1117 // 4. has_trivial_stab  true if the stabiliser of the values already in MAP is
1118 //                      trivial, false otherwise.
1119 // 5. rank              The current number of distinct values in MAP.
1120 // 6. max_results       The maximum number of results to find.
1121 // 7. hint              The desired number of distinct points in a (full)
1122 //                      homomorphism.
1123 // 8. count             The number of homomorphisms found so far.
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)1124 static void find_digraph_homos(uint16_t        depth,
1125                                uint16_t        pos,
1126                                uint16_t        rep_depth,
1127                                bool            has_trivial_stab,
1128                                uint16_t        rank,
1129                                uint64_t const  max_results,
1130                                uint64_t const  hint,
1131                                uint64_t* const count) {
1132 #ifdef DIGRAPHS_ENABLE_STATS
1133   update_stats(STATS, depth);
1134 #endif
1135   if (depth == DIGRAPH1->nr_vertices) {
1136     // we've assigned every position in <MAP>
1137     if (hint != UNDEFINED && rank != hint) {
1138 #ifdef DIGRAPHS_ENABLE_STATS
1139       STATS->nr_dead_branches++;
1140 #endif
1141       return;
1142     }
1143     external_order_map_digraph(DIGRAPH1);
1144     Obj ret = HOOK(
1145         USER_PARAM, MAX(DIGRAPH1->nr_vertices, DIGRAPH2->nr_vertices), MAP);
1146     internal_order_map_digraph(DIGRAPH1);
1147     (*count)++;
1148     if (*count >= max_results || ret == True) {
1149       longjmp(OUTOFHERE, 1);
1150     }
1151     return;
1152   }
1153 
1154   uint16_t next = 0;          // the next position to fill
1155   uint16_t min  = UNDEFINED;  // the minimum number of candidates for MAP[next]
1156   uint16_t i;
1157 
1158   if (depth > 0) {  // this is not the first call of the function
1159     copy_bit_array(
1160         MAP_UNDEFINED[depth], MAP_UNDEFINED[depth - 1], DIGRAPH1->nr_vertices);
1161     FOR_SET_BITS(MAP_UNDEFINED[depth], DIGRAPH1->nr_vertices, i) {
1162       uint16_t const n = digraph_homo_update_conditions(depth, pos, i);
1163       if (n == 0) {
1164 #ifdef DIGRAPHS_ENABLE_STATS
1165         STATS->nr_dead_branches++;
1166 #endif
1167         pop_conditions(CONDITIONS, depth);
1168         return;
1169       }
1170       STORE_MIN(min, next, n, i);
1171     }
1172     END_FOR_SET_BITS
1173   } else {  // depth == 0
1174     for (i = 0; i < DIGRAPH1->nr_vertices; ++i) {
1175       size_t const n = size_conditions(CONDITIONS, i);
1176       STORE_MIN_BREAK(min, next, n, i);
1177     }
1178   }
1179 
1180   BitArray* possible = BIT_ARRAY_BUFFER[depth];
1181 
1182   if (rank < hint) {
1183     copy_bit_array(
1184         possible, get_conditions(CONDITIONS, next), DIGRAPH2->nr_vertices);
1185     complement_bit_arrays(possible, VALS, DIGRAPH2->nr_vertices);
1186     intersect_bit_arrays(possible, REPS[rep_depth], DIGRAPH2->nr_vertices);
1187     FOR_SET_BITS(possible, DIGRAPH2->nr_vertices, i) {
1188       MAP[next] = i;
1189       set_bit_array(VALS, i, true);
1190       set_bit_array(MAP_UNDEFINED[depth], next, false);
1191       if (!has_trivial_stab) {
1192         find_digraph_homos(depth + 1,
1193                            next,
1194                            rep_depth + 1,
1195                            compute_stabs_and_orbit_reps(DIGRAPH1->nr_vertices,
1196                                                         DIGRAPH2->nr_vertices,
1197                                                         rep_depth + 1,
1198                                                         depth,
1199                                                         i,
1200                                                         false),
1201                            rank + 1,
1202                            max_results,
1203                            hint,
1204                            count);
1205       } else {
1206         find_digraph_homos(depth + 1,
1207                            next,
1208                            rep_depth,
1209                            true,
1210                            rank + 1,
1211                            max_results,
1212                            hint,
1213                            count);
1214       }
1215       MAP[next] = UNDEFINED;
1216       set_bit_array(VALS, i, false);
1217       set_bit_array(MAP_UNDEFINED[depth], next, true);
1218     }
1219     END_FOR_SET_BITS
1220   }
1221   copy_bit_array(
1222       possible, get_conditions(CONDITIONS, next), DIGRAPH2->nr_vertices);
1223   intersect_bit_arrays(possible, VALS, DIGRAPH2->nr_vertices);
1224   FOR_SET_BITS(possible, DIGRAPH2->nr_vertices, i) {
1225     MAP[next] = i;
1226     set_bit_array(MAP_UNDEFINED[depth], next, false);
1227     find_digraph_homos(depth + 1,
1228                        next,
1229                        rep_depth,
1230                        has_trivial_stab,
1231                        rank,
1232                        max_results,
1233                        hint,
1234                        count);
1235     MAP[next] = UNDEFINED;
1236     set_bit_array(MAP_UNDEFINED[depth], next, true);
1237   }
1238   END_FOR_SET_BITS
1239   pop_conditions(CONDITIONS, depth);
1240 }
1241 
1242 // Helper for the main recursive monomorphism of digraphs function.
1243 static ALWAYS_INLINE uint16_t
digraph_mono_update_conditions(uint16_t const depth,uint16_t const last_defined,uint16_t const vertex)1244 digraph_mono_update_conditions(uint16_t const depth,
1245                                uint16_t const last_defined,
1246                                uint16_t const vertex) {
1247   push_conditions(CONDITIONS, depth, vertex, NULL);
1248   if (is_adjacent_digraph(DIGRAPH1, last_defined, vertex)) {
1249     intersect_bit_arrays(get_conditions(CONDITIONS, vertex),
1250                          DIGRAPH2->out_neighbours[MAP[last_defined]],
1251                          DIGRAPH2->nr_vertices);
1252   }
1253   if (is_adjacent_digraph(DIGRAPH1, vertex, last_defined)) {
1254     intersect_bit_arrays(get_conditions(CONDITIONS, vertex),
1255                          DIGRAPH2->in_neighbours[MAP[last_defined]],
1256                          DIGRAPH2->nr_vertices);
1257   }
1258   store_size_conditions(CONDITIONS, vertex);
1259 
1260   return size_conditions(CONDITIONS, vertex);
1261 }
1262 
1263 // The main recursive function for monomorphisms of digraphs.
1264 //
1265 // The arguments are:
1266 //
1267 // 1. depth             The current depth of the search (i.e. number of
1268 //                      positions in MAP that are assigned).
1269 // 2. pos               The last position in MAP that was filled
1270 // 3. rep_depth         The depth of the stabilisers of points in MAP already.
1271 //                      This is different than depth because we stop increasing
1272 //                      the rep_depth when we first encounter a trivial
1273 //                      stabiliser.
1274 // 4. has_trivial_stab  true if the stabiliser of the values already in MAP is
1275 //                      trivial, false otherwise.
1276 // 5. max_results       The maximum number of results to find.
1277 // 6. count             The number of embeddings found so far.
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)1278 static void find_digraph_monos(uint16_t        depth,
1279                                uint16_t        pos,
1280                                uint16_t        rep_depth,
1281                                bool            has_trivial_stab,
1282                                uint64_t const  max_results,
1283                                uint64_t* const count) {
1284 #ifdef DIGRAPHS_ENABLE_STATS
1285   update_stats(STATS, depth);
1286 #endif
1287   if (depth == DIGRAPH1->nr_vertices) {
1288     // we've assigned every position in <MAP>
1289     external_order_map_digraph(DIGRAPH1);
1290     Obj ret = HOOK(
1291         USER_PARAM, MAX(DIGRAPH1->nr_vertices, DIGRAPH2->nr_vertices), MAP);
1292     internal_order_map_digraph(DIGRAPH1);
1293     (*count)++;
1294     if (*count >= max_results || ret == True) {
1295       longjmp(OUTOFHERE, 1);
1296     }
1297     return;
1298   }
1299 
1300   uint16_t next = 0;          // the next position to fill
1301   uint16_t min  = UNDEFINED;  // the minimum number of candidates for MAP[next]
1302   uint16_t i;
1303 
1304   if (depth > 0) {  // this is not the first call of the function
1305     copy_bit_array(
1306         MAP_UNDEFINED[depth], MAP_UNDEFINED[depth - 1], DIGRAPH1->nr_vertices);
1307     FOR_SET_BITS(MAP_UNDEFINED[depth], DIGRAPH1->nr_vertices, i) {
1308       size_t const n = digraph_mono_update_conditions(depth, pos, i);
1309       if (n == 0) {
1310 #ifdef DIGRAPHS_ENABLE_STATS
1311         STATS->nr_dead_branches++;
1312 #endif
1313         pop_conditions(CONDITIONS, depth);
1314         return;
1315       }
1316       STORE_MIN(min, next, n, i);
1317     }
1318     END_FOR_SET_BITS
1319   } else {
1320     for (i = 0; i < DIGRAPH1->nr_vertices; i++) {
1321       size_t const n = size_conditions(CONDITIONS, i);
1322       STORE_MIN_BREAK(min, next, n, i);
1323     }
1324   }
1325 
1326   BitArray* possible = BIT_ARRAY_BUFFER[depth];
1327   copy_bit_array(
1328       possible, get_conditions(CONDITIONS, next), DIGRAPH2->nr_vertices);
1329   intersect_bit_arrays(possible, REPS[rep_depth], DIGRAPH2->nr_vertices);
1330   complement_bit_arrays(possible, VALS, DIGRAPH2->nr_vertices);
1331   FOR_SET_BITS(possible, DIGRAPH2->nr_vertices, i) {
1332     MAP[next] = i;
1333     set_bit_array(VALS, i, true);
1334     set_bit_array(MAP_UNDEFINED[depth], next, false);
1335     if (!has_trivial_stab) {
1336       find_digraph_monos(depth + 1,
1337                          next,
1338                          rep_depth + 1,
1339                          compute_stabs_and_orbit_reps(DIGRAPH1->nr_vertices,
1340                                                       DIGRAPH2->nr_vertices,
1341                                                       rep_depth + 1,
1342                                                       depth,
1343                                                       i,
1344                                                       false),
1345                          max_results,
1346                          count);
1347     } else {
1348       find_digraph_monos(depth + 1, next, rep_depth, true, max_results, count);
1349     }
1350     MAP[next] = UNDEFINED;
1351     set_bit_array(VALS, i, false);
1352     set_bit_array(MAP_UNDEFINED[depth], next, true);
1353   }
1354   END_FOR_SET_BITS
1355   pop_conditions(CONDITIONS, depth);
1356 }
1357 
1358 // Helper for the main recursive embedding digraphs function.
1359 static ALWAYS_INLINE uint16_t
digraph_embed_update_conditions(uint16_t const depth,uint16_t const last_def,uint16_t const vertex)1360 digraph_embed_update_conditions(uint16_t const depth,
1361                                 uint16_t const last_def,
1362                                 uint16_t const vertex) {
1363   push_conditions(CONDITIONS, depth, vertex, NULL);
1364   if (is_adjacent_digraph(DIGRAPH1, last_def, vertex)) {
1365     intersect_bit_arrays(get_conditions(CONDITIONS, vertex),
1366                          DIGRAPH2->out_neighbours[MAP[last_def]],
1367                          DIGRAPH2->nr_vertices);
1368   } else {
1369     complement_bit_arrays(get_conditions(CONDITIONS, vertex),
1370                           DIGRAPH2->out_neighbours[MAP[last_def]],
1371                           DIGRAPH2->nr_vertices);
1372   }
1373   if (is_adjacent_digraph(DIGRAPH1, vertex, last_def)) {
1374     intersect_bit_arrays(get_conditions(CONDITIONS, vertex),
1375                          DIGRAPH2->in_neighbours[MAP[last_def]],
1376                          DIGRAPH2->nr_vertices);
1377   } else {
1378     complement_bit_arrays(get_conditions(CONDITIONS, vertex),
1379                           DIGRAPH2->in_neighbours[MAP[last_def]],
1380                           DIGRAPH2->nr_vertices);
1381   }
1382   store_size_conditions(CONDITIONS, vertex);
1383   return size_conditions(CONDITIONS, vertex);
1384 }
1385 
1386 // The main recursive function for embeddings of digraphs.
1387 //
1388 // The arguments are:
1389 //
1390 // 1. depth             The current depth of the search (i.e. number of
1391 //                      positions in MAP that are assigned).
1392 // 2. pos               The last position in MAP that was filled
1393 // 3. rep_depth         The depth of the stabilisers of points in MAP already.
1394 //                      This is different than depth because we stop increasing
1395 //                      the rep_depth when we first encounter a trivial
1396 //                      stabiliser.
1397 // 4. has_trivial_stab  true if the stabiliser of the values already in MAP is
1398 //                      trivial, false otherwise.
1399 // 5. max_results       The maximum number of results to find.
1400 // 6. count             The number of embeddings found so far.
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)1401 static void find_digraph_embeddings(uint16_t        depth,
1402                                     uint16_t        pos,
1403                                     uint16_t        rep_depth,
1404                                     bool            has_trivial_stab,
1405                                     uint64_t const  max_results,
1406                                     uint64_t* const count) {
1407 #ifdef DIGRAPHS_ENABLE_STATS
1408   update_stats(STATS, depth);
1409 #endif
1410   if (depth == DIGRAPH1->nr_vertices) {
1411     // we've assigned every position in <MAP>
1412     external_order_map_digraph(DIGRAPH1);
1413     Obj ret = HOOK(
1414         USER_PARAM, MAX(DIGRAPH1->nr_vertices, DIGRAPH2->nr_vertices), MAP);
1415     internal_order_map_digraph(DIGRAPH1);
1416     (*count)++;
1417     if (*count >= max_results || ret == True) {
1418       longjmp(OUTOFHERE, 1);
1419     }
1420     return;
1421   }
1422 
1423   uint16_t next = 0;          // the next position to fill
1424   uint16_t min  = UNDEFINED;  // the minimum number of candidates for MAP[next]
1425   uint16_t i;
1426 
1427   if (depth > 0) {  // this is not the first call of the function
1428     copy_bit_array(
1429         MAP_UNDEFINED[depth], MAP_UNDEFINED[depth - 1], DIGRAPH1->nr_vertices);
1430     FOR_SET_BITS(MAP_UNDEFINED[depth], DIGRAPH1->nr_vertices, i) {
1431       size_t const n = digraph_embed_update_conditions(depth, pos, i);
1432       if (n == 0) {
1433 #ifdef DIGRAPHS_ENABLE_STATS
1434         STATS->nr_dead_branches++;
1435 #endif
1436         pop_conditions(CONDITIONS, depth);
1437         return;
1438       }
1439       STORE_MIN(min, next, n, i);
1440     }
1441     END_FOR_SET_BITS
1442   } else {
1443     for (i = 0; i < DIGRAPH1->nr_vertices; i++) {
1444       size_t const n = size_conditions(CONDITIONS, i);
1445       STORE_MIN_BREAK(min, next, n, i);
1446     }
1447   }
1448 
1449   BitArray* possible = get_conditions(CONDITIONS, next);
1450   copy_bit_array(
1451       possible, get_conditions(CONDITIONS, next), DIGRAPH2->nr_vertices);
1452   intersect_bit_arrays(possible, REPS[rep_depth], DIGRAPH2->nr_vertices);
1453   complement_bit_arrays(possible, VALS, DIGRAPH2->nr_vertices);
1454 
1455   FOR_SET_BITS(possible, DIGRAPH2->nr_vertices, i) {
1456     MAP[next] = i;
1457     set_bit_array(VALS, i, true);
1458     set_bit_array(MAP_UNDEFINED[depth], next, false);
1459     if (!has_trivial_stab) {
1460       find_digraph_embeddings(
1461           depth + 1,
1462           next,
1463           rep_depth + 1,
1464           compute_stabs_and_orbit_reps(DIGRAPH1->nr_vertices,
1465                                        DIGRAPH2->nr_vertices,
1466                                        rep_depth + 1,
1467                                        depth,
1468                                        i,
1469                                        false),
1470           max_results,
1471           count);
1472     } else {
1473       find_digraph_embeddings(
1474           depth + 1, next, rep_depth, true, max_results, count);
1475     }
1476     MAP[next] = UNDEFINED;
1477     set_bit_array(VALS, i, false);
1478     set_bit_array(MAP_UNDEFINED[depth], next, true);
1479   }
1480   END_FOR_SET_BITS
1481   pop_conditions(CONDITIONS, depth);
1482 }
1483 
1484 ////////////////////////////////////////////////////////////////////////////////
1485 // The next function should be called before find_digraph_homos.  This function
1486 // simulates the first steps of the recursion using the information in
1487 // partial_map_obj (if any) so that the values specified in partial_map_obj are
1488 // installed in MAP first.
1489 ////////////////////////////////////////////////////////////////////////////////
1490 
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)1491 static void init_partial_map_and_find_digraph_homos(Obj partial_map_obj,
1492                                                     uint64_t const  max_results,
1493                                                     uint64_t const  hint,
1494                                                     uint64_t* const count,
1495                                                     Obj injective_obj) {
1496   uint16_t depth                = 0;
1497   uint16_t rep_depth            = 0;
1498   uint16_t last_defined         = UNDEFINED;
1499   bool     last_stab_is_trivial = (STAB_GENS[0]->size == 0 ? true : false);
1500   uint16_t rank                 = 0;
1501   uint16_t next                 = UNDEFINED;
1502 
1503   if (partial_map_obj != Fail) {
1504     for (next = 0; next < LEN_LIST(partial_map_obj); ++next) {
1505       if (ISB_LIST(partial_map_obj, next + 1)) {
1506         if (depth > 0) {
1507           DIGRAPHS_ASSERT(last_defined != UNDEFINED);
1508           copy_bit_array(MAP_UNDEFINED[depth],
1509                          MAP_UNDEFINED[depth - 1],
1510                          DIGRAPH1->nr_vertices);
1511           uint16_t i;  // variable for the next for-loop
1512           FOR_SET_BITS(MAP_UNDEFINED[depth], DIGRAPH1->nr_vertices, i) {
1513             uint16_t n;
1514             if (INT_INTOBJ(injective_obj) == 0) {
1515               n = digraph_homo_update_conditions(depth, last_defined, i);
1516             } else if (INT_INTOBJ(injective_obj) == 1) {
1517               n = digraph_mono_update_conditions(depth, last_defined, i);
1518             } else {
1519               DIGRAPHS_ASSERT(INT_INTOBJ(injective_obj) == 2);
1520               n = digraph_embed_update_conditions(depth, last_defined, i);
1521             }
1522             if (n == 0) {
1523               return;
1524             }
1525           }
1526           END_FOR_SET_BITS
1527         }
1528         uint16_t val = INT_INTOBJ(ELM_LIST(partial_map_obj, next + 1)) - 1;
1529         next         = (ORDERED ? INVERSE_ORDER[next] : next);
1530         MAP[next]    = val;
1531         if (!get_bit_array(VALS, MAP[next])) {
1532           rank++;
1533           if (rank > hint) {
1534             return;
1535           }
1536         }
1537         set_bit_array(VALS, MAP[next], true);
1538         set_bit_array(MAP_UNDEFINED[depth], next, false);
1539         if (!last_stab_is_trivial) {
1540           last_stab_is_trivial =
1541               compute_stabs_and_orbit_reps(DIGRAPH1->nr_vertices,
1542                                            DIGRAPH2->nr_vertices,
1543                                            rep_depth + 1,
1544                                            depth,
1545                                            MAP[next],
1546                                            false);
1547           rep_depth++;
1548         }
1549         depth++;
1550         last_defined = next;
1551         next         = (ORDERED ? ORDER[next] : next);
1552       }
1553     }
1554   }
1555   if (INT_INTOBJ(injective_obj) == 0) {
1556     find_digraph_homos(depth,
1557                        last_defined,
1558                        rep_depth,
1559                        last_stab_is_trivial,
1560                        rank,
1561                        max_results,
1562                        hint,
1563                        count);
1564   } else if (INT_INTOBJ(injective_obj) == 1) {
1565     find_digraph_monos(depth,
1566                        last_defined,
1567                        rep_depth,
1568                        last_stab_is_trivial,
1569                        max_results,
1570                        count);
1571   } else {
1572     DIGRAPHS_ASSERT(INT_INTOBJ(injective_obj) == 2);
1573     find_digraph_embeddings(depth,
1574                             last_defined,
1575                             rep_depth,
1576                             last_stab_is_trivial,
1577                             max_results,
1578                             count);
1579   }
1580   return;
1581 }
1582 
1583 ////////////////////////////////////////////////////////////////////////////////
1584 // 7. The GAP-level function (and helpers)
1585 ////////////////////////////////////////////////////////////////////////////////
1586 
1587 // Initialises the data structures required by the recursive functions for
1588 // finding homomorphisms. If true is returned everything was initialised ok, if
1589 // false is returned, then the arguments already imply that there can be no
1590 // homomorphisms.
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)1591 static bool init_data_from_args(Obj digraph1_obj,
1592                                 Obj digraph2_obj,
1593                                 Obj hook_obj,
1594                                 Obj user_param_obj,
1595                                 Obj max_results_obj,
1596                                 Obj hint_obj,
1597                                 Obj injective_obj,
1598                                 Obj image_obj,
1599                                 Obj partial_map_obj,
1600                                 Obj colors1_obj,
1601                                 Obj colors2_obj,
1602                                 Obj order_obj,
1603                                 Obj aut_grp_obj) {
1604   static bool is_initialized = false;  // did we call this method before?
1605   if (!is_initialized) {
1606     // srand(time(0));
1607     is_initialized = true;
1608 #ifdef DIGRAPHS_ENABLE_STATS
1609     STATS = malloc(sizeof(HomoStats));
1610 #endif
1611 
1612     DIGRAPH1 = new_digraph(MAXVERTS);
1613     DIGRAPH2 = new_digraph(MAXVERTS);
1614 
1615     GRAPH1 = new_graph(MAXVERTS);
1616     GRAPH2 = new_graph(MAXVERTS);
1617 
1618     IMAGE_RESTRICT = new_bit_array(MAXVERTS);
1619     ORB_LOOKUP     = new_bit_array(MAXVERTS);
1620     REPS           = malloc(MAXVERTS * sizeof(BitArray*));
1621     for (uint16_t i = 0; i < MAXVERTS; i++) {
1622       BLISS_GRAPH[i]      = bliss_digraphs_new(i);
1623       REPS[i]             = new_bit_array(MAXVERTS);
1624       BIT_ARRAY_BUFFER[i] = new_bit_array(MAXVERTS);
1625       MAP_UNDEFINED[i]    = new_bit_array(MAXVERTS);
1626       STAB_GENS[i]        = new_perm_coll(MAXVERTS, MAXVERTS);
1627     }
1628     VALS          = new_bit_array(MAXVERTS);
1629     CONDITIONS    = new_conditions(MAXVERTS, MAXVERTS);
1630     SCHREIER_SIMS = new_schreier_sims();
1631   }
1632 #ifdef DIGRAPHS_ENABLE_STATS
1633   clear_stats(STATS);
1634 #endif
1635 
1636   uint16_t nr1 = DigraphNrVertices(digraph1_obj);
1637   uint16_t nr2 = DigraphNrVertices(digraph2_obj);
1638 
1639   init_bit_array(MAP_UNDEFINED[0], true, nr1);
1640   init_bit_array(VALS, false, nr2);
1641 
1642   if (IS_LIST(order_obj)) {
1643     ORDERED = true;
1644     for (uint16_t i = 0; i < nr1; i++) {
1645       ORDER[i]                = INT_INTOBJ(ELM_LIST(order_obj, i + 1)) - 1;
1646       INVERSE_ORDER[ORDER[i]] = i;
1647     }
1648   } else {
1649     ORDERED = false;
1650   }
1651 
1652   bool is_undirected;
1653   if (CALL_1ARGS(IsSymmetricDigraph, digraph1_obj) == True
1654       && CALL_1ARGS(IsSymmetricDigraph, digraph2_obj) == True) {
1655     init_graph_from_digraph_obj(GRAPH1, digraph1_obj, ORDERED);
1656     init_graph_from_digraph_obj(GRAPH2, digraph2_obj, false);
1657     is_undirected = true;
1658   } else {
1659     init_digraph_from_digraph_obj(DIGRAPH1, digraph1_obj, ORDERED);
1660     init_digraph_from_digraph_obj(DIGRAPH2, digraph2_obj, false);
1661     is_undirected = false;
1662   }
1663 
1664   if (hook_obj != Fail) {
1665     GAP_FUNC = hook_obj;
1666     HOOK     = homo_hook_gap;
1667   } else {
1668     HOOK = homo_hook_collect;
1669   }
1670   USER_PARAM = user_param_obj;
1671 
1672   clear_conditions(CONDITIONS, nr1, nr2);
1673 
1674   // IMAGE_RESTRICT is a pointer to a BitArray of possible image values for the
1675   // homomorphisms, it is also used by orbit_reps so that orbit reps are chosen
1676   // from among the restricted values of the image . . .
1677   set_bit_array_from_gap_list(IMAGE_RESTRICT, image_obj);
1678   if (INT_INTOBJ(injective_obj) > 0
1679       && size_bit_array(IMAGE_RESTRICT, nr1) < nr1) {
1680     // homomorphisms should be injective (by injective_obj) but are not since
1681     // the image is too restricted.
1682     return false;
1683   }
1684 
1685   init_bit_array(BIT_ARRAY_BUFFER[0], false, nr2);
1686   if (partial_map_obj != Fail) {
1687     for (uint16_t i = 1; i <= LEN_LIST(partial_map_obj); i++) {
1688       if (ISB_LIST(partial_map_obj, i)) {
1689         Obj o = ELM_LIST(partial_map_obj, i);
1690         DIGRAPHS_ASSERT(IS_INTOBJ(o));
1691         DIGRAPHS_ASSERT(INT_INTOBJ(o) > 0);
1692         DIGRAPHS_ASSERT(INT_INTOBJ(o) <= nr2);
1693         if (INT_INTOBJ(injective_obj) > 0) {
1694           if (get_bit_array(BIT_ARRAY_BUFFER[0], INT_INTOBJ(o) - 1)) {
1695             // partial_map_obj should be injective, but is not
1696             return false;
1697           }
1698           set_bit_array(BIT_ARRAY_BUFFER[0], INT_INTOBJ(o) - 1, true);
1699         }
1700         // The only value that vertex `i` can have is `o`!
1701         if (ORDERED) {
1702           init_bit_array(
1703               get_conditions(CONDITIONS, INVERSE_ORDER[i - 1]), false, nr2);
1704           set_bit_array_from_gap_int(
1705               get_conditions(CONDITIONS, INVERSE_ORDER[i - 1]), o);
1706         } else {
1707           init_bit_array(get_conditions(CONDITIONS, i - 1), false, nr2);
1708           set_bit_array_from_gap_int(get_conditions(CONDITIONS, i - 1), o);
1709         }
1710       }
1711       // Intersect everything in the first row of the conditions with <image>,
1712       intersect_bit_arrays(
1713           get_conditions(CONDITIONS, i - 1), IMAGE_RESTRICT, nr2);
1714     }
1715   }
1716 
1717   init_bit_array(BIT_ARRAY_BUFFER[0], false, nr2);
1718   if (is_undirected) {
1719     for (uint16_t i = 0; i < nr2; i++) {
1720       if (is_adjacent_graph(GRAPH2, i, i)) {
1721         set_bit_array(BIT_ARRAY_BUFFER[0], i, true);
1722       }
1723     }
1724     // Loops in digraph1 can only MAP to loops in digraph2
1725     for (uint16_t i = 0; i < nr1; i++) {
1726       if (is_adjacent_graph(GRAPH1, i, i)) {
1727         intersect_bit_arrays(
1728             get_conditions(CONDITIONS, i), BIT_ARRAY_BUFFER[0], nr2);
1729       }
1730     }
1731   } else {
1732     for (uint16_t i = 0; i < nr2; i++) {
1733       if (is_adjacent_digraph(DIGRAPH2, i, i)) {
1734         set_bit_array(BIT_ARRAY_BUFFER[0], i, true);
1735       }
1736     }
1737     // Loops in digraph1 can only MAP to loops in digraph2
1738     for (uint16_t i = 0; i < nr1; i++) {
1739       if (is_adjacent_digraph(DIGRAPH1, i, i)) {
1740         intersect_bit_arrays(
1741             get_conditions(CONDITIONS, i), BIT_ARRAY_BUFFER[0], nr2);
1742       }
1743     }
1744   }
1745 
1746   // Process the vertex colours . . .
1747   uint16_t* colors;
1748   if (colors1_obj != Fail && colors2_obj != Fail) {
1749     DIGRAPHS_ASSERT(IS_LIST(colors1_obj));
1750     DIGRAPHS_ASSERT(IS_LIST(colors2_obj));
1751     DIGRAPHS_ASSERT(LEN_LIST(colors1_obj) == nr1);
1752     DIGRAPHS_ASSERT(LEN_LIST(colors2_obj) == nr2);
1753 
1754     for (uint16_t i = 1; i <= LEN_LIST(colors2_obj); i++) {
1755       DIGRAPHS_ASSERT(ISB_LIST(colors2_obj, i));
1756       DIGRAPHS_ASSERT(IS_INTOBJ(ELM_LIST(colors2_obj, i)));
1757       COLORS2[i - 1] = INT_INTOBJ(ELM_LIST(colors2_obj, i)) - 1;
1758     }
1759     for (uint16_t i = 1; i <= LEN_LIST(colors1_obj); i++) {
1760       init_bit_array(BIT_ARRAY_BUFFER[0], false, nr2);
1761       DIGRAPHS_ASSERT(ISB_LIST(colors1_obj, i));
1762       DIGRAPHS_ASSERT(IS_INTOBJ(ELM_LIST(colors1_obj, i)));
1763       for (uint16_t j = 1; j <= LEN_LIST(colors2_obj); j++) {
1764         if (INT_INTOBJ(ELM_LIST(colors1_obj, i))
1765             == INT_INTOBJ(ELM_LIST(colors2_obj, j))) {
1766           set_bit_array(BIT_ARRAY_BUFFER[0], j - 1, true);
1767         }
1768       }
1769       if (ORDERED) {
1770         intersect_bit_arrays(get_conditions(CONDITIONS, INVERSE_ORDER[i - 1]),
1771                              BIT_ARRAY_BUFFER[0],
1772                              nr2);
1773       } else {
1774         intersect_bit_arrays(
1775             get_conditions(CONDITIONS, i - 1), BIT_ARRAY_BUFFER[0], nr2);
1776       }
1777       // can only map vertices of color i to vertices of color i
1778     }
1779     colors = COLORS2;
1780   } else {
1781     colors = NULL;
1782   }
1783 
1784   // Ensure that the sizes of all conditions are known before we start and
1785   // define the MAP
1786 
1787   for (uint16_t i = 0; i < nr1; i++) {
1788     store_size_conditions(CONDITIONS, i);
1789     MAP[i] = UNDEFINED;
1790   }
1791   for (uint16_t i = nr1; i < MAX(nr1, nr2); i++) {
1792     MAP[i] = i;
1793   }
1794 
1795   // Get generators of the automorphism group of the second (di)graph, and the
1796   // orbit reps
1797   PERM_DEGREE = nr2;
1798   if (aut_grp_obj == Fail) {
1799     if (colors == NULL) {
1800       get_automorphism_group_from_gap(digraph2_obj, STAB_GENS[0]);
1801     } else if (is_undirected) {
1802 #ifdef DIGRAPHS_WITH_INCLUDED_BLISS
1803       automorphisms_graph(
1804           GRAPH2, colors, STAB_GENS[0], BLISS_GRAPH[PERM_DEGREE]);
1805 #else
1806       automorphisms_graph(GRAPH2, colors, STAB_GENS[0]);
1807 #endif
1808     } else {
1809 #ifdef DIGRAPHS_WITH_INCLUDED_BLISS
1810       automorphisms_digraph(
1811           DIGRAPH2, colors, STAB_GENS[0], BLISS_GRAPH[3 * PERM_DEGREE]);
1812 #else
1813       automorphisms_digraph(DIGRAPH2, colors, STAB_GENS[0]);
1814 #endif
1815     }
1816   } else {
1817     set_automorphisms(aut_grp_obj, STAB_GENS[0]);
1818   }
1819 
1820   compute_stabs_and_orbit_reps(nr1, nr2, 0, 0, UNDEFINED, true);
1821   return true;
1822 }
1823 
1824 // The next function is the main function for accessing the homomorphisms code.
1825 //
1826 // The arguments are:
1827 //
1828 // 1. digraph1_obj      the source/domain of the homomorphisms sought
1829 // 2. digraph2_obj      the range/codomain of the homomorphisms sought
1830 // 3. hook_obj          apply this function to every homomorphism found, or Fail
1831 //                      for no function.
1832 // 4. user_param_obj    GAP variable that can be used in the hook_obj, must be a
1833 //                      plist if hook_obj is Fail.
1834 // 5. max_results_obj   the maximum number of homomorphisms to find
1835 // 6. hint_obj          the rank of any homomorphisms found
1836 // 7. injective_obj     should be 0 for non-injective, 1 for injective, 2 for
1837 //                      embedding, or (for backwards compatibility, true for
1838 //                      injective, or false for non-injective).
1839 // 8. image_obj         a list of vertices in digraph2_obj from which the images
1840 //                      of any homomorphism will be taken.
1841 // 9. partial_map_obj   a partial map from digraph1_obj to digraph2_obj, only
1842 //                      homomorphisms extending this partial map are found (if
1843 //                      any). Can also be fail to indicate no partial mapping is
1844 //                      defined.
1845 // 10. colors1_obj      a coloring of digraph1_obj
1846 // 11. colors2_obj      a coloring of digraph2_obj, only homomorphisms that
1847 //                      respect these colorings are found.
1848 // 12. order_obj        an optional argument that can specify the order the
1849 //                      vertices in digraph1_obj should be installed in the
1850 //                      Graph or Digraph used in the recursive search. This does
1851 //                      not directly specify which order vertices are visited in
1852 //                      the search, but can have a large impact on the run time
1853 //                      of the function. It seems in many cases to be a good
1854 //                      idea for this to be the DigraphWelshPowellOrder, i.e.
1855 //                      vertices ordered from highest to lowest degree.
1856 // 13. aut_grp_obj      an optional argument that can specifiy the
1857 //                      automorphisms of the graph that will be used in the
1858 //                      recursive search. If not given, the full automorphism
1859 //                      group will be used.
1860 
FuncHomomorphismDigraphsFinder(Obj self,Obj args)1861 Obj FuncHomomorphismDigraphsFinder(Obj self, Obj args) {
1862   if (LEN_PLIST(args) != 11 && LEN_PLIST(args) != 12 && LEN_PLIST(args) != 13) {
1863     ErrorQuit(
1864         "there must be 11 or 12 arguments, found %d,", LEN_PLIST(args), 0L);
1865   }
1866   Obj digraph1_obj    = ELM_PLIST(args, 1);
1867   Obj digraph2_obj    = ELM_PLIST(args, 2);
1868   Obj hook_obj        = ELM_PLIST(args, 3);
1869   Obj user_param_obj  = ELM_PLIST(args, 4);
1870   Obj max_results_obj = ELM_PLIST(args, 5);
1871   Obj hint_obj        = ELM_PLIST(args, 6);
1872   Obj injective_obj   = ELM_PLIST(args, 7);
1873   Obj image_obj       = ELM_PLIST(args, 8);
1874   Obj partial_map_obj = ELM_PLIST(args, 9);
1875   Obj colors1_obj     = ELM_PLIST(args, 10);
1876   Obj colors2_obj     = ELM_PLIST(args, 11);
1877   Obj order_obj       = Fail;
1878   Obj aut_grp_obj     = Fail;
1879   if (LEN_PLIST(args) == 12) {
1880     if (IS_LIST(ELM_PLIST(args, 12))) {
1881       order_obj = ELM_PLIST(args, 12);
1882     } else {
1883       aut_grp_obj = ELM_PLIST(args, 12);
1884     }
1885   }
1886 
1887   if (LEN_PLIST(args) == 13) {
1888     order_obj   = ELM_PLIST(args, 12);
1889     aut_grp_obj = ELM_PLIST(args, 13);
1890   }
1891 
1892   // Validate the arguments
1893   if (CALL_1ARGS(IsDigraph, digraph1_obj) != True) {
1894     ErrorQuit("the 1st argument <digraph1> must be a digraph, not %s,",
1895               (Int) TNAM_OBJ(digraph1_obj),
1896               0L);
1897   } else if (DigraphNrVertices(digraph1_obj) > MAXVERTS) {
1898     ErrorQuit("the 1st argument <digraph1> must have at most 512 vertices, "
1899               "found %d,",
1900               DigraphNrVertices(digraph1_obj),
1901               0L);
1902   }
1903   if (CALL_1ARGS(IsDigraph, digraph2_obj) != True) {
1904     ErrorQuit("the 2nd argument <digraph2> must be a digraph, not %s,",
1905               (Int) TNAM_OBJ(digraph2_obj),
1906               0L);
1907   } else if (DigraphNrVertices(digraph2_obj) > MAXVERTS) {
1908     ErrorQuit("the 2nd argument <digraph2> must have at most 512 vertices, "
1909               "found %d,",
1910               DigraphNrVertices(digraph2_obj),
1911               0L);
1912   }
1913   if (hook_obj == Fail) {
1914     if (!IS_LIST(user_param_obj) || !IS_MUTABLE_OBJ(user_param_obj)) {
1915       ErrorQuit("the 3rd argument <hook> is fail and so the 4th argument must "
1916                 "be a mutable list, not %s,",
1917                 (Int) TNAM_OBJ(user_param_obj),
1918                 0L);
1919     }
1920   } else if (!IS_FUNC(hook_obj) || NARG_FUNC(hook_obj) != 2) {
1921     ErrorQuit(
1922         "the 3rd argument <hook> must be a function with 2 arguments,", 0L, 0L);
1923   }
1924   if (!IS_INTOBJ(max_results_obj) && max_results_obj != Infinity) {
1925     ErrorQuit("the 5th argument <max_results> must be an integer "
1926               "or infinity, not %s,",
1927               (Int) TNAM_OBJ(max_results_obj),
1928               0L);
1929   }
1930   if (!IS_INTOBJ(hint_obj) && hint_obj != Fail) {
1931     ErrorQuit("the 6th argument <hint> must be an integer "
1932               "or fail, not %s,",
1933               (Int) TNAM_OBJ(hint_obj),
1934               0L);
1935   } else if (IS_INTOBJ(hint_obj) && INT_INTOBJ(hint_obj) <= 0) {
1936     ErrorQuit("the 6th argument <hint> must be a positive integer, "
1937               "not %d,",
1938               INT_INTOBJ(hint_obj),
1939               0L);
1940   }
1941 
1942   if (!IS_INTOBJ(injective_obj) && injective_obj != True
1943       && injective_obj != False) {
1944     ErrorQuit("the 7th argument <injective> must be an integer "
1945               "or true or false, not %s,",
1946               (Int) TNAM_OBJ(injective_obj),
1947               0L);
1948   } else if (IS_INTOBJ(injective_obj)) {
1949     if (INT_INTOBJ(injective_obj) < 0 || INT_INTOBJ(injective_obj) > 2) {
1950       ErrorQuit("the 7th argument <injective> must 0, 1, or 2, not %d,",
1951                 INT_INTOBJ(injective_obj),
1952                 0L);
1953     }
1954   } else if (injective_obj == True) {
1955     injective_obj = INTOBJ_INT(1);
1956   } else {
1957     injective_obj = INTOBJ_INT(0);
1958   }
1959 
1960   if (!IS_LIST(image_obj) && image_obj != Fail) {
1961     ErrorQuit("the 8th argument <image> must be a list or fail, not %s,",
1962               (Int) TNAM_OBJ(image_obj),
1963               0L);
1964   } else if (IS_LIST(image_obj)) {
1965     for (Int i = 1; i <= LEN_LIST(image_obj); ++i) {
1966       if (!ISB_LIST(image_obj, i)) {
1967         ErrorQuit("the 8th argument <image> must be a dense list,", 0L, 0L);
1968       } else if (!IS_POS_INT(ELM_LIST(image_obj, i))) {
1969         ErrorQuit("the 8th argument <image> must only contain positive "
1970                   "integers, but found %s in position %d,",
1971                   (Int) TNAM_OBJ(ELM_LIST(image_obj, i)),
1972                   i);
1973       } else if (INT_INTOBJ(ELM_LIST(image_obj, i))
1974                  > DigraphNrVertices(digraph2_obj)) {
1975         ErrorQuit("in the 8th argument <image> position %d is out of range, "
1976                   "must be in the range [1, %d],",
1977                   i,
1978                   DigraphNrVertices(digraph2_obj));
1979       } else if (INT_INTOBJ(
1980                      POS_LIST(image_obj, ELM_LIST(image_obj, i), INTOBJ_INT(0)))
1981                  < i) {
1982         ErrorQuit(
1983             "in the 8th argument <image> position %d is a duplicate,", i, 0L);
1984       }
1985     }
1986   }
1987   if (!IS_LIST(partial_map_obj) && partial_map_obj != Fail) {
1988     ErrorQuit("the 9th argument <partial_map> must be a list or fail, not %s,",
1989               (Int) TNAM_OBJ(partial_map_obj),
1990               0L);
1991   } else if (IS_LIST(partial_map_obj)) {
1992     if (LEN_LIST(partial_map_obj) > DigraphNrVertices(digraph1_obj)) {
1993       ErrorQuit("the 9th argument <partial_map> is too long, must be at most "
1994                 "%d, found %d,",
1995                 DigraphNrVertices(digraph1_obj),
1996                 LEN_LIST(partial_map_obj));
1997     }
1998     for (Int i = 1; i <= LEN_LIST(partial_map_obj); ++i) {
1999       if (ISB_LIST(partial_map_obj, i)
2000           && POS_LIST(image_obj, ELM_LIST(partial_map_obj, i), INTOBJ_INT(0))
2001                  == Fail) {
2002         ErrorQuit("in the 9th argument <partial_map> the value %d in position "
2003                   "%d does not belong to the 7th argument <image>,",
2004                   INT_INTOBJ(ELM_LIST(partial_map_obj, i)),
2005                   i);
2006       }
2007     }
2008   }
2009 
2010   if (!IS_LIST(colors1_obj) && colors1_obj != Fail) {
2011     ErrorQuit("the 10th argument <colors1> must be a list or fail, not %s,",
2012               (Int) TNAM_OBJ(colors1_obj),
2013               0L);
2014   } else if (IS_LIST(colors1_obj)) {
2015     colors1_obj = CALL_2ARGS(DIGRAPHS_ValidateVertexColouring,
2016                              INTOBJ_INT(DigraphNrVertices(digraph1_obj)),
2017                              colors1_obj);
2018   }
2019   if (!IS_LIST(colors2_obj) && colors2_obj != Fail) {
2020     ErrorQuit("the 11th argument <colors2> must be a list or fail, not %s,",
2021               (Int) TNAM_OBJ(colors2_obj),
2022               0L);
2023   } else if (IS_LIST(colors2_obj)) {
2024     colors2_obj = CALL_2ARGS(DIGRAPHS_ValidateVertexColouring,
2025                              INTOBJ_INT(DigraphNrVertices(digraph2_obj)),
2026                              colors2_obj);
2027   }
2028   if ((IS_LIST(colors1_obj) && !IS_LIST(colors2_obj))
2029       || (colors1_obj == Fail && colors2_obj != Fail)) {
2030     ErrorQuit("the 10th and 11th arguments <colors1> and <colors2> must both "
2031               "be lists or both be fail,",
2032               0L,
2033               0L);
2034   }
2035   if (!IS_LIST(order_obj) && order_obj != Fail) {
2036     ErrorQuit("the 12th argument <order> must be a list or fail, not %s,",
2037               (Int) TNAM_OBJ(order_obj),
2038               0L);
2039   } else if (IS_LIST(order_obj)) {
2040     if (LEN_LIST(order_obj) != DigraphNrVertices(digraph1_obj)) {
2041       ErrorQuit(
2042           "the 12th argument <order> must be a list of length %d, not %d,",
2043           DigraphNrVertices(digraph1_obj),
2044           LEN_LIST(order_obj));
2045     }
2046     for (Int i = 1; i <= LEN_LIST(order_obj); ++i) {
2047       if (!ISB_LIST(order_obj, i)) {
2048         ErrorQuit("the 12th argument <order> must be a dense list, but "
2049                   "position %d is not bound,",
2050                   i,
2051                   0L);
2052       } else if (!IS_INTOBJ(ELM_LIST(order_obj, i))) {
2053         ErrorQuit("the 12th argument <order> must consist of integers, but "
2054                   "found %s in position %d,",
2055                   (Int) TNAM_OBJ(order_obj),
2056                   i);
2057       } else if (INT_INTOBJ(ELM_LIST(order_obj, i)) <= 0
2058                  || INT_INTOBJ(ELM_LIST(order_obj, i))
2059                         > DigraphNrVertices(digraph1_obj)) {
2060         ErrorQuit("the 12th argument <order> must consist of integers, in the "
2061                   "range [1, %d] but found %d,",
2062                   DigraphNrVertices(digraph1_obj),
2063                   INT_INTOBJ(ELM_LIST(order_obj, i)));
2064       } else if (INT_INTOBJ(
2065                      POS_LIST(order_obj, ELM_LIST(order_obj, i), INTOBJ_INT(0)))
2066                  < i) {
2067         ErrorQuit("the 12th argument <order> must be duplicate-free, but "
2068                   "the value %d in position %d is a duplicate,",
2069                   INT_INTOBJ(ELM_LIST(order_obj, i)),
2070                   i);
2071       }
2072     }
2073   }
2074   if (aut_grp_obj != Fail) {
2075     if (CALL_1ARGS(IsPermGroup, aut_grp_obj) != True) {
2076       ErrorQuit(
2077           "the 12th or 13th argument <aut_grp> must be a permutation group "
2078           "or fail, not %s,",
2079           (Int) TNAM_OBJ(aut_grp_obj),
2080           0L);
2081     }
2082     Obj gens = CALL_1ARGS(GeneratorsOfGroup, aut_grp_obj);
2083     DIGRAPHS_ASSERT(IS_LIST(gens));
2084     DIGRAPHS_ASSERT(LEN_LIST(gens) > 0);
2085     UInt lmp = INT_INTOBJ(CALL_1ARGS(LargestMovedPointPerms, gens));
2086     if (lmp > 0 && LEN_LIST(gens) >= lmp) {
2087       ErrorQuit("expected at most %d generators in the 12th or 13th argument "
2088                 "but got %d,",
2089                 lmp - 1,
2090                 LEN_LIST(gens));
2091     }
2092     for (UInt i = 1; i <= LEN_LIST(gens); ++i) {
2093       if (colors2_obj != Fail) {
2094         if (CALL_3ARGS(IsDigraphAutomorphism,
2095                        digraph2_obj,
2096                        ELM_LIST(gens, i),
2097                        colors2_obj)
2098             != True) {
2099           ErrorQuit("expected group of automorphisms, but found a "
2100                     "non-automorphism in position %d of the group generators,",
2101                     i,
2102                     0L);
2103         }
2104       } else {
2105         if (CALL_2ARGS(IsDigraphAutomorphism, digraph2_obj, ELM_LIST(gens, i))
2106             != True) {
2107           ErrorQuit("expected group of automorphisms, but found a "
2108                     "non-automorphism in position %d of the group generators,",
2109                     i,
2110                     0L);
2111         }
2112       }
2113     }
2114   }
2115 
2116   // Some conditions that immediately rule out there being any homomorphisms.
2117   if (((INT_INTOBJ(injective_obj) == 1 || INT_INTOBJ(injective_obj) == 2)
2118        && ((hint_obj != Fail
2119             && INT_INTOBJ(hint_obj) != DigraphNrVertices(digraph1_obj))
2120            || DigraphNrVertices(digraph1_obj)
2121                   > DigraphNrVertices(digraph2_obj)))
2122       || (hint_obj != Fail
2123           && (INT_INTOBJ(hint_obj) > DigraphNrVertices(digraph1_obj)
2124               || INT_INTOBJ(hint_obj) > DigraphNrVertices(digraph2_obj)))
2125       || LEN_LIST(image_obj) == 0
2126       || (hint_obj != Fail && INT_INTOBJ(hint_obj) > LEN_LIST(image_obj))) {
2127     // Can't print stats here because they are not initialised, also why would
2128     // we want to, nothing has actually happened yet!
2129     return user_param_obj;
2130   }
2131 
2132   // Initialise all of the global variables that are used in the recursion.
2133   // Returns false if the arguments somehow rule out there being any
2134   // homomorphisms (i.e. if injective_obj indicates that the homomorphisms
2135   // should be injective, and image_obj is too small).
2136   if (!init_data_from_args(digraph1_obj,
2137                            digraph2_obj,
2138                            hook_obj,
2139                            user_param_obj,
2140                            max_results_obj,
2141                            hint_obj,
2142                            injective_obj,
2143                            image_obj,
2144                            partial_map_obj,
2145                            colors1_obj,
2146                            colors2_obj,
2147                            order_obj,
2148                            aut_grp_obj)) {
2149 #ifdef DIGRAPHS_ENABLE_STATS
2150     print_stats(STATS);
2151 #endif
2152     return user_param_obj;
2153   }
2154 
2155   uint64_t max_results =
2156       (max_results_obj == Infinity ? SMALLINTLIMIT
2157                                    : INT_INTOBJ(max_results_obj));
2158   uint16_t hint  = (IS_INTOBJ(hint_obj) ? INT_INTOBJ(hint_obj) : UNDEFINED);
2159   uint64_t count = 0;
2160 
2161   // go!
2162   if (setjmp(OUTOFHERE) == 0) {
2163     if (CALL_1ARGS(IsSymmetricDigraph, digraph1_obj) == True
2164         && CALL_1ARGS(IsSymmetricDigraph, digraph2_obj) == True) {
2165       init_partial_map_and_find_graph_homos(
2166           partial_map_obj, max_results, hint, &count, injective_obj);
2167     } else {
2168       init_partial_map_and_find_digraph_homos(
2169           partial_map_obj, max_results, hint, &count, injective_obj);
2170     }
2171   }
2172 #ifdef DIGRAPHS_ENABLE_STATS
2173   print_stats(STATS);
2174 #endif
2175   return user_param_obj;
2176 }
2177