1 /*****************************************************************************
2 * *
3 * Vertex-invariants source file for nauty 2.7. *
4 * *
5 * Copyright (1989-2013) Brendan McKay. All rights reserved. *
6 * Subject to waivers and disclaimers in nauty.h. *
7 * *
8 * CHANGE HISTORY *
9 * 13-Mar-90 : initial release for version 1.5 *
10 * 10-Nov-90 : changes for version 1.6 : *
11 * - added dummy routine nautinv_null() *
12 * 27-Aug-92 : renamed to version 1.7, no changes to this file *
13 * 5-Jun-93 : renamed to version 1.7+, no changes to this file *
14 * 18-Aug-93 : renamed to version 1.8, no changes to this file *
15 * 17-Sep-93 : changes for version 1.9 : *
16 * - added invariant routine adjacencies() *
17 * 20-Jun-96 : changes for version 2.0 : *
18 * - added invariants cellfano() and cellfano2() *
19 * 11-Jul-96 - added dynamic allocation *
20 * 21-Oct-98 - use shortish in place of short for BIGNAUTY *
21 * 9-Jan-00 - added nautinv_check() *
22 * 12-Feb-00 - minor code formatting *
23 * 16-Nov-00 - made changes listed in nauty.h *
24 * 22-Apr-01 : changes for version 2.1 : *
25 * - made all large dynamic memory external to routines *
26 * - added nautinv_freedyn() to free all such memory *
27 * - include nautinv.h rather than naututil.h *
28 * - removed nautinv_null() *
29 * - added code to facilitate compilation into Magma *
30 * - removed EXTDEFS *
31 * 12-Jul-01 - use invararg in distances() *
32 * - fixed comments in ind and cliq routines *
33 * 21-Nov-01 : use NAUTYREQUIRED in nautinv_check() *
34 * 10-Dec-06 : remove BIGNAUTY *
35 * 10-Nov-09 : remove types shortish and permutation *
36 * 23-Nov-09 : add refinvar() *
37 * 12-Jun-10 : fixed identical errors in cellcliq() and cellind() *
38 * 15-Jan-12 : add TLS_ATTR attributes *
39 * 23-Aug-12 : fix getbigcells(), thanks to Fatih Demirkale *
40 * 23-Jan-13 : add some parens to satisfy icc *
41 * *
42 *****************************************************************************/
43
44 #define ONE_WORD_SETS
45 #include "nautinv.h"
46
47 #if MAXM==1
48 #define M 1
49 #else
50 #define M m
51 #endif
52
53 #define MASH(l,i) ((((l) ^ 056345) + (i)) & 077777)
54 /* : expression whose long value depends only on long l and int/long i.
55 Anything goes, preferably non-commutative. */
56
57 #define CLEANUP(l) ((int)((l) % 077777))
58 /* : expression whose value depends on long l and is less than 077777
59 when converted to int then short. Anything goes. */
60
61 #define ACCUM(x,y) x = (((x) + (y)) & 077777)
62 /* : must be commutative. */
63
64 static const int fuzz1[] = {037541,061532,005257,026416};
65 static const int fuzz2[] = {006532,070236,035523,062437};
66
67 #define FUZZ1(x) ((x) ^ fuzz1[(x)&3])
68 #define FUZZ2(x) ((x) ^ fuzz2[(x)&3])
69
70 #define MAXCLIQUE 10 /* max clique size for cliques() and maxindset() */
71
72 #if MAXN
73 static TLS_ATTR int workshort[MAXN+2];
74 static TLS_ATTR int vv[MAXN],ww[MAXN];
75 static TLS_ATTR int workperm[MAXN];
76 static TLS_ATTR int bucket[MAXN+2];
77 static TLS_ATTR int count[MAXN];
78 static TLS_ATTR set workset[MAXM];
79 static TLS_ATTR set w01[MAXM],w02[MAXM],w03[MAXM],w12[MAXM],w13[MAXM],w23[MAXM];
80 static TLS_ATTR set pt0[MAXM],pt1[MAXM],pt2[MAXM];
81 static TLS_ATTR set wss[MAXCLIQUE-1][MAXM];
82 static TLS_ATTR set ws1[MAXM],ws2[MAXM];
83 #else
84 DYNALLSTAT(int,workshort,workshort_sz);
85 DYNALLSTAT(int,vv,vv_sz);
86 DYNALLSTAT(int,ww,ww_sz);
87 DYNALLSTAT(int,workperm,workperm_sz);
88 DYNALLSTAT(int,bucket,bucket_sz);
89 DYNALLSTAT(int,count,count_sz);
90 DYNALLSTAT(set,ws1,ws1_sz);
91 DYNALLSTAT(set,ws2,ws2_sz);
92 DYNALLSTAT(set,workset,workset_sz);
93 DYNALLSTAT(set,w01,w01_sz);
94 DYNALLSTAT(set,w02,w02_sz);
95 DYNALLSTAT(set,w03,w03_sz);
96 DYNALLSTAT(set,w12,w12_sz);
97 DYNALLSTAT(set,w13,w13_sz);
98 DYNALLSTAT(set,w23,w23_sz);
99 DYNALLSTAT(set,pt0,pt0_sz);
100 DYNALLSTAT(set,pt1,pt1_sz);
101 DYNALLSTAT(set,pt2,pt2_sz);
102 DYNALLSTAT(set,wss,wss_sz);
103 #endif
104
105 /* aproto: header new_nauty_protos.h */
106
107 /*****************************************************************************
108 * *
109 * This file contains a number of procedures which compute vertex-invariants *
110 * for stronger partition refinement. Since entirely different *
111 * vertex-invariants seem to work better for different types of graph, we *
112 * cannot do more than give a small collection of representative examples. *
113 * Any serious computations with difficult graphs may well need to use *
114 * specially-written vertex-invariants. The use of vertex-invariants *
115 * procedures is supported by nauty from version 1.5 onwards, via the *
116 * options userinvarproc, mininvarlevel, maxinvarlevel and invararg. *
117 * The meaning of these fields in detail are as follows: *
118 * userinvarproc is the address of the vertex-invariant procedure. If *
119 * no vertex-invariants is required, this field should *
120 * have the value NULL. *
121 * maxinvarlevel The absolute value of this is the maximum level in the *
122 * search tree at which the vertex-invariant will be *
123 * computed. The root of the tree is at level 1, so the *
124 * vertex-invariant will not be invoked at all if *
125 * maxinvarlevel==0. Negative values of maxinvarlevel *
126 * request nauty to not compute the vertex-invariant at *
127 * a level greater than that of the earliest node (if any) *
128 * on the path to the first leaf of the search tree at *
129 * which the vertex-invariant refines the partition. *
130 * mininvarlevel The absolute value of this is the minimum level in the *
131 * search tree at which the vertex-invariant will be *
132 * computed. The root of the tree is at level 1, so there *
133 * is no effective limit if mininvarlevel is -1, 0 or 1. *
134 * Negative values of mininvarlevel request nauty to not *
135 * compute the vertex-invariant at a level less than *
136 * that of the earliest node (if any) on the path to the *
137 * first leaf of the search tree at which the *
138 * vertex-invariant refines the partition. *
139 * invararg is passed to the vertex-invariant procedure via the *
140 * argument of the same name. It can be used by the *
141 * procedure for any purpose. *
142 * Note that negative values of maxinvarlevel and mininvarlevel make the *
143 * canonical labelling invalid, but can speed automorphism group finding. *
144 * Nauty already knows this and takes their absolute values. *
145 * *
146 * A vertex-invariant must be declared thus: *
147 * void invarproc(g,lab,ptn,level,numcells,tvpos,invar,invararg,digraph,m,n) *
148 * All of these arguments must be treated as read-only except for invar. *
149 * g : the graph, exactly as passed to nauty() *
150 * lab,ptn : the current partition nest (see nauty.h for the format) *
151 * level : the level of this node in the search tree. *
152 * numcells : the number of cells in the partition at this node. *
153 * tvpos : the index in (lab,ptn) of one cell in the partition. *
154 * If level <= 1, the cell will be the first fragment of the *
155 * first active cell (as provided by the initial call to nauty), *
156 * or the first cell, if there were no active cells. *
157 * If level > 1, the cell will be the singleton cell which was *
158 * created to make this node of the search tree from its parent. *
159 * invararg : a copy of options.invararg *
160 * digraph : a copy of options.digraph *
161 * m,n : size parameters as passed to nauty() *
162 * invar : an array to return the answer in. The procedure must put in *
163 * each invar[i] (0 <= i < n) an invariant of the 6-tuple *
164 * (<vertex i>,g,<the partition nest to this level>,level, *
165 * invararg,digraph) *
166 * Note that invar[] is declared as an int array. Since the *
167 * absolute value of the invariant is irrelevant, only the *
168 * comparative values, any short, int or long value can be *
169 * assigned to the entries of invar[] without fear. However, *
170 * you should assign a value less than 077777 to ensure machine- *
171 * independence of the canonical labelling. *
172 * *
173 * The refinement procedure has already been called before the invariant *
174 * procedure is called. That means that the partition is equitable if *
175 * digraph==FALSE. *
176 * *
177 *****************************************************************************/
178
179 /*****************************************************************************
180 * *
181 * twopaths() assigns to each vertex v the sum of the weights of each vertex *
182 * which can be reached from v along a walk of length two (including itself *
183 * usually). The weight of each vertex w is defined as the ordinal number *
184 * of the cell containing w, starting at 1 for the first cell. *
185 * *
186 *****************************************************************************/
187
188 void
twopaths(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)189 twopaths(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
190 int *invar, int invararg, boolean digraph, int m, int n)
191 {
192 int i,v,w;
193 int wt;
194 set *gv,*gw;
195
196 #if !MAXN
197 DYNALLOC1(set,workset,workset_sz,m,"twopaths");
198 DYNALLOC1(int,workshort,workshort_sz,n+2,"twopaths");
199 #endif
200
201 wt = 1;
202 for (i = 0; i < n; ++i)
203 {
204 workshort[lab[i]] = wt;
205 if (ptn[i] <= level) ++wt;
206 }
207
208 for (v = 0, gv = (set*)g; v < n; ++v, gv += M)
209 {
210 EMPTYSET(workset,m);
211 w = -1;
212 while ((w = nextelement(gv,M,w)) >= 0)
213 {
214 gw = GRAPHROW(g,w,m);
215 for (i = M; --i >= 0;) UNION(workset[i],gw[i]);
216 }
217 wt = 0;
218 w = -1;
219 while ((w = nextelement(workset,M,w)) >= 0) ACCUM(wt,workshort[w]);
220 invar[v] = wt;
221 }
222 }
223
224 /*****************************************************************************
225 * *
226 * quadruples() assigns to each vertex v a value depending on the set of *
227 * weights w(v,v1,v2,v3), where w(v,v1,v2,v3) depends on the number of *
228 * vertices adjacent to an odd number of {v,v1,v2,v3}, and to the cells *
229 * that v,v1,v2,v3 belong to. {v,v1,v2,v3} are permitted to range over all *
230 * distinct 4-tuples which contain at least one member in the cell tvpos. *
231 * *
232 *****************************************************************************/
233
234 void
quadruples(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)235 quadruples(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
236 int *invar, int invararg, boolean digraph, int m, int n)
237 {
238 int i,pc;
239 setword sw;
240 set *gw;
241 int wt;
242 int v,iv,v1,v2,v3;
243 set *gv;
244 long wv,wv1,wv2,wv3;
245
246 #if !MAXN
247 DYNALLOC1(int,workshort,workshort_sz,n+2,"quadruples");
248 DYNALLOC1(set,ws1,ws1_sz,m,"quadruples");
249 DYNALLOC1(set,workset,workset_sz,m,"quadruples");
250 #endif
251
252 for (i = n; --i >= 0;) invar[i] = 0;
253
254 wt = 1;
255 for (i = 0; i < n; ++i)
256 {
257 workshort[lab[i]] = FUZZ2(wt);
258 if (ptn[i] <= level) ++wt;
259 }
260
261 iv = tvpos - 1;
262 do
263 {
264 v = lab[++iv];
265 gv = GRAPHROW(g,v,m);
266 wv = workshort[v];
267 for (v1 = 0; v1 < n-2; ++v1)
268 {
269 wv1 = workshort[v1];
270 if (wv1 == wv && v1 <= v) continue;
271 wv1 += wv;
272 gw = GRAPHROW(g,v1,m);
273 for (i = M; --i >= 0;) workset[i] = gv[i] ^ gw[i];
274 for (v2 = v1+1; v2 < n-1; ++v2)
275 {
276 wv2 = workshort[v2];
277 if (wv2 == wv && v2 <= v) continue;
278 wv2 += wv1;
279 gw = GRAPHROW(g,v2,m);
280 for (i = M; --i >= 0;) ws1[i] = workset[i] ^ gw[i];
281 for (v3 = v2+1; v3 < n; ++v3)
282 {
283 wv3 = workshort[v3];
284 if (wv3 == wv && v3 <= v) continue;
285 wv3 += wv2;
286 gw = GRAPHROW(g,v3,m);
287 pc = 0;
288 for (i = M; --i >= 0;)
289 if ((sw = ws1[i] ^ gw[i]) != 0) pc += POPCOUNT(sw);
290 wt = (FUZZ1(pc)+wv3) & 077777;
291 wt = FUZZ2(wt);
292 ACCUM(invar[v],wt);
293 ACCUM(invar[v1],wt);
294 ACCUM(invar[v2],wt);
295 ACCUM(invar[v3],wt);
296 }
297 }
298 }
299 }
300 while (ptn[iv] > level);
301 }
302
303 /*****************************************************************************
304 * *
305 * triples() assigns to each vertex v a value depending on the set of *
306 * weights w(v,v1,v2), where w(v,v1,v2) depends on the number of vertices *
307 * adjacent to an odd number of {v,v1,v2}, and to the cells that *
308 * v,v1,v2 belong to. {v,v1,v2} are permitted to range over all distinct *
309 * triples which contain at least one member in the cell tvpos. *
310 * *
311 *****************************************************************************/
312
313 void
triples(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)314 triples(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
315 int *invar, int invararg, boolean digraph, int m, int n)
316 {
317 int i,pc;
318 setword sw;
319 set *gw;
320 int wt;
321 int v,iv,v1,v2;
322 set *gv;
323 long wv,wv1,wv2;
324
325 #if !MAXN
326 DYNALLOC1(set,workset,workset_sz,m,"triples");
327 DYNALLOC1(int,workshort,workshort_sz,n+2,"triples");
328 #endif
329
330 for (i = n; --i >= 0;) invar[i] = 0;
331
332 wt = 1;
333 for (i = 0; i < n; ++i)
334 {
335 workshort[lab[i]] = FUZZ1(wt);
336 if (ptn[i] <= level) ++wt;
337 }
338
339 iv = tvpos - 1;
340 do
341 {
342 v = lab[++iv];
343 gv = GRAPHROW(g,v,m);
344 wv = workshort[v];
345 for (v1 = 0; v1 < n-1; ++v1)
346 {
347 wv1 = workshort[v1];
348 if (wv1 == wv && v1 <= v) continue;
349 wv1 += wv;
350 gw = GRAPHROW(g,v1,m);
351 for (i = M; --i >= 0;) workset[i] = gv[i] ^ gw[i];
352 for (v2 = v1+1; v2 < n; ++v2)
353 {
354 wv2 = workshort[v2];
355 if (wv2 == wv && v2 <= v) continue;
356 wv2 += wv1;
357 gw = GRAPHROW(g,v2,m);
358 pc = 0;
359 for (i = M; --i >= 0;)
360 if ((sw = workset[i] ^ gw[i]) != 0) pc += POPCOUNT(sw);
361 wt = (FUZZ1(pc)+wv2) & 077777;
362 wt = FUZZ2(wt);
363 ACCUM(invar[v],wt);
364 ACCUM(invar[v1],wt);
365 ACCUM(invar[v2],wt);
366 }
367 }
368 }
369 while (ptn[iv] > level);
370 }
371
372 /*****************************************************************************
373 * *
374 * adjtriang() assigns to each vertex v a value depending on the numbers *
375 * of common neighbours between each pair {v1,v2} of neighbours of v, and *
376 * which cells v1 and v2 lie in. The vertices v1 and v2 must be adjacent *
377 * if invararg == 0 and not adjacent if invararg == 1. *
378 * *
379 *****************************************************************************/
380
381 void
adjtriang(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)382 adjtriang(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
383 int *invar, int invararg, boolean digraph, int m, int n)
384 {
385 int j,pc;
386 setword sw;
387 set *gi;
388 int wt;
389 int i,v1,v2;
390 boolean v1v2;
391 set *gv1,*gv2;
392
393 #if !MAXN
394 DYNALLOC1(set,workset,workset_sz,m,"adjtriang");
395 DYNALLOC1(int,workshort,workshort_sz,n+2,"adjtriang");
396 #endif
397
398 for (i = n; --i >= 0;) invar[i] = 0;
399
400 wt = 1;
401 for (i = 0; i < n; ++i)
402 {
403 workshort[lab[i]] = FUZZ1(wt);
404 if (ptn[i] <= level) ++wt;
405 }
406
407 for (v1 = 0, gv1 = g; v1 < n; ++v1, gv1 += M)
408 {
409 for (v2 = (digraph ? 0 : v1+1); v2 < n; ++v2)
410 {
411 if (v2 == v1) continue;
412 v1v2 = (ISELEMENT(gv1,v2) != 0);
413 if ((invararg == 0 && !v1v2)
414 || (invararg == 1 && v1v2)) continue;
415 wt = workshort[v1];
416 ACCUM(wt,workshort[v2]);
417 ACCUM(wt,v1v2);
418
419 gv2 = GRAPHROW(g,v2,m);
420 for (i = M; --i >= 0;) workset[i] = gv1[i] & gv2[i];
421 i = -1;
422 while ((i = nextelement(workset,M,i)) >= 0)
423 {
424 pc = 0;
425 gi = GRAPHROW(g,i,m);
426 for (j = M; --j >= 0;)
427 if ((sw = workset[j] & gi[j]) != 0) pc += POPCOUNT(sw);
428 pc = (pc + wt) & 077777;
429 ACCUM(invar[i],pc);
430 }
431 }
432 }
433 }
434
435 /*****************************************************************************
436 * *
437 * getbigcells(ptn,level,minsize,bigcells,cellstart,cellsize,n) is an *
438 * auxiliary procedure to make a list of all the large cells in the current *
439 * partition. On entry, ptn, level and n have their usual meanings, *
440 * while minsize is the smallest size of an interesting cell. On return, *
441 * bigcells is the number of cells of size at least minsize, cellstart[0...] *
442 * contains their starting positions in ptn, and cellsize[0...] contain *
443 * their sizes. These two arrays are in increasing order of cell size, *
444 * then position. *
445 * *
446 *****************************************************************************/
447
448 void
getbigcells(int * ptn,int level,int minsize,int * bigcells,int * cellstart,int * cellsize,int n)449 getbigcells(int *ptn, int level, int minsize, int *bigcells,
450 int *cellstart, int *cellsize, int n)
451 {
452 int cell1,cell2,j;
453 int si,st;
454 int bc,i,h;
455
456 bc = 0;
457 for (cell1 = 0; cell1 < n; cell1 = cell2 + 1)
458 {
459 for (cell2 = cell1; ptn[cell2] > level; ++cell2) {}
460
461 if (cell2 >= cell1 + minsize - 1)
462 {
463 cellstart[bc] = cell1;
464 cellsize[bc] = cell2 - cell1 + 1;
465 ++bc;
466 }
467 }
468 *bigcells = bc;
469
470 j = bc / 3;
471 h = 1;
472 do
473 h = 3 * h + 1;
474 while (h < j);
475
476 do /* shell sort */
477 {
478 for (i = h; i < bc; ++i)
479 {
480 st = cellstart[i];
481 si = cellsize[i];
482 for (j = i; cellsize[j-h] > si ||
483 (cellsize[j-h] == si && cellstart[j-h] > st); )
484 {
485 cellsize[j] = cellsize[j-h];
486 cellstart[j] = cellstart[j-h];
487 if ((j -= h) < h) break;
488 }
489 cellsize[j] = si;
490 cellstart[j] = st;
491 }
492 h /= 3;
493 }
494 while (h > 0);
495 }
496
497 /*****************************************************************************
498 * *
499 * celltrips() assigns to each vertex v a value depending on the set of *
500 * weights w(v,v1,v2), where w(v,v1,v2) depends on the number of vertices *
501 * adjacent to an odd number of {v,v1,v2}. {v,v1,v2} are constrained to *
502 * belong to the same cell. We try the cells in increasing order of size, *
503 * and stop as soon as any cell splits. *
504 * *
505 *****************************************************************************/
506
507 void
celltrips(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)508 celltrips(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
509 int *invar, int invararg, boolean digraph, int m, int n)
510 {
511 int i,pc;
512 setword sw;
513 set *gw;
514 int wt;
515 int v,iv,v1,iv1,v2,iv2;
516 int icell,bigcells,cell1,cell2;
517 int *cellstart,*cellsize;
518 set *gv;
519
520 #if !MAXN
521 DYNALLOC1(set,workset,workset_sz,m,"celltrips");
522 DYNALLOC1(int,workshort,workshort_sz,n+2,"celltrips");
523 #endif
524
525 for (i = n; --i >= 0;) invar[i] = 0;
526
527 cellstart = workshort;
528 cellsize = workshort + (n/2);
529 getbigcells(ptn,level,3,&bigcells,cellstart,cellsize,n);
530
531 for (icell = 0; icell < bigcells; ++icell)
532 {
533 cell1 = cellstart[icell];
534 cell2 = cell1 + cellsize[icell] - 1;
535 for (iv = cell1; iv <= cell2 - 2; ++iv)
536 {
537 v = lab[iv];
538 gv = GRAPHROW(g,v,m);
539 for (iv1 = iv + 1; iv1 <= cell2 - 1; ++iv1)
540 {
541 v1 = lab[iv1];
542 gw = GRAPHROW(g,v1,m);
543 for (i = M; --i >= 0;) workset[i] = gv[i] ^ gw[i];
544 for (iv2 = iv1 + 1; iv2 <= cell2; ++iv2)
545 {
546 v2 = lab[iv2];
547 gw = GRAPHROW(g,v2,m);
548 pc = 0;
549 for (i = M; --i >= 0;)
550 if ((sw = workset[i] ^ gw[i]) != 0)
551 pc += POPCOUNT(sw);
552 wt = FUZZ1(pc);
553 ACCUM(invar[v],wt);
554 ACCUM(invar[v1],wt);
555 ACCUM(invar[v2],wt);
556 }
557 }
558 }
559 wt = invar[lab[cell1]];
560 for (i = cell1 + 1; i <= cell2; ++i)
561 if (invar[lab[i]] != wt) return;
562 }
563 }
564
565 /*****************************************************************************
566 * *
567 * cellquads() assigns to each vertex v a value depending on the set of *
568 * weights w(v,v1,v2,v3), where w(v,v1,v2,v3) depends on the number of *
569 * vertices adjacent to an odd number of {v,v1,v2,v3}. {v,v1,v2,v3} are *
570 * constrained to belong to the same cell. We try the cells in increasing *
571 * order of size, and stop as soon as any cell splits. *
572 * *
573 *****************************************************************************/
574
575 void
cellquads(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)576 cellquads(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
577 int *invar, int invararg, boolean digraph, int m, int n)
578 {
579 int i,pc;
580 setword sw;
581 set *gw;
582 int wt;
583 int v,iv,v1,iv1,v2,iv2,v3,iv3;
584 int icell,bigcells,cell1,cell2;
585 int *cellstart,*cellsize;
586 set *gv;
587
588 #if !MAXN
589 DYNALLOC1(set,workset,workset_sz,m,"cellquads");
590 DYNALLOC1(int,workshort,workshort_sz,n+2,"cellquads");
591 DYNALLOC1(set,ws1,ws1_sz,m,"cellquads");
592 #endif
593
594 for (i = n; --i >= 0;) invar[i] = 0;
595
596 cellstart = workshort;
597 cellsize = workshort + (n/2);
598 getbigcells(ptn,level,4,&bigcells,cellstart,cellsize,n);
599
600 for (icell = 0; icell < bigcells; ++icell)
601 {
602 cell1 = cellstart[icell];
603 cell2 = cell1 + cellsize[icell] - 1;
604 for (iv = cell1; iv <= cell2 - 3; ++iv)
605 {
606 v = lab[iv];
607 gv = GRAPHROW(g,v,m);
608 for (iv1 = iv + 1; iv1 <= cell2 - 2; ++iv1)
609 {
610 v1 = lab[iv1];
611 gw = GRAPHROW(g,v1,m);
612 for (i = M; --i >= 0;) workset[i] = gv[i] ^ gw[i];
613 for (iv2 = iv1 + 1; iv2 <= cell2 - 1; ++iv2)
614 {
615 v2 = lab[iv2];
616 gw = GRAPHROW(g,v2,m);
617 for (i = M; --i >= 0;) ws1[i] = workset[i] ^ gw[i];
618 for (iv3 = iv2 + 1; iv3 <= cell2; ++iv3)
619 {
620 v3 = lab[iv3];
621 gw = GRAPHROW(g,v3,m);
622 pc = 0;
623 for (i = M; --i >= 0;)
624 if ((sw = ws1[i] ^ gw[i]) != 0)
625 pc += POPCOUNT(sw);
626 wt = FUZZ1(pc);
627 ACCUM(invar[v],wt);
628 ACCUM(invar[v1],wt);
629 ACCUM(invar[v2],wt);
630 ACCUM(invar[v3],wt);
631 }
632 }
633 }
634 }
635 wt = invar[lab[cell1]];
636 for (i = cell1 + 1; i <= cell2; ++i)
637 if (invar[lab[i]] != wt) return;
638 }
639 }
640
641 /*****************************************************************************
642 * *
643 * cellquins() assigns to each vertex v a value depending on the set of *
644 * weights w(v,v1,v2,v3,v4), where w(v,v1,v2,v3,v4) depends on the number *
645 * of vertices adjacent to an odd number of {v,v1,v2,v3,v4}. *
646 * {v,v1,v2,v3,v4} are constrained to belong to the same cell. We try the *
647 * cells in increasing order of size, and stop as soon as any cell splits. *
648 * *
649 *****************************************************************************/
650
651 void
cellquins(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)652 cellquins(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
653 int *invar, int invararg, boolean digraph, int m, int n)
654 {
655 int i,pc;
656 setword sw;
657 set *gw;
658 int wt;
659 int v,iv,v1,iv1,v2,iv2,v3,iv3,v4,iv4;
660 int icell,bigcells,cell1,cell2;
661 int *cellstart,*cellsize;
662 set *gv;
663
664 #if !MAXN
665 DYNALLOC1(set,workset,workset_sz,m,"cellquins");
666 DYNALLOC1(int,workshort,workshort_sz,n+2,"cellquins");
667 DYNALLOC1(set,ws1,ws1_sz,m,"cellquins");
668 DYNALLOC1(set,ws2,ws2_sz,m,"cellquins");
669 #endif
670
671 for (i = n; --i >= 0;) invar[i] = 0;
672
673 cellstart = workshort;
674 cellsize = workshort + (n/2);
675 getbigcells(ptn,level,5,&bigcells,cellstart,cellsize,n);
676
677 for (icell = 0; icell < bigcells; ++icell)
678 {
679 cell1 = cellstart[icell];
680 cell2 = cell1 + cellsize[icell] - 1;
681 for (iv = cell1; iv <= cell2 - 4; ++iv)
682 {
683 v = lab[iv];
684 gv = GRAPHROW(g,v,m);
685 for (iv1 = iv + 1; iv1 <= cell2 - 3; ++iv1)
686 {
687 v1 = lab[iv1];
688 gw = GRAPHROW(g,v1,m);
689 for (i = M; --i >= 0;) workset[i] = gv[i] ^ gw[i];
690 for (iv2 = iv1 + 1; iv2 <= cell2 - 2; ++iv2)
691 {
692 v2 = lab[iv2];
693 gw = GRAPHROW(g,v2,m);
694 for (i = M; --i >= 0;) ws1[i] = workset[i] ^ gw[i];
695 for (iv3 = iv2 + 1; iv3 <= cell2 - 1; ++iv3)
696 {
697 v3 = lab[iv3];
698 gw = GRAPHROW(g,v3,m);
699 for (i = M; --i >= 0;) ws2[i] = ws1[i] ^ gw[i];
700 for (iv4 = iv3 + 1; iv4 <= cell2; ++iv4)
701 {
702 v4 = lab[iv4];
703 gw = GRAPHROW(g,v4,m);
704 pc = 0;
705 for (i = M; --i >= 0;)
706 if ((sw = ws2[i] ^ gw[i]) != 0)
707 pc += POPCOUNT(sw);
708 wt = FUZZ1(pc);
709 ACCUM(invar[v],wt);
710 ACCUM(invar[v1],wt);
711 ACCUM(invar[v2],wt);
712 ACCUM(invar[v3],wt);
713 ACCUM(invar[v4],wt);
714 }
715 }
716 }
717 }
718 }
719 wt = invar[lab[cell1]];
720 for (i = cell1 + 1; i <= cell2; ++i)
721 if (invar[lab[i]] != wt) return;
722 }
723 }
724
725 /*****************************************************************************
726 * *
727 * uniqinter(s1,s2,m) returns the number in both sets if it is unique, *
728 * or -1 if there is none or it is not unique. *
729 *****************************************************************************/
730
731 static int
uniqinter(set * s1,set * s2,int m)732 uniqinter(set *s1, set *s2, int m)
733 {
734 int i,j;
735 setword w;
736
737 for (i = 0; i < M; ++i)
738 {
739 if ((w = s1[i] & s2[i]) != 0)
740 {
741 j = FIRSTBITNZ(w);
742 if (w != BITT[j]) return -1;
743 j += TIMESWORDSIZE(i);
744 for (++i; i < M; ++i)
745 if (s1[i] & s2[i]) return -1;
746 return j;
747 }
748 }
749 return -1;
750 }
751
752 /*****************************************************************************
753 * *
754 * cellfano2() assigns to each vertex v a value depending on the set of *
755 * weights w(v,v1,v2,v3), where w(v,v1,v2,v3) depends on the number of *
756 * fano-plane analogues containing {v,v1,v2,v3}. {v,v1,v2,v3} are *
757 * constrained to belong to the same cell and being independent and *
758 * non-collinear. We try the cells in increasing order of size, and stop *
759 * as soon as any cell splits. *
760 * *
761 *****************************************************************************/
762
763 void
cellfano2(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)764 cellfano2(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
765 int *invar, int invararg, boolean digraph, int m, int n)
766 {
767 int i,pc;
768 setword sw;
769 int wt;
770 int v0,v1,v2,v3,iv0,iv1,iv2,iv3;
771 int icell,bigcells,cell1,cell2;
772 int *cellstart,*cellsize;
773 int nw,x01,x02,x03,x12,x13,x23;
774 int pnt0,pnt1,pnt2;
775 set *gv0,*gv1,*gv2,*gv3;
776 set *gp0,*gp1,*gp2;
777
778 #if !MAXN
779 DYNALLOC1(int,workshort,workshort_sz,n+2,"cellfano2");
780 DYNALLOC1(int,vv,vv_sz,n,"cellfano2");
781 DYNALLOC1(int,ww,ww_sz,n,"cellfano2");
782 #endif
783
784 for (i = n; --i >= 0;) invar[i] = 0;
785
786 cellstart = workshort;
787 cellsize = workshort + (n/2);
788 getbigcells(ptn,level,4,&bigcells,cellstart,cellsize,n);
789
790 for (icell = 0; icell < bigcells; ++icell)
791 {
792 cell1 = cellstart[icell];
793 cell2 = cell1 + cellsize[icell] - 1;
794 for (iv0 = cell1; iv0 <= cell2 - 3; ++iv0)
795 {
796 v0 = lab[iv0];
797 gv0 = GRAPHROW(g,v0,m);
798 nw = 0;
799 for (iv1 = iv0 + 1; iv1 <= cell2; ++iv1)
800 {
801 v1 = lab[iv1];
802 if (ISELEMENT(gv0,v1)) continue;
803 if ((x01 = uniqinter(gv0,GRAPHROW(g,v1,m),m)) < 0) continue;
804 vv[nw] = v1;
805 ww[nw] = x01;
806 ++nw;
807 }
808
809 for (iv1 = 0; iv1 < nw-2; ++iv1)
810 {
811 v1 = vv[iv1];
812 gv1 = GRAPHROW(g,v1,m);
813 x01 = ww[iv1];
814
815 for (iv2 = iv1 + 1; iv2 < nw-1; ++iv2)
816 {
817 x02 = ww[iv2];
818 if (x02 == x01) continue;
819 v2 = vv[iv2];
820 if (ISELEMENT(gv1,v2)) continue;
821 gv2 = GRAPHROW(g,v2,m);
822 if ((x12 = uniqinter(gv1,gv2,m)) < 0) continue;
823
824 for (iv3 = iv2 + 1; iv3 < nw; ++iv3)
825 {
826 x03 = ww[iv3];
827 if (x03 == x01 || x03 == x02) continue;
828 v3 = vv[iv3];
829 if (ISELEMENT(gv1,v3) || ISELEMENT(gv2,v3))
830 continue;
831 gv3 = GRAPHROW(g,v3,m);
832 if ((x13 = uniqinter(gv1,gv3,m)) < 0) continue;
833 if ((x23 = uniqinter(gv2,gv3,m)) < 0
834 || x23 == x13) continue;
835
836 if ((pnt0 = uniqinter(GRAPHROW(g,x01,m),
837 GRAPHROW(g,x23,m),m)) < 0)
838 continue;
839 if ((pnt1 = uniqinter(GRAPHROW(g,x02,m),
840 GRAPHROW(g,x13,m),m)) < 0)
841 continue;
842 if ((pnt2 = uniqinter(GRAPHROW(g,x03,m),
843 GRAPHROW(g,x12,m),m)) < 0)
844 continue;
845
846 gp0 = GRAPHROW(g,pnt0,m);
847 gp1 = GRAPHROW(g,pnt1,m);
848 gp2 = GRAPHROW(g,pnt2,m);
849
850 pc = 0;
851 for (i = M; --i >= 0;)
852 {
853 sw = gp0[i] & gp1[i] & gp2[i];
854 if (sw) pc += POPCOUNT(sw);
855 }
856 wt = FUZZ1(pc);
857 ACCUM(invar[v0],wt);
858 ACCUM(invar[v1],wt);
859 ACCUM(invar[v2],wt);
860 ACCUM(invar[v3],wt);
861 }
862 }
863 }
864 }
865 wt = invar[lab[cell1]];
866 for (i = cell1 + 1; i <= cell2; ++i)
867 if (invar[lab[i]] != wt) return;
868 }
869 }
870
871 /*****************************************************************************
872 * *
873 * setnbhd(g,m,n,w,wn) is an auxiliary routine that sets wn to the union *
874 * of the neighbours of the vertices in w. *
875 * *
876 *****************************************************************************/
877
878 void
setnbhd(graph * g,int m,int n,set * w,set * wn)879 setnbhd(graph *g, int m, int n, set *w, set *wn)
880 {
881 int i,j;
882 set *gi;
883
884 i = nextelement(w,M,-1);
885 if (i < 0)
886 {
887 EMPTYSET(wn,M);
888 return;
889 }
890
891 gi = GRAPHROW(g,i,M);
892 for (j = M; --j >= 0;) wn[j] = gi[j];
893
894 while ((i = nextelement(w,M,i)) >= 0)
895 {
896 gi = GRAPHROW(g,i,M);
897 for (j = M; --j >= 0;) wn[j] |= gi[j];
898 }
899 }
900
901 /*****************************************************************************
902 * *
903 * cellfano() assigns to each vertex v a value depending on the set of *
904 * weights w(v,v1,v2,v3), where w(v,v1,v2,v3) depends on the number of *
905 * fano-plane analogues containing {v,v1,v2,v3}. {v,v1,v2,v3} are *
906 * constrained to belong to the same cell and being independent. We try *
907 * the cells in increasing order of size, and stop as soon as any cell *
908 * splits. *
909 * *
910 *****************************************************************************/
911
912 void
cellfano(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)913 cellfano(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
914 int *invar, int invararg, boolean digraph, int m, int n)
915 {
916 int i,pc;
917 setword sw;
918 int wt;
919 int v0,v1,v2,v3,iv0,iv1,iv2,iv3;
920 int icell,bigcells,cell1,cell2;
921 int *cellstart,*cellsize;
922 set *gv0,*gv1,*gv2,*gv3;
923
924 #if !MAXN
925 DYNALLOC1(int,workshort,workshort_sz,n+2,"cellfano");
926 DYNALLOC1(set,w01,w01_sz,m,"cellfano");
927 DYNALLOC1(set,w02,w02_sz,m,"cellfano");
928 DYNALLOC1(set,w03,w03_sz,m,"cellfano");
929 DYNALLOC1(set,w12,w12_sz,m,"cellfano");
930 DYNALLOC1(set,w13,w13_sz,m,"cellfano");
931 DYNALLOC1(set,w23,w23_sz,m,"cellfano");
932 DYNALLOC1(set,pt0,pt0_sz,m,"cellfano");
933 DYNALLOC1(set,pt1,pt1_sz,m,"cellfano");
934 DYNALLOC1(set,pt2,pt2_sz,m,"cellfano");
935 DYNALLOC1(set,workset,workset_sz,m,"cellfano");
936 #else
937 #endif
938
939 for (i = n; --i >= 0;) invar[i] = 0;
940
941 cellstart = workshort;
942 cellsize = workshort + (n/2);
943 getbigcells(ptn,level,4,&bigcells,cellstart,cellsize,n);
944
945 for (icell = 0; icell < bigcells; ++icell)
946 {
947 cell1 = cellstart[icell];
948 cell2 = cell1 + cellsize[icell] - 1;
949 for (iv0 = cell1; iv0 <= cell2 - 3; ++iv0)
950 {
951 v0 = lab[iv0];
952 gv0 = GRAPHROW(g,v0,m);
953 for (iv1 = iv0 + 1; iv1 <= cell2 - 2; ++iv1)
954 {
955 v1 = lab[iv1];
956 if (ISELEMENT(gv0,v1)) continue;
957 gv1 = GRAPHROW(g,v1,m);
958 for (i = M; --i >= 0;) workset[i] = gv0[i] & gv1[i];
959 setnbhd(g,m,n,workset,w01);
960
961 for (iv2 = iv1 + 1; iv2 <= cell2 - 1; ++iv2)
962 {
963 v2 = lab[iv2];
964 if (ISELEMENT(gv0,v2) || ISELEMENT(gv1,v2))
965 continue;
966 gv2 = GRAPHROW(g,v2,m);
967 for (i = M; --i >= 0;) workset[i] = gv0[i] & gv2[i];
968 setnbhd(g,m,n,workset,w02);
969 for (i = M; --i >= 0;) workset[i] = gv1[i] & gv2[i];
970 setnbhd(g,m,n,workset,w12);
971
972 for (iv3 = iv2 + 1; iv3 <= cell2; ++iv3)
973 {
974 v3 = lab[iv3];
975 if (ISELEMENT(gv0,v3) || ISELEMENT(gv1,v3) ||
976 ISELEMENT(gv2,v3))
977 continue;
978 gv3 = GRAPHROW(g,v3,m);
979 for (i = M; --i >= 0;) workset[i] = gv0[i] & gv3[i];
980 setnbhd(g,m,n,workset,w03);
981 for (i = M; --i >= 0;) workset[i] = gv1[i] & gv3[i];
982 setnbhd(g,m,n,workset,w13);
983 for (i = M; --i >= 0;) workset[i] = gv2[i] & gv3[i];
984 setnbhd(g,m,n,workset,w23);
985
986 for (i = M; --i >= 0;) workset[i] = w01[i] & w23[i];
987 setnbhd(g,m,n,workset,pt0);
988 for (i = M; --i >= 0;) workset[i] = w03[i] & w12[i];
989 setnbhd(g,m,n,workset,pt1);
990 for (i = M; --i >= 0;) workset[i] = w02[i] & w13[i];
991 setnbhd(g,m,n,workset,pt2);
992 pc = 0;
993 for (i = M; --i >= 0;)
994 {
995 sw = pt0[i] & pt1[i] & pt2[i];
996 if (sw) pc += POPCOUNT(sw);
997 }
998 wt = FUZZ1(pc);
999 ACCUM(invar[v0],wt);
1000 ACCUM(invar[v1],wt);
1001 ACCUM(invar[v2],wt);
1002 ACCUM(invar[v3],wt);
1003 }
1004 }
1005 }
1006 }
1007 wt = invar[lab[cell1]];
1008 for (i = cell1 + 1; i <= cell2; ++i)
1009 if (invar[lab[i]] != wt) return;
1010 }
1011 }
1012
1013 /*****************************************************************************
1014 * *
1015 * distances() assigns to each vertex v a value depending on the number of *
1016 * vertices at each distance from v, and what cells they lie in. *
1017 * If we find any cell which is split in this manner, we don't try any *
1018 * further cells. *
1019 * *
1020 *****************************************************************************/
1021
1022 void
distances(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)1023 distances(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
1024 int *invar, int invararg, boolean digraph, int m, int n)
1025 {
1026 int i;
1027 set *gw;
1028 int wt;
1029 int d,dlim,cell1,cell2,iv,v,w;
1030 boolean success;
1031
1032 #if !MAXN
1033 DYNALLOC1(set,workset,workset_sz,m,"distances");
1034 DYNALLOC1(int,workshort,workshort_sz,n+2,"distances");
1035 DYNALLOC1(set,ws1,ws1_sz,m,"distances");
1036 DYNALLOC1(set,ws2,ws2_sz,m,"distances");
1037 #endif
1038
1039 for (i = n; --i >= 0;) invar[i] = 0;
1040
1041 wt = 1;
1042 for (i = 0; i < n; ++i)
1043 {
1044 workshort[lab[i]] = FUZZ1(wt);
1045 if (ptn[i] <= level) ++wt;
1046 }
1047
1048 if (invararg > n || invararg == 0) dlim = n;
1049 else dlim = invararg+1;
1050
1051 success = FALSE;
1052 for (cell1 = 0; cell1 < n; cell1 = cell2 + 1)
1053 {
1054 for (cell2 = cell1; ptn[cell2] > level; ++cell2) {}
1055 if (cell2 == cell1) continue;
1056
1057 for (iv = cell1; iv <= cell2; ++iv)
1058 {
1059 v = lab[iv];
1060 EMPTYSET(ws1,m);
1061 ADDELEMENT(ws1,v);
1062 EMPTYSET(ws2,m);
1063 ADDELEMENT(ws2,v);
1064 for (d = 1; d < dlim; ++d)
1065 {
1066 EMPTYSET(workset,m);
1067 wt = 0;
1068 w = -1;
1069 while ((w = nextelement(ws2,M,w)) >= 0)
1070 {
1071 gw = GRAPHROW(g,w,m);
1072 ACCUM(wt,workshort[w]);
1073 for (i = M; --i >= 0;) workset[i] |= gw[i];
1074 }
1075 if (wt == 0) break;
1076 ACCUM(wt,d);
1077 wt = FUZZ2(wt);
1078 ACCUM(invar[v],wt);
1079 for (i = M; --i >= 0;)
1080 {
1081 ws2[i] = workset[i] & ~ws1[i];
1082 ws1[i] |= ws2[i];
1083 }
1084 }
1085 if (invar[v] != invar[lab[cell1]]) success = TRUE;
1086 }
1087 if (success) break;
1088 }
1089 }
1090
1091 /*****************************************************************************
1092 * *
1093 * indsets() assigns to each vertex v a value depending on which cells the *
1094 * vertices which join v in an independent set lie in. The size of the *
1095 * independent sets which are used is the smallest of invararg and MAXCLIQUE.*
1096 * *
1097 *****************************************************************************/
1098
1099 void
indsets(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)1100 indsets(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
1101 int *invar, int invararg, boolean digraph, int m, int n)
1102 {
1103 int i;
1104 int wt;
1105 set *gv;
1106 int ss,setsize;
1107 int v[MAXCLIQUE];
1108 long wv[MAXCLIQUE];
1109 set *s0,*s1;
1110
1111 #if !MAXN
1112 DYNALLOC1(int,workshort,workshort_sz,n+2,"indsets");
1113 DYNALLOC2(set,wss,wss_sz,m,MAXCLIQUE-1,"indsets");
1114 #endif
1115
1116 for (i = n; --i >= 0;) invar[i] = 0;
1117
1118 if (invararg <= 1 || digraph) return;
1119
1120 if (invararg > MAXCLIQUE) setsize = MAXCLIQUE;
1121 else setsize = invararg;
1122
1123 wt = 1;
1124 for (i = 0; i < n; ++i)
1125 {
1126 workshort[lab[i]] = FUZZ2(wt);
1127 if (ptn[i] <= level) ++wt;
1128 }
1129
1130 for (v[0] = 0; v[0] < n; ++v[0])
1131 {
1132 wv[0] = workshort[v[0]];
1133 s0 = (set*)wss;
1134 EMPTYSET(s0,m);
1135 for (i = v[0]+1; i < n; ++i) ADDELEMENT(s0,i);
1136 gv = GRAPHROW(g,v[0],m);
1137 for (i = M; --i >= 0;) s0[i] &= ~gv[i];
1138 ss = 1;
1139 v[1] = v[0];
1140 while (ss > 0)
1141 {
1142 if (ss == setsize)
1143 {
1144 wt = FUZZ1(wv[ss-1]);
1145 for (i = ss; --i >= 0;) ACCUM(invar[v[i]],wt);
1146 --ss;
1147 }
1148 else if ((v[ss] = nextelement((set*)wss+M*(ss-1),M,v[ss])) < 0)
1149 --ss;
1150 else
1151 {
1152 wv[ss] = wv[ss-1] + workshort[v[ss]];
1153 ++ss;
1154 if (ss < setsize)
1155 {
1156 gv = GRAPHROW(g,v[ss-1],m);
1157 s1 = (set*)wss + M*(ss-2);
1158 for (i = M; --i >= 0;) s1[i+M] = s1[i] & ~gv[i];
1159 v[ss] = v[ss-1];
1160 }
1161 }
1162 }
1163 }
1164 }
1165
1166 /*****************************************************************************
1167 * *
1168 * cliques() assigns to each vertex v a value depending on which cells the *
1169 * vertices which join v in a clique lie in. The size of the cliques used *
1170 * is the smallest of invararg and MAXCLIQUE. *
1171 * *
1172 *****************************************************************************/
1173
1174 void
cliques(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)1175 cliques(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
1176 int *invar, int invararg, boolean digraph, int m, int n)
1177 {
1178 int i;
1179 int wt;
1180 set *gv;
1181 int ss,setsize;
1182 int v[MAXCLIQUE];
1183 long wv[MAXCLIQUE];
1184 set *ns;
1185
1186 #if !MAXN
1187 DYNALLOC1(int,workshort,workshort_sz,n+2,"cliques");
1188 DYNALLOC2(set,wss,wss_sz,m,MAXCLIQUE-1,"cliques");
1189 #else
1190 set wss[MAXCLIQUE-1][MAXM];
1191 #endif
1192
1193 for (i = n; --i >= 0;) invar[i] = 0;
1194
1195 if (invararg <= 1 || digraph) return;
1196
1197 if (invararg > MAXCLIQUE) setsize = MAXCLIQUE;
1198 else setsize = invararg;
1199
1200 wt = 1;
1201 for (i = 0; i < n; ++i)
1202 {
1203 workshort[lab[i]] = FUZZ2(wt);
1204 if (ptn[i] <= level) ++wt;
1205 }
1206
1207 for (v[0] = 0; v[0] < n; ++v[0])
1208 {
1209 wv[0] = workshort[v[0]];
1210 gv = GRAPHROW(g,v[0],m);
1211 ns = (set*)wss;
1212 for (i = M; --i >= 0;) ns[i] = gv[i];
1213 ss = 1;
1214 v[1] = v[0];
1215 while (ss > 0)
1216 {
1217 if (ss == setsize)
1218 {
1219 wt = FUZZ1(wv[ss-1]);
1220 for (i = ss; --i >= 0;) ACCUM(invar[v[i]],wt);
1221 --ss;
1222 }
1223 else if ((v[ss] = nextelement((set*)wss+M*(ss-1),M,v[ss])) < 0)
1224 --ss;
1225 else
1226 {
1227 wv[ss] = wv[ss-1] + workshort[v[ss]];
1228 ++ss;
1229 if (ss < setsize)
1230 {
1231 gv = GRAPHROW(g,v[ss-1],m);
1232 ns = (set*)wss + M*(ss-2);
1233 for (i = M; --i >= 0;) ns[i+M] = ns[i] & gv[i];
1234 v[ss] = v[ss-1];
1235 }
1236 }
1237 }
1238 }
1239 }
1240
1241 /*****************************************************************************
1242 * *
1243 * cellcliq() assigns to each vertex v a value depending on the number of *
1244 * cliques which v lies in and which lie in the same cell as v. The size *
1245 * of clique counted is the smallest of invararg and MAXCLIQUE. We try the *
1246 * cells in increasing order of size and stop as soon as any cell splits. *
1247 * *
1248 *****************************************************************************/
1249
1250 void
cellcliq(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)1251 cellcliq(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
1252 int *invar, int invararg, boolean digraph, int m, int n)
1253 {
1254 int i;
1255 int wt;
1256 set *gv;
1257 int ss,setsize;
1258 int v[MAXCLIQUE];
1259 set *ns;
1260 int *cellstart,*cellsize;
1261 int iv,icell,bigcells,cell1,cell2;
1262 int pc;
1263 setword sw;
1264
1265 #if !MAXN
1266 DYNALLOC1(set,workset,workset_sz,m,"cellcliq");
1267 DYNALLOC1(int,workshort,workshort_sz,n+2,"cellcliq");
1268 DYNALLOC2(set,wss,wss_sz,m,MAXCLIQUE-1,"cellcliq");
1269 #endif
1270
1271 for (i = n; --i >= 0;) invar[i] = 0;
1272
1273 if (invararg <= 1 || digraph) return;
1274
1275 if (invararg > MAXCLIQUE) setsize = MAXCLIQUE;
1276 else setsize = invararg;
1277
1278 cellstart = workshort;
1279 cellsize = workshort + (n/2);
1280 getbigcells(ptn,level,setsize > 6 ? setsize : 6,&bigcells,
1281 cellstart,cellsize,n);
1282
1283 for (icell = 0; icell < bigcells; ++icell)
1284 {
1285 cell1 = cellstart[icell];
1286 cell2 = cell1 + cellsize[icell] - 1;
1287
1288 EMPTYSET(workset,m);
1289 for (iv = cell1; iv <= cell2; ++iv) ADDELEMENT(workset,lab[iv]);
1290
1291 for (iv = cell1; iv <= cell2; ++iv)
1292 {
1293 v[0] = lab[iv];
1294 gv = GRAPHROW(g,v[0],m);
1295 ns = (set*)wss;
1296 pc = 0;
1297
1298 for (i = M; --i >= 0;)
1299 {
1300 ns[i] = gv[i] & workset[i];
1301 if ((sw = ns[i]) != 0) pc += POPCOUNT(sw);
1302 }
1303 if (pc <= 1 || pc >= cellsize[icell] - 2) continue;
1304
1305 ss = 1;
1306 v[1] = v[0];
1307 while (ss > 0)
1308 {
1309 if (ss == setsize)
1310 {
1311 for (i = ss; --i >= 0;) ++invar[v[i]];
1312 --ss;
1313 }
1314 else if ((v[ss]
1315 = nextelement((set*)wss+M*(ss-1),M,v[ss])) < 0)
1316 --ss;
1317 else
1318 {
1319 ++ss;
1320 if (ss < setsize)
1321 {
1322 gv = GRAPHROW(g,v[ss-1],m);
1323 ns = (set*)wss + M*(ss-2);
1324 for (i = M; --i >= 0;) ns[i+M] = ns[i] & gv[i];
1325 v[ss] = v[ss-1];
1326 }
1327 }
1328 }
1329 }
1330 wt = invar[lab[cell1]];
1331 for (iv = cell1 + 1; iv <= cell2; ++iv)
1332 if (invar[lab[iv]] != wt) return;
1333 }
1334 }
1335
1336 /*****************************************************************************
1337 * *
1338 * cellind() assigns to each vertex v a value depending on the number of *
1339 * independent sets which v lies in and which lie in the same cell as v. *
1340 * The size of clique counted is the smallest of invararg and MAXCLIQUE. *
1341 * We try the cells in increasing order of size and stop as soon as any *
1342 * cell splits. *
1343 * *
1344 *****************************************************************************/
1345
1346 void
cellind(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)1347 cellind(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
1348 int *invar, int invararg, boolean digraph, int m, int n)
1349 {
1350 int i;
1351 int wt;
1352 set *gv;
1353 int ss,setsize;
1354 int v[MAXCLIQUE];
1355 set *ns;
1356 int *cellstart,*cellsize;
1357 int iv,icell,bigcells,cell1,cell2;
1358 int pc;
1359 setword sw;
1360
1361 #if !MAXN
1362 DYNALLOC1(set,workset,workset_sz,m,"cellind");
1363 DYNALLOC1(int,workshort,workshort_sz,n+2,"cellind");
1364 DYNALLOC2(set,wss,wss_sz,m,MAXCLIQUE-1,"cellind");
1365 #endif
1366
1367 for (i = n; --i >= 0;) invar[i] = 0;
1368
1369 if (invararg <= 1 || digraph) return;
1370
1371 if (invararg > MAXCLIQUE) setsize = MAXCLIQUE;
1372 else setsize = invararg;
1373
1374 cellstart = workshort;
1375 cellsize = workshort + (n/2);
1376 getbigcells(ptn,level,setsize > 6 ? setsize : 6,&bigcells,
1377 cellstart,cellsize,n);
1378
1379 for (icell = 0; icell < bigcells; ++icell)
1380 {
1381 cell1 = cellstart[icell];
1382 cell2 = cell1 + cellsize[icell] - 1;
1383
1384 EMPTYSET(workset,m);
1385 for (iv = cell1; iv <= cell2; ++iv) ADDELEMENT(workset,lab[iv]);
1386
1387 for (iv = cell1; iv <= cell2; ++iv)
1388 {
1389 v[0] = lab[iv];
1390 gv = GRAPHROW(g,v[0],m);
1391 ns = (set*)wss;
1392 pc = 0;
1393
1394 for (i = M; --i >= 0;)
1395 {
1396 ns[i] = ~gv[i] & workset[i];
1397 if ((sw = ns[i]) != 0) pc += POPCOUNT(sw);
1398 }
1399 if (pc <= 1 || pc >= cellsize[icell] - 2) continue;
1400
1401 ss = 1;
1402 v[1] = v[0];
1403 while (ss > 0)
1404 {
1405 if (ss == setsize)
1406 {
1407 for (i = ss; --i >= 0;) ++invar[v[i]];
1408 --ss;
1409 }
1410 else if ((v[ss]
1411 = nextelement((set*)wss+M*(ss-1),M,v[ss])) < 0)
1412 --ss;
1413 else
1414 {
1415 ++ss;
1416 if (ss < setsize)
1417 {
1418 gv = GRAPHROW(g,v[ss-1],m);
1419 ns = (set*)wss + M*(ss-2);
1420 for (i = M; --i >= 0;) ns[i+M] = ns[i] & ~gv[i];
1421 v[ss] = v[ss-1];
1422 }
1423 }
1424 }
1425 }
1426 wt = invar[lab[cell1]];
1427 for (iv = cell1 + 1; iv <= cell2; ++iv)
1428 if (invar[lab[iv]] != wt) return;
1429 }
1430 }
1431
1432 /*****************************************************************************
1433 * *
1434 * adjacencies() assigns to each vertex v a code depending on which cells *
1435 * it is joined to and from, and how many times. It is intended to provide *
1436 * better partitioning that the normal refinement routine for digraphs. *
1437 * It will not help with undirected graphs in nauty at all. *
1438 * *
1439 *****************************************************************************/
1440
1441 void
adjacencies(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)1442 adjacencies(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
1443 int *invar, int invararg, boolean digraph, int m, int n)
1444 {
1445 int i,v,w;
1446 int vwt,wwt;
1447 set *gv;
1448
1449 #if !MAXN
1450 DYNALLOC1(int,workshort,workshort_sz,n+2,"adjacencies");
1451 #endif
1452
1453 vwt = 1;
1454 for (i = 0; i < n; ++i)
1455 {
1456 workshort[lab[i]] = vwt;
1457 if (ptn[i] <= level) ++vwt;
1458 invar[i] = 0;
1459 }
1460
1461 for (v = 0, gv = (set*)g; v < n; ++v, gv += M)
1462 {
1463 vwt = FUZZ1(workshort[v]);
1464 wwt = 0;
1465 w = -1;
1466 while ((w = nextelement(gv,M,w)) >= 0)
1467 {
1468 ACCUM(wwt,FUZZ2(workshort[w]));
1469 ACCUM(invar[w],vwt);
1470 }
1471 ACCUM(invar[v],wwt);
1472 }
1473 }
1474
1475 /*****************************************************************************
1476 * *
1477 * nautinv_check() checks that this file is compiled compatibly with the *
1478 * given parameters. If not, call exit(1). *
1479 * *
1480 *****************************************************************************/
1481
1482 void
nautinv_check(int wordsize,int m,int n,int version)1483 nautinv_check(int wordsize, int m, int n, int version)
1484 {
1485 if (wordsize != WORDSIZE)
1486 {
1487 fprintf(ERRFILE,"Error: WORDSIZE mismatch in nautinv.c\n");
1488 exit(1);
1489 }
1490
1491 #if MAXN
1492 if (m > MAXM)
1493 {
1494 fprintf(ERRFILE,"Error: MAXM inadequate in nautinv.c\n");
1495 exit(1);
1496 }
1497
1498 if (n > MAXN)
1499 {
1500 fprintf(ERRFILE,"Error: MAXN inadequate in nautinv.c\n");
1501 exit(1);
1502 }
1503 #endif
1504
1505 if (version < NAUTYREQUIRED)
1506 {
1507 fprintf(ERRFILE,"Error: nautinv.c version mismatch\n");
1508 exit(1);
1509 }
1510 }
1511
1512 /*****************************************************************************
1513 * *
1514 * nautinv_freedyn() - free the dynamic memory in this module *
1515 * *
1516 *****************************************************************************/
1517
1518 void
nautinv_freedyn(void)1519 nautinv_freedyn(void)
1520 {
1521 #if !MAXN
1522 DYNFREE(workset,workset_sz);
1523 DYNFREE(workshort,workshort_sz);
1524 DYNFREE(ws1,ws1_sz);
1525 DYNFREE(ws2,ws2_sz);
1526 DYNFREE(vv,vv_sz);
1527 DYNFREE(ww,ww_sz);
1528 DYNFREE(w01,w01_sz);
1529 DYNFREE(w02,w02_sz);
1530 DYNFREE(w03,w03_sz);
1531 DYNFREE(w12,w12_sz);
1532 DYNFREE(w13,w13_sz);
1533 DYNFREE(w23,w23_sz);
1534 DYNFREE(pt0,pt0_sz);
1535 DYNFREE(pt1,pt1_sz);
1536 DYNFREE(pt2,pt2_sz);
1537 DYNFREE(wss,wss_sz);
1538 #endif
1539 }
1540
1541 /*===================================================================*/
1542
1543
1544 /*****************************************************************************
1545 * *
1546 * semirefine(g,lab,ptn,level,numcells,strength,active,m,n) performs a *
1547 * refinement operation on the partition at the specified level of the *
1548 * partition nest (lab,ptn). *numcells is assumed to contain the number of *
1549 * cells on input, and is updated. The initial set of active cells (alpha *
1550 * in the paper) is specified in the set active. Precisely, x is in active *
1551 * iff the cell starting at index x in lab is active. *
1552 * *code is set to a value which depends on the fine detail of the *
1553 * algorithm, but which is independent of the labelling of the graph. *
1554 * *
1555 *****************************************************************************/
1556
1557 static int
semirefine(graph * g,int * lab,int * ptn,int level,int * numcells,int strength,set * active,int m,int n)1558 semirefine(graph *g, int *lab, int *ptn, int level, int *numcells,
1559 int strength, set *active, int m, int n)
1560 {
1561 int i,c1,c2,labc1;
1562 setword x;
1563 set *set1,*set2;
1564 int split1,split2,cell1,cell2;
1565 int cnt,bmin,bmax;
1566 long longcode;
1567 set *gptr;
1568 int maxcell,maxpos,hint;
1569
1570 #if !MAXN
1571 DYNALLOC1(int,workperm,workperm_sz,n,"refine");
1572 DYNALLOC1(set,workset,workset_sz,m,"refine");
1573 DYNALLOC1(int,bucket,bucket_sz,n+2,"refine");
1574 DYNALLOC1(int,count,count_sz,n,"refine");
1575 #endif
1576
1577 longcode = *numcells;
1578 split1 = -1;
1579 hint = 0;
1580 while (*numcells < n && ((split1 = hint, ISELEMENT(active,split1))
1581 || (split1 = nextelement(active,M,split1)) >= 0
1582 || (split1 = nextelement(active,M,-1)) >= 0))
1583 {
1584 DELELEMENT(active,split1);
1585 for (split2 = split1; ptn[split2] > level; ++split2) {}
1586 longcode = MASH(longcode,split1+split2);
1587 if (split1 == split2) /* trivial splitting cell */
1588 {
1589 gptr = GRAPHROW(g,lab[split1],M);
1590 for (cell1 = 0; cell1 < n; cell1 = cell2 + 1)
1591 {
1592 for (cell2 = cell1; ptn[cell2] > level; ++cell2) {}
1593 if (cell1 == cell2) continue;
1594 c1 = cell1;
1595 c2 = cell2;
1596 while (c1 <= c2)
1597 {
1598 labc1 = lab[c1];
1599 if (ISELEMENT(gptr,labc1))
1600 ++c1;
1601 else
1602 {
1603 lab[c1] = lab[c2];
1604 lab[c2] = labc1;
1605 --c2;
1606 }
1607 }
1608 if (c2 >= cell1 && c1 <= cell2)
1609 {
1610 ptn[c2] = level;
1611 longcode = MASH(longcode,FUZZ1(c2));
1612 ++*numcells;
1613 if (ISELEMENT(active,cell1) || c2-cell1 >= cell2-c1)
1614 {
1615 ADDELEMENT(active,c1);
1616 if (c1 == cell2) hint = c1;
1617 }
1618 else
1619 {
1620 ADDELEMENT(active,cell1);
1621 if (c2 == cell1) hint = cell1;
1622 }
1623 }
1624 }
1625 }
1626 else /* nontrivial splitting cell */
1627 {
1628 EMPTYSET(workset,m);
1629 for (i = split1; i <= split2; ++i)
1630 ADDELEMENT(workset,lab[i]);
1631 longcode = MASH(longcode,FUZZ2(split2-split1+1));
1632
1633 for (cell1 = 0; cell1 < n; cell1 = cell2 + 1)
1634 {
1635 for (cell2 = cell1; ptn[cell2] > level; ++cell2) {}
1636 if (cell1 == cell2) continue;
1637 i = cell1;
1638 set1 = workset;
1639 set2 = GRAPHROW(g,lab[i],m);
1640 cnt = 0;
1641 for (c1 = m; --c1 >= 0;)
1642 if ((x = ((*set1++) & (*set2++))) != 0)
1643 cnt += POPCOUNT(x);
1644
1645 count[i] = bmin = bmax = cnt;
1646 bucket[cnt] = 1;
1647 while (++i <= cell2)
1648 {
1649 set1 = workset;
1650 set2 = GRAPHROW(g,lab[i],m);
1651 cnt = 0;
1652 for (c1 = m; --c1 >= 0;)
1653 if ((x = ((*set1++) & (*set2++))) != 0)
1654 cnt += POPCOUNT(x);
1655
1656 while (bmin > cnt) bucket[--bmin] = 0;
1657 while (bmax < cnt) bucket[++bmax] = 0;
1658 ++bucket[cnt];
1659 count[i] = cnt;
1660 }
1661 if (bmin == bmax)
1662 {
1663 longcode = MASH(longcode,FUZZ1(bmin+cell1));
1664 continue;
1665 }
1666 c1 = cell1;
1667 maxcell = -1;
1668 for (i = bmin; i <= bmax; ++i)
1669 if (bucket[i])
1670 {
1671 c2 = c1 + bucket[i];
1672 bucket[i] = c1;
1673 longcode = MASH(longcode,i+c1);
1674 if (c2-c1 > maxcell)
1675 {
1676 maxcell = c2-c1;
1677 maxpos = c1;
1678 }
1679 if (c1 != cell1)
1680 {
1681 ADDELEMENT(active,c1);
1682 if (c2-c1 == 1) hint = c1;
1683 ++*numcells;
1684 }
1685 if (c2 <= cell2) ptn[c2-1] = level;
1686 c1 = c2;
1687 }
1688 for (i = cell1; i <= cell2; ++i)
1689 workperm[bucket[count[i]]++] = lab[i];
1690 for (i = cell1; i <= cell2; ++i) lab[i] = workperm[i];
1691 if (!ISELEMENT(active,cell1))
1692 {
1693 ADDELEMENT(active,cell1);
1694 DELELEMENT(active,maxpos); /* check maxpos is alwas defined */
1695 }
1696 }
1697 }
1698 if (--strength == 0) break; /* negative is fine! */
1699 }
1700
1701 longcode = MASH(longcode,FUZZ2(*numcells));
1702 return CLEANUP(longcode);
1703 }
1704
1705 void
refinvar(graph * g,int * lab,int * ptn,int level,int numcells,int tvpos,int * invar,int invararg,boolean digraph,int m,int n)1706 refinvar(graph *g, int *lab, int *ptn, int level, int numcells, int tvpos,
1707 int *invar, int invararg, boolean digraph, int m, int n)
1708 {
1709 int i,j;
1710 int wt;
1711 int icell,bigcells,cell1,cell2;
1712 int *cellstart,*cellsize;
1713 int newnumcells;
1714
1715 #if !MAXN
1716 DYNALLOC1(int,workshort,workshort_sz,n+2,"refinvar");
1717 DYNALLOC1(int,vv,vv_sz,n,"refinvar");
1718 DYNALLOC1(int,ww,ww_sz,n,"refinvar");
1719 DYNALLOC1(set,ws1,ws1_sz,n,"refinvar");
1720 #endif
1721
1722 for (i = n; --i >= 0;) invar[i] = 0;
1723
1724 cellstart = workshort;
1725 cellsize = workshort + (n/2);
1726 getbigcells(ptn,level,2,&bigcells,cellstart,cellsize,n);
1727
1728 for (icell = 0; icell < bigcells; ++icell)
1729 {
1730 cell1 = cellstart[icell];
1731 cell2 = cell1 + cellsize[icell] - 1;
1732 for (i = cell1; i <= cell2; ++i)
1733 {
1734 for (j = 0; j < n; ++j)
1735 {
1736 vv[j] = lab[j];
1737 ww[j] = ptn[j];
1738 }
1739 newnumcells = numcells + 1;
1740 ww[cell1] = level;
1741 EMPTYSET(ws1,m);
1742 ADDELEMENT(ws1,cell1);
1743 vv[i] = lab[cell1];
1744 vv[cell1] = lab[i];
1745 invar[lab[i]] = semirefine(g,vv,ww,level,&newnumcells,
1746 invararg,ws1,m,n);
1747 }
1748 wt = invar[lab[cell1]];
1749 for (i = cell1 + 1; i <= cell2; ++i)
1750 if (invar[lab[i]] != wt) return;
1751 }
1752 }
1753