latest version

1 files changed, 717 insertions, 0 deletions

diff --git a/graph-checker/nauty/naugraph.c b/graph-checker/nauty/naugraph.c
new file mode 100644
index 0000000..a880dae
--- /dev/null
+++ b/graph-checker/nauty/naugraph.c
@@ -0,0 +1,717 @@
+/*****************************************************************************
+*                                                                            *
+*  Graph-specific auxiliary source file for version 2.8 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()       *
+*        6-Apr-21 : increase work space in 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<rw>,permset(),workperm<rw>                      *
+*                                                                            *
+*****************************************************************************/
+
+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<rw>                                  *
+*                                                                            *
+*****************************************************************************/
+
+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<w>,bit<r>,nextelement(),bucket<w>,workperm<w>  *
+*                                                                            *
+*****************************************************************************/
+
+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<r>,workperm<rw>,workset<rw>,bucket<rw>              *
+*                                                                            *
+*****************************************************************************/
+
+static int
+bestcell(const graph *g, const int *lab, const 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
+}