@*Intro. Johan de Ruiter presented a beautiful puzzle on 14 March 2018, based on the first 32 digits of~$\pi$. It's a special case of the following self-referential problem: Given a directed graph, find all vertex labelings such that each vertex is labeled with the number of distinct labels on its successors. In Johan's puzzle, some of the labels are given, and we're supposed to find the others. He also presented the digraph in terms of a $10\times10$ array, with each cell pointing either north, south, east, or west; its successors are the cells in that direction. I've written this program so that it could be applied to fairly arbitrary digraphs, if I decide to make it more general. The program uses bitmaps in interesting ways, not complicated. @d N (0<<4) @d S (1<<4) @d E (2<<4) @d W (3<<4) @d debug 1 /* for optional verbose printing */ @d verts 100 /* vertices in the digraph */ @d maxd 9 /* maximum out-degree in the digraph; must be less than 16 */ @d bitmax (1<<(maxd+1)) @d infinity (unsigned long long)(-1) @d o mems++ @d oo mems+=2 @d ooo mems+=3 @c #include #include long mems; char johan[10][10]={@| {S+3,W+1,E+4,W+0,S+1,W+0,S+5,S+0,S+9,S+0},@| {E+0,S+0,W+2,S+6,S+0,E+0,S+0,W+0,E+0,S+5},@| {E+0,S+0,E+0,E+0,S+0,S+0,E+3,S+5,W+8,W+9},@| {E+0,E+0,S+0,N+0,S+0,E+0,W+0,S+0,W+7,W+0},@| {E+9,E+0,S+3,S+0,S+0,S+0,W+0,W+0,S+0,W+0},@| {E+0,E+0,E+0,W+0,S+0,E+0,S+0,E+2,S+0,S+3},@| {E+0,E+8,S+0,N+0,S+0,S+0,N+0,W+0,N+0,W+0},@| {N+4,E+6,S+2,N+6,S+0,E+0,S+0,W+0,S+0,N+0},@| {N+4,E+0,E+0,E+0,S+0,W+0,W+3,W+3,W+0,N+0},@| {E+0,E+8,N+0,W+3,N+0,N+2,W+0,W+7,N+9,N+5}}; int nu[bitmax],gnu[bitmax],un[bitmax]; /* $\nu k$, $2^{\nu k}$, and $\rho[k]$ */ @; @; main() { register int a,d,g,i,j,k,l,q,t,u,v,x; register unsigned long long p; @; @; @; @; @; @; } @ @= for (o,gnu[0]=1,k=0;k>1],nu[k+1]=nu[k]+1,gnu[k]=gnu[k>>1],gnu[k+1]=gnu[k]<<1; for (k=1;k<=maxd;k++) o,un[1<= for (i=0;i<10;i++) for (j=0;j<10;j++) { v=inx(i,j); sprintf(name[v],"%02d", v); known[v]=(johan[i][j]&0xf? johan[i][j]&0xf: -1); d=0; switch (johan[i][j]>>4) { case N>>4: for (k=0;k>4: for (k=9;k>i;k--) newarc(k,j);@+break; case E>>4: for (k=9;k>j;k--) newarc(i,k);@+break; case W>>4: for (k=0;kmaxd) { fprintf(stderr,"The outdegree of %s should be at most %d, not %d!\n", name[v],maxd,d); exit(-1); } o,deg[v]=d; if (d<=1) known[v]=d; /* we can consider this label prespecified */ } @ The set of possible labels for vertex |v| is kept in |bits[v]|, a (|maxd|+1)-bit number. It's either a single bit (if |v|'s label was prespecified) or $1+2+\cdots+2^d$ (if |v| has degree~$d$ and wasn't given a label). @= for (i=0;i<10;i++) for (j=0;j<10;j++) { o,v=inx(i,j),l=johan[i][j]&0xf; if (known[v]>=0) o,bits[v]=1<= int llink[verts+1],rlink[verts+1]; int deg[verts],arcs[verts],scra[verts],bits[verts],isactive[verts],known[verts]; unsigned long long size[verts]; char name[verts][8]; /* each vertex name is assumed to be at most seven characters */ int tip[2*verts*verts],next[2*verts*verts]; int arcptr=0; /* this many entries of |tip| and |next| are in use */ @ @= for (v=0;v= void printvert(int v,FILE*stream) { register int b,d; fprintf(stream,"%s(", name[v]); for (b=bits[v],d=0;(1<fprintf(stream,")"); } @ @= void printact(void) { register int v; for (v=rlink[active];v!=active;v=rlink[v]) { if (llink[v]!=active) fprintf(stderr," "); printvert(v,stderr); } fprintf(stderr,"\n"); } @*The stability test. Here's the fun routine that motivated me to write this program. The total number of solutions $(x_1,\ldots,x_d)$ to |v|'s stability problem is at most |size[v]|. But of course we hope to cut this number way down. The nicest part of the following code is its calculation of |goal| bits, to rule out impossible partial solutions. Once again I follow Algorithm 7.2.2B. @= b1: mems+=4,w[0]=v,wb[0]=bits[v],wbp[0]=0; for (o,a=arcs[v],d=0;a;o,a=next[a],d++) mems+=5,w[d+1]=tip[a],wb[d+1]=bits[tip[a]],wbp[d+1]=0; for (o,k=d,g=bits[v];k;k--) o,goal[k]=g,g=(g|(g>>1))&((1<d) @; o,x=wb[l]&-wb[l]; /* the lowest bit */ b3: oo,s[l]=s[l-1]|x; if (oo,gnu[s[l]]&goal[l]) { o,move[l++]=x; goto b2; } b4:@+for (x<<=1;o,x<=wb[l];x<<=1) if (x&wb[l]) goto b3; b5:@+if (--l) { o,x=move[l]; goto b4; } @; @ @= { if (debug) printsol(d); for (k=1;k= void printsol(int d) { register int k; fprintf(stderr," %s->", name[w[0]]); for (k=1;k<=d;k++) fprintf(stderr,"%d", un[move[k]]); fprintf(stderr,"\n"); } @ If there were no solutions, we've been given an impossible problem. @= for (k=0;k<=d;k++) if (oo,wbp[k]!=wb[k]) { if (wbp[k]==0) { fprintf(stderr,"Contradiction reached while testing stability of %s!\n", name[w[0]]); exit(-666); } o,u=w[k]; oo,bits[u]=wbp[k]; if (debug) { fprintf(stderr," now "); printvert(u,stderr); fprintf(stderr,"\n"); } @; for (o,a=scra[u];a;o,a=next[a]) { o,u=tip[a]; ooo,size[u]=(size[u]/nu[wb[k]])*nu[wbp[k]]; @; } } @ @= if (o,!isactive[u]) { mems+=5,t=llink[active],rlink[t]=llink[active]=u,llink[u]=t,rlink[u]=active; o,isactive[u]=roundno+1; } @ @= int w[maxd+1],wb[maxd+1],wbp[maxd+1],goal[maxd+1],move[maxd+1],s[maxd+1]; @*The main loop. Hurray: We're ready to put everything together. Besides our desire to choose an active item of minimum size, we want to keep cycling through the array. So we choose each item at most once per round.'' The value of |isactive[v]| tells in which round |v| will become activated. @= while (o,rlink[active]!=active) { for (o,u=rlink[active],q=roundno+2;u!=active;o,u=rlink[u]) if (o,isactive[u] "); for (a=arcs[v];a;a=next[a]) printvert(tip[a],stderr); fprintf(stderr,"\n"); } @; } @ @= fprintf(stderr,"Stability achieved after %d tests, %d rounds, %ld mems.\n", tests,roundno,mems); for (v=0;v= int tests,roundno; @*Index.