gb_basic.w

来自「模拟器提供了一个简单易用的平台」· W 代码 · 共 1,575 行 · 第 1/5 页

    panic(no_room); /* out of memory before we're even started */
  sprintf(new_graph->id,"perms(%ld,%ld,%ld,%ld,%ld,%lu,%d)",
    n0,n1,n2,n3,n4,max_inv,directed?1:0);
  strcpy(new_graph->util_types,"VVZZZZZZZZZZZZ"); /* hash table will be used */
}

@ After multiplication by $(1-z^{k+s})/(1-z^k)$, the coefficients of the
power series will be nonnegative, because they are the coefficients of
a $z$-multinomial coefficient.

@<Multiply the power series coefficients by
       $\prod_{1\le k\le n_j}(1-z^{s+k})/(1-z^k)$@>=
for (k=1;k<=nn[j];k++) {@+register long ii;
  for (i=max_inv,ii=i-k-s;ii>=0;ii--,i--) coef[i]-=coef[ii];
  for (i=k,ii=0;i<=max_inv;i++,ii++) {
    coef[i]+=coef[ii];
    if (coef[i]>1000000000) panic(very_bad_specs+1); /* way too big */
  }
}

@ As we generate the permutations, we maintain a table $(y_1,\ldots,y_n)$,
where $y_k$ is the number of
inversions whose first element is the $k$th element of the multiset.
For example, if the multiset is $\{0,0,1,2\}$ and the current permutation is
$(2,0,1,0)$, the inversion table is $(y_1,y_2,y_3,y_4)=(0,0,1,3)$. Clearly
$0\le y_k<k$, and $y_k\le y_{k-1}$ when the $k$th element of the multiset
is the same as the $(k-1)$st element. These conditions are necessary
and sufficient to define a valid inversion table. We will generate
permutations in lexicographic order of their inversion tables.

For convenience, we set up another array~|z|, which holds the
initial inversion-free permutation.

@<Name the permutations and create the arcs or edges@>=
{@+register long *xtab,*ytab,*ztab; /* permutations and their inversions */
  long m=0; /* current number of inversions */
  @<Initialize |xtab|, |ytab|, and |ztab|@>;
  v=new_graph->vertices;
  while (1) {
    @<Assign a symbolic name for $(x_1,\ldots,x_n)$ to vertex~|v|@>;
    @<Create arcs or edges from previous permutations to~|v|@>;
    v++;
    @<Advance to the next perm; |goto last| if there are no more solutions@>;
  }
  last:@+if (v!=new_graph->vertices+new_graph->n)
    panic(impossible); /* can't happen */
  gb_free(working_storage);
}

@ @<Initialize |xtab|, |ytab|, and |ztab|@>=
xtab=gb_typed_alloc(3*n+3,long,working_storage);
if (gb_trouble_code) { /* can't allocate |xtab| */
 gb_recycle(new_graph);@+panic(no_room+2);@+}
ytab=xtab+(n+1);
ztab=ytab+(n+1);
for (j=0,k=1,s=nn[0];;k++) {
  xtab[k]=ztab[k]=j; /* |ytab[k]=0| */
  if (k==s) {
    if (++j>d) break;
    else s+=nn[j];
  }
}

@ Here is the heart of the permutation logic. We find the largest~$k$
such that $y_k$ can legitimately be increased by~1. When we encounter
a~$k$ for which $y_k$ cannot be increased, we set $y_k=0$ and adjust
the $x$'s accordingly. If no $y_k$ can be increased, we are done.

@<Advance to the next perm...@>=
for (k=n;k;k--) {
  if (m<max_inv && ytab[k]<k-1)
    if (ytab[k]<ytab[k-1] || ztab[k]>ztab[k-1]) goto move;
  if (ytab[k]) {
    for (j=k-ytab[k];j<k;j++) xtab[j]=xtab[j+1];
    m-=ytab[k];
    ytab[k]=0;
    xtab[k]=ztab[k];
  }
}
goto last;
move: j=k-ytab[k]; /* the current location of the $k$th element, $z_k$ */
xtab[j]=xtab[j-1];@+xtab[j-1]=ztab[k];
ytab[k]++;@+m++;

@ A permutation is encoded as a sequence of nonblank characters,
using an abbreviated copy of the |imap| code from {\sc GB\_\,IO} and omitting
the characters that need to be quoted within strings. If the
number of distinct elements in the multiset is at most~62, only digits
and letters will appear in the vertex name.

@<Private variables@>=
static char *short_imap="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ\
abcdefghijklmnopqrstuvwxyz_^~&@@,;.:?!%#$+-*/|<=>()[]{}`'";

@ @<Assign a symbolic name for $(x_1,\ldots,x_n)...@>=
{@+register char *p; register long *q;
  for (p=&buffer[n-1],q=&xtab[n];q>xtab;p--,q--) *p=short_imap[*q];
  v->name=gb_save_string(buffer);
  hash_in(v); /* enter |v->name| into the hash table
                  (via utility fields |u,v|) */
}

