gb_econ.w
上传用户:rrhhcc
上传日期:2015-12-11
资源大小:54129k
文件大小:27k
源码类别:

通讯编程

开发平台:

Visual C++

  1. % This file is part of the Stanford GraphBase (c) Stanford University 1993
  2. @i boilerplate.w %<< legal stuff: PLEASE READ IT BEFORE MAKING ANY CHANGES!
  3. @i gb_types.w
  4. deftitle{GB_,ECON}
  5. prerequisites{GB_,GRAPH}{GB_,IO}
  6. @* Introduction. This GraphBase module contains the |econ| subroutine,
  7. which creates a family of directed graphs related to the flow of money
  8. between industries.  An example of the use of this procedure can be
  9. found in the demo program {sc ECON_,ORDER}.
  10. @(gb_econ.h@>=
  11. extern Graph *econ();
  12. @ The subroutine call |econ(n,omit,threshold,seed)|
  13. constructs a directed graph based on the information in .{econ.dat}.
  14. Each vertex of the graph corresponds to one of 81 sectors of the U.S.
  15. economy. The data values come from the year 1985; they were derived from
  16. tables published in {sl Survey of Current Business/ bf70} (1990), 41--56.
  17. If |omit=threshold=0|, the directed graph is a ``circulation'';
  18. that is, each arc has an associated |flow| value, and
  19. the sum of arc flows leaving each vertex is equal to the
  20. sum of arc flows entering. This sum is called the ``total commodity output''
  21. for the sector in question. The flow in an arc from sector $j$~to
  22. sector~$k$ is the amount of the commodity made by sector~$j$ that was
  23. used by sector~$k$, rounded to millions of dollars at producers' prices.
  24. For example, the total commodity output of the sector called .{Apparel}
  25. is 54031, meaning that the total cost of making all kinds of apparel in
  26. 1985 was about 54 billion dollars. There is an arc from .{Apparel} to
  27. itself with a flow of 9259, meaning that 9.259 billion dollars' worth
  28. of apparel went from one group within the apparel industry to another.
  29. There also is an arc of flow~44 from .{Apparel} to .{Household}
  30. .{furniture}, indicating that some 44 million dollars' worth of apparel
  31. went into the making of household furniture. By looking at all
  32. arcs that leave the .{Apparel} vertex, you can see where all that
  33. new apparel went; by looking at all arcs that enter .{Apparel}, you can
  34. see what ingredients the apparel industry needed to make~it.
  35. One vertex, called .{Users}, represents people like you and me, the
  36. non-industrial end users of everything. The arc from .{Apparel} to
  37. .{Users} has flow 42172; this is the ``total final demand'' for
  38. apparel, the amount that didn't flow into other sectors of the economy
  39. before it reached people like us. The arc from .{Users} to .{Apparel}
  40. has flow 19409, which is called the ``value added'' by users; it
  41. represents wages and salaries paid to support the manufacturing
  42. process. The sum of total final demand over all sectors, which also
  43. equals the sum of value added over all sectors, is conventionally
  44. called the Gross National Product (GNP). In 1985 the GNP was 3999362,
  45. nearly 4 trillion dollars, according to .{econ.dat}. (The sum of all
  46. arc flows coming out of all vertices was 7198680; this sum
  47. overestimates the total economic activity, because it counts some
  48. items more than once---statistics are recorded whenever an item
  49. passes a statistics gatherer. Economists try to adjust the data so that
  50. they avoid double-counting as much as possible.)
  51. Speaking of economists, there is another special vertex called
  52. .{Adjustments}, included by economists so that GNP is measured
  53. more accurately. This vertex takes account of such things as changes in
  54. the value of inventories, and imported materials that cannot be obtained
  55. within the U.S., as well as work done for the government and for foreign
  56. concerns. In 1985, these adjustments accounted for about 11% of the GNP.
  57. Incidentally, some of the ``total final demand'' arcs are negative.
  58. For example, the arc from .{Petroleum} .{and} .{natural} .{gas}
  59. .{production} to .{Users} has flow $-27032$. This might seem strange
  60. at first, but it makes sense when imports are considered, because
  61. crude oil and natural gas go more to other industries than to end users.
  62. Total final demand does not mean total user demand.
  63. @d flow a.I /* utility field |a| specifies the flow in an arc */
  64. @ If |omit=1|, the .{Users} vertex is omitted from the digraph; in
  65. particular, this will eliminate all arcs of negative flow. If
  66. |omit=2|, the .{Adjustments} vertex is also omitted, thereby leaving
  67. 79~sectors with arcs showing inter-industry flow. (The graph is no
  68. longer a ``circulation,'' of course, when |omit>0|.)  If .{Users} and
  69. .{Adjustments} are not omitted, .{Users} is the last vertex of the
  70. graph, and .{Adjustments} is next-to-last.
  71. If |threshold=0|, the digraph has an arc for every nonzero |flow|.
  72. But if |threshold>0|, the digraph becomes more sparse;
  73. there is then an arc from $j$ to~$k$ if and
  74. only if the amount of commodity $j$ used by sector~$k$ exceeds
  75. |threshold/65536| times the total input of sector~$k$.  (The total
  76. input figure always includes value added, even if |omit>0|.)
  77. Thus the arcs go to each sector from
  78. that sector's main suppliers. When |n=79|, |omit=2|, and
  79. |threshold=0|, the digraph has 4602 arcs out of a possible
  80. $79times79=6241$. Raising |threshold| to 1 decreases the number of
  81. arcs to 4473; raising it to 6000 leaves only~72 arcs.
  82. The |len| field in each arc is~1.
  83. The constructed graph will have $min(n,81-|omit|)$ vertices. If |n| is less
  84. than |81-omit|, the |n| vertices will be selected by repeatedly combining
  85. related sectors. For example, two of the 81 original sectors are called
  86. `.{Paper} .{products,} .{except} .{containers}' and
  87. `.{Paperboard} .{containers} .{and} .{boxes}'; these might be combined
  88. into a sector called `.{Paper} .{products}'. There is a binary tree
  89. with 79 leaves, which describes a fixed hierarchical breakdown of the
  90. 79 non-special sectors. This tree is
  91. pruned, if necessary, by replacing pairs of leaves by their parent node,
  92. which becomes a new leaf; pruning continues
  93. until just |n| leaves remain. Although pruning is a bottom-up process, its
  94. effect can also be obtained from the top down if we imagine ``growing''
  95. the tree, starting out with a whole economy as a single sector and
  96. repeatedly subdividing a sector into two parts. For example,
  97. if |omit=2| and |n=2|, the two sectors will
  98. be called .{Goods} and .{Services}. If |n=3|, .{Goods} might be
  99. subdivided into .{Natural} .{Resources} and .{Manufacturing}; or
  100. .{Services} might be subdivided into .{Indirect} .{Services} and
  101. .{Direct} .{Services}.
  102. If |seed=0|, the binary tree is pruned in such a way that the |n|
  103. resulting sectors are as equal as possible with respect to total
  104. input and output, while respecting the tree structure. If |seed>0|,
  105. the pruning is carried out at random, in such a way that all |n|-leaf
  106. subtrees of the original tree are obtained with approximately equal
  107. probability (depending on |seed| in a machine-independent fashion).
  108. Any |seed| value from 1 to $2^{31}-1=2147483647$ is permissible.
  109. As usual in GraphBase routines, you can set |n=0| to get the default
  110. situation where |n| has its maximum value. For example, either
  111. |econ(0,0,0,0)| or |econ(81,0,0,0)| produces the full graph;
  112. |econ(0,2,0,0)| or |econ(79,2,0,0)| produces the full graph except 
  113. for the two special vertices.
  114. @d MAX_N 81 /* maximum number of vertices in constructed graph */
  115. @d NORM_N MAX_N-2 /* the number of normal SIC sectors */
  116. @d ADJ_SEC MAX_N-1 /* code number for the .{Adjustments} sector */
  117. @ The U.S. Bureau of Economic Analysis and the U.S. Bureau of the Census have
  118. assigned code numbers 1--79 to the individual sectors for which
  119. statistics are given in .{econ.dat}. These sector numbers are
  120. traditionally called Standard Industrial Classification (SIC) codes.
  121. If for some reason you want to know the SIC codes for
  122. all sectors represented by vertex |v| of a graph generated by |econ|,
  123. you can access them via a list of |Arc| nodes starting at the utility
  124. field |v->SIC_codes|.
  125. This list is linked by |next| fields in the usual way, and each
  126. SIC code appears in the |len| field; the |tip| field is unused.
  127. The special vertex .{Adjustments} is given code number~80; it is
  128. actually a composite of six different SIC categories, numbered 80--86 in their
  129. published tables.
  130. For example, if |n=80| and |omit=1|, each list will have length~1.
  131. Hence |v->SIC_codes->next| will equal |NULL| for each~|v|, and
  132. |v->SIC_codes->len| will be |v|'s SIC code, a number between 1 and~80.
  133. The special vertex .{Users} has no SIC code; it is the only vertex
  134. whose |SIC_codes| field will be null in the graph returned by |econ|.
  135. @d SIC_codes z.A /* utility field |z| leads to the SIC codes for a vertex */
  136. @ The total output of each sector, which also equals the total input of that
  137. sector, is placed in utility field |sector_total| of the corresponding vertex.
  138. @d sector_total y.I /* utility field |y| holds the total flow in and out */
  139. @(gb_econ.h@>=
  140. #define flow @tquad@> a.I
  141.    /* definitions of utility fields in the header file */
  142. #define SIC_codes @tquad@> z.A
  143. #define sector_total @tquad@> y.I
  144. @ If the |econ| routine encounters a problem, it returns |NULL|
  145. (.{NULL}), after putting a nonzero number into the external variable
  146. |panic_code|. This code number identifies the type of failure.
  147. Otherwise |econ| returns a pointer to the newly created graph, which
  148. will be represented with the data structures explained in {sc GB_,GRAPH}.
  149. (The external variable |panic_code| is itself defined in
  150. {sc GB_,GRAPH}.)
  151. @d panic(c) @+{@+panic_code=c;@+gb_trouble_code=0;@+return NULL;@+}
  152. @ The CEE/ file .{gb_econ.c} has the following overall shape:
  153. @p
  154. #include "gb_io.h" /* we will use the {sc GB_,IO} routines for input */
  155. #include "gb_flip.h"
  156.  /* we will use the {sc GB_,FLIP} routines for random numbers */
  157. #include "gb_graph.h"
  158.  /* and of course we'll use the {sc GB_,GRAPH} data structures */
  159. @h@#
  160. @<Type declarations@>@;
  161. @<Private variables@>@;
  162. @#
  163. Graph *econ(n,omit,threshold,seed)
  164.   unsigned long n; /* number of vertices desired */
  165.   unsigned long omit; /* number of special vertices to omit */
  166.   unsigned long threshold; /* minimum per-64K-age in arcs leading in */
  167.   long seed; /* random number seed */
  168. {@+@<Local variables@>@;@#
  169.   gb_init_rand(seed);
  170.   init_area(working_storage);
  171.   @<Check the parameters and adjust them for defaults@>;
  172.   @<Set up a graph with |n| vertices@>;
  173.   @<Read .{econ.dat} and note the binary tree structure@>;
  174.   @<Determine the |n| sectors to use in the graph@>;
  175.   @<Put the appropriate arcs into the graph@>;
  176.   if (gb_close()!=0)
  177.     panic(late_data_fault);
  178.      /* something's wrong with |"econ.dat"|; see |io_errors| */
  179.   gb_free(working_storage);
  180.   if (gb_trouble_code) {
  181.     gb_recycle(new_graph);
  182.     panic(alloc_fault); /* oops, we ran out of memory somewhere back there */
  183.   }
  184.   return new_graph;
  185. }
  186. @ @<Local var...@>=
  187. Graph *new_graph; /* the graph constructed by |econ| */
  188. register long j,k; /* all-purpose indices */
  189. Area working_storage; /* tables needed while |econ| does its thinking */
  190. @ @<Check the param...@>=
  191. if (omit>2) omit=2;
  192. if (n==0 || n>MAX_N-omit) n=MAX_N-omit;
  193. else if (n+omit<3) omit=3-n; /* we need at least one normal sector */
  194. if (threshold>65536) threshold=65536;
  195. @ @<Set up a graph with |n| vertices@>=
  196. new_graph=gb_new_graph(n);
  197. if (new_graph==NULL)
  198.   panic(no_room); /* out of memory before we're even started */
  199. sprintf(new_graph->id,"econ(%lu,%lu,%lu,%ld)",n,omit,threshold,seed);
  200. strcpy(new_graph->util_types,"ZZZZIAIZZZZZZZ");
  201. @* The economic tree.
  202. As we read in the data, we construct a sequential list of nodes,
  203. each of which represents either a micro-sector of the economy (one of
  204. the basic SIC sectors) or a macro-sector (which is the union of two subnodes).
  205. In more technical terms, the nodes form an extended binary tree,
  206. whose external nodes correspond to micro-sectors and whose internal nodes
  207. correspond to macro-sectors. The nodes of the tree appear in preorder.
  208. Subsequently we will do a variety of operations on this binary tree,
  209. proceeding either top-down (from the beginning of the list to the end)
  210. or bottom-up (from the end to the beginning).
  211. Each node is a rather large record, because we will store a complete
  212. vector of sector output data in each node.
  213. @<Type declarations@>=
  214. typedef struct node_struct { /* records for micro- and macro-sectors */
  215.   struct node_struct *rchild; /* pointer to right child of macro-sector */
  216.   char title[44]; /* |"Sector name"| */
  217.   long table[MAX_N+1]; /* outputs from this sector */
  218.   unsigned long total; /* total input to this sector ($=$ total output) */
  219.   long thresh; /* |flow| must exceed |thresh| in arcs to this sector */
  220.   long SIC; /* SIC code number; initially zero in macro-sectors */
  221.   long tag; /* 1 if this node will be a vertex in the graph */
  222.   struct node_struct *link; /* next smallest unexplored sector */
  223.   Arc *SIC_list; /* first item on list of SIC codes */
  224. } node;
  225. @ When we read the given data in preorder, we'll need a stack to remember
  226. what nodes still need to have their |rchild| pointer filled in.
  227. (There is a no need for an \{lchild} pointer, because the left child
  228. always follows its parent immediately in preorder.)
  229. @<Private v...@>=
  230. static node *stack[NORM_N+NORM_N];
  231. static node **stack_ptr; /* current position in |stack| */
  232. static node *node_block; /* array of nodes, specifies the tree in preorder */
  233. static node *node_index[MAX_N+1]; /* which node has a given SIC code */
  234. @ @<Local v...@>=
  235. register node *p,*pl,*pr; /* current node and its children */
  236. register node *q; /* register for list manipulation */
  237. @ @<Read .{econ.dat} and note the binary tree structure@>=
  238. node_block=gb_typed_alloc(2*MAX_N-3,node,working_storage);
  239. if (gb_trouble_code) panic(no_room+1); /* no room to copy the data */
  240. if (gb_open("econ.dat")!=0)
  241.   panic(early_data_fault);
  242.    /* couldn't open |"econ.dat"| using GraphBase conventions */
  243. @<Read and store the sector names and SIC numbers@>;
  244. for (k=1; k<=MAX_N; k++)
  245.  @<Read and store the output coefficients for sector |k|@>;
  246. @ The first part of .{econ.dat} specifies the nodes of the binary
  247. tree in preorder. Each line contains a node name
  248. followed by a colon, and the colon is followed by the SIC number if
  249. that node is a leaf.
  250. The tree is uniquely specified in this way,
  251. because of the nature of preorder. (Think of Polish prefix notation,
  252. in which a formula like `${+}x{+}xx$' means `${+}(x,{+}(x,x))$'; the
  253. parentheses in Polish notation are redundant.)
  254. The two special sector names don't appear in the file; we manufacture
  255. them ourselves.
  256. The program here is careful not to clobber itself in the
  257. presence of arbitrarily garbled data.
  258. @<Read and store the sector names...@>=
  259. stack_ptr=stack;
  260. for (p=node_block; p<node_block+NORM_N+NORM_N-1; p++) {@+register long c;
  261.   gb_string(p->title,':');
  262.   if (strlen(p->title)>43) panic(syntax_error); /* sector name too long */
  263.   if (gb_char()!=':') panic(syntax_error+1); /* missing colon */
  264.   p->SIC=c=gb_number(10);
  265.   if (c==0) /* macro-sector */
  266.     *stack_ptr++=p; /* left child is |p+1|, we'll know |rchild| later */
  267.   else { /* micro-sector; |p+1| will be somebody's right child */
  268.     node_index[c]=p;
  269.     if (stack_ptr>stack) (*--stack_ptr)->rchild=p+1;
  270.   }
  271.   if (gb_char()!='n') panic(syntax_error+2); /* garbage on the line */
  272.   gb_newline();
  273. }
  274. if (stack_ptr!=stack) panic(syntax_error+3); /* tree malformed */
  275. for (k=NORM_N;k;k--) if (node_index[k]==0)
  276.   panic(syntax_error+4); /* SIC code not mentioned in the tree */
  277. strcpy(p->title,"Adjustments");@+p->SIC=ADJ_SEC;@+node_index[ADJ_SEC]=p;
  278. strcpy((p+1)->title,"Users");@+node_index[MAX_N]=p+1;
  279. @ The remaining part of .{econ.dat} is an $81times80$ matrix in which
  280. the $k$th row contains the outputs of sector~$k$ to all sectors except
  281. .{Users}. Each row consists of a blank line followed by 8 data lines;
  282. each data line contains 10 numbers separated by commas.
  283. Zeroes are represented by |""| instead of by |"0"|.
  284. For example, the data line $$hbox{tt
  285. 8490,2182,42,467,,,,,,}$$ follows the initial blank line; it means
  286. that sector~1 output 8490 million dollars to itself, $2182M to
  287. sector~2, dots, $0M to sector~10.
  288. @<Read and store the output...@>=
  289. {@+register long s=0; /* row sum */
  290.   register long x; /* entry read from .{econ.dat} */
  291.   if (gb_char()!='n') panic(syntax_error+5);
  292.    /* blank line missing between rows */
  293.   gb_newline();
  294.   p=node_index[k];
  295.   for (j=1;j<MAX_N;j++) {
  296.     p->table[j]=x=gb_number(10);@+s+=x;
  297.     node_index[j]->total+=x;
  298.     if ((j%10)==0) {
  299.       if (gb_char()!='n') panic(syntax_error+6);
  300.        /* out of synch in input file */
  301.       gb_newline();
  302.     }@+else if (gb_char()!=',') panic(syntax_error+7);
  303.      /* missing comma after entry */
  304.   }
  305.   p->table[MAX_N]=s; /* sum of |table[1]| through |table[80]| */
  306. }
  307. @* Growing a subtree.
  308. Once all the data appears in |node_block|, we want to extract from it
  309. and combine~it as specified by parameters |n|, |omit|, and |seed|.
  310. This amalgamation process effectively prunes the tree; it can also be
  311. regarded as a procedure that grows a subtree of the full economic tree.
  312. @<Determine the |n| sectors to use in the graph@>=
  313. {@+long l=n+omit-2; /* the number of leaves in the desired subtree */
  314.   if (l==NORM_N) @<Choose all sectors@>@;
  315.   else if (seed) @<Grow a random subtree with |l| leaves@>@;
  316.   else @<Grow a subtree with |l| leaves by subdividing largest sectors first@>;
  317. }
  318. @ The chosen leaves of our subtree are identified by having their
  319. |tag| field set to~1.
  320. @<Choose all sectors@>=
  321. for (k=NORM_N;k;k--) node_index[k]->tag=1;
  322. @ To grow the |l|-leaf subtree when |seed=0|, we first pass over the
  323. tree bottom-up to compute the total input (and output) of each macro-sector;
  324. then we proceed from the top down to subdivide sectors in decreasing
  325. order of their total input. This process provides a good introduction to the
  326. bottom-up and top-down tree methods we will be using in several other
  327. parts of the program.
  328. The |special| node is used here for two purposes: It is the head of a
  329. linked list of unexplored nodes, sorted by decreasing order of
  330. their |total| fields; and it appears at the end of that list, because
  331. |special->total=0|.
  332. @<Grow a subtree with |l| leaves by subdividing largest sectors first@>=
  333. {@+register node *special=node_index[MAX_N];
  334.      /* the .{Users} node at the end of |node_block| */
  335.   for (p=node_index[ADJ_SEC]-1;p>=node_block;p--) /* bottom up */
  336.     if (p->rchild)
  337.       p->total=(p+1)->total+p->rchild->total;
  338.   special->link=node_block;@+node_block->link=special; /* start at the root */
  339.   k=1; /* |k| is the number of nodes we have tagged or put onto the list */
  340.   while (k<l) @<If the first node on the list is a leaf, delete it and tag it;
  341.                 otherwise replace it by its two children@>;
  342.   for (p=special->link;p!=special;p=p->link)
  343.     p->tag=1; /* tag everything on the list */
  344. }
  345. @ @<If the first node on the list is a leaf,...@>=
  346. {
  347.   p=special->link; /* remove |p|, the node with greatest |total| */
  348.   special->link=p->link;
  349.   if (p->rchild==0) p->tag=1; /* |p| is a leaf */
  350.   else {
  351.     pl=p+1;@+pr=p->rchild;
  352.     for (q=special;q->link->total>pl->total;q=q->link) ;
  353.     pl->link=q->link;@+q->link=pl; /* insert left child in its proper place */
  354.     for (q=special;q->link->total>pr->total;q=q->link) ;
  355.     pr->link=q->link;@+q->link=pr; /* insert right child in its proper place */
  356.     k++;
  357.   }
  358. }
  359. @ We can obtain a uniformly distributed |l|-leaf subtree of a given tree
  360. by choosing the root when |l=1| or by using the following idea when |l>1|:
  361. Suppose the given tree~$T$ has subtrees $T_0$ and $T_1$. Then it has
  362. $T(l)$ subtrees with |l|~leaves, where $T(l)=sum_k T_0(k)T_1(l-k)$.
  363. We choose a random number $r$ between 0 and $T(l)-1$, and we find the
  364. smallest $m$ such that $sum_{kle m}T_0(k)T_1(l-k)>r$. Then we
  365. proceed recursively to
  366. compute a random $m$-leaf subtree of~$T_0$ and a random $(l-m)$-leaf
  367. subtree of~$T_1$.
  368. A difficulty arises when $T(l)$ is $2^{31}$ or more. But then we can replace
  369. $T_0(k)$ and $T_1(l-k)$ in the formulas above by $lceil T_0(k)/d_0rceil$
  370. and $lceil T_1(k)/d_1rceil$, respectively, where $d_0$ and $d_1$ are
  371. arbitrary constants; this yields smaller values
  372. $T(l)$ that define approximately the same distribution of~$k$.
  373. The program here computes the $T(l)$ values bottom-up, then grows a
  374. random tree top-down. If node~|p| is not a leaf, its |table[0]| field
  375. will be set to the number of leaves below it; and its |table[l]| field
  376. will be set to $T(l)$, for |1<=l<=table[0]|.
  377. The data in .{econ.dat} is sufficiently simple that most of the $T(l)$
  378. values are less than $2^{31}$. We need to scale them
  379. down to avoid overflow only at the root node of the tree; this
  380. case is handled separately.
  381. We set the |tag| field of a node equal to the number of leaves to be
  382. grown in the subtree rooted at that node. This convention is consistent
  383. with our previous stipulation that |tag=1| should characterize the
  384. nodes that are chosen to be vertices.
  385. @<Grow a random subtree with |l| leaves@>=
  386. {
  387.   node_block->tag=l;
  388.   for (p=node_index[ADJ_SEC]-1;p>node_block;p--) /* bottom up, except root */
  389.     if (p->rchild) @<Compute the $T(l)$ values for subtree |p|@>;
  390.   for (p=node_block;p<node_index[ADJ_SEC];p++) /* top down, from root */
  391.     if (p->tag>1) {
  392.       l=p->tag;
  393.       pl=p+1;@+pr=p->rchild;
  394.       if (pl->rchild==NULL) {
  395.         pl->tag=1;@+pr->tag=l-1;
  396.       }@+else if (pr->rchild==NULL) {
  397.         pl->tag=l-1;@+pr->tag=1;
  398.       }@+else @<Stochastically determine the number of leaves to grow in
  399.                     each of |p|'s children@>;
  400.     }
  401. }
  402. @ Here we are essentially multiplying two generating functions.
  403. Suppose $f(z)=sum_l T(l)z^l$; then we are computing $f_p(z)=
  404. z+f_{pl}(z)f_{pr}(z)$.
  405. @<Compute the $T(l)$ values for subtree |p|@>=
  406. {
  407.   pl=p+1;@+pr=p->rchild;
  408.   p->table[1]=p->table[2]=1; /* $T(1)$ and $T(2)$ are always 1 */
  409.   if (pl->rchild==0) { /* left child is a leaf */
  410.     if (pr->rchild==0) p->table[0]=2; /* and so is the right child */
  411.     else { /* no, it isn't */
  412.       for (k=2;k<=pr->table[0];k++) p->table[1+k]=pr->table[k];
  413.       p->table[0]=pr->table[0]+1;
  414.     }
  415.   }@+else if (pr->rchild==0) { /* right child is a leaf */
  416.     for (k=2;k<=pl->table[0];k++) p->table[1+k]=pl->table[k];
  417.     p->table[0]=pl->table[0]+1;
  418.   }@+else { /* neither child is a leaf */
  419.     @<Set |p->table[2]|, |p->table[3]|, dots to convolution of
  420.       |pl| and |pr| table entries@>;
  421.     p->table[0]=pl->table[0]+pr->table[0];
  422.   }
  423. }
  424. @ @<Set |p->table[2]|, |p->table[3]|, dots to convolution...@>=
  425. p->table[2]=0;
  426. for (j=pl->table[0];j;j--) {@+register long t=pl->table[j];
  427.   for (k=pr->table[0];k;k--)
  428.     p->table[j+k]+=t*pr->table[k];
  429. }
  430. @ @<Stochastically determine the number of leaves to grow...@>=
  431. {@+register long ss,rr;
  432.   j=0; /* we will set |j=1| if scaling is necessary at the root */
  433.   if (p==node_block) {
  434.     ss=0;
  435.     if (l>29 && l<67) {
  436.       j=1; /* more than $2^{31}$ possibilities exist */
  437.       for (k=(l>pr->table[0]? l-pr->table[0]: 1);k<=pl->table[0] && k<l;k++)
  438.         ss+=((pl->table[k]+0x3ff)>>10)*pr->table[l-k];
  439.               /* scale with $d_0=1024$, $d_1=1$ */
  440.     }@+else
  441.       for (k=(l>pr->table[0]? l-pr->table[0]: 1);k<=pl->table[0] && k<l;k++)
  442.         ss+=pl->table[k]*pr->table[l-k];
  443.   }@+else ss=p->table[l];
  444.   rr=gb_unif_rand(ss);
  445.   if (j)
  446.     for (ss=0,k=(l>pr->table[0]? l-pr->table[0]: 1);ss<=rr;k++)
  447.       ss+=((pl->table[k]+0x3ff)>>10)*pr->table[l-k];
  448.   else for (ss=0,k=(l>pr->table[0]? l-pr->table[0]: 1);ss<=rr;k++)
  449.       ss+=pl->table[k]*pr->table[l-k];
  450.   pl->tag=k-1;@+pr->tag=l-k+1;
  451. }
  452. @* Arcs.
  453. In the general case, we have to combine some of the basic micro-sectors
  454. into macro-sectors by adding together the appropriate input/output
  455. coefficients. This is a bottom-up pruning process.
  456. Suppose |p| is being formed as the union of |pl| and~|pr|.
  457. Then the arcs leading out of |p| are obtaining by summing the numbers
  458. on arcs leading out of |pl| and~|pr|; the arcs leading into |p| are
  459. obtained by summing the numbers on arcs leading into |pl| and~|pr|;
  460. the arcs from |p| to itself are obtained by summing the four numbers
  461. on arcs leading from |pl| or~|pr| to |pl| or~|pr|.
  462. We maintain the |node_index| table so that its non-|NULL| entries
  463. contain all the currently active nodes. When |pl| and~|pr| are
  464. being pruned in favor of~|p|, node |p|~inherits |pl|'s place in
  465. |node_index|; |pr|'s former place becomes~|NULL|.
  466. @<Put the appropriate arcs into the graph@>=
  467. @<Prune the sectors that are used in macro-sectors, and form
  468.   the lists of SIC sector codes@>;
  469. @<Make the special nodes invisible if they are omitted, visible otherwise@>;
  470. @<Compute individual thresholds for each chosen sector@>;
  471. {@+register Vertex *v=new_graph->vertices+n;
  472.   for (k=MAX_N;k;k--)
  473.     if ((p=node_index[k])!=NULL) {
  474.       vert_index[k]=--v;
  475.       v->name=gb_save_string(p->title);
  476.       v->SIC_codes=p->SIC_list;
  477.       v->sector_total=p->total;
  478.     }@+else vert_index[k]=NULL;
  479.   if (v!=new_graph->vertices)
  480.     panic(impossible); /* bug in algorithm; this can't happen */
  481.   for (j=MAX_N;j;j--)
  482.     if ((p=node_index[j])!=NULL) {@+register Vertex *u=vert_index[j];
  483.       for (k=MAX_N;k;k--)
  484.         if ((v=vert_index[k])!=NULL)
  485.           if (p->table[k]!=0 && p->table[k]>node_index[k]->thresh) {
  486.             gb_new_arc(u,v,1L);
  487.             u->arcs->flow=p->table[k];
  488.           }
  489.       }
  490. }
  491. @ @<Private v...@>=
  492. static Vertex *vert_index[MAX_N+1]; /* the vertex assigned to an SIC code */
  493. @ The theory underlying this step is the following, for integers
  494. $a,b,c,d$ with $b,d>0$:
  495. $$ {aover b}>{cover d} qquadiffqquad
  496.   a>biggllfloor{bover d}biggrrfloor,c +
  497.        biggllfloor{(bbmod d)cover d}biggrrfloor,.$$
  498. In our case, $b=hbox{|p->total|}$ and $c=thresholdle d=65536=2^{16}$, hence
  499. the multiplications cannot overflow. (But they can come awfully darn close.)
  500. @<Compute individual thresholds for each chosen sector@>=
  501. for (k=MAX_N;k;k--)
  502.   if ((p=node_index[k])!=NULL) {
  503.     if (threshold==0) p->thresh=-99999999;
  504.     else p->thresh=((p->total>>16)*threshold)+
  505.       (((p->total&0xffff)*threshold)>>16);
  506.   }
  507. @ @<Prune the sectors that are used in macro-sectors, and form
  508.   the lists of SIC sector codes@>=
  509. for (p=node_index[ADJ_SEC];p>=node_block;p--) { /* bottom up */
  510.   if (p->SIC) { /* original leaf */
  511.     p->SIC_list=gb_virgin_arc();
  512.     p->SIC_list->len=p->SIC;
  513.   }@+else {
  514.     pl=p+1;@+pr=p->rchild;
  515.     if (p->tag==0) p->tag=pl->tag+pr->tag;
  516.     if (p->tag<=1) @<Replace |pl| and |pr| by their union, |p|@>;
  517.   }
  518. }
  519.     
  520. @ @<Replace |pl| and |pr| by their union, |p|@>=
  521. {@+register Arc *a=pl->SIC_list;
  522.   register long jj=pl->SIC, kk=pr->SIC;
  523.   p->SIC_list=a;
  524.   while (a->next) a=a->next;
  525.   a->next=pr->SIC_list;
  526.   for (k=MAX_N;k;k--)
  527.     if ((q=node_index[k])!=NULL) {
  528.       if (q!=pl && q!=pr) q->table[jj]+=q->table[kk];
  529.       p->table[k]=pl->table[k]+pr->table[k];
  530.     }
  531.   p->total=pl->total+pr->total;
  532.   p->SIC=jj;
  533.   p->table[jj]+=p->table[kk];
  534.   node_index[jj]=p;
  535.   node_index[kk]=NULL;
  536. }
  537. @ If the .{Users} vertex is not omitted, we need to compute each
  538. sector's total final demand, which is calculated so that the row sums
  539. and column sums of the input/output coefficients come out equal. We've
  540. already computed the column sum, |p->total|; we've also computed
  541. |p->table[1]+@[@thbox{$cdots$}@>@]+p->table[ADJ_SEC]|, and put it into
  542. |p->table[MAX_N]|. So now we want to replace |p->table[MAX_N]| by
  543. |p->total-p->table[MAX_N]|. As remarked earlier, this quantity might
  544. be negative.
  545. In the special node |p| for the .{Users} vertex, the preliminary
  546. processing has made |p->total=0|; moreover, |p->table[MAX_N]| is the
  547. sum of value added, or GNP.  We want to switch those fields.
  548. We don't have to set the |tag| fields to 1 in the special nodes, because
  549. the remaining parts of the arc-generation algorithm don't look at those fields.
  550. @<Make the special nodes invisible if they are omitted, visible otherwise@>=
  551. if (omit==2) node_index[ADJ_SEC]=node_index[MAX_N]=NULL;
  552. else if (omit==1) node_index[MAX_N]=NULL;
  553. else {
  554.   for (k=ADJ_SEC;k;k--)
  555.     if ((p=node_index[k])!=NULL) p->table[MAX_N]=p->total-p->table[MAX_N];
  556.   p=node_index[MAX_N]; /* the special node */
  557.   p->total=p->table[MAX_N];
  558.   p->table[MAX_N]=0;
  559. }
  560.       
  561. @* Index. As usual, we close with an index that
  562. shows where the identifiers of \{gb_econ} are defined and used.