1 /* gutil2.c: Some more graph utilities. */
2
3 #include "gtools.h"
4 #include "gutils.h"
5
6 /**************************************************************************/
7
8 int
loopcount(graph * g,int m,int n)9 loopcount(graph *g, int m, int n)
10 /* Number of loops */
11 {
12 set *gi;
13 int i,nl;
14
15 nl = 0;
16 for (i = 0, gi = g; i < n; ++i, gi += m)
17 if (ISELEMENT(gi,i)) ++nl;
18
19 return nl;
20 }
21
22 /**************************************************************************/
23
24 long
pathcount1(graph * g,int start,setword body,setword last)25 pathcount1(graph *g, int start, setword body, setword last)
26 /* Number of paths in g starting at start, lying within body and
27 ending in last. {start} and last should be disjoint subsets of body. */
28 {
29 long count;
30 setword gs,w;
31 int i;
32
33 gs = g[start];
34 w = gs & last;
35 count = POPCOUNT(w);
36
37 body &= ~bit[start];
38 w = gs & body;
39 while (w)
40 {
41 TAKEBIT(i,w);
42 count += pathcount1(g,i,body,last&~bit[i]);
43 }
44
45 return count;
46 }
47
48 /**************************************************************************/
49
50 long
cyclecount1(graph * g,int n)51 cyclecount1(graph *g, int n)
52 /* The total number of cycles in g (assumed no loops), m=1 only */
53 {
54 setword body,nbhd;
55 long total;
56 int i,j;
57
58 body = ALLMASK(n);
59 total = 0;
60
61 for (i = 0; i < n-2; ++i)
62 {
63 body ^= bit[i];
64 nbhd = g[i] & body;
65 while (nbhd)
66 {
67 TAKEBIT(j,nbhd);
68 total += pathcount1(g,j,body,nbhd);
69 }
70 }
71
72 return total;
73 }
74
75 /**************************************************************************/
76
77 long
cyclecount(graph * g,int m,int n)78 cyclecount(graph *g, int m, int n)
79 /* The total number of cycles in g (assumed no loops) */
80 {
81 if (n == 0) return 0;
82 if (m == 1) return cyclecount1(g,n);
83
84 gt_abort(">E cycle counting is only implemented for n <= WORDSIZE\n");
85 return 0;
86 }
87
88 /**************************************************************************/
89
90 long
indpathcount1(graph * g,int start,setword body,setword last)91 indpathcount1(graph *g, int start, setword body, setword last)
92 /* Number of induced paths in g starting at start, extravertices within
93 * body and ending in last.
94 * {start}, body and last should be disjoint. */
95 {
96 long count;
97 setword gs,w;
98 int i;
99
100 gs = g[start];
101 w = gs & last;
102 count = POPCOUNT(w);
103
104 w = gs & body;
105 while (w)
106 {
107 TAKEBIT(i,w);
108 count += indpathcount1(g,i,body&~gs,last&~bit[i]&~gs);
109 }
110
111 return count;
112 }
113
114 /**************************************************************************/
115
116 long
indcyclecount1(graph * g,int n)117 indcyclecount1(graph *g, int n)
118 /* The total number of induced cycles in g (assumed no loops), m=1 only */
119 {
120 setword body,last,cni;
121 long total;
122 int i,j;
123
124 body = ALLMASK(n);
125 total = 0;
126
127 for (i = 0; i < n-2; ++i)
128 {
129 body ^= bit[i];
130 last = g[i] & body;
131 cni = g[i] | bit[i];
132 while (last)
133 {
134 TAKEBIT(j,last);
135 total += indpathcount1(g,j,body&~cni,last);
136 }
137 }
138
139 return total;
140 }
141
142 /**************************************************************************/
143
144 long
indcyclecount(graph * g,int m,int n)145 indcyclecount(graph *g, int m, int n)
146 /* The total number of induced cycles in g (assumed no loops) */
147 {
148 if (n == 0) return 0;
149 if (m == 1) return indcyclecount1(g,n);
150
151 gt_abort(
152 ">E induced cycle counting is only implemented for n <= WORDSIZE\n");
153 return 0;
154 }
155
156 /**************************************************************************/
157
158 long
numtriangles1(graph * g,int n)159 numtriangles1(graph *g, int n)
160 /* The number of triangles in g; undirected only */
161 {
162 int i,j;
163 setword gi,w;
164 long total;
165
166 total = 0;
167 for (i = 0; i < n-2; ++i)
168 {
169 gi = g[i] & BITMASK(i);
170 while (gi)
171 {
172 TAKEBIT(j,gi);
173 w = g[j] & gi;
174 if (w) total += POPCOUNT(w);
175 }
176 }
177
178 return total;
179 }
180
181 /**************************************************************************/
182
183 long
numtriangles(graph * g,int m,int n)184 numtriangles(graph *g, int m, int n)
185 /* The number of triangles in g; undirected only */
186 {
187 int i,j,k,kw;
188 setword *gi,*gj,w;
189 long total;
190
191 if (m == 1) return numtriangles1(g,n);
192
193 total = 0;
194 for (i = 0, gi = g; i < n-2; ++i, gi += m)
195 for (j = i; (j = nextelement(gi,m,j)) > 0; )
196 {
197 gj = GRAPHROW(g,j,m);
198 kw = SETWD(j);
199 w = gi[kw] & gj[kw] & BITMASK(SETBT(j));
200 if (w) total += POPCOUNT(w);
201 for (k = kw+1; k < m; ++k)
202 {
203 w = gi[k] & gj[k];
204 if (w) total += POPCOUNT(w);
205 }
206 }
207
208 return total;
209 }
210
211 /**************************************************************************/
212
213 long
numdirtriangles(graph * g,int m,int n)214 numdirtriangles(graph *g, int m, int n)
215 /* The number of directed triangles in g */
216 {
217 long total;
218 int i,j,k;
219 set *gi,*gj;
220
221 total = 0;
222 for (i = 0, gi = g; i < n-2; ++i, gi += m)
223 for (j = i; (j = nextelement(gi,m,j)) >= 0;)
224 {
225 gj = GRAPHROW(g,j,m);
226 for (k = i; (k = nextelement(gj,m,k)) >= 0; )
227 if (k != j && ISELEMENT(GRAPHROW(g,k,m),i)) ++total;
228 }
229
230 return total;
231 }
232
233 /**************************************************************************/
234
235 void
commonnbrs(graph * g,int * minadj,int * maxadj,int * minnon,int * maxnon,int m,int n)236 commonnbrs(graph *g, int *minadj, int *maxadj, int *minnon, int *maxnon,
237 int m, int n)
238 /* Count the common neighbours of pairs of vertices, and give the minimum
239 and maximum for adjacent and non-adjacent vertices. Undirected only.
240 Null minimums are n+1 and null maximums are -1.
241 */
242 {
243 int j,k;
244 int mina,maxa,minn,maxn;
245 int cn;
246 set *gi,*gj;
247 setword w;
248
249 if (n == 0)
250 {
251 *minadj = *maxadj = *minnon = *maxnon = 0;
252 return;
253 }
254
255 mina = minn = n+1;
256 maxa = maxn = -1;
257
258 for (j = 0, gj = g; j < n; ++j, gj += m)
259 for (gi = g; gi != gj; gi += m)
260 {
261 cn = 0;
262 for (k = 0; k < m; ++k)
263 {
264 w = gi[k] & gj[k];
265 if (w) cn += POPCOUNT(w);
266 }
267
268 if (ISELEMENT(gi,j))
269 {
270 if (cn < mina) mina = cn;
271 if (cn > maxa) maxa = cn;
272 }
273 else
274 {
275 if (cn < minn) minn = cn;
276 if (cn > maxn) maxn = cn;
277 }
278 }
279
280 *minadj = mina;
281 *maxadj = maxa;
282 *minnon = minn;
283 *maxnon = maxn;
284 }
285
286 /**************************************************************************/
287
288 void
delete1(graph * g,graph * h,int v,int n)289 delete1(graph *g, graph *h, int v, int n)
290 /* Delete vertex v from g, result in h */
291 {
292 setword mask1,mask2,gi;
293 int i;
294
295 mask1 = ALLMASK(v);
296 mask2 = BITMASK(v);
297
298 for (i = 0; i < v; ++i)
299 {
300 gi = g[i];
301 h[i] = (gi & mask1) | ((gi & mask2) << 1);
302 }
303 for (i = v; i < n-1; ++i)
304 {
305 gi = g[i+1];
306 h[i] = (gi & mask1) | ((gi & mask2) << 1);
307 }
308 }
309
310 /**************************************************************************/
311
312 void
contract1(graph * g,graph * h,int v,int w,int n)313 contract1(graph *g, graph *h, int v, int w, int n)
314 /* Contract distinct vertices v and w (not necessarily adjacent)
315 with result in h. No loops are created. */
316 {
317 int x,y;
318 setword bitx,bity,mask1,mask2;
319 int i;
320
321 if (w < v)
322 {
323 x = w;
324 y = v;
325 }
326 else
327 {
328 x = v;
329 y = w;
330 }
331
332 bitx = bit[x];
333 bity = bit[y];
334 mask1 = ALLMASK(y);
335 mask2 = BITMASK(y);
336
337 for (i = 0; i < n; ++i)
338 if (g[i] & bity)
339 h[i] = (g[i] & mask1) | bitx | ((g[i] & mask2) << 1);
340 else
341 h[i] = (g[i] & mask1) | ((g[i] & mask2) << 1);
342
343 h[x] |= h[y];
344 for (i = y+1; i < n; ++i) h[i-1] = h[i];
345 h[x] &= ~bitx;
346 }
347
348 /**************************************************************************/
349
350 static TLS_ATTR int knm[18][16]; /* knm[n,m] = conncontent(K_n - m*K_2) */
351 static TLS_ATTR boolean knm_computed = FALSE;
352
353 int
conncontent(graph * g,int m,int n)354 conncontent(graph *g, int m, int n)
355 /* number of connected spanning subgraphs with an even number
356 of edges minus the number with an odd number of edges */
357 {
358 graph h[WORDSIZE];
359 setword gj;
360 int i,j,v1,v2,x,y;
361 int minv,mindeg,deg,goodv;
362 long ne;
363
364 if (m > 1) ABORT("conncontent only implemented for m=1");
365
366 /* First handle tiny graphs */
367
368 if (n <= 3)
369 {
370 if (n == 1) return 1;
371 if (n == 2) return (g[0] ? -1 : 0);
372 if (!g[0] || !g[1] || !g[2]) return 0; /* disconnected */
373 if (g[0]^g[1]^g[2]) return 1; /* path */
374 return 2; /* triangle */
375 }
376
377 /* Now compute
378 ne = number of edges
379 mindeg = minimum degree
380 minv = a vertex of minimum degree
381 goodv = a vertex with a clique neighbourhood (-1 if none)
382 */
383
384 mindeg = n;
385 ne = 0;
386 goodv = -1;
387 for (j = 0; j < n; ++j)
388 {
389 gj = g[j];
390 deg = POPCOUNT(gj);
391 ne += deg;
392 if (deg < mindeg)
393 {
394 mindeg = deg;
395 minv = j;
396 if (deg == 1) goodv = j;
397 }
398 if (deg >= 3 && deg <= 4 && goodv < 0)
399 {
400 while (gj)
401 {
402 TAKEBIT(i,gj);
403 if (gj & ~g[i]) break;
404 }
405 if (!gj) goodv = j;
406 }
407 }
408 ne /= 2;
409
410 /* Cases of isolated vertex or tree */
411
412 if (mindeg == 0) return 0;
413
414 #if 0
415 if (mindeg == 1 && ne == n-1)
416 {
417 if (isconnected1(g,n)) return ((n&1) ? 1 : -1);
418 else return 0;
419 }
420 #endif
421
422 /* Cases of clique and near-clique */
423
424 if (mindeg == n-1)
425 {
426 j = -1;
427 for (i = 2; i < n; ++i) j *= -i;
428 return j;
429 }
430
431 if (mindeg == n-2 && n < 16)
432 {
433 if (!knm_computed)
434 {
435 knm_computed = TRUE;
436 knm[1][0] = 1;
437 for (i = 2; i < 16; ++i)
438 {
439 knm[i][0] = -knm[i-1][0] * (i-1);
440 for (j = 1; j+j <= i; ++j)
441 knm[i][j] = knm[i][j-1] + knm[i-1][j-1];
442 }
443 }
444 return knm[n][(n*n-n)/2-ne];
445 }
446
447 /* Case of vertex with clique neighbourhood */
448
449 if (goodv >= 0)
450 {
451 delete1(g,h,goodv,n);
452 return -POPCOUNT(g[goodv]) * conncontent(h,m,n-1);
453 }
454
455 /* Case of minimum degree 2 */
456
457 if (mindeg == 2)
458 {
459 x = FIRSTBITNZ(g[minv]);
460 y = FIRSTBITNZ(g[minv]^bit[x]);
461 if (x > minv) --x;
462 if (y > minv) --y;
463 delete1(g,h,minv,n);
464 v1 = conncontent(h,m,n-1);
465 if (h[x] & bit[y]) return -2*v1; /* adjacent neighbours */
466
467 h[x] |= bit[y];
468 h[y] |= bit[x];
469 v2 = conncontent(h,m,n-1);
470 return -v1 - v2;
471 }
472
473 /* Case of more than 2/3 dense but not complete */
474
475 if (3*ne > n*n-n)
476 {
477 j = FIRSTBITNZ(g[minv] ^ bit[minv] ^ ALLMASK(n)); /* non-neighbour */
478
479 g[minv] ^= bit[j];
480 g[j] ^= bit[minv];
481 v1 = conncontent(g,m,n);
482 g[minv] ^= bit[j];
483 g[j] ^= bit[minv];
484
485 contract1(g,h,minv,j,n);
486 v2 = conncontent(h,m,n-1);
487
488 return v1 + v2;
489 }
490
491 /* All remaining cases */
492
493 j = FIRSTBITNZ(g[minv]); /* neighbour */
494
495 g[minv] ^= bit[j];
496 g[j] ^= bit[minv];
497 v1 = conncontent(g,m,n);
498 g[minv] ^= bit[j];
499 g[j] ^= bit[minv];
500
501 contract1(g,h,minv,j,n);
502 v2 = conncontent(h,m,n-1);
503
504 return v1 - v2;
505 }
506
507 boolean
stronglyconnected(graph * g,int m,int n)508 stronglyconnected(graph *g, int m, int n)
509 /* test if digraph g is strongly connected */
510 {
511 int sp,v,vc;
512 int numvis;
513 set *gv;
514 #if MAXN
515 int num[MAXN],lowlink[MAXN],stack[MAXN];
516 #else
517 DYNALLSTAT(int,num,num_sz);
518 DYNALLSTAT(int,lowlink,lowlink_sz);
519 DYNALLSTAT(int,stack,stack_sz);
520 #endif
521
522 #if !MAXN
523 DYNALLOC1(int,num,num_sz,n,"stronglyconnected");
524 DYNALLOC1(int,lowlink,lowlink_sz,n,"stronglyconnected");
525 DYNALLOC1(int,stack,stack_sz,n,"stronglyconnected");
526 #endif
527
528 if (n == 0) return FALSE;
529
530 num[0] = 0;
531 for (v = 1; v < n; ++v) num[v] = -1;
532 lowlink[0] = 0;
533 numvis = 1;
534 sp = 0;
535 v = 0;
536 vc = -1;
537 gv = (set*)g;
538
539 for (;;)
540 {
541 vc = nextelement(gv,m,vc);
542 if (vc < 0)
543 {
544 if (sp == 0) break;
545 if (lowlink[v] == num[v]) return FALSE;
546 vc = v;
547 v = stack[--sp];
548 gv = GRAPHROW(g,v,m);
549 if (lowlink[vc] < lowlink[v]) lowlink[v] = lowlink[vc];
550 }
551 else if (num[vc] < 0)
552 {
553 stack[++sp] = vc;
554 v = vc;
555 gv = GRAPHROW(g,v,m);
556 vc = -1;
557 lowlink[v] = num[v] = numvis++;
558 }
559 else if (vc != v)
560 {
561 if (num[vc] < lowlink[v]) lowlink[v] = num[vc];
562 }
563 }
564
565 return numvis == n;
566 }
567