% (c) 2005 D E Knuth; use and enjoy, but please don't blame me for mistakes... \datethis \def\title{ZEILBERGER} \def\[#1]{[\hbox{$\mkern1mu\thickmuskip=\thinmuskip#1\mkern1mu$}]} % Iverson @*Introduction. This little program implements Zeilberger's bijection between $n$-node forests with pruning order~$s$ and $n$-node binary trees with Strahler number~$s$. [Doron Zeilberger, A bijection from ordered trees to binary trees that sends the pruning order to the Strahler number,'' {\sl Discrete Mathematics\/ \bf82} (1990), 89--92.] As Zeilberger says in his paper, First, definitions!'' I won't define forests (by which I always mean {\it ordered forests}, as in {\sl The Art of Computer Programming\/}), nor need I define binary trees. But I should explain the concepts of pruning order and Strahler number, since they are less familiar. A {\it filament\/} of a forest is a sequence of one or more nodes $x_1$, \dots,~$x_k$, where $x_{j+1}$ is the only child of $x_j$ for $1\le j int l[n+2], r[n+1]; /* leftmost child and right sibling */ int ll[n+1], rr[n+1]; /* links of the binary tree */ int lll[n+1], rrr[n+1]; /* links of the inverse forest */ int q[n+1], s[n+1]; /* data needed by Zeilberger's algorithm */ int serial; /* total number of cases checked */ int count[10]; /* individual counts by pruning order */ @@; @# main() { register int j,k,y,p; printf("Checking all forests with %d nodes...\n",n); q[0]=0, s[0]=1000000; /* see below */ for (k=1;k; @; @; @ } for (p=1;count[p];p++) printf("Altogether %d cases with pruning order %d.\n",count[p],p); } @ @= for (j=1;!l[j];j++) r[j]=0, l[j]=j+1; if (j>n) break; for (k=0,y=l[j];r[y];k=y,y=r[y]); if (k>0) r[k]=0;@+else l[j]=0; r[y]=r[j], r[j]=y; @* The main algorithm. The nodes in the forest are represented by positive integers, according to their rank in preorder. For each node~$x>0$, we let$l[x]$be its leftmost child (namely$x+1$), if it has one, but$l[x]=0$when$x$is childless. Similarly,$r[x]$is$x$'s right sibling, or$r[x]=0$when$x$is rightmost in its family. Zeilberger's method implicitly begins by determining the pruning order of the given forest and its principal subforests. We can carry this out by defining three numbers$p[x]$,$q[x]$, and$s[x]$for each node~$x$; here$p[x]$is the pruning order of the subtree rooted at~$x$,$q[x]$is the maximum pruning order of$x$and all of its right siblings, and$s[x]$is the number of nodes in which that maximum occurs. (Zeilberger called$s[l[x]]$the number of skewers'' of~$x$.) Formally, we can compute$p[x]$,$q[x]$, and$s[x]$via the following recursive definitions, if we set$q[0]=0$and$s[0]=\infty$:$\openup1\jot\eqalign{ p[x]&=q[l[x]]+\[s[l[x]]>1];\cr q[x]&=\max(p[x],q[r[x]]);\cr s[x]&=\cases{1,&ifp[x]>q[r[x]]$;\cr s[r[x]]+1,&if$p[x]=q[r[x]]$;\cr s[r[x]],&if$p[x]= void label(register int x) { register int p,qr; if (l[x]) label(l[x]); if (r[x]) label(r[x]); p=q[l[x]]+(s[l[x]]>1); qr=q[r[x]]; if (p>qr) q[x]=p, s[x]=1; else q[x]=qr,s[x]=s[r[x]]+(p==qr); } @ Zeilberger's bijection is recursively defined too, and we will implement it by writing another recursive procedure. But before we do so, let's examine how that procedure will be invoked in the program. @= for (k=1;k<=n;k++) ll[k]=l[k], rr[k]=r[k]; /* clone the forest */ label(1); /* compute all the $q$'s and $s$'s */ zeil(1); /* transform the binary tree */ @ When subroutine |zeil(x)| is called, |x| is a node of a forest, represented via left-child and right-sibling links |ll| and |rr|. The mission of |zeil(x)| is to transform that forest in such a way that |x| will be a node of the final binary tree, while |ll[x]| and |rr[x]| will be the binary subtrees that are obtained by applying |zeil| to two subforests. Zeilberger's paper considered three cases, depending on the relative sizes of |q[x]|, |q[ll[x]]|, and |q[rr[x]]|. In Case 1, no change is needed; in Case 2, some of |x|'s children are promoted to be siblings of~|x|; in Case 3, some of |x|'s children swap places with all of the right siblings. @ @= void zeil(register int x) { register int k,ql,qr,p,y,yy,z,zz,ss; k=q[x], ql=q[ll[x]], qr=q[rr[x]]; if (ql==k) @@; else if (qr; /* otherwise we do Zeilberger's Case 1, which involves no action */ if (ll[x]) zeil(ll[x]); if (rr[x]) zeil(rr[x]); } @ Case 2, the tricky case, detaches all but the first of |x|'s skewers.'' When right links change, we must update the |q| and |s| numbers in each subnode of~|x| before applying |zeil| to the new subforests. In particular, we must recompute all the |q|'s and |s|'s for the children of~|x| that lie between the first and second skewer, because those nodes will no longer lie in the shadow of the second skewer. This recomputation is a bit tricky because it must be done bottom-up, while our links go downward. An auxiliary recursive procedure could be introduced at this point, in order to achieve bottom-up behavior. But there's a more efficient (and more fun) way to do the job, namely to reverse the links temporarily as we descend the chain, and to reverse them again as we go back up. @= { for (y=ll[x],ss=s[y];s[y]==ss;yy=y,y=rr[yy]) s[y]=1; /* at this point node |yy| is the first skewer of |x| */ for (z=yy,ss--; s[rr[y]]==ss && q[rr[y]]==ql;@| yy=z, z=y, y=rr[z], rr[z]=yy); /* reverse links */ @t@>@; /* now node |y| is the second skewer of |x| */ @t@>@; /* its left sibling, |z|, will become the last child of |x| */ if (z==yy) rr[z]=0; /* easy case: the skewers were adjacent */ else@+ for (zz=0;rr[z]!=zz; yy=zz,zz=z,z=rr[zz],rr[zz]=yy) { p=q[ll[z]]+(s[ll[z]]>1); /* we will recompute |q[z]| and |s[z]| */ if (p>q[zz]) q[z]=p, s[z]=1; else q[z]=q[zz], s[z]=s[zz]+(p==q[zz]); } for (zz=x,z=rr[zz];z;zz=z,z=rr[zz]) q[z]=ql, s[z]=ss; rr[zz]=y; /* |y| and its right siblings become |x|'s final siblings */ } @ In Case 3, |x| has only one skewer. But the transformation in this case demotes siblings to children, where they might become additional skewers. It also might promote some children to siblings. @= { if (qr==k) ss=s[rr[x]]+1;@+else ss=1; for (y=ll[x];q[y]==k;yy=y,y=rr[yy]) s[y]=ss; rr[yy]=rr[x], rr[x]=y; } @* Checking the Strahler number. If our implementation of Zeilberger's transformation is correct, it will have set |q[x]| to the Strahler number of the binary subtree rooted at |x| with respect to the |ll| and |rr| links, for every node~|x|. Therefore we want to check this condition. And we might as well do the checking by brute force, so that the evidence is convincing. @ @= int strahler(register int x) { register int sl,sr,s; if (ll[x]) sl=strahler(ll[x]); else sl=0; if (rr[x]) sr=strahler(rr[x]); else sr=0; s=(sl>sr? sl: sl= count[strahler(1)]++; serial++; @* The inverse algorithm. The evidence of correctness is mounting. But our argument in favor of Zeilberger's algorithm is still not compelling, because there are lots of ways to convert a forest with pruning order~$s$ into a binary tree with Strahler number~$s$. For example, we need only compute the pruning order, then choose our favorite binary tree that has the desired Strahler number. A bijection must do more than this: It must not destroy information. We must be able to go back from each binary tree to the original forest that produced it. Thus, our implementation of Zeilberger's procedure is incomplete until we have also implemented its inverse, |unzeil|. Our |unzeil| algorithm will recreate the forest in new arrays |lll| and |rrr|, just to emphasize that no cheating is going on. @= for (k=1;k<=n;k++) lll[k]=ll[k], rrr[k]=rr[k], q[k]=s[k]=0; /* clone the binary tree */ unzeil(1); /* attempt to recreate the original forest */ for (k=1;k<=n;k++) if (lll[k]!=l[k] || rrr[k]!=r[k]) fprintf(stderr,"Rejection at node %d of case %d!\n",k,serial); @ The |unzeil| procedure also computes the |q| and |s| values at each node of the reconstructed forest. @= void unzeil(register int x) { register int ql,qr,p,y,yy,z,zz,ss; if (rrr[x]) unzeil(rrr[x]); if (lll[x]) unzeil(lll[x]); ql=q[lll[x]], qr=q[rrr[x]]; if (ql>qr) @@; else if (ql==qr && s[lll[x]]==1) @; /* otherwise we undo Zeilberger's Case 1, which involves no action */ p=ql+(s[lll[x]]>1); if (p>qr) q[x]=p, s[x]=1; else q[x]=qr, s[x]=s[rrr[x]]+(p==qr); } @ @= { if (s[lll[x]]>1) { for (y=lll[x];s[rrr[y]]==s[y];y=rrr[y]) s[y]=1; s[y]=1; }@+else@+ for (y=lll[x];q[rrr[y]]==q[y];y=rrr[y]); yy=rrr[y], rrr[y]=rrr[x], rrr[x]=yy; /* swap siblings with children */ qr=q[yy]; /* the subsequent program assumes that |qr=q[rrr[x]]| */ } @ We have saved the other tricky case for last: Again we do a double-reversal in order to recompute any |q|'s and |s|'s that have been obliterated. In this case the recomputation affects the near siblings of~|x|. @= { for (z=rrr[x],zz=x;s[rrr[z]]==s[z]&&q[rrr[z]]==q[z];@| yy=zz,zz=z,z=rrr[zz],rrr[zz]=yy); if (zz==x) yy=rrr[x]=0; else@+for (yy=0;zz!=x;y=zz,zz=rrr[y],rrr[y]=yy,yy=y) { p=q[lll[zz]]+(s[lll[zz]]>1); /* we will recompute |q[zz]| and |s[zz]| */ if (p>q[yy]) q[zz]=p, s[zz]=1; else q[zz]=q[yy], s[zz]=s[yy]=(p==q[yy]); } /* |x|'s former siblings are now undone */ qr=q[yy],ss=s[z]+1; /* at this point |yy=rrr[x]| */ for (yy=lll[x],y=rrr[yy];y;yy=y,y=rrr[yy]) { s[yy]=ss; if (q[yy]==ql) ss--; } s[yy]=ss, rrr[yy]=z; /* unpromote |x|'s former children */ } @* Index.