\datethis @i gb_types.w @*Intro. This is a simple program to find all maximal cliques of a graph, when the graph has at most 64 vertices and is given by its adjacency matrix. The algorithm, due to Moody and Hollis, is described in Appendix 5 of the book {\sl Mathematical Taxonomy\/} by Jardine and Sibson (1971). I'm writing this program hastily today because it's a nice example of bitwise manipulation in graph theory, probably suitable for Section 7.1.3 of {\sl The Art of Computer Programming}. In brief, we have vectors $\rho_v$ and $\delta_v$ for each vertex~$v$, where $\rho_v$ is $v$'s row in the adjacency matrix and $\delta_v$ is all 1s except for 0 in position~$v$. (Position~$v$ in $\rho_v$ is always~1; in other words, we assume that each vertex is adjacent to itself.) One can easily show that a set of vertices~$q$ is a clique if and only if $q$ is the bitwise intersection of $(v{\in}q?\,\rho_v{:}\,\delta_v)$ over all vertices~$v$. To find all cliques, we could compute all $2^n$ such bitwise intersections, discarding duplicates. To find all maximal cliques, we could compute all $2^n$ such bitwise intersections and discard any clique $q$ that is found to be contained in another. To speed the process up, we can do those discards much more cleverly. The input graph is specified by an ASCII file in Stanford GraphBase format. The name of that file, say \.{foo.gb}', is the one-and-only command-line parameter. I've instrumented this program to see how many memory references it makes (|mems|) and how many words of workspace it needs (|space|), exclusive of input-output. (Note: I consider a clique'' to be any complete subgraph of a graph. Many of the older books on graph theory define it to be a {\it maximal\/} complete subgraph; but that terminology is now dying out. The earlier definition is less desirable, because for example we'd like a clique of~$G$ to be equivalent to an independent set of $\overline G$.) Beware: Memory overflow is not checked. This is not intended to be a robust program; I wrote it only as a experiment. @d size 100000 /* size reserved for the workspace */ @d o mems++ @d oo mems+=2 @c #include #include "gb_graph.h" #include "gb_save.h" unsigned long long rho[64]; /* rows of the adjacency matrix */ unsigned long long work[size]; unsigned int mems; int space; char table[64]; main(int argc, char* argv[]) { register int i,j,k,l,m,n,p,q,r; register Graph *g; register Arc *a; unsigned long long u,v,w; @; @; @; printf("(The computation took %u mems, using %d words of workspace.)\n", mems,space); } @ @= if (argc!=2) { fprintf(stderr,"Usage: %s inputgraph.gb\n",argv[0]); exit(-1); } g=restore_graph(argv[1]); if (!g) { fprintf(stderr,"I can't input the graph %s (panic code %ld)!\n", argv[1],panic_code); exit(-2); } n=g->n; if (n>64) { fprintf(stderr,"Sorry, that graph has %d vertices; ",n); fprintf(stderr,"I can't handle more than 64!\n"); exit(-3); } @; @ @= for (j=0;jvertices+j)->arcs;a;a=a->next) w|=1LL<<(a->tip-g->vertices); rho[j]=w; } @ @d deBruijn 0x03f79d71b4ca8b09 /* the least de Bruijn cycle of length 64 */ @= printf("Graph %s has %d maximal cliques:\n",g->id,m); @; for (k=1;k<=m;k++) { for (w=work[k];w;w=w^v) { v=w&-w; u=v*deBruijn; j=table[u>>58]; printf(" %s",(g->vertices+j)->name); } printf("\n"); } @ The ruler-function calculation in the previous section isn't part of the inner loop; so I could use a slower, brute-force scheme that simply examines each bit, one bit at a time. But I'm using the de Bruijn method anyway, in order to get experience with it. @= for (j=0,v=1;v;j++,v<<=1) { u=v*deBruijn; table[u>>58]=j; } @*The algorithm. The |work| area holds all maximal elements of the $2^i$ bitwise-ands of either $\rho_v$ or $\delta_v$ for the first $i$ vertices. There are $m$ of them. A vertex that's adjacent to all other vertices is just carried along in all entries, so we needn't bother with it. @= w=1LL<<(n-1); oo,work[1]=work[size-1]=(w<<1)-1; /* $n$ 1s */ m=1, space=4; for (i=0;i; @; } @ To visualize the current situation, suppose bit~$v$ is the leftmost. Then we want to rearrange |work[1]| thru |work[m]| so that they have the form \vcenter{\halign{#\thinspace&#\hfil\cr 0&\alpha_1\cr &\vdots\cr 0&\alpha_k\cr 1&\beta_1\cr &\vdots\cr 1&\beta_{m-k}\cr}} This loop is like the partitioning step of radix-exchange. There always is at least one element $u$ with $u\AND v\ne0$. For convenience (and speed) we keep |work[0]=0|. @= j=1, k=m; while (1) { while ((o,work[j]&v)==0) j++; while ((o,work[k]&v)!=0) k--; if (j>k) break; oo,u=work[j], work[j]=work[k], work[k]=u; j++,k--; } @ Now comes the interesting part. At this point $j=k+1$. Entries 1 thru $k$ of the workspace should simply be carried over to the next round. For if $u=0\,\alpha_i$ in the notation above, we have $u\AND\rho_i\subseteq u\AND\delta_i=u$; and $u$~(which is currently maximal) won't be contained in any other entries. Thus we focus our attention on the remaining entries in the current workspace, namely |work[j]| thru |work[m]|, which need to be `split.'' Let $u$ be the current entry; we want to split it into $u'=u\AND\rho_v$ and $u''=u\AND\delta_v$. The first one, $u'=(1\beta)\AND\rho_v$, must be checked against all other first entries before it is accepted for the new round. (It can't be contained in a second entry, because it has 1~in position~$v$.) If it is contained in a first entry already generated, we drop it. But if it contains one of those entries, we might need to drop more than one of them. The second entry, $u''=(1\beta)\AND\delta_v=0\beta$, is fairly easy to deal with. It can be contained in some first entry $(1\beta')\AND\rho_v$ only if $\beta\subseteq\beta'$; hence $\beta=\beta'$ and $1\beta\subseteq\rho_v$. If $1\beta\not\subseteq\rho_v$, we accept $0\beta$ if and only if it isn't contained in any held-over entry $0\alpha$. In the following steps, positions $p$ thru |size-2| of the workspace hold the first entries $u'$ that have been tentatively accepted. Positions $k+1$ thru $l$ hold the second entries $u'$ that have definitely been accepted. Notice that the resulting algorithm avoids linked memory, so it ought to be cache-friendly. @= for (l=k,p=size-1;j<=m;j++) { o,u=work[j],q=size-2; w=u&rho[i]; /* $w=u'$; we've already fetched |rho[i]| from memory */ if (u!=w) { for (;q>=p;q--) { if ((o,w&work[q])==w) goto second_entry; if ((w&work[q])==work[q]) goto absorb; } o,work[--p]=w; /* accept $u'$, tentatively */ if (space; if (u==w) continue; second_entry: w=u&~v; /* $w=u''$ */ for (q=1;q<=k;q++) if ((o,w&work[q])==w) goto done_with_u; o,work[++l]=w; /* accept $u''$ */ done_with_u: continue; } for (m=l;p= o, r=p, work[p-1]=0; while (1) { while (o,(w|work[q])!=w) q--; while (o,(w|work[r])==w) r++; if (q