gb_rand.w

模拟器提供了一个简单易用的平台

  if (dist_to) {
    while (nn<n) nn+=nn, kk--;
    to_table=walker(n,nn,dist_to,new_graph);
  }
  if (gb_trouble_code) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* oops, we ran out of memory somewhere back there */
  }
}

@ @<Private...@>=
typedef struct {
  long prob; /* a probability, multiplied by $2^{31}$ and translated */
  long inx; /* index that might be selected */
} @[magic_entry@];

@ Once the magic tables have been set up, we can generate
nonuniform vertices by using the following code:

@<Generate a random vertex |u|...@>=
{@+register magic_entry *magic;
  register long uu=gb_next_rand(); /* uniform random number */
  k=uu>>kk;
  magic=from_table+k;
  if (uu<=magic->prob) u=new_graph->vertices+k;
  else u=new_graph->vertices+magic->inx;
}

@ @<Generate a random vertex |v|...@>=
{@+register magic_entry *magic;
  register long uu=gb_next_rand(); /* uniform random number */
  k=uu>>kk;
  magic=to_table+k;
  if (uu<=magic->prob) v=new_graph->vertices+k;
  else v=new_graph->vertices+magic->inx;
}

@ So all we have to do is set up those magic tables. If |uu| is a uniform
random integer between 0 and $2^{31}-1$, the index |k=uu>>kk| is a
uniform random integer between 0
and |nn-1|, because of the relation between |nn| and |kk|. Once |k| is
computed, the code above selects vertex~|k| with probability
|(p+1-(k<<kk))|/$2^{31}$, where |p=magic->prob| and |magic| is the $k$th
element of the magic table; otherwise the code selects
vertex |magic->inx|. The trick is to set things up so that each vertex
is selected with the proper overall probability.

Let's imagine that the given distribution vector has length |nn|,
instead of~|n|, by extending it if necessary with zeroes. Then the
average entry among these |nn| integers is exactly $t=2^{30}/|nn|$.
If some entry, say entry~|i|, exceeds |t|, there must be another entry
that's less than |t|, say entry~|j|. We can set the $j$th entry
of the magic table so that its |prob| field selects vertex~$j$ with the
correct probability, and so that its |inx| field equals~|i|. Then
we are selecting vertex~|i| with a certain residual probability; so we
subtract that residual from |i|'s present probability, and repeat the
process with vertex~|j| eliminated. The average of the remaining entries
is still~|t|, so we can repeat this procedure until all remaining entries
are exactly equal to~|t|. The rest is easy.

During the calculation, we maintain two linked lists of
|(prob,inx)| pairs. The |hi| list contains entries with |prob>t|,
and the |lo| list contains the rest. During this part of the computation
we call these list elements `nodes', and we use the field names
|key| and~|j| instead of |prob| and |inx|.

@<Private...@>=
typedef struct node_struct {
  long key; /* a numeric quantity */
  struct node_struct *link; /* the next node on the list */
  long j; /* a vertex number to be selected with probability $|key|/2^{30}$ */
} node;
static Area temp_nodes; /* nodes will be allocated in this area */
static node *base_node; /* beginning of a block of nodes */

@ @<Internal...@>=
static magic_entry *walker(n,nn,dist,g)
  long n; /* length of |dist| vector */
  long nn; /* $2^{\lceil\mskip1mu\lg n\rceil}$ */
  register long *dist;
    /* start of distribution table, which sums to $2^{30}$ */
  Graph *g; /* tables will be allocated for this graph's vertices */
{@+magic_entry *table; /* this will be the magic table we compute */
  long t; /* average |key| value */
  node *hi=NULL, *lo=NULL; /* nodes not yet included in magic table */
  register node *p, *q; /* pointer variables for list manipulation */
  base_node=gb_typed_alloc(nn,node,temp_nodes);
  table=gb_typed_alloc(nn,magic_entry,g->aux_data);
  if (!gb_trouble_code) {
    @<Initialize the |hi| and |lo| lists@>;
    while (hi) @<Remove a |lo| element and match it with a |hi| element;
        deduct the residual probability from that |hi|~element@>;
    while (lo) @<Remove a |lo| element of |key| value |t|@>;
  }
  gb_free(temp_nodes);
  return table; /* if |gb_trouble_code| is nonzero, the table is empty */
}

@ @<Initialize the |hi| and |lo| lists@>=
t=0x40000000/nn; /* this division is exact */
p=base_node;
while (nn>n) {
  p->key=0;
  p->link=lo;
  p->j=--nn;
  lo=p++;
}
for (dist=dist+n-1; n>0; dist--,p++) {
  p->key=*dist;
  p->j=--n;
  if (*dist>t)
    p->link=hi,@, hi=p;
  else p->link=lo,@, lo=p;
}

@ When we change the scale factor from $2^{30}$ to $2^{31}$, we need to
be careful lest integer overflow occur. The introduction of register |x| into
this code removes the risk.

@<Remove a |lo| element and match it with a |hi| element...@>=
{@+register magic_entry *r; register long x;
  p=hi,@, hi=p->link;
  q=lo,@, lo=q->link;
  r=table+q->j;
  x=t*q->j+q->key-1;
  r->prob=x+x+1;
  r->inx=p->j;
  /* we have just given |q->key| units of probability to vertex |q->j|,
     and |t-q->key| units to vertex |p->j| */
  if ((p->key-=t-q->key)>t)
    p->link=hi,@, hi=p;
  else p->link=lo,@, lo=p;
}

@ When all remaining entries have the average probability, the
|inx| component need not be set, because it will never be used.

@<Remove a |lo| element of |key| value |t|@>=
{@+register magic_entry *r; register long x;
  q=lo, lo=q->link;
  r=table+q->j;
  x=t*q->j+t-1;
  r->prob=x+x+1;
  /* that's |t| units of probability for vertex |q->j| */
}

@*Random bipartite graphs. The procedure call
$$\hbox{|random_bigraph(n1,n2,m,multi,dist1,dist2,min_len,max_len,seed)|}$$
is designed to produce a pseudo-random bipartite graph
with |n1| vertices in one part and |n2| in the other, having |m| edges.
The remaining parameters |multi|, |dist1|, |dist2|, |min_len|, |max_len|,
and |seed| have the same meaning as the analogous parameters of |random_graph|.

In fact, |random_bigraph| does its work by reducing its parameters
to a special case of |random_graph|. Almost all that needs to be done is
to pad |dist1| with |n2| trailing zeroes and |dist2| with |n1| leading
zeroes. The only slightly tricky part occurs when |dist1| and/or |dist2| are
null, since non-null distribution vectors summing exactly to $2^{30}$ must then
be fabricated.

@<External f...@>=
Graph *random_bigraph(n1,n2,m,multi,dist1,dist2,min_len,max_len,seed)
  unsigned long n1,n2; /* number of vertices desired in each part */
  unsigned long m; /* number of edges desired */
  long multi; /* allow duplicate edges? */
  long *dist1, *dist2; /* distribution of edge endpoints */
  long min_len,max_len; /* bounds on random lengths */
  long seed; /* random number seed */
{@+unsigned long n=n1+n2; /* total number of vertices */
  Area new_dists;
  long *dist_from, *dist_to;
  Graph *new_graph;
  init_area(new_dists);
  if (n1==0 || n2==0) panic(bad_specs); /* illegal options */
  if (min_len>max_len) panic(very_bad_specs); /* what are you trying to do? */
  if (((unsigned long)(max_len))-((unsigned long)(min_len))>=
      ((unsigned long)0x80000000)) panic(bad_specs+1); /* too much range */
  dist_from=gb_typed_alloc(n,long,new_dists);
  dist_to=gb_typed_alloc(n,long,new_dists);
  if (gb_trouble_code) {
    gb_free(new_dists);
    panic(no_room+2); /* no room for auxiliary distribution tables */
  }
  @<Compute the entries of |dist_from| and |dist_to|@>;
  new_graph=random_graph(n,m,multi,0L,0L,
                dist_from,dist_to,min_len,max_len,seed);
  sprintf(new_graph->id,"random_bigraph(%lu,%lu,%lu,%d,%s,%s,%ld,%ld,%ld)",@|
    n1,n2,m,multi>0?1:multi<0?-1:0,dist_code(dist1),dist_code(dist2),@|
    min_len,max_len,seed);
  mark_bipartite(new_graph,n1);
  gb_free(new_dists);
  return new_graph;
}

@ The relevant identity we need here is the replicative law for the
floor function:
$$\left\lfloor x\over n\right\rfloor+\left\lfloor x+1\over n\right\rfloor
+ \cdots + \left\lfloor x+n-1\over n\right\rfloor = \lfloor x\rfloor\,.$$

@<Compute the entries...@>=
{@+register long *p, *q; /* traversers of the dists */
  register long k; /* vertex count */
  p=dist1; q=dist_from;
  if (p)
    while (p<dist1+n1) *q++=*p++;
  else for (k=0; k<n1; k++) *q++=(0x40000000+k)/n1;
  p=dist2; q=dist_to+n1;
  if (p)
    while (p<dist2+n2) *q++=*p++;
  else for (k=0; k<n2; k++) *q++=(0x40000000+k)/n2;
}

@* Random lengths. The subroutine call
$$\hbox{|random_lengths(g,directed,min_len,max_len,dist,seed)|}$$
takes an existing graph and assigns new lengths to
each of its arcs. If |dist=NULL|, the lengths will be uniformly distributed
between |min_len| and |max_len| inclusive; otherwise |dist|
should be a probability distribution vector of length |max_len-min_len+1|,
like those in |random_graph|.

If |directed=0|, pairs of arcs $u\to v$ and $v\to u$ will be regarded as
a single edge, both arcs receiving the same length.

The procedure returns a nonzero value if something goes wrong; in that
case, graph |g| will not have been changed.

Alias tables for generating nonuniform random lengths will survive
in |g->aux_data|.

@<External f...@>=
long random_lengths(g,directed,min_len,max_len,dist,seed)
  Graph *g; /* graph whose lengths will be randomized */
  long directed; /* is it directed? */
  long min_len,max_len; /* bounds on random lengths */
  long *dist; /* distribution of lengths */
  long seed; /* random number seed */
{@+register Vertex *u,*v; /* current vertices of interest */
  register Arc *a; /* current arc of interest */
  long nn=1, kk=31; /* variables for nonuniform generation */
  magic_entry *dist_table; /* alias table for nonuniform generation */
  if (g==NULL) return missing_operand; /* where is |g|? */
  gb_init_rand(seed);
  if (min_len>max_len) return very_bad_specs; /* what are you trying to do? */
  if (((unsigned long)(max_len))-((unsigned long)(min_len))>=
      ((unsigned long)0x80000000)) return bad_specs; /* too much range */
  @<Check |dist| for validity, and set up the |dist_table|@>;
  sprintf(buffer,",%d,%ld,%ld,%s,%ld)",directed?1:0,@|
     min_len,max_len,dist_code(dist),seed);
  make_compound_id(g,"random_lengths(",g,buffer);
  @<Run through all arcs and assign new lengths@>;
  return 0;
}

@ @<Private dec...@>=
static char buffer[]="1,-1000000001,-1000000000,dist,1000000000)";

@ @<Check |dist| for validity...@>=
if (dist) {@+register long acc; /* sum of probabilities */
  register long *p; /* pointer to current probability of interest */
  register long n=max_len-min_len+1;
  for (acc=0,p=dist; p<dist+n; p++) {
    if (*p<0) return -1; /* negative probability */
    if (*p>0x40000000-acc) return 1; /* probability too high */
    acc+=*p;
  }
  if (acc!=0x40000000) return 2; /* probabilities don't sum to 1 */
  while (nn<n) nn+=nn,kk--;
  dist_table=walker(n,nn,dist,g);
  if (gb_trouble_code) {
    gb_trouble_code=0;
    return alloc_fault; /* not enough room to generate the magic tables */
  }
}

@ @<Run through all arcs and assign new lengths@>=
for (u=g->vertices;u<g->vertices+g->n;u++)
  for (a=u->arcs;a;a=a->next) {
    v=a->tip;
    if (directed==0 && u>v) a->len=(a-1)->len;
    else {@+register long len; /* a random length */
      if (dist==0) len=rand_len;
      else {@+long uu=gb_next_rand();
        long k=uu>>kk;
        magic_entry *magic=dist_table+k;
        if (uu<=magic->prob) len=min_len+k;
        else len=min_len+magic->inx;
      }
      a->len=len;
      if (directed==0 && u==v && a->next==a+1) (++a)->len=len;
    }
  }

@* Index. Here is a list that shows where the identifiers of this program are
defined and used.