summaryrefslogtreecommitdiff
path: root/nauty/gutil2.c
diff options
context:
space:
mode:
Diffstat (limited to 'nauty/gutil2.c')
-rw-r--r--nauty/gutil2.c762
1 files changed, 762 insertions, 0 deletions
diff --git a/nauty/gutil2.c b/nauty/gutil2.c
new file mode 100644
index 0000000..87af574
--- /dev/null
+++ b/nauty/gutil2.c
@@ -0,0 +1,762 @@
+/* gutil2.c: Some more graph utilities. */
+
+#include "gtools.h"
+#include "gutils.h"
+
+/**************************************************************************/
+
+long
+digoncount(graph *g, int m, int n)
+/* Number of digons (cycles of length 2). Useful for digraphs. */
+{
+ int i,j;
+ set *gi;
+ setword w;
+ long ans;
+
+ ans = 0;
+
+ if (m == 1)
+ {
+ for (i = 0; i < n; ++i)
+ {
+ w = g[i] & BITMASK(i);
+ while (w)
+ {
+ TAKEBIT(j,w);
+ if ((g[j] & bit[i])) ++ans;
+ }
+ }
+ }
+ else
+ {
+ for (i = 0, gi = g; i < n; ++i, gi += m)
+ {
+ for (j = i; (j = nextelement(gi,m,j)) > 0; )
+ if (ISELEMENT(g+j*m,i)) ++ans;
+ }
+ }
+
+ return ans;
+}
+
+/**************************************************************************/
+
+int
+loopcount(graph *g, int m, int n)
+/* Number of loops */
+{
+ set *gi;
+ int i,nl;
+
+ nl = 0;
+ for (i = 0, gi = g; i < n; ++i, gi += m)
+ if (ISELEMENT(gi,i)) ++nl;
+
+ return nl;
+}
+
+/**************************************************************************/
+
+long
+pathcount1(graph *g, int start, setword body, setword last)
+/* Number of paths in g starting at start, lying within body and
+ ending in last. {start} and last should be disjoint subsets of body. */
+{
+ long count;
+ setword gs,w;
+ int i;
+
+ gs = g[start];
+ w = gs & last;
+ count = POPCOUNT(w);
+
+ body &= ~bit[start];
+ w = gs & body;
+ while (w)
+ {
+ TAKEBIT(i,w);
+ count += pathcount1(g,i,body,last&~bit[i]);
+ }
+
+ return count;
+}
+
+/**************************************************************************/
+
+long
+cyclecount1(graph *g, int n)
+/* The total number of cycles in g (assumed no loops), m=1 only */
+{
+ setword body,nbhd;
+ long total;
+ int i,j;
+
+ body = ALLMASK(n);
+ total = 0;
+
+ for (i = 0; i < n-2; ++i)
+ {
+ body ^= bit[i];
+ nbhd = g[i] & body;
+ while (nbhd)
+ {
+ TAKEBIT(j,nbhd);
+ total += pathcount1(g,j,body,nbhd);
+ }
+ }
+
+ return total;
+}
+
+/**************************************************************************/
+
+long
+cyclecount(graph *g, int m, int n)
+/* The total number of cycles in g (assumed no loops) */
+{
+ if (n == 0) return 0;
+ if (m == 1) return cyclecount1(g,n);
+
+ gt_abort(">E cycle counting is only implemented for n <= WORDSIZE\n");
+ return 0;
+}
+
+/**************************************************************************/
+
+long
+indpathcount1(graph *g, int start, setword body, setword last)
+/* Number of induced paths in g starting at start, extravertices within
+ * body and ending in last.
+ * {start}, body and last should be disjoint. */
+{
+ long count;
+ setword gs,w;
+ int i;
+
+ gs = g[start];
+ w = gs & last;
+ count = POPCOUNT(w);
+
+ w = gs & body;
+ while (w)
+ {
+ TAKEBIT(i,w);
+ count += indpathcount1(g,i,body&~gs,last&~bit[i]&~gs);
+ }
+
+ return count;
+}
+
+/**************************************************************************/
+
+long
+indcyclecount1(graph *g, int n)
+/* The total number of induced cycles in g (assumed no loops), m=1 only */
+{
+ setword body,last,cni;
+ long total;
+ int i,j;
+
+ body = ALLMASK(n);
+ total = 0;
+
+ for (i = 0; i < n-2; ++i)
+ {
+ body ^= bit[i];
+ last = g[i] & body;
+ cni = g[i] | bit[i];
+ while (last)
+ {
+ TAKEBIT(j,last);
+ total += indpathcount1(g,j,body&~cni,last);
+ }
+ }
+
+ return total;
+}
+
+/**************************************************************************/
+
+long
+indcyclecount(graph *g, int m, int n)
+/* The total number of induced cycles in g (assumed no loops) */
+{
+ if (n == 0) return 0;
+ if (m == 1) return indcyclecount1(g,n);
+
+ gt_abort(
+ ">E induced cycle counting is only implemented for n <= WORDSIZE\n");
+ return 0;
+}
+
+/**************************************************************************/
+
+long
+numind3sets1(graph *g, int n)
+/* The number of triangles in g; undirected only */
+{
+ int i,j;
+ setword gi,w;
+ long total;
+
+ total = 0;
+ for (i = 2; i < n; ++i)
+ {
+ gi = ALLMASK(i) & ~g[i];
+ while (gi)
+ {
+ TAKEBIT(j,gi);
+ w = gi & ~g[j];
+ total += POPCOUNT(w);
+ }
+ }
+
+ return total;
+}
+
+/**************************************************************************/
+
+long
+numind3sets(graph *g, int m, int n)
+/* The number of independent 3-sets in g; undirected only */
+{
+ if (m == 1) return numind3sets1(g,n);
+
+ gt_abort(">E numind3sets is only implemented for n <= WORDSIZE\n");
+ return 0;
+}
+
+/**************************************************************************/
+
+long
+numtriangles1(graph *g, int n)
+/* The number of triangles in g; undirected only */
+{
+ int i,j;
+ setword gi,w;
+ long total;
+
+ total = 0;
+ for (i = 0; i < n-2; ++i)
+ {
+ gi = g[i] & BITMASK(i);
+ while (gi)
+ {
+ TAKEBIT(j,gi);
+ w = g[j] & gi;
+ total += POPCOUNT(w);
+ }
+ }
+
+ return total;
+}
+
+/**************************************************************************/
+
+long
+numtriangles(graph *g, int m, int n)
+/* The number of triangles in g; undirected only */
+{
+ int i,j,k,kw;
+ setword *gi,*gj,w;
+ long total;
+
+ if (m == 1) return numtriangles1(g,n);
+
+ total = 0;
+ for (i = 0, gi = g; i < n-2; ++i, gi += m)
+ for (j = i; (j = nextelement(gi,m,j)) > 0; )
+ {
+ gj = GRAPHROW(g,j,m);
+ kw = SETWD(j);
+ w = gi[kw] & gj[kw] & BITMASK(SETBT(j));
+ if (w) total += POPCOUNT(w);
+ for (k = kw+1; k < m; ++k)
+ {
+ w = gi[k] & gj[k];
+ total += POPCOUNT(w);
+ }
+ }
+
+ return total;
+}
+
+/**************************************************************************/
+
+long
+numdirtriangles1(graph *g, int n)
+{
+ long total;
+ int i,j,k;
+ setword biti,bm,wi,wj;
+
+ total = 0;
+ for (i = 0; i < n; ++i)
+ {
+ biti = bit[i];
+ bm = BITMASK(i);
+ wi = g[i] & bm;
+ while (wi)
+ {
+ TAKEBIT(j,wi);
+ wj = g[j] & bm;
+ while (wj)
+ {
+ TAKEBIT(k,wj);
+ if ((g[k] & biti)) ++total;
+ }
+ }
+ }
+
+ return total;
+}
+
+/**************************************************************************/
+
+long
+numdirtriangles(graph *g, int m, int n)
+/* The number of directed triangles in g */
+{
+ long total;
+ int i,j,k;
+ set *gi,*gj;
+
+ if (m == 1) return numdirtriangles1(g,n);
+
+ total = 0;
+ for (i = 0, gi = g; i < n-2; ++i, gi += m)
+ for (j = i; (j = nextelement(gi,m,j)) >= 0;)
+ {
+ gj = GRAPHROW(g,j,m);
+ for (k = i; (k = nextelement(gj,m,k)) >= 0; )
+ if (k != j && ISELEMENT(GRAPHROW(g,k,m),i)) ++total;
+ }
+
+ return total;
+}
+
+/**************************************************************************/
+
+void
+commonnbrs(graph *g, int *minadj, int *maxadj, int *minnon, int *maxnon,
+ int m, int n)
+/* Count the common neighbours of pairs of vertices, and give the minimum
+ and maximum for adjacent and non-adjacent vertices. Undirected only.
+ Null minimums are n+1 and null maximums are -1.
+*/
+{
+ int j,k;
+ int mina,maxa,minn,maxn;
+ int cn;
+ set *gi,*gj;
+ setword w;
+
+ if (n == 0)
+ {
+ *minadj = *maxadj = *minnon = *maxnon = 0;
+ return;
+ }
+
+ mina = minn = n+1;
+ maxa = maxn = -1;
+
+ for (j = 0, gj = g; j < n; ++j, gj += m)
+ for (gi = g; gi != gj; gi += m)
+ {
+ cn = 0;
+ for (k = 0; k < m; ++k)
+ {
+ w = gi[k] & gj[k];
+ if (w) cn += POPCOUNT(w);
+ }
+
+ if (ISELEMENT(gi,j))
+ {
+ if (cn < mina) mina = cn;
+ if (cn > maxa) maxa = cn;
+ }
+ else
+ {
+ if (cn < minn) minn = cn;
+ if (cn > maxn) maxn = cn;
+ }
+ }
+
+ *minadj = mina;
+ *maxadj = maxa;
+ *minnon = minn;
+ *maxnon = maxn;
+}
+
+/**************************************************************************/
+
+void
+delete1(graph *g, graph *h, int v, int n)
+/* Delete vertex v from g, result in h */
+{
+ setword mask1,mask2,gi;
+ int i;
+
+ mask1 = ALLMASK(v);
+ mask2 = BITMASK(v);
+
+ for (i = 0; i < v; ++i)
+ {
+ gi = g[i];
+ h[i] = (gi & mask1) | ((gi & mask2) << 1);
+ }
+ for (i = v; i < n-1; ++i)
+ {
+ gi = g[i+1];
+ h[i] = (gi & mask1) | ((gi & mask2) << 1);
+ }
+}
+
+/**************************************************************************/
+
+void
+contract1(graph *g, graph *h, int v, int w, int n)
+/* Contract distinct vertices v and w (not necessarily adjacent)
+ with result in h. No loops are created. */
+{
+ int x,y;
+ setword bitx,bity,mask1,mask2;
+ int i;
+
+ if (w < v)
+ {
+ x = w;
+ y = v;
+ }
+ else
+ {
+ x = v;
+ y = w;
+ }
+
+ bitx = bit[x];
+ bity = bit[y];
+ mask1 = ALLMASK(y);
+ mask2 = BITMASK(y);
+
+ for (i = 0; i < n; ++i)
+ if (g[i] & bity)
+ h[i] = (g[i] & mask1) | bitx | ((g[i] & mask2) << 1);
+ else
+ h[i] = (g[i] & mask1) | ((g[i] & mask2) << 1);
+
+ h[x] |= h[y];
+ for (i = y+1; i < n; ++i) h[i-1] = h[i];
+ h[x] &= ~bitx;
+}
+
+/**************************************************************************/
+
+static TLS_ATTR int knm[18][16]; /* knm[n,m] = conncontent(K_n - m*K_2) */
+static TLS_ATTR boolean knm_computed = FALSE;
+
+int
+conncontent(graph *g, int m, int n)
+/* number of connected spanning subgraphs with an even number
+ of edges minus the number with an odd number of edges */
+{
+ graph h[WORDSIZE];
+ setword gj;
+ int i,j,v1,v2,x,y;
+ int minv,mindeg,deg,goodv;
+ long ne;
+
+ if (m > 1) ABORT("conncontent only implemented for m=1");
+
+ /* First handle tiny graphs */
+
+ if (n <= 3)
+ {
+ if (n == 1) return 1;
+ if (n == 2) return (g[0] ? -1 : 0);
+ if (!g[0] || !g[1] || !g[2]) return 0; /* disconnected */
+ if (g[0]^g[1]^g[2]) return 1; /* path */
+ return 2; /* triangle */
+ }
+
+ /* Now compute
+ ne = number of edges
+ mindeg = minimum degree
+ minv = a vertex of minimum degree
+ goodv = a vertex with a clique neighbourhood (-1 if none)
+ */
+
+ mindeg = n;
+ ne = 0;
+ goodv = -1;
+ for (j = 0; j < n; ++j)
+ {
+ gj = g[j];
+ deg = POPCOUNT(gj);
+ ne += deg;
+ if (deg < mindeg)
+ {
+ mindeg = deg;
+ minv = j;
+ if (deg == 1) goodv = j;
+ }
+ if (deg >= 3 && deg <= 4 && goodv < 0)
+ {
+ while (gj)
+ {
+ TAKEBIT(i,gj);
+ if (gj & ~g[i]) break;
+ }
+ if (!gj) goodv = j;
+ }
+ }
+ ne /= 2;
+
+/* Cases of isolated vertex or tree */
+
+ if (mindeg == 0) return 0;
+
+#if 0
+ if (mindeg == 1 && ne == n-1)
+ {
+ if (isconnected1(g,n)) return ((n&1) ? 1 : -1);
+ else return 0;
+ }
+#endif
+
+/* Cases of clique and near-clique */
+
+ if (mindeg == n-1)
+ {
+ j = -1;
+ for (i = 2; i < n; ++i) j *= -i;
+ return j;
+ }
+
+ if (mindeg == n-2 && n < 16)
+ {
+ if (!knm_computed)
+ {
+ knm_computed = TRUE;
+ knm[1][0] = 1;
+ for (i = 2; i < 16; ++i)
+ {
+ knm[i][0] = -knm[i-1][0] * (i-1);
+ for (j = 1; j+j <= i; ++j)
+ knm[i][j] = knm[i][j-1] + knm[i-1][j-1];
+ }
+ }
+ return knm[n][(n*n-n)/2-ne];
+ }
+
+/* Case of vertex with clique neighbourhood */
+
+ if (goodv >= 0)
+ {
+ delete1(g,h,goodv,n);
+ return -POPCOUNT(g[goodv]) * conncontent(h,m,n-1);
+ }
+
+/* Case of minimum degree 2 */
+
+ if (mindeg == 2)
+ {
+ x = FIRSTBITNZ(g[minv]);
+ y = FIRSTBITNZ(g[minv]^bit[x]);
+ if (x > minv) --x;
+ if (y > minv) --y;
+ delete1(g,h,minv,n);
+ v1 = conncontent(h,m,n-1);
+ if (h[x] & bit[y]) return -2*v1; /* adjacent neighbours */
+
+ h[x] |= bit[y];
+ h[y] |= bit[x];
+ v2 = conncontent(h,m,n-1);
+ return -v1 - v2;
+ }
+
+/* Case of more than 2/3 dense but not complete */
+
+ if (3*ne > n*n-n)
+ {
+ j = FIRSTBITNZ(g[minv] ^ bit[minv] ^ ALLMASK(n)); /* non-neighbour */
+
+ g[minv] ^= bit[j];
+ g[j] ^= bit[minv];
+ v1 = conncontent(g,m,n);
+ g[minv] ^= bit[j];
+ g[j] ^= bit[minv];
+
+ contract1(g,h,minv,j,n);
+ v2 = conncontent(h,m,n-1);
+
+ return v1 + v2;
+ }
+
+/* All remaining cases */
+
+ j = FIRSTBITNZ(g[minv]); /* neighbour */
+
+ g[minv] ^= bit[j];
+ g[j] ^= bit[minv];
+ v1 = conncontent(g,m,n);
+ g[minv] ^= bit[j];
+ g[j] ^= bit[minv];
+
+ contract1(g,h,minv,j,n);
+ v2 = conncontent(h,m,n-1);
+
+ return v1 - v2;
+}
+
+boolean
+stronglyconnected(graph *g, int m, int n)
+/* test if digraph g is strongly connected */
+{
+ int sp,v,vc;
+ int numvis;
+ set *gv;
+#if MAXN
+ int num[MAXN],lowlink[MAXN],stack[MAXN];
+#else
+ DYNALLSTAT(int,num,num_sz);
+ DYNALLSTAT(int,lowlink,lowlink_sz);
+ DYNALLSTAT(int,stack,stack_sz);
+#endif
+
+#if !MAXN
+ DYNALLOC1(int,num,num_sz,n,"stronglyconnected");
+ DYNALLOC1(int,lowlink,lowlink_sz,n,"stronglyconnected");
+ DYNALLOC1(int,stack,stack_sz,n,"stronglyconnected");
+#endif
+
+ if (n == 0) return FALSE;
+
+ num[0] = 0;
+ for (v = 1; v < n; ++v) num[v] = -1;
+ lowlink[0] = 0;
+ numvis = 1;
+ sp = 0;
+ v = 0;
+ vc = -1;
+ gv = (set*)g;
+
+ for (;;)
+ {
+ vc = nextelement(gv,m,vc);
+ if (vc < 0)
+ {
+ if (sp == 0) break;
+ if (lowlink[v] == num[v]) return FALSE;
+ vc = v;
+ v = stack[--sp];
+ gv = GRAPHROW(g,v,m);
+ if (lowlink[vc] < lowlink[v]) lowlink[v] = lowlink[vc];
+ }
+ else if (num[vc] < 0)
+ {
+ stack[++sp] = vc;
+ v = vc;
+ gv = GRAPHROW(g,v,m);
+ vc = -1;
+ lowlink[v] = num[v] = numvis++;
+ }
+ else if (vc != v)
+ {
+ if (num[vc] < lowlink[v]) lowlink[v] = num[vc];
+ }
+ }
+
+ return numvis == n;
+}
+
+/**************************************************************************/
+
+long
+numsquares(graph *g, int m, int n)
+/* Number of 4-cycles. Undirected graphs only, loops ok. */
+{
+ setword w,bitij;
+ int i,j,k;
+ graph *gi,*gj;
+ unsigned long total,t;
+ boolean iloop,jloop;
+
+ total = 0;
+
+ if (m == 1)
+ {
+ for (j = 1; j < n; ++j)
+ for (i = 0; i < j; ++i)
+ {
+ w = (g[i] & g[j]) & ~(bit[i] | bit[j]);
+ t = POPCOUNT(w);
+ total += t*(t-1)/2;
+ }
+ return (long)(total / 2);
+ }
+ else
+ {
+ for (j = 1, gj = g + m; j < n; ++j, gj += m)
+ {
+ jloop = ISELEMENT(gj,j);
+ if (jloop) DELELEMENT(gj,j);
+ for (i = 0, gi = g; i < j; ++i, gi += m)
+ {
+ iloop = ISELEMENT(gi,i);
+ if (iloop) DELELEMENT(gi,i);
+ t = 0;
+ for (k = 0; k < m; ++k)
+ t += POPCOUNT(gi[k]&gj[k]);
+ total += t*(t-1)/2;
+ if (iloop) ADDELEMENT(gi,i);
+ }
+ if (jloop) ADDELEMENT(gj,j);
+ }
+ return (long)(total / 2);
+ }
+}
+
+/**************************************************************************/
+
+long
+numdiamonds(graph *g, int m, int n)
+/* Number of diamonds (squares with diagonal). Undirected only, no loops. */
+{
+ int i,j,k;
+ setword w;
+ long ans,c;
+ set *gi,*gj;
+
+ ans = 0;
+
+ if (m == 1)
+ {
+ for (i = 0; i < n; ++i)
+ {
+ w = g[i] & BITMASK(i);
+ while (w)
+ {
+ TAKEBIT(j,w);
+ c = POPCOUNT(g[i]&g[j]);
+ ans += c*(c-1)/2;
+ }
+ }
+ }
+ else
+ {
+ for (i = 0, gi = g; i < n; ++i, gi += m)
+ {
+ for (j = i; (j = nextelement(gi,m,j)) >= 0; )
+ {
+ gj = g + m*j;
+ c = 0;
+ for (k = 0; k < m; ++k) c += POPCOUNT(gi[k]&gj[k]);
+ ans += c*(c-1)/2;
+ }
+ }
+ }
+
+ return ans;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-17 07:27:23 +0000