@ Since we are generating the vertices in lexicographic order of their
inversions, it is easy to identify all adjacent vertices that
precede the current setting of $(x_1,\ldots,x_n)$. We locate them
via their symbolic names.

@<Create arcs or edges from previous permutations to~|v|@>=
for (j=1;j<n;j++)
  if (xtab[j]>xtab[j+1]) {@+register Vertex *u;
                         /* previous vertex adjacent to |v| */
    buffer[j-1]=short_imap[xtab[j+1]];@+buffer[j]=short_imap[xtab[j]];
    u=hash_out(buffer);
    if (u==NULL) panic(impossible+2); /* can't happen */
    if (directed) gb_new_arc(u,v,1L);
    else gb_new_edge(u,v,1L);
    buffer[j-1]=short_imap[xtab[j]];@+buffer[j]=short_imap[xtab[j+1]];
  }

@* Partition graphs. The subroutine call
|parts(n,max_parts,max_size,directed)|
creates a graph whose vertices represent the different ways to partition
the integer~|n| into at most |max_parts| parts, where each part is at most
|max_size|. Two partitions are adjacent in the graph if
one can be obtained from the other by combining two parts.
Each arc has length~1.

For example, the partitions of~5 are
$$5,\quad 4+1,\quad 3+2,\quad 3+1+1,\quad 2+2+1,
       \quad 2+1+1+1,\quad 1+1+1+1+1.$$
Here 5 is adjacent to $4+1$ and to $3+2$; $4+1$ is adjacent also to
$3+1+1$ and to $2+2+1$; $3+2$ is adjacent also to $3+1+1$ and to $2+2+1$; etc.
If |max_size| is 3, the partitions 5 and $4+1$ would not be included in
the graph. If |max_parts| is 3, the partitions $2+1+1+1$ and $1+1+1+1+1$
would not be included.

If |max_parts| or |max_size| are zero, they are reset to be equal to~|n|
so that they make no restriction on the partitions.

If |directed| is nonzero, the graph will contain only directed arcs from
partitions to their neighbors that have exactly one more part.

The special case when we want to generate all $p(n)$ partitions of the
integer~$n$ can be obtained by calling |all_parts(n,directed)|.

@(gb_basic.h@>=
#define all_parts(n,directed) @[parts((long)(n),0L,0L,(long)(directed))@]

@ The program for |parts| is very similar in structure to the program
for |perms| already considered.

@<Basic subroutines@>=
Graph *parts(n,max_parts,max_size,directed)
  unsigned long n; /* the number being partitioned */
  unsigned long max_parts; /* maximum number of parts */
  unsigned long max_size; /* maximum size of each part */ 
  long directed; /* should the graph be directed? */
{@+@<Vanilla local variables@>@;@#
  if (max_parts==0 || max_parts>n) max_parts=n;
  if (max_size==0 || max_size>n) max_size=n;
  if (max_parts>MAX_D) panic(bad_specs); /* too many parts allowed */
  @<Create a graph with one vertex for each partition@>;
  @<Name the partitions and create the arcs or edges@>;
  if (gb_trouble_code) {
    gb_recycle(new_graph);
    panic(alloc_fault);
     /* doggone it, we ran out of memory somewhere back there */
  }
  return new_graph;
}

@ The number of vertices is the coefficient of $z^n$
in the $z$-binomial coefficient ${m+p\choose m}_z$, where $m=|max_parts|$
and $p=|max_size|$. This coefficient is calculated as in the |perms| routine.

@<Create a graph with one vertex for each partition@>=
{@+long nverts; /* the number of vertices */
  register long *coef=gb_typed_alloc(n+1,long,working_storage);
  if (gb_trouble_code) panic(no_room+1); /* can't allocate |coef| array */
  coef[0]=1;
  for (k=1;k<=max_parts;k++) {
    for (j=n,i=n-k-max_size;i>=0;i--,j--) coef[j]-=coef[i];
    for (j=k,i=0;j<=n;i++,j++) {
      coef[j]+=coef[i];
      if (coef[j]>1000000000) panic(very_bad_specs); /* way too big */
    }
  }
  nverts=coef[n];
  gb_free(working_storage); /* recycle the |coef| array */
  new_graph=gb_new_graph(nverts);
  if (new_graph==NULL)
    panic(no_room); /* out of memory before we're even started */
  sprintf(new_graph->id,"parts(%lu,%lu,%lu,%d)",
    n,max_parts,max_size,directed?1:0);
  strcpy(new_graph->util_types,"VVZZZZZZZZZZZZ"); /* hash table will be used */
}

@ As we generate the partitions, we maintain
the numbers $\sigma_j=n-(x_1+\cdots+x_{j-1})=x_j+x_{j+1}+\cdots\,$,
somewhat as we did in the |simplex| routine. We set $x_0=|max_size|$
and $y_j=|max_parts|+1-j$; then when values $(x_1,\ldots,x_{j-1})$ are
given, the conditions
$$\sigma_j/y_j\le x_j\le \sigma_j,\qquad x_j\le x_{j-1}$$
characterize the legal values of~$x_j$.

@<Name the partitions and create the arcs or edges@>=
v=new_graph->vertices;
xx[0]=max_size;@+sig[1]=n;
for (k=max_parts,s=1;k>0;k--,s++) yy[k]=s;
if (max_size*max_parts>=n) {
  k=1;@+xx[1]=(n-1)/max_parts+1; /* $\lceil n/|max_parts|\rceil$ */
  while (1) {
    @<Complete the partial solution $(x_1,\ldots,x_k)$@>;
    @<Assign the name $x_1+\cdots+x_d$ to vertex~|v|@>;
    @<Create arcs or edges from |v| to previous partitions@>;
    v++;
    @<Advance to the next partial solution $(x_1,\ldots,x_k)$, where |k| is
       as large as possible; |goto last| if there are no more solutions@>;
  }
}
last:@+if (v!=new_graph->vertices+new_graph->n)
  panic(impossible); /* can't happen */

@ @<Complete the partial solution $(x_1,\ldots,x_k)$@>=
for (s=sig[k]-xx[k],k++;s;k++) {
  sig[k]=s;
  xx[k]=(s-1)/yy[k]+1;
  s-=xx[k];
}
d=k-1; /* the smallest part is $x_d$ */

@ Here we seek the largest $k$ such that $x_k$ can be increased without
violating the necessary and sufficient conditions stated earlier.

@<Advance to the next partial solution $(x_1,\ldots,x_k)$...@>=
if (d==1) goto last;
for (k=d-1;;k--) {
  if (xx[k]<sig[k] && xx[k]<xx[k-1]) break;
  if (k==1) goto last;
}
xx[k]++;

@ @<Assign the name $x_1+...@>=
{@+register char *p=buffer; /* string pointer */
  for (k=1;k<=d;k++) {
    sprintf(p,"+%ld",xx[k]);
    while (*p) p++;
  }
  v->name=gb_save_string(&buffer[1]); /* omit |buffer[0]|, which is |'+'| */
  hash_in(v); /* enter |v->name| into the hash table
                  (via utility fields |u,v|) */
}

@ Since we are generating the partitions in lexicographic order of their
parts, it is reasonably easy to identify all adjacent vertices that
precede the current setting of $(x_1,\ldots,x_d)$, by splitting
$x_j$ into two parts when $x_j\ne x_{j+1}$. We locate previous partitions
via their symbolic names.

@<Create arcs or edges from |v| to previous partitions@>=
if (d<max_parts) {
  xx[d+1]=0;
  for (j=1;j<=d;j++) {
    if (xx[j]!=xx[j+1]) {@+long a,b;
      for (b=xx[j]/2,a=xx[j]-b;b;a++,b--)
        @<Generate a subpartition $(n_1,\ldots,n_{d+1})$ by
              splitting $x_j$ into $a+b$, and make that subpartition
              adjacent to~|v|@>;
    }
    nn[j]=xx[j];
  }
}

@ The values of $(x_1,\ldots,x_{j-1})$ have already been copied into
$(n_1,\ldots,n_{j-1})$. Our job is to copy the smaller parts
$(x_{j+1},\ldots,x_d)$ while
inserting $a$ and $b$ in their proper places, knowing that $a\ge b$.

@<Generate a subpartition $(n_1,\ldots,n_{d+1})$...@>=
{@+register Vertex *u; /* previous vertex adjacent to |v| */
  register char *p=buffer;
  for (k=j+1;xx[k]>a;k++) nn[k-1]=xx[k];
  nn[k-1]=a;
  for (;xx[k]>b;k++) nn[k]=xx[k];
  nn[k]=b;
  for (;k<=d;k++) nn[k+1]=xx[k];
  for (k=1;k<=d+1;k++) {
    sprintf(p,"+%ld",nn[k]);
    while (*p) p++;
  }
  u=hash_out(&buffer[1]);
  if (u==NULL) panic(impossible+2); /* can't happen */
  if (directed) gb_new_arc(v,u,1L);
  else gb_new_edge(v,u,1L);
}

@* Binary tree graphs. The subroutine call
|binary(n,max_height,directed)|
creates a graph whose vertices represent the binary trees with $n$ internal
nodes and with all leaves at a distance that is at most |max_height|
from the root. Two binary trees are adjacent in the graph if
one can be obtained from the other by a single application of the
associative law for binary operations, i.e., by replacing some subtree
of the form $(\alpha\cdot\beta)\cdot\gamma$ by the subtree $\alpha\cdot
(\beta\cdot\gamma)$. (This transformation on binary trees is often
called a ``rotation.'') If the |directed| parameter is nonzero, the
directed arcs go from a tree containing $(\alpha\cdot\beta)\cdot\gamma$
to a tree containing $\alpha\cdot(\beta\cdot\gamma)$ in its place; otherwise
the graph is undirected. Each arc has length~1.

For example, the binary trees with three internal nodes form a circuit of
length~5. They are
$$\mathcode`.="2201 % \cdot
(a.b).(c.d),\quad a.(b.(c.d)),\quad a.((b.c).d),\quad (a.(b.c)).d,\quad
((a.b).c).d,$$