/***************************************************************************** * * * Graph-specific auxiliary source file for version 2.7 of nauty. * * * * Copyright (1984-2019) Brendan McKay. All rights reserved. * * Subject to waivers and disclaimers in nauty.h. * * * * CHANGE HISTORY * * 16-Nov-00 : initial creation out of nautil.c * * 22-Apr-01 : added aproto line for Magma * * EXTDEFS is no longer required * * removed dynamic allocation from refine1() * * 21-Nov-01 : use NAUTYREQUIRED in naugraph_check() * * 23-Nov-06 : add targetcell(); make bestcell() local * * 10-Dec-06 : remove BIGNAUTY * * 10-Nov-09 : remove shortish and permutation types * * 23-May-10 : add densenauty() * * 15-Jan-12 : add TLS_ATTR attributes * * 23-Jan-13 : add some parens to make icc happy * * 15-Oct-19 : fix default size of dnwork[] to match densenauty() * * * *****************************************************************************/ #define ONE_WORD_SETS #include "nauty.h" /* macros for hash-codes: */ #define MASH(l,i) ((((l) ^ 065435) + (i)) & 077777) /* : expression whose long value depends only on long l and int/long i. Anything goes, preferably non-commutative. */ #define CLEANUP(l) ((int)((l) % 077777)) /* : expression whose value depends on long l and is less than 077777 when converted to int then short. Anything goes. */ #if MAXM==1 #define M 1 #else #define M m #endif /* aproto: header new_nauty_protos.h */ dispatchvec dispatch_graph = {isautom,testcanlab,updatecan,refine,refine1,cheapautom,targetcell, naugraph_freedyn,naugraph_check,NULL,NULL}; #if !MAXN DYNALLSTAT(set,workset,workset_sz); DYNALLSTAT(int,workperm,workperm_sz); DYNALLSTAT(int,bucket,bucket_sz); DYNALLSTAT(set,dnwork,dnwork_sz); #else static TLS_ATTR set workset[MAXM]; /* used for scratch work */ static TLS_ATTR int workperm[MAXN]; static TLS_ATTR int bucket[MAXN+2]; static TLS_ATTR set dnwork[2*500*MAXM]; #endif /***************************************************************************** * * * isautom(g,perm,digraph,m,n) = TRUE iff perm is an automorphism of g * * (i.e., g^perm = g). Symmetry is assumed unless digraph = TRUE. * * * *****************************************************************************/ boolean isautom(graph *g, int *perm, boolean digraph, int m, int n) { set *pg; int pos; set *pgp; int posp,i; for (pg = g, i = 0; i < n; pg += M, ++i) { pgp = GRAPHROW(g,perm[i],M); pos = (digraph ? -1 : i); while ((pos = nextelement(pg,M,pos)) >= 0) { posp = perm[pos]; if (!ISELEMENT(pgp,posp)) return FALSE; } } return TRUE; } /***************************************************************************** * * * testcanlab(g,canong,lab,samerows,m,n) compares g^lab to canong, * * using an ordering which is immaterial since it's only used here. The * * value returned is -1,0,1 if g^lab <,=,> canong. *samerows is set to * * the number of rows (0..n) of canong which are the same as those of g^lab. * * * * GLOBALS ACCESSED: workset,permset(),workperm * * * *****************************************************************************/ int testcanlab(graph *g, graph *canong, int *lab, int *samerows, int m, int n) { int i,j; set *ph; #if !MAXN DYNALLOC1(int,workperm,workperm_sz,n,"testcanlab"); DYNALLOC1(set,workset,workset_sz,m,"testcanlab"); #endif for (i = 0; i < n; ++i) workperm[lab[i]] = i; for (i = 0, ph = canong; i < n; ++i, ph += M) { permset(GRAPHROW(g,lab[i],M),workset,M,workperm); for (j = 0; j < M; ++j) if (workset[j] < ph[j]) { *samerows = i; return -1; } else if (workset[j] > ph[j]) { *samerows = i; return 1; } } *samerows = n; return 0; } /***************************************************************************** * * * updatecan(g,canong,lab,samerows,m,n) sets canong = g^lab, assuming * * the first samerows of canong are ok already. * * * * GLOBALS ACCESSED: permset(),workperm * * * *****************************************************************************/ void updatecan(graph *g, graph *canong, int *lab, int samerows, int m, int n) { int i; set *ph; #if !MAXN DYNALLOC1(int,workperm,workperm_sz,n,"updatecan"); #endif for (i = 0; i < n; ++i) workperm[lab[i]] = i; for (i = samerows, ph = GRAPHROW(canong,samerows,M); i < n; ++i, ph += M) permset(GRAPHROW(g,lab[i],M),ph,M,workperm); } /***************************************************************************** * * * refine(g,lab,ptn,level,numcells,count,active,code,m,n) performs a * * refinement operation on the partition at the specified level of the * * partition nest (lab,ptn). *numcells is assumed to contain the number of * * cells on input, and is updated. The initial set of active cells (alpha * * in the paper) is specified in the set active. Precisely, x is in active * * iff the cell starting at index x in lab is active. * * The resulting partition is equitable if active is correct (see the paper * * and the Guide). * * *code is set to a value which depends on the fine detail of the * * algorithm, but which is independent of the labelling of the graph. * * count is used for work space. * * * * GLOBALS ACCESSED: workset,bit,nextelement(),bucket,workperm * * * *****************************************************************************/ void refine(graph *g, int *lab, int *ptn, int level, int *numcells, int *count, set *active, int *code, int m, int n) { #if MAXM==1 refine1(g,lab,ptn,level,numcells,count,active,code,m,n); } #else int i,c1,c2,labc1; setword x; set *set1,*set2; int split1,split2,cell1,cell2; int cnt,bmin,bmax; long longcode; set *gptr; int maxcell,maxpos,hint; #if !MAXN DYNALLOC1(int,workperm,workperm_sz,n,"refine"); DYNALLOC1(set,workset,workset_sz,m,"refine"); DYNALLOC1(int,bucket,bucket_sz,n+2,"refine"); #endif longcode = *numcells; split1 = -1; hint = 0; while (*numcells < n && ((split1 = hint, ISELEMENT(active,split1)) || (split1 = nextelement(active,M,split1)) >= 0 || (split1 = nextelement(active,M,-1)) >= 0)) { DELELEMENT(active,split1); for (split2 = split1; ptn[split2] > level; ++split2) {} longcode = MASH(longcode,split1+split2); if (split1 == split2) /* trivial splitting cell */ { gptr = GRAPHROW(g,lab[split1],M); for (cell1 = 0; cell1 < n; cell1 = cell2 + 1) { for (cell2 = cell1; ptn[cell2] > level; ++cell2) {} if (cell1 == cell2) continue; c1 = cell1; c2 = cell2; while (c1 <= c2) { labc1 = lab[c1]; if (ISELEMENT(gptr,labc1)) ++c1; else { lab[c1] = lab[c2]; lab[c2] = labc1; --c2; } } if (c2 >= cell1 && c1 <= cell2) { ptn[c2] = level; longcode = MASH(longcode,c2); ++*numcells; if (ISELEMENT(active,cell1) || c2-cell1 >= cell2-c1) { ADDELEMENT(active,c1); if (c1 == cell2) hint = c1; } else { ADDELEMENT(active,cell1); if (c2 == cell1) hint = cell1; } } } } else /* nontrivial splitting cell */ { EMPTYSET(workset,m); for (i = split1; i <= split2; ++i) ADDELEMENT(workset,lab[i]); longcode = MASH(longcode,split2-split1+1); for (cell1 = 0; cell1 < n; cell1 = cell2 + 1) { for (cell2 = cell1; ptn[cell2] > level; ++cell2) {} if (cell1 == cell2) continue; i = cell1; set1 = workset; set2 = GRAPHROW(g,lab[i],m); cnt = 0; for (c1 = m; --c1 >= 0;) if ((x = ((*set1++) & (*set2++))) != 0) cnt += POPCOUNT(x); count[i] = bmin = bmax = cnt; bucket[cnt] = 1; while (++i <= cell2) { set1 = workset; set2 = GRAPHROW(g,lab[i],m); cnt = 0; for (c1 = m; --c1 >= 0;) if ((x = ((*set1++) & (*set2++))) != 0) cnt += POPCOUNT(x); while (bmin > cnt) bucket[--bmin] = 0; while (bmax < cnt) bucket[++bmax] = 0; ++bucket[cnt]; count[i] = cnt; } if (bmin == bmax) { longcode = MASH(longcode,bmin+cell1); continue; } c1 = cell1; maxcell = -1; for (i = bmin; i <= bmax; ++i) if (bucket[i]) { c2 = c1 + bucket[i]; bucket[i] = c1; longcode = MASH(longcode,i+c1); if (c2-c1 > maxcell) { maxcell = c2-c1; maxpos = c1; } if (c1 != cell1) { ADDELEMENT(active,c1); if (c2-c1 == 1) hint = c1; ++*numcells; } if (c2 <= cell2) ptn[c2-1] = level; c1 = c2; } for (i = cell1; i <= cell2; ++i) workperm[bucket[count[i]]++] = lab[i]; for (i = cell1; i <= cell2; ++i) lab[i] = workperm[i]; if (!ISELEMENT(active,cell1)) { ADDELEMENT(active,cell1); DELELEMENT(active,maxpos); } } } } longcode = MASH(longcode,*numcells); *code = CLEANUP(longcode); } #endif /* else case of MAXM==1 */ /***************************************************************************** * * * refine1(g,lab,ptn,level,numcells,count,active,code,m,n) is the same as * * refine(g,lab,ptn,level,numcells,count,active,code,m,n), except that * * m==1 is assumed for greater efficiency. The results are identical in all * * respects. See refine (above) for the specs. * * * *****************************************************************************/ void refine1(graph *g, int *lab, int *ptn, int level, int *numcells, int *count, set *active, int *code, int m, int n) { int i,c1,c2,labc1; setword x; int split1,split2,cell1,cell2; int cnt,bmin,bmax; long longcode; set *gptr,workset0; int maxcell,maxpos,hint; #if !MAXN DYNALLOC1(int,workperm,workperm_sz,n,"refine1"); DYNALLOC1(int,bucket,bucket_sz,n+2,"refine1"); #endif longcode = *numcells; split1 = -1; hint = 0; while (*numcells < n && ((split1 = hint, ISELEMENT1(active,split1)) || (split1 = nextelement(active,1,split1)) >= 0 || (split1 = nextelement(active,1,-1)) >= 0)) { DELELEMENT1(active,split1); for (split2 = split1; ptn[split2] > level; ++split2) {} longcode = MASH(longcode,split1+split2); if (split1 == split2) /* trivial splitting cell */ { gptr = GRAPHROW(g,lab[split1],1); for (cell1 = 0; cell1 < n; cell1 = cell2 + 1) { for (cell2 = cell1; ptn[cell2] > level; ++cell2) {} if (cell1 == cell2) continue; c1 = cell1; c2 = cell2; while (c1 <= c2) { labc1 = lab[c1]; if (ISELEMENT1(gptr,labc1)) ++c1; else { lab[c1] = lab[c2]; lab[c2] = labc1; --c2; } } if (c2 >= cell1 && c1 <= cell2) { ptn[c2] = level; longcode = MASH(longcode,c2); ++*numcells; if (ISELEMENT1(active,cell1) || c2-cell1 >= cell2-c1) { ADDELEMENT1(active,c1); if (c1 == cell2) hint = c1; } else { ADDELEMENT1(active,cell1); if (c2 == cell1) hint = cell1; } } } } else /* nontrivial splitting cell */ { workset0 = 0; for (i = split1; i <= split2; ++i) ADDELEMENT1(&workset0,lab[i]); longcode = MASH(longcode,split2-split1+1); for (cell1 = 0; cell1 < n; cell1 = cell2 + 1) { for (cell2 = cell1; ptn[cell2] > level; ++cell2) {} if (cell1 == cell2) continue; i = cell1; if ((x = workset0 & g[lab[i]]) != 0) cnt = POPCOUNT(x); else cnt = 0; count[i] = bmin = bmax = cnt; bucket[cnt] = 1; while (++i <= cell2) { if ((x = workset0 & g[lab[i]]) != 0) cnt = POPCOUNT(x); else cnt = 0; while (bmin > cnt) bucket[--bmin] = 0; while (bmax < cnt) bucket[++bmax] = 0; ++bucket[cnt]; count[i] = cnt; } if (bmin == bmax) { longcode = MASH(longcode,bmin+cell1); continue; } c1 = cell1; maxcell = -1; for (i = bmin; i <= bmax; ++i) if (bucket[i]) { c2 = c1 + bucket[i]; bucket[i] = c1; longcode = MASH(longcode,i+c1); if (c2-c1 > maxcell) { maxcell = c2-c1; maxpos = c1; } if (c1 != cell1) { ADDELEMENT1(active,c1); if (c2-c1 == 1) hint = c1; ++*numcells; } if (c2 <= cell2) ptn[c2-1] = level; c1 = c2; } for (i = cell1; i <= cell2; ++i) workperm[bucket[count[i]]++] = lab[i]; for (i = cell1; i <= cell2; ++i) lab[i] = workperm[i]; if (!ISELEMENT1(active,cell1)) { ADDELEMENT1(active,cell1); DELELEMENT1(active,maxpos); } } } } longcode = MASH(longcode,*numcells); *code = CLEANUP(longcode); } /***************************************************************************** * * * cheapautom(ptn,level,digraph,n) returns TRUE if the partition at the * * specified level in the partition nest (lab,ptn) {lab is not needed here} * * satisfies a simple sufficient condition for its cells to be the orbits of * * some subgroup of the automorphism group. Otherwise it returns FALSE. * * It always returns FALSE if digraph!=FALSE. * * * * nauty assumes that this function will always return TRUE for any * * partition finer than one for which it returns TRUE. * * * *****************************************************************************/ boolean cheapautom(int *ptn, int level, boolean digraph, int n) { int i,k,nnt; if (digraph) return FALSE; k = n; nnt = 0; for (i = 0; i < n; ++i) { --k; if (ptn[i] > level) { ++nnt; while (ptn[++i] > level) {} } } return (k <= nnt + 1 || k <= 4); } /***************************************************************************** * * * bestcell(g,lab,ptn,level,tc_level,m,n) returns the index in lab of the * * start of the "best non-singleton cell" for fixing. If there is no * * non-singleton cell it returns n. * * This implementation finds the first cell which is non-trivially joined * * to the greatest number of other cells. * * * * GLOBALS ACCESSED: bit,workperm,workset,bucket * * * *****************************************************************************/ static int bestcell(graph *g, int *lab, int *ptn, int level, int tc_level, int m, int n) { int i; set *gp; setword setword1,setword2; int v1,v2,nnt; #if !MAXN DYNALLOC1(int,workperm,workperm_sz,n,"bestcell"); DYNALLOC1(set,workset,workset_sz,m,"bestcell"); DYNALLOC1(int,bucket,bucket_sz,n+2,"bestcell"); #endif /* find non-singleton cells: put starts in workperm[0..nnt-1] */ i = nnt = 0; while (i < n) { if (ptn[i] > level) { workperm[nnt++] = i; while (ptn[i] > level) ++i; } ++i; } if (nnt == 0) return n; /* set bucket[i] to # non-trivial neighbours of n.s. cell i */ for (i = nnt; --i >= 0;) bucket[i] = 0; for (v2 = 1; v2 < nnt; ++v2) { EMPTYSET(workset,m); i = workperm[v2] - 1; do { ++i; ADDELEMENT(workset,lab[i]); } while (ptn[i] > level); for (v1 = 0; v1 < v2; ++v1) { gp = GRAPHROW(g,lab[workperm[v1]],m); #if MAXM==1 setword1 = *workset & *gp; setword2 = *workset & ~*gp; #else setword1 = setword2 = 0; for (i = m; --i >= 0;) { setword1 |= workset[i] & gp[i]; setword2 |= workset[i] & ~gp[i]; } #endif if (setword1 != 0 && setword2 != 0) { ++bucket[v1]; ++bucket[v2]; } } } /* find first greatest bucket value */ v1 = 0; v2 = bucket[0]; for (i = 1; i < nnt; ++i) if (bucket[i] > v2) { v1 = i; v2 = bucket[i]; } return (int)workperm[v1]; } /***************************************************************************** * * * targetcell(g,lab,ptn,level,tc_level,digraph,hint,m,n) returns the index * * in lab of the next cell to split. * * hint is a suggestion for the answer, which is obeyed if it is valid. * * Otherwise we use bestcell() up to tc_level and the first non-trivial * * cell after that. * * * *****************************************************************************/ int targetcell(graph *g, int *lab, int *ptn, int level, int tc_level, boolean digraph, int hint, int m, int n) { int i; if (hint >= 0 && ptn[hint] > level && (hint == 0 || ptn[hint-1] <= level)) return hint; else if (level <= tc_level) return bestcell(g,lab,ptn,level,tc_level,m,n); else { for (i = 0; i < n && ptn[i] <= level; ++i) {} return (i == n ? 0 : i); } } /***************************************************************************** * * * densenauty(g,lab,ptn,orbits,&options,&stats,m,n,h) * * is a slightly simplified interface to nauty(). It allocates enough * * workspace for 500 automorphisms and checks that the densegraph dispatch * * vector is in use. * * * *****************************************************************************/ void densenauty(graph *g, int *lab, int *ptn, int *orbits, optionblk *options, statsblk *stats, int m, int n, graph *h) { if (options->dispatch != &dispatch_graph) { fprintf(ERRFILE,"Error: densenauty() needs standard options block\n"); exit(1); } #if !MAXN /* Don't increase 2*500*m in the next line unless you also increase the default declaration of dnwork[] earlier. */ DYNALLOC1(set,dnwork,dnwork_sz,2*500*m,"densenauty malloc"); #endif nauty(g,lab,ptn,NULL,orbits,options,stats,dnwork,2*500*m,m,n,h); } /***************************************************************************** * * * naugraph_check() checks that this file is compiled compatibly with the * * given parameters. If not, call exit(1). * * * *****************************************************************************/ void naugraph_check(int wordsize, int m, int n, int version) { if (wordsize != WORDSIZE) { fprintf(ERRFILE,"Error: WORDSIZE mismatch in naugraph.c\n"); exit(1); } #if MAXN if (m > MAXM) { fprintf(ERRFILE,"Error: MAXM inadequate in naugraph.c\n"); exit(1); } if (n > MAXN) { fprintf(ERRFILE,"Error: MAXN inadequate in naugraph.c\n"); exit(1); } #endif if (version < NAUTYREQUIRED) { fprintf(ERRFILE,"Error: naugraph.c version mismatch\n"); exit(1); } } /***************************************************************************** * * * naugraph_freedyn() - free the dynamic memory in this module * * * *****************************************************************************/ void naugraph_freedyn(void) { #if !MAXN DYNFREE(workset,workset_sz); DYNFREE(workperm,workperm_sz); DYNFREE(bucket,bucket_sz); DYNFREE(dnwork,dnwork_sz); #endif }