Moreover, the optimum flow is unique; it uses all red edges to
    their full capacity, and no green edges.  
    
      PROOF: Note that an optimum flow must use all the "W" sink edges
      and all the "W" source edges.

      Now, let "e" be any long green edge. There is at least one
      directed cut "U:V" of "G - {s,t}" that contains "e" but no red
      edges.  We can choose this cut to be a `threshold cut' i.e.
      such that any vertex of "U" is smaller than any of "V".  Let "k"
      be the number of red paths contained in "G[U + {s,t}]".  (This
      number is well-defined because no red path crosses the cut.)
      There are exactly "k" source edges of "G" entering "U", and "k"
      sink edges leaving it, which must be used by the flow.  All
      other edges in the frontier of "U" are green and directed
      inwards.  Therefore, none of these green edges may carry any
      flow.

      We have thus proved that any optimum flow must be confined to
      the red and short green edges.  Let "H" be the subgraph 
      with those edges.  For each short green edge "e",
      there is at least one threshold cut "U:V" of "H" that includes
      "e". This cut must include exactly one edge from each red path,
      directed from "U" to "V"; all other edges in the cut are
      directed from "V" to "U".  The "W" units of flow must use all
      these "W" red edges, and cannot use any of the green.  Thus, the
      flow is confined to the red edges exclusively. The claim then
      follows. QED.


PROCEDURE MakeCoins(seed: CARDINAL): Random.T =
  (*
    Creates a random number generator, primed with the given "seed".

    This routine initializes the random number generator so that each
    seed generates a different but fixed sequence.  Unfortunately, the
    random-number generator in release 3.5 does not provide a method
    to explicitly set the seed.  So we merely discarding the first ""
    elements of the Random.Default's fixed sequence.
  *)
  BEGIN
    WITH
      rnd = NEW(Random.Default).init(fixed := TRUE)
    DO
      Wr.PutText(stderr, "warming up the random number generator...\n");
      FOR i := 1 TO 200*seed DO
        EVAL rnd.integer()
      END;
      RETURN rnd
    END
  END MakeCoins;
  

  PROCEDURE PickGreenEdge(pLo, pHi: CARDINAL): ARRAY [0..1] OF Vertex =
    (* 
      Selects a pair of vertices "u,v" suitable for a green edge.
    *)
    VAR u, v: INTEGER;
    BEGIN
      WITH
        NV = vHi - vLo + 1,
        p = coins.integer(pLo, pHi)
      DO
        <* ASSERT 2*pHi < NV*(NV-1) *>
        <* ASSERT pLo <= pHi *>
        
        (* Map pair index to vertex pair: *)
        WITH 
          NR = NV DIV 2,
          NC = NV - (1 - NV MOD 2)
        DO
          u := p DIV NC;
          v := p MOD NC;
          IF u > v THEN
            u := NV - 1 - u;
            v := NC - v;
            <* ASSERT u < v *>
          END;
        END;
        
        (* Fiddle 
      END
                    
    END PickGreenEdge
  
  
PROCEDURE MakeNestGraph(
    NV: Count;           (* Total vertex count. *)
    M: Count;            (* Maximum flow. *)
    L: Count;            (* Mean flow length. *)
    NVX: Count;          (* Number of vertices in skeleton graph. *)
    NEX: Count;          (* Number of edges in skeleton graph. *)
    MX: Count;           (* Max flow for skeleton graph; must divide "M". *)
    LX: Count;           (* Mean flow length for skeleton graph; must divide "M". *)
    rnd: Random.T;       (* Randomness source. *)
    VAR e: Edges;        (* OUT: edge vector. *)
    VAR s: VertexNum;       (* OUT: source vertex. *)
    VAR t: VertexNum;       (* OUT: sink vertex. *)
  ) =
  (*
    The "Nest" graph is essentially a "Rete" graph, the `skeleton',
    whose red edges have been replaced by copies of another
    other "Rete" graph, the `atom'.
    
    The skeleton graph has max flow "MX" and path length "LX".  The
    corresponding values for the atom graph are "MY = M DIV MX" and
    "LY = L DIV LX".  (The divisions must be exact.) The final graph
    has therefore max flow "M = MY * MX" and path length "L = LY *
    LX".
    
    The number of red edges "R" will be "RY*RX = MY*LY +
    MX*LX".  The skeleton graph will have a total of "NVX" vertices
    and "NEX" edges. Thus, it will have "NVX - 2 - MX*(LX-1)" 
    green vertices (not on any red path), and "NEX - MX*LX" green edges.
    Each instance of the atom will have "NV - 
    
    Excess edges, which will be green, are divided
    equally between the atoms, subject to
    natural bounds, with the remainder assigned to the skeleton.
  *)
  CONST MaxVertices = 256*256;
  BEGIN
    <* ASSERT M > 0 *>
    <* ASSERT MX > 0 *>
    <* ASSERT L > 0 *>
    <* ASSERT LX > 0 *>
    <* ASSERT NV >= 2 *>
    <* ASSERT NUMBER(e) > 0 *>
    Require("M", M, 1, (MaxVertices - 2) DIV 2);
    IF L = 1 THEN Require("M", M, 1, 1) END;
    IF LX = 1 THEN Require("MX", MX, 1, 1) END;
    Require("M MOD MX", M MOD MX, 0, 0);
    Require("L MOD LX", L MOD LX, 0, 0);
    WITH
      NE = NUMBER(e),
      
      (* Outer (skeleton) graph: *)
      RX = MX * LX,
      minEX = RX,
      minVX = MX * (LX - 1) + 2,
      
      (* Inner (atom) graph: *)
      MY = M DIV MX,
      LY = L DIV LX,
      RY = MY * LY,
      minEY = RY,
      minVY = MY * (LY - 1) + 2,

      minV = RX * (minVY - 2) + minVX,
      maxV = MaxVertices
    DO
      <* ASSERT RX > 0 *>
      <* ASSERT RY > 0 *>
      Require("NV", NV, minV, maxV);

      WITH

        (* Distribution of vertices: *)
        NVY = (NV - minVX) DIV RX + 2,

        NVX = NV - (NVY - 2) * RX,
        
        (* Distribution of edges: *)
        maxEY = MY * MIN(LY, 2) + Choose2(NVY-2),
        NEY = MIN(MAX(NE DIV RX, minEY), maxEY),

        maxEX = MX * MIN(LX, 2) + Choose2(NVX-2),
        NEX = RX + (NE - RX * NEY),
        
        minE = RX * minEY,
        maxE = RX * maxEY + Choose2(NVX-2)
      DO
        <* ASSERT NEX >= 0 *>
        Require("NE", NE, minE, maxE);
        
        <* ASSERT NVY >= 2 *>
        
        WITH
          eX = NEW(REF Edges, NEX)^
        DO
          (* Create skeleton in "eX": *)
          s := 0;
          t := NVX-1;
          MakeReteSubGraph(
            s := s, t := t, 
            M := MX, (L-1) := (LX-1),
            vLo := 1, vHY := NVX-2,
            eLo := 0, eHY := NEX-1,
            rnd := rnd,
            e := eX
          );
          (* Create expanded graph in "e": *)
          ek := 0;
          vk := NVX;
          FOR iX := 0 TO LAST(eX) DO
            WITH uv = eX[iX], u = uv[0], v = uv[1] DO
              IF u > v THEN 
                (* Green edge, just copy: *)
                e[ek] := uv; INC(ek)
              ELSIF u < v THEN
                (* Red edge, replace by an instance of atom: *)
                MakeReteSubGraph(
                  s := u, t := v, 
                  M := MY, L := LY,
                  vLo := vk, vHY := vk + NVY-3,
                  eLo := ek, eHY := ek + NEY-1,
                  rnd := rnd,
                  e := e
                );
                ek := ek + NEY
              ELSE
                (* Self-loop *)
                <* ASSERT FALSE *>
              END
            END
          END;
          <* ASSERT ek = NE *>
        END
      END
    END
  END MakeNestGraph;
  


PROCEDURE AddSpinSubgraph(
    s, t: CARDINAL;          (* The source and sink vertices. *)
    W: CARDINAL;             (* Maximum flow. *)
    L: CARDINAL;             (* Mean flow length. *)
    vLo, vHi: VertexNum;     (* New vertices (excl. "s" and "t") are "[vlo..vHi-1]". *)
    rnd: Random.T;           (* Randomness source *)
    (*OUT*) VAR e: Edges;    (* Edge vector. *)
    (*OUT*) VAR m: Marks;    (* Set to TRUE for solution edges. *)
    (*OUT*) VAR c: Costs;    (* Edge capacities (all 1 for now). *)
  ) =
  (*
    Makes a "Spin" subgraph with maximum flow "W" and mean flow
    length "L".  The width "W" must be at least 2,
    and the length "L" must be an odd integer, "L = 2*S+1". 
    
    There are "2^S" maximum flows from "s" to "t".  Ignoring eddies,
    there are "W*(3*S+1)" edges that take essential part in any of
    these flows (the `red' edges); these edges will have their "m"
    bits set to TRUE.  The remaining edges are painted `green'.
    
    Besides the source "s" and sink "t", the red edges are incident to
    "W*(L-1)" `red' vertices, arranged as "L-1" layers of "W" vertices
    each. Let "red[i,k]" denote the "k"th red vertex in layer "i", in
    circular order.
    
    The red edges comprise:
      "W" `entry' edges connecting "s" to the vertices of layer 1;
      "W" `exit' edges connecting the vertices of layer "L-1" to the sink "t";
      "2*W*(L-2)" `switch' edges, connecting every other 
    
    The switch edges are organized into "L-2" layers of "2*W" edges
    each. The "i"th edge layer lies between vertex layers "i" and
    "i+1".  More precisely, the "k"th edge in layer "i" goes
    from "red[i, (k DIV 2) MOD W]" to "red[i+1, ((k+1) DIV 2) MOD w]", for
    "k" in "[0..2*W-1]".  That edge is said to be `even' or `odd'
    according to the parity of "k".  Note that two consecutive switch edges
    in a layer have either the same origin or the same destination. 
    
    The green vertices and green edges, if any, are divided into one or more 
    blocks.  Each block is attached to only one red vertex, and contains
    a random collection of edges, with arbitrary oriantations.
    
    The optimum flow is obviously "W". A trivial cut at "s" or "t" proves the 
    upper bound.  The flow consisting of the entry, exit, and even switch edges
    demonstrates the lower bound.
    
    Moreover, the optimum flow is unique; it uses all red edges to
    their full capacity, and no green edges.
    
      PROOF: Note that an optimum flow must use all the "W" sink edges
      and all the "W" source edges.

      Now, let "e" be any long green edge. There is at least one
      directed cut "U:V" of "G - {s,t}" that contains "e" but no red
      edges.  We can choose this cut to be a `threshold cut' i.e.
      such that any vertex of "U" is smaller than any of "V".  Let "k"
      be the number of red paths contained in "G[U + {s,t}]".  (This
      number is well-defined because no red path crosses the cut.)
      There are exactly "k" source edges of "G" entering "U", and "k"
      sink edges leaving it, which must be used by the flow.  All
      other edges in the frontier of "U" are green and directed
      inwards.  Therefore, none of these green edges may carry any
      flow.

      We have thus proved that any optimum flow must be confined to
      the red and short green edges.  Let "H" be the subgraph 
      with those edges.  For each short green edge "e",
      there is at least one threshold cut "U:V" of "H" that includes
      "e". This cut must include exactly one edge from each red path,
      directed from "U" to "V"; all other edges in the cut are
      directed from "V" to "U".  The "W" units of flow must use all
      these "W" red edges, and cannot use any of the green.  Thus, the
      flow is confined to the red edges exclusively. The claim then
      follows. QED.
    
    If "L > 1", the procedure requires "NV >= W*(L-1) + 2", 
    "NE >= W*L", and "NE <= 2*W + (NV-2:2)", where "(n:k)" denotes
    the binomial coefficient: "n!/k!/(n-k)" (or 0 if "k>n").
    
    If "L = 1", then "W" must be 1, and "NE" must be at
    least "W" and at most "W + (NV-2:2)". The red subgraph will be
    either empty or a single edge connecting "s" and "t".
  *)
  
  TYPE PathRange = ARRAY [0..1] OF VertexNum;
  
  PROCEDURE SelectPathNodes(): REF ARRAY OF PathRange =
    (*
      Returns a vector "pr" such that "pr[k][0]" is the 
      second vertex of path "k", and "pr[k][1]" is its
      next-to-last vertex.  
      
      The vertices are consecutive an increasing along each path. 
      The remaining vertices are distributed more or less uniformly
      between the paths, with a few before the first path and a few
      after the last path.
    *)
    BEGIN
      <* ASSERT W > 0 *>
      <* ASSERT L > 0 *>
      WITH
        NV = vHi - vLo,
        rpr = NEW(REF ARRAY OF PathRange, W), pr = rpr^,
        NGV = NV - W * (L - 1)
      DO
        (* Enough vertices for all red paths? *)
        <* ASSERT NV >= W * (L - 1) *>   
          
        FOR k := 0 TO W-1 DO
          WITH g = ((0 + 3*k) * NGV) DIV (0 + 3*(W-1)) DO
            pr[k][0] := vLo + g + k * (L-1);
            IF k > 0 THEN <* ASSERT pr[k][0] > pr[k-1][1] *> END;
            pr[k][1] := pr[k][0] + L - 2;
            <* ASSERT pr[k][1] < vHi *>
          END
        END;
        RETURN rpr
      END
    END SelectPathNodes;
    
  TYPE PathIndex = INTEGER;
  
  PROCEDURE AssignNodesToPaths(READONLY pr: ARRAY OF PathRange): REF ARRAY OF PathIndex =
    (* 
      Returns a vector "px" such that "px[i]" is "k" if
      vertex "vLo + i" belongs to red path "k"; or "-1" if 
      it doesn't belong to any red path.
    *)
    BEGIN
      WITH
        NV = vHi - vLo,
        rpx = NEW(REF ARRAY OF PathIndex, NV), px = rpx^
      DO
        FOR v := 0 TO NV-1 DO px[v] := -1 END;
        FOR k := 0 TO W-1 DO
          FOR v := pr[k][0] TO pr[k][1] DO
            px[v - vLo] := k
          END
        END;
        RETURN rpx
      END 
    END AssignNodesToPaths;
    
  VAR ek: CARDINAL;
  
  PROCEDURE AddRedEdges(
      READONLY pr: ARRAY OF PathRange;
      <*UNUSED*> READONLY px: ARRAY OF PathIndex;
    ) =
    (* 
      Adds the red edges.
    *)
    VAR ctXRE, ctIRE: Count := 0;
    BEGIN
      <* ASSERT L > 0 *>
      <* ASSERT W > 0 *>
      IF L = 1 THEN 
        IF W = 1 THEN 
          e[ek] := Edge{s,t}; W[ek] := TRUE;
          INC(ek); INC(ctXRE)
        ELSE
          <* ASSERT W = 0 *>
        END
      ELSE
        <* ASSERT vLo < vHi *>
        <* ASSERT eLo < eHi *>
        WITH
          NE = eHi - eLo,
          NV = vHi - vLo,
          (* Split "NE" among the red, short green, and long green edges: *)
          NRE = W * L,        (* Total number of red edges *)
          NXRE = 2 * W,       (* External red edges, attached to "s" and/or "t" *)
          NIRE = NRE - NXRE   (* Internal red edges. *) 
        DO
          (* Enough space in "e" for all red edges? *)
          <* ASSERT NE >= NRE *>    
          
          (* Enough vertices for the red edges? *)
          <* ASSERT NV >= W * (L - 1) *>

          FOR k := 0 TO W-1 DO
            WITH u = pr[k][0], v = pr[k][1] DO
              e[ek] := Edge{s, u}; W[ek] := TRUE;
              INC(ek); INC(ctXRE);
              FOR t := u TO v - 1 DO
                e[ek] := Edge{t, t+1}; m[ek] := TRUE;
                INC(ek); INC(ctIRE)
              END;
              e[ek] := Edge{v, t}; m[ek] := TRUE;
              INC(ek); INC(ctXRE);
            END
          END;
          <* ASSERT ctXRE = NXRE *>
          <* ASSERT ctIRE = NIRE *>
        END;
      END;
    END AddRedEdges;

  PROCEDURE AddGreenEdges(
      READONLY pr: ARRAY OF PathRange;
      READONLY px: ARRAY OF PathIndex;
    ) =
    (* 
      Adds the green edges. Biased towards `perverse shortcuts'.
    *)
    VAR ctGE: Count := 0;   (* Number of green edges already added. *)
        remGE: Count;       (* Number of remaining chances to pick a green edge. *)
    BEGIN
      <* ASSERT L > 0 *>
      <* ASSERT W > 0 *>
      (* Enumerate the decreasing pairs in order of decreasing desirability.*)
      (* Select each pair with probability "1 - 0.25*(ctGE/NGE)". *)
      WITH
        NE = eHi - eLo,
        NV = vHi - vLo,
        (* Split "NE" among the red, short green, and long green edges: *)
        NRE = W * L,        (* Required number of red edges *)
        NGE = NE - NRE,     (* Required number of green edges. *)
        MaxGE = Choose2(NV) (* Max possible number of green edges. *)
      DO
        (* Enough vertices for all green edges? *)
        <* ASSERT NGE <= MaxGE *>
        
        remGE := MaxGE;
        (* First, all long RR edges:  *)
        FOR t := 1 TO W-1 DO
          FOR d := 0 TO 2*(L-1)-1 DO
            (* Enumerate long RR edges that create an "s-t" path of length "d+3": *)
            FOR r := 0 TO d DO
              FOR ku := 0 TO W-1-t DO
                WITH kv = ku + t DO
                  IF ctGE >= NGE THEN RETURN END;
                  WITH 
                    v = pr[kv][0] + r,
                    u = pr[ku][1] - (d - r)
                  DO
                    IF v <= pr[kv][1] AND u >= pr[ku][0] THEN
                      (* Found another long RR edge: *)
                      <* ASSERT u < v *>
                      WITH,
                        prob = 0.75d0 + 0.25d0 * FLOAT(NGE-ctGE, LONG)/FLOAT(remGE, LONG)
                      DO
                        IF rnd.longreal(0.0d0, 1.0d0) <= prob THEN 
                          e[ek] := Edge{v, u}; m[ek] := FALSE;
                          INC(ek); INC(ctGE);
                        END
                      END;
                      DEC(remGE)
                    END
                  END
                END
              END
            END
          END;
        END;

        (* Enumerate other decreasing edges, except the ones above:  *)
        FOR t := 1 TO NV - 1 DO
          FOR u := vLo TO vHi - 1 - t DO
            IF ctGE >= NGE THEN RETURN END;
            WITH v = u + t DO
              <* ASSERT u < v *>
              <* ASSERT v < vHi *>
              IF px[u - vLo] = -1 
              OR px[v - vLo] = -1 
              OR px[u - vLo] = px[v - vLo] THEN
                (* Found another non-long or non-RR edge *)
                <* ASSERT u < v *>
                WITH
                  prob = 0.25d0 + 0.75d0 * FLOAT(NGE-ctGE, LONG)/FLOAT(remGE, LONG)
                DO
                  IF rnd.longreal(0.0d0, 1.0d0) <= prob THEN 
                    e[ek] := Edge{v, u}; INC(ek); INC(ctGE);
                  END
                END;
                DEC(remGE)
              END
            END
          END
        END;
        <* ASSERT ctGE = NGE *>
        <* ASSERT remGE = 0 *>
      END     
    END AddGreenEdges;

  BEGIN (* AddSpinSubGraph *)
    <* ASSERT W > 0 *>
    <* ASSERT L > 0 *>
    <* ASSERT eLo < eHi *>
    IF L = 1 THEN Require("W", W, 1, 1) END;
    WITH
      NV = vHi - vLo,
      NE = NUMBER(e),
      minV = W * (L - 1) + 2,
      maxV = MaxVertices,
      minE = W * L,
      maxE = MIN(L, 2) * W + Choose2(NV-2)
    DO
      Require("NV", NV, minV, maxV);
      Require("NE", NE, minE, maxE);
      WITH
        pr = SelectPathNodes()^,
        px = AssignNodesToPaths(pr)^
      DO
        ek := 0;
        AddRedEdges(pr, px);
        AddGreenEdges(pr, px);
        <* ASSERT ek = NE-1 *>
      END
    END;
  END AddSpinSubGraph;