summaryrefslogtreecommitdiff
path: root/support/graphbase/gb_basic.w
blob: eff1ee68fde012149d91b70cc3f33269046e3d47 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063
2064
2065
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112
2113
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
2143
2144
2145
2146
2147
2148
2149
2150
2151
2152
2153
2154
2155
2156
2157
2158
2159
2160
2161
2162
2163
2164
2165
2166
2167
2168
2169
2170
2171
2172
2173
2174
2175
2176
2177
2178
2179
2180
2181
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191
2192
2193
2194
2195
2196
2197
2198
2199
2200
2201
2202
2203
2204
2205
2206
2207
2208
2209
2210
2211
2212
2213
2214
2215
2216
2217
2218
2219
2220
2221
2222
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238
2239
2240
2241
2242
2243
2244
2245
2246
2247
2248
2249
2250
2251
2252
2253
2254
2255
2256
2257
2258
2259
2260
2261
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274
2275
2276
2277
2278
2279
2280
2281
2282
2283
2284
2285
2286
2287
2288
2289
2290
2291
2292
2293
2294
2295
2296
2297
2298
2299
2300
2301
2302
2303
2304
2305
2306
2307
2308
2309
2310
2311
2312
2313
2314
2315
2316
2317
2318
2319
2320
2321
2322
2323
2324
2325
2326
2327
2328
2329
2330
2331
2332
2333
2334
2335
2336
2337
2338
2339
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349
2350
2351
2352
2353
2354
2355
2356
2357
2358
2359
2360
2361
2362
2363
2364
2365
2366
2367
2368
2369
2370
2371
2372
2373
2374
2375
2376
2377
2378
2379
2380
2381
2382
2383
2384
2385
2386
2387
2388
2389
2390
2391
2392
2393
2394
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404
2405
2406
2407
2408
2409
2410
2411
2412
2413
2414
2415
% This file is part of the Stanford GraphBase (c) Stanford University 1992
\def\title{GB\_\thinspace BASIC}
@i boilerplate.w %<< legal stuff: PLEASE READ IT BEFORE MAKING ANY CHANGES!

\prerequisite{GB\_\thinspace GRAPH}
@* Introduction. This GraphBase module contains six subroutines that generate
standard graphs of various types, together with six routines that combine or
transform existing graphs.

Simple examples of the use of these routines can be found in the
demonstration programs |queen| and |queen_wrap|.

@f Graph int /* |gb_graph| defines the |Graph| type and a few others */
@f Vertex int
@f Arc int
@f Area int

@<gb_basic.h@>=
extern Graph *board(); /* moves on generalized chessboards */
extern Graph *simplex(); /* generalized triangular configurations */
extern Graph *subsets(); /* patterns of subset intersection */
extern Graph *perms(); /* permutations of a multiset */
extern Graph *parts(); /* partitions of an integer */
extern Graph *binary(); /* binary trees */
@#
extern Graph *complement(); /* the complement of a graph */
extern Graph *gunion(); /* the union of two graphs */
extern Graph *intersection(); /* the intersection of two graphs */
extern Graph *lines(); /* the line graph of a graph */
extern Graph *product(); /* the product of two graphs */
extern Graph *induced(); /* a graph induced from another */

@ The \Cee\ file \.{gb\_basic.c} has the following overall shape:

@p
#include "gb_graph.h" /* we use the |gb_graph| data structures */
@#
@<Private variables@>@;
@<Basic subroutines@>@;
@<Applications of basic subroutines@>@;

@ Several of the programs below allocate arrays that will be freed again
before the routine is finished.

@<Private variables@>=
static Area working_storage;

@ If a graph-generating subroutine encounters a problem, it returns |NULL|
(that is, \.{NULL}), after putting a code number into the external variable
|panic_code|. This code number identifies the type of failure.
Otherwise the routine returns a pointer to the newly created graph, which
will be represented with the data structures explained in |gb_graph|.
(The external variable |@!panic_code| is itself defined in |gb_graph|.)

@d panic(c) 
  {@+panic_code=c;
    gb_free(working_storage);
    gb_alloc_trouble=0;
    return NULL;
  }

@ The names of vertices are sometimes formed from the names of other
vertices, or from potentially long sequences of numbers. We assemble
them in the |buffer| array, which is sufficiently long that the
vast majority of applications will be unconstrained by size limitations.
The programs do always make sure that |BUF_SIZE| is not exceeded, but
they assume that it is rather large.

@d BUF_SIZE 4096

@<Private v...@>=
static char buffer[BUF_SIZE];

@*Grids and game boards. The subroutine call
`|board(n1,n2,n3,n4,piece,wrap,directed)|'
constructs a graph based on the moves of generalized chesspieces on a
generalized rectangular board. Each vertex of the graph corresponds to a
position on the board. Each arc of the graph corresponds to a move from
one position to another.

The first parameters, |n1| through |n4|, specify the size of the board.
If, for example, a two-dimensional board with $n_1$ rows and $n_2$ columns
is desired, you set $|n1|=n_1$, $|n2|=n_2$, and $|n3|=0$; the resulting
graph will have $n_1n_2$ vertices. If you want a three-dimensional
board with $n_3$ layers, set $|n3|=n_3$ and $n_4=0$. If you want
a four-dimensional board, put the number of 4th coordinates in~|n4|.
If you want a $d$-dimensional board with $2^d$ positions, set |n1=2|
and |n2=-d|.

In general, the |board| subroutine determines the dimensions by scanning the
sequence |(n1,n2,n3,n4,0)=@t$(n_1,n_2,n_3,n_4,0)$@>| from left to right
until coming to the first nonpositive parameter $n_{k+1}$. If $k=0$
(i.e., if |n1<=0|), the default size $8\times8$ will be used; this is
an ordinary chessboard with 8~rows and 8~columns. Otherwise if $n_{k+1}=0$,
the board will have $k$~dimensions $n_1$, \dots,~$n_k$. Otherwise
we must have $n_{k+1}<0$; in this case the board will have $d=\vert n_{k+1}
\vert$ dimensions, chosen as the first $d$ elements of the infinite
periodic sequence $(n_1,\ldots,n_k,n_1,\ldots,n_k,n_1,\ldots\,)$.
For example, the specification |(n1,n2,n3,n4)=(2,3,5,-7)| is about as
tricky as you can get. It produces a seven-dimensional board with
dimensions $(n_1,\ldots,n_7)=(2,3,5,2,3,5,2)$, hence a graph with
$2\cdot3\cdot5\cdot2\cdot3\cdot5\cdot2=1800$ vertices.

The |piece| parameter specifies the legal moves of a generalized chesspiece.
If |piece>0|, a move from position~|u| to position~|v| is considered legal
if and only if the Euclidean distance between points |u| and~|v| is
equal to $\sqrt{\vphantom1\smash{|piece|}}$.
For example, if |piece=1| and if we have a
two-dimensional board, the legal moves from $(x,y)$ are to $(x,y\pm1)$ and
$(x\pm1,y)$; these are the moves of a so-called wazir, the only moves that
a king and a rook can both make. If |piece=2|, the legal moves from $(x,y)$
are to $(x\pm1,y\pm1)$; these are the four moves that a king and a bishop
can both make. (A piece that can make only these moves was called a ``fers''
in ancient Muslim chess.) If |piece=5|, the legal moves are those of a
knight, from $(x,y)$ to $(x\pm1,y\pm2)$ or to $(x\pm2,y\pm1)$. If |piece=3|,
there are no legal moves on a two-dimensional board, but moves from
$(x,y,z)$ to $(x\pm1,y\pm1,z\pm1)$ would be legal in three dimensions.
If |piece=0|, it is changed to the default value |piece=1|.

If the value of |piece| is negative, arbitrary multiples of the basic moves
for $\vert|piece|\vert$ are permitted. For example, |piece=-1| defines the
moves of a rook, from $(x,y)$ to $(x\pm a,y)$ or to $(x,y\pm a)$ for all
$a>0$; |piece=-2| defines the moves of a bishop, from $(x,y)$ to
$(x\pm a,y\pm a)$. The literature of ``fairy chess'' assigns standard names
to the following |piece| values: $\rm wazir=1$, $\rm fers=2$, $\rm dabbaba=4$,
$\rm knight=5$, $\rm alfil=8$, $\rm camel=10$, $\rm zebra=13$, $\rm giraffe
=17$, $\rm fiveleaper=25$, $\hbox{root-50-leaper}=50$, etc.; $\rm rook=-1$,
$\rm bishop=-2$, $\rm unicorn=-3$, $\rm dabbabarider=-4$, $\rm nightrider=-5$,
$\rm alfilrider=-8$, $\rm camelrider=-10$, etc.

To generate a board with the moves of a king, you can use the |gunion|
subroutine below to take the union of boards with |piece=1| and
|piece=2|. Similarly, you can get queen moves by taking the union of
boards with |piece=-1| and |piece=-2|.

If |piece>0|, all arcs of the graph will have length~1. If |piece<0|, the
length of each arc will be the number of multiples of a basic move that
produced the arc.

@ If the |wrap| parameter is nonzero, it specifies a subset of coordinates
in which values are computed modulo the corresponding size.
For example, the coordinates $(x,y)$ for vertices on a two-dimensional
board are restricted to the range $0\le x<n_1$, $0\le y<n_2$; when
|wrap=0|, a move from $(x,y)$ to $(x+\delta_1,y+\delta_2)$ is
therefore legal only if $0\le x+\delta_1<n_1$ and $0\le
y+\delta_2<n_2$. But when |wrap=1|, the $x$~coordinates are allowed to
``wrap around''; the move would then be made to $((x+\delta_1)\bmod
n_1,y+\delta_2)$, provided that $0\le y+\delta_2<n_2$. Setting
|wrap=1| effectively makes the board into a cylinder instead of a
rectangle. Similarly, the $y$~coordinates are allowed to wrap around
when |wrap=2|. Both $x$ and $y$ coordinates are treated modulo their
corresponding sizes when |wrap=3|; the board is then effectively a
torus.  In general, coordinates $k_1$, $k_2$, \dots~ will wrap around
when $|wrap|=2^{k_1-1}+2^{k_2-1}+\cdots\,$. Setting |wrap=-1| causes all
coordinates to be computed modulo their size.

The graph constructed by |board| will be undirected unless |directed!=0|.
Directed |board| graphs will be acyclic when |wrap=0|, but they may
have cycles when |wrap!=0|. Precise rules defining the directed arcs
are given below.

Several important special cases are worth noting: To get the complete graph
on |n| vertices, you can say |board(n,0,0,0,-1,0,0)|. To get the
transitive tournament on |n| vertices, i.e., the directed graph
with arcs from |u| to |v| when |u<v|, you can say |board(n,0,0,0,-1,0,1)|.
To get the empty graph on |n| vertices, you can say |board(n,0,0,0,2,0,0)|.
To get a circuit (undirected) or a cycle (directed) of length~|n|,
you can say |board(n,0,0,0,1,1,0)| and |board(n,0,0,0,1,1,1)|,
respectively.

@(gb_basic.h@>=
#define complete(n) @[board(n,0,0,0,-1,0,0)@]
#define transitive(n) @[board(n,0,0,0,-1,0,1)@]
#define empty(n) @[board(n,0,0,0,2,0,0)@]
#define circuit(n) @[board(n,0,0,0,1,1,0)@]
#define cycle(n) @[board(n,0,0,0,1,1,1)@]


@ @<Basic subroutines@>=
Graph *board(n1,n2,n3,n4,piece,wrap,directed)
  int n1,n2,n3,n4; /* size of board desired */
  int piece; /* type of moves desired */
  long wrap; /* mask for coordinate positions that wrap around */
  int directed; /* should the graph be directed? */
{@+@<Vanilla local variables@>@;
  long n; /* total number of vertices */
  int p; /* $\vert|piece|\vert$ */
  int l; /* length of current arc */
  @<Normalize the board-size parameters@>;
  @<Set up a graph with |n| vertices@>;
  @<Insert arcs or edges for all legal moves@>;
  if (gb_alloc_trouble) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* alas, we ran out of memory somewhere back there */
  }
  return new_graph;
}

@ Most of the subroutines in |gb_basic| use the following local
variables.

@<Vanilla local variables@>=
Graph *new_graph; /* the graph being constructed */
register int i,j,k; /* all-purpose indices */
register int d; /* the number of dimensions */
register Vertex *v; /* the current vertex of interest */
register long s; /* accumulator */

@ Several arrays will facilitate the calculations that |board| needs to make:
The number of distinct values in coordinate position~$k$ will be |nn[k]|;
this coordinate position will wrap around if and only if |wr[k]!=0|.
The current moves under consideration will be from $(x_1,\ldots,x_d)$
to $(x_1+\delta_1,\ldots, x_k+\delta_k)$, where $\delta_k$ is stored
in |del[k]|. An auxiliary array |sig| holds the sums
$\sigma_k=\delta_1^2+\cdots+\delta_{k-1}^2$.  Additional arrays |xx|
and |yy| hold coordinates of vertices before and after a move is made.

Some of these arrays are also used for other purposes by other programs
besides |board|; we will meet those programs later.

We limit the number of dimensions to 91 or less. This is hardly a limitation,
since the number of vertices would be astronomical even if the dimensionality
were only half this big. But some of our later programs will be able
to make good use of 40 or 50 dimensions and perhaps more; the number 91
is an upper limit imposed by the number of standard printable characters
(see the convention for vertex names in the |perms| routine).

@d MAX_D 91

@<Private...@>=
static int nn[MAX_D+1]; /* component sizes */
static int wr[MAX_D+1]; /* does this component wrap around? */
static int del[MAX_D+1]; /* displacements for the current move */
static int sig[MAX_D+2]; /* partial sums of squares of displacements */
static int xx[MAX_D+1], yy[MAX_D+1]; /* coordinate values */

@ @<Normalize the board-size parameters@>=
if (piece==0) piece=1;
if (n1<=0) {@+n1=n2=8;@+n3=0;@+}
nn[1]=n1;
if (n2<=0) {@+k=2;@+d=-n2;@+n3=n4=0;@+}
else {
  nn[2]=n2;
  if (n3<=0) {@+k=3;@+d=-n3;@+n4=0;@+}
  else {
    nn[3]=n3;
    if (n4<=0) {@+k=4;@+d=-n4;@+}
    else {@+nn[4]=n4;@+d=4;@+goto done;@+}
  }
}
if (d==0) {@+d=k-1;@+goto done;@+}
@<Compute component sizes periodically for |d| dimensions@>;
done: /* now |nn[1]| through |nn[d]| are set up */

@ At this point, |nn[1]| through |nn[k-1]| are the component sizes
that should be replicated periodically. In unusual cases, the number
of dimensions might not be as large as the number of specifications.

@<Compute component sizes periodically...@>=
if (d>MAX_D) panic(bad_specs); /* too many dimensions */
for (j=1; k<=d; j++,k++) nn[k]=nn[j];

@ We want to make the subroutine idiot-proof, so we use floating-point
arithmetic to make sure that boards with more than a billion cells have
not been specified.

@d MAX_NNN 1000000000.0

@<Set up a graph with |n| vertices@>=
{@+float nnn; /* approximate size */
  for (n=1,nnn=1.0,j=1; j<=d; j++) {
    nnn *= (float)nn[j];
    if (nnn>MAX_NNN) panic(very_bad_specs); /* way too big */
    n *= nn[j]; /* this multiplication cannot cause integer overflow */
  }
  new_graph=gb_new_graph(n);
  if (new_graph==NULL)
    panic(no_room); /* out of memory before we're even started */
  sprintf(new_graph->id,"board(%d,%d,%d,%d,%d,%ld,%d)",
    n1,n2,n3,n4,piece,wrap,directed?1:0);
  strcpy(new_graph->format,"ZZZIIIZZZZZZZZ");
  @<Give names to the vertices@>;
}

@ The symbolic name of a board position like $(3,1)$ will be the string
`\.{3.1}'. The first three coordinates are also stored as integers, in
utility fields |x.i|, |y.i|, and |z.i|, because immediate access to
those values will be helpful in certain applications. (The coordinates can,
of course, always be recovered in a slower fashion from the vertex name,
via |sscanf|.)

The process of assigning coordinate values and names is equivalent to
adding unity in a mixed-radix number system. Vertex $(x_1,\ldots,x_d)$
will be in position $x_1n_2\ldots n_d+\cdots+x_{d-1}n_d+x_d$ relative
to the first vertex of the new graph; therefore it is also possible to
deduce the coordinates of a vertex from its address.

@<Give names...@>=
{@+register char *p; /* string pointer */
  nn[0]=xx[0]=xx[1]=xx[2]=xx[3]=0;
  for (k=4;k<=d;k++) xx[k]=0;
  for (v=new_graph->vertices;;v++) {
    p=buffer;
    for (k=1;k<=d;k++) {
      sprintf(p,".%d",xx[k]);
      while (*p) p++;
    }
    v->name=gb_save_string(&buffer[1]); /* omit |buffer[0]|, which is |'.'| */
    v->x.i=xx[1];@+v->y.i=xx[2];@+v->z.i=xx[3];
    for (k=d;xx[k]+1==nn[k];k--) xx[k]=0;
    if (k==0) break; /* a ``carry'' has occurred all the way to the left */
    xx[k]++; /* increase coordinate |k| */
  }
}
      
@ Now we come to a slightly tricky part of the routine, the move generator.
Let $p=\vert|piece|\vert$. The outer loop of this procedure runs through all
solutions of the equation $\delta_1^2+\cdots+\delta_d^2=p$, where the
$\delta$'s are nonnegative integers. Within that loop, we attach signs
to the $\delta$'s, but always leaving $\delta_k$ positive if $\delta_1=
\cdots=\delta_{k-1}=0$. For every such vector~$\delta$, we generate moves
from |v| to $v+\delta$ for every vertex |v|. When |directed=0|,
we use |gb_new_edge| instead of |gb_new_arc|, so that the reverse arc
from $v+\delta$ to~|v| is also generated.

@<Insert arcs or edges for all legal moves@>=
@<Initialize the |wr|, |sig|, and |del| tables@>;
p=piece;
if (p<0) p=-p;
while (1) {
  @<Advance to the next nonnegative |del| vector, or |break| if done@>;
  while (1) {
    @<Generate moves for the current |del| vector@>;
    @<Advance to the next signed |del| vector, or restore |del|
          to nonnegative values and |break|@>;
  }
}

@ The \Cee\ language does not define |>>| unambiguously. If |w| is negative,
the assignment `|w>>=1|' here should keep |w| negative.
(However, this technicality doesn't matter except in highly unusual cases
when there are more than 32 dimensions.)
@^system dependencies@>

@<Initialize the |wr|, |sig|, and |del| tables@>=
{@+register long w=wrap;
  for (k=1;k<=d;k++,w>>=1) {
    wr[k]=w&1;
    del[k]=sig[k]=0;
  }
  sig[0]=del[0]=sig[d+1]=0;
}

@ @<Advance to the next nonnegative |del|...@>=
for (k=d;sig[k]+(del[k]+1)*(del[k]+1)>p;k--) del[k]=0;
if (k==0) break;
del[k]++;
sig[k+1]=sig[k]+del[k]*del[k];
for (k++;k<=d;k++) sig[k+1]=sig[k];
if (sig[d+1]<p) continue;

@ @<Advance to the next signed |del| vector, or restore |del|
          to nonnegative values and |break|@>=
for (k=d;del[k]<=0;k--) del[k]=-del[k];
if (sig[k]==0) break; /* all but |del[k]| were negative or zero */
del[k]=-del[k]; /* some entry preceding |del[k]| is positive */

@ We use the mixed-radix addition technique again when generating moves.

@<Generate moves for the current |del| vector@>=
for (k=1;k<=d;k++) xx[k]=0;
for (v=new_graph->vertices;;v++) {
  @<Generate moves from |v| corresponding to |del|@>;
  for (k=d;xx[k]+1==nn[k];k--) xx[k]=0;
  if (k==0) break; /* a ``carry'' has occurred all the way to the left */
  xx[k]++; /* increase coordinate |k| */
}

@ The legal moves when |piece| is negative are derived as follows, in
the presence of possible wraparound: Starting at $(x_1,\ldots,x_d)$, we
move to $(x_1+\delta_1,\ldots,x_d+\delta_d)$, $(x_1+2\delta_1,\ldots,
x_d+2\delta_d)$,~\dots, until either coming to a position with a nonwrapped
coordinate out of range or coming back to the original point.

A subtle technicality should be noted: When coordinates are wrapped and
|piece>0|, self-loops are possible---for example, in |board(1,0,0,0,1,1,1)|.
But self-loops never arise when |piece<0|.

@<Generate moves from |v|...@>=
for (k=1;k<=d;k++) yy[k]=xx[k]+del[k];
for (l=1;;l++) {
  @<Correct for wraparound, or |goto no_more| if off the board@>;
  if (piece<0) @<Go to |no_more| if |yy=xx|@>;
  @<Record a legal move from |xx| to |yy|@>;
  if (piece>0) goto no_more;
  for (k=1;k<=d;k++) yy[k]+=del[k];
}
no_more:

@ @<Go to |no_more|...@>=
{
  for (k=1;k<=d;k++) if (yy[k]!=xx[k]) goto unequal;
  goto no_more;
  unequal:;
}

@ @<Correct for wraparound, or |goto no_more| if off the board@>=
for (k=1;k<=d;k++) {
  if (yy[k]<0) {
    if (!wr[k]) goto no_more;
    do yy[k]+=nn[k];@+ while (yy[k]<0);
  } else if (yy[k]>=nn[k]) {
    if (!wr[k]) goto no_more;
    do yy[k]-=nn[k];@+ while (yy[k]>=nn[k]);
  }
}

@ @<Record a legal move from |xx| to |yy|@>=
for (k=2,j=yy[1];k<=d;k++) j=nn[k]*j+yy[k];
if (directed) gb_new_arc(v,new_graph->vertices+j,l);
else gb_new_edge(v,new_graph->vertices+j,l);

@* Generalized triangular boards. The subroutine call
`|simplex(n,n0,n1,n2,n3,n4,directed)|' creates a graph based on
generalized triangular or tetrahedral configurations. Such graphs are
similar in spirit to the game boards created by |board|, but they
pertain to nonrectangular grids like those in ``Chinese checkers.'' As
with |board| in the case |piece=1|, the vertices represent board positions,
and the arcs run from board positions to their nearest neighbors. Each arc has
length~1.{\tolerance=1000\par}

More formallly, the vertices can be defined as sequences of nonnegative
integers $(x_0,x_1,\ldots,x_d)$ whose sum is~|n|, where two sequences
are considered adjacent if and only if they differ by $\pm1$ in exactly
two components---equivalently, if the Euclidean distance between them
is~$\sqrt2$. When $d=2$, for example, the vertices can be visualized
as a triangular array
$$\vcenter{\halign{&\hbox to 2em{\hss$#$\hss}\cr
&&&(0,0,3)\cr
&&(0,1,2)&&(1,0,2)\cr
&(0,2,1)&&(1,1,1)&&(2,0,1)\cr
(0,3,0)&&(1,2,0)&&(2,1,0)&&(3,0,0)\cr}}$$
containing $(n+1)(n+2)/2$ elements, illustrated here when $n=3$; each vertex of
the array has up to 6 neighbors. When $d=3$ the vertices form a tetrahedral
array, a stack of triangular layers, and they can have as many as 12
neighbors. In general, a vertex in a $d$-simplicial array will have up to
$d(d+1)$ neighbors.

If the |directed| parameter is nonzero, arcs run only form vertices to neighbors
that are lexicographically greater---for example, downward or to the right
in the triangular array shown. The directed graph is therefore acyclic,
and a vertex of a $d$-simplicial array has out-degree at most $d(d+1)/2$.

@ The first parameter, |n|, specifies the sum of the coordinates
$(x_0,x_1,\ldots,x_d)$. The following parameters |n0| through |n4| specify
upper bounds on those coordinates, and they also specify the dimensionality~|d|.

If, for example, |n0|, |n1|, and |n2| are positive while |n3=0|, the
value of~|d| will be~2 and the coordinates will be constrained to
satisfy $0\le x_0\le|n0|$, $0\le x_1\le|n1|$, $0\le x_2\le|n2|$. These
upper bounds essentially lop off the corners of the triangular array.
We obtain a hexagonal board with $6m$ boundary cells by asking for
|simplex(3m,2m,2m,2m,0,0,0)|. We obtain the diamond-shaped board used
in the game of Hex [Martin Gardner, {\sl The Scientific American
Book of Mathematical Puzzles {\char`\&} Diversions\/} (Simon {\char`\&}
Schuster, 1959), Chapter~8] by calling |simplex(20,10,20,10,0,0,0)|.

In general, |simplex| determines |d| and upper bounds $(n_0,n_1,\ldots,n_d)$
in the following way: Let the first nonpositive entry of the sequence
|(n0,n1,n2,n3,n4,0)|$\null=(n_0,n_1,n_2,n_3,n_4,0)$ be~$n_k$. If $k>0$
and $n_k=0$, the value of~$d$ will be $k-1$ and the coordinates will be
bounded by the given numbers $(n_0,\ldots,n_d)$. If $k>0$ and $n_k<0$,
the value of~$d$ will be $\vert n_k\vert$ and the coordinates will be
bounded by the first $d+1$ elements of the infinite periodic sequence
$(n_0,\ldots,n_{k-1},n_0,\ldots,n_{k-1},n_0,\ldots\,)$. If $k=0$ and
$n_0<0$, the value of~$d$ will be $\vert n_0\vert$ and the coordinates
will be unbounded; equivalently, we may set $n_0=\cdots=n_d=n$. In
this case the number of vertices will be $n+d\choose d$. Finally,
if $k=0$ and $n_0=0$, we have the default case of a triangular array
with $3n$ boundary cells, exactly as if $n_0=-2$.

For example, the specification |n0=3|, |n1=-5| will produce all vertices
$(x_0,x_1,\ldots,x_5)$ such that $x_0+x_1+\cdots+x_5=n$ and $0\le x_j\le3$.
The specification |n0=1|, |n1=-d| will essentially produce all $n$-element
subsets of the $(d+1)$-element set $\{0,1,\ldots,d\}$, because we can
regard an element~$k$ as being present in the set if $x_k=1$, absent
if $x_k=0$. In that case two subsets are adjacent if and only if
they have exactly $n-1$ elements in common. 

@ @<Basic subroutines@>=
Graph *simplex(n,n0,n1,n2,n3,n4,directed)
  unsigned n; /* the constant sum of all coordinates */
  int n0,n1,n2,n3,n4; /* constraints on coordinates */
  int directed; /* should the graph be directed? */
{@+@<Vanilla local variables@>@;
  @<Normalize the simplex parameters@>;
  @<Create a graph with one vertex for each point@>;
  @<Name the points and create the arcs or edges@>;
  if (gb_alloc_trouble) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* darn, we ran out of memory somewhere back there */
  }
  return new_graph;
}

@ @<Normalize the simplex parameters@>=
if (n0==0) n0=-2;
if (n0<0) {@+k=2;@+nn[0]=n;@+d=-n0;@+n1=n2=n3=n4=0;@+}
else {
  if (n0>n) n0=n;
  nn[0]=n0;
  if (n1<=0) {@+k=2;@+d=-n1;@+n2=n3=n4=0;@+}
  else {
    if (n1>n) n1=n;
    nn[1]=n1;
    if (n2<=0) {@+k=3;@+d=-n2;@+n3=n4=0;@+}
    else {
      if (n2>n) n2=n;
      nn[2]=n2;
      if (n3<=0) {@+k=4;@+d=-n3;@+n4=0;@+}
      else {
        if (n3>n) n3=n;
        nn[3]=n3;
        if (n4<=0) {@+k=5;@+d=-n4;@+}
        else {@+if (n4>n) n4=n;
          nn[4]=n4;@+d=4;@+goto done;@+}
      }
    }
  }
}
if (d==0) {@+d=k-2;@+goto done;@+}
nn[k-1]=nn[0];
@<Compute component sizes periodically...@>;
done: /* now |nn[0]| through |nn[d]| are set up */

@ @<Create a graph with one vertex for each point@>=
@<Determine the number of feasible $(x_0,\ldots,x_d)$, and allocate the graph@>;
sprintf(new_graph->id,"simplex(%u,%d,%d,%d,%d,%d,%d)",
    n,n0,n1,n2,n3,n4,directed?1:0);
strcpy(new_graph->format,"VVZIIIZZZZZZZZ"); /* hash table will be used */

@ We determine the number of vertices by determining the coefficient of~$z^n$
in the power series
$$(1+z+\cdots+z^{n_0})(1+z+\cdots+z^{n_1})\ldots(1+z+\cdots+z^{n_d}).$$

@<Determine the number of feasible $(x_0,\ldots,x_d)$...@>=
{@+long nverts; /* the number of vertices */
  register long *coef=gb_alloc_type(n+1,@[long@],working_storage);
  if (gb_alloc_trouble) panic(no_room+1); /* can't allocate |coef| array */
  for (k=0;k<=nn[0];k++) coef[k]=1;
   /* now |coef| represents the coefficients of $1+z+\cdots+z^{n_0}$ */
  for (j=1;j<=d;j++)
    @<Multiply the power series coefficients by $1+z+\cdots+z^{n_j}$@>;
  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 */
}

@ There's a neat way to multiply by $1+z+\cdots+z^{n_j}$: We multiply
first by $1-z^{n_j+1}$, then sum the coefficients.

We want to detect impossibly large specifications without risking
integer overflow. It is easy to do this because multiplication is being
done via addition.

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

@ As we generate the vertices, it proves convenient to precompute an
array containing the numbers $y_j=n_j+\cdots+n_d$, representing the
largest possible sum of $x_j+\cdots+x_d$. We also want to maintain
the numbers $\sigma_j=n-(x_0+\cdots+x_{j-1})=x_j+\cdots+x_d$. The
conditions
$$0\le x_j\le n_j, \qquad \sigma_j-y_{j+1}\le x_j\le \sigma_j$$
are ``necessary and sufficient,'' in the sense that we can find at least
one way to complete a partial solution $(x_0,\ldots,x_k)$ to a full
solution $(x_0,\ldots,x_d)$ if and only if the conditions hold for
all $j\le k$.

There is at least one solution if and only if $n\le y_0$.

We enter the name string into a hash table, using the |hash_in|
routine of |gb_graph|, because there is no simple way to compute the
location of a vertex from its coordinates.

@<Name the points and create the arcs or edges@>=
v=new_graph->vertices;
yy[d+1]=0;@+sig[0]=n;
for (k=d;k>=0;k--) yy[k]=yy[k+1]+nn[k];
if (yy[0]>=n) {
  k=0;@+xx[0]=(yy[1]>=n? 0: n-yy[1]);
  while (1) {
    @<Complete the partial solution $(x_0,\ldots,x_k)$@>;
    @<Assign a symbolic name for $(x_0,\ldots,x_d)$ to vertex~|v|@>;
    hash_in(v); /* enter |v->name| into the hash table
                  (via utility fields |u,v|) */
    @<Create arcs or edges from previous points to~|v|@>;
    v++;
    @<Advance to the next partial solution $(x_0,\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_0,\ldots,x_k)$@>=
for (s=sig[k]-xx[k],k++;k<=d;s-=xx[k],k++) {
  sig[k]=s;
  if (s<=yy[k+1]) xx[k]=0;
  else xx[k]=s-yy[k+1];
}
if (s!=0) panic(impossible+1) /* can't happen */

@ 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_0,\ldots,x_k)$...@>=
for (k=d-1;;k--) {
  if (xx[k]<sig[k] && xx[k]<nn[k]) break;
  if (k==0) goto last;
}
xx[k]++;

@ As in the |board| routine, we represent the sequence of coordinates
$(2,0,1)$ by the string `\.{2.0.1}'.
The string won't exceed |BUF_SIZE|, because the ratio |BUF_SIZE/MAX_D| is
plenty big.

The first three coordinate values, $(x_0,x_1,x_2)$, are placed into
utility fields |x|, |y|, and |z|, so that they can be accessed immediately
if an application needs them.

@<Assign a symbolic name for $(x_0,\ldots,x_d)$ to vertex~|v|@>=
{@+register char *p=buffer; /* string pointer */
  for (k=0;k<=d;k++) {
    sprintf(p,".%d",xx[k]);
    while (*p) p++;
  }
  v->name=gb_save_string(&buffer[1]); /* omit |buffer[0]|, which is |'.'| */
  v->x.i=xx[0];@+v->y.i=xx[1];@+v->z.i=xx[2];
}

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

@<Create arcs or edges from previous points to~|v|@>=
for (j=0;j<d;j++)
  if (xx[j]) {@+register Vertex *u; /* previous vertex adjacent to |v| */
    xx[j]--;
    for (k=j+1;k<=d;k++)
      if (xx[k]<nn[k]) {@+register char *p=buffer; /* string pointer */
        xx[k]++;
        for (i=0;i<=d;i++) {
          sprintf(p,".%d",xx[i]);
          while (*p) p++;
        }
        u=hash_out(&buffer[1]);
        if (u==NULL) panic(impossible+2); /* can't happen */
        if (directed) gb_new_arc(u,v,1);
        else gb_new_edge(u,v,1);
        xx[k]--;
      }
    xx[j]++;
  }

@* Subset graphs. The subroutine call
`|subsets(n,n0,n1,n2,n3,n4,size_bits,directed)|'
creates a graph having the same vertices as
|simplex(n,n0,n1,n2,n3,n4,directed)| but with a quite different notion
of adjacency. In this we interpret a solution $(x_0,x_1,\ldots,x_d)$ to
the conditions $x_0+x_1+\cdots+x_d=n$ and $0\le x_j\le n_j$ not as a
position on a game board but as a submultiset of the multiset
$\{n_0\cdot0,n_1\cdot 1,\ldots,n_d\cdot d\}$, having $x_j$ elements
equal to~$j$. (If each $n_j=1$, the multiset is a set; this is an
important special case.) Two vertices are adjacent if and only if
their intersection has a cardinality that matches one of the bits in
|size_bits|, which is an unsigned integer. Each arc has length~1.

For example, suppose $n=3$ and |(n0,n1,n2,n3)=(2,2,2,0)|. Then the vertices
are the 3-element submultisets of $\{0,0,1,1,2,2\}$, namely
$$\{0,0,1\},\quad \{0,0,2\},\quad \{0,1,2\},\quad
\{0,2,2\},\quad \{1,1,2\},\quad \{1,2,2\},$$
which are represented by the respective vectors
$$(2,1,0),\quad (2,0,1),\quad (1,1,1),\quad
(1,0,2),\quad (0,2,1),\quad (0,1,2).$$
The intersection of multisets represented by $(x_0,x_1,\ldots,x_d)$ and
$(y_0,y_1,\ldots,y_d)$ is $$\bigl(\min(x_0,y_0),\min(x_1,y_1),\ldots,
\min(x_d,y_d)\bigr);$$ each element occurs as often as it occurs
in both multisets being intersected. If now |size_bits=3|, the
multisets will be considered adjacent whenever their
intersection contains exactly 0 or~1 elements, because $3=2^0+2^1$.
The vertices adjacent to $\{0,0,1\}$ will, for example, be
$\{0,2,2\}$ and $\{1,2,2\}$. In this case, every pair of submultisets
has a nonempty intersection, so the same graph would be obtained
if |size_bits=2|.

If |directed| is nonzero, the graph will have directed arcs, from |u|
to~|v| only if $u\le v$. Notice that the graph will have self-loops if
and only if the binary representation of |size_bits| contains the term
$2^n$, in which case there will be a loop from every vertex to itself.
(In an undirected graph, such loops are represented by two arcs.)

We define a macro |disjoint_subsets(n,k)| for the case
of $n\choose k$ vertices, adjacent if and only if they represent
disjoint $k$-subsets of an $n$-set.
One important special case is the Petersen graph, whose vertices
are the 2-element subsets of $\{0,1,2,3,4\}$, adjacent when they
are disjoint. This graph is remarkable because it contains 10 vertices,
each of degree~3, but it has no circuits of length less than~5.

@(gb_basic.h@>=
#define disjoint_subsets(n,k) @[subsets(k,1,1-n,0,0,0,1,0)@]
#define petersen() @[disjoint_subsets(5,2)@]

@ @<Basic subroutines@>=
Graph *subsets(n,n0,n1,n2,n3,n4,size_bits,directed)
  unsigned n; /* the number of elements in the multiset */
  int n0,n1,n2,n3,n4; /* multiplicities of elements */
  unsigned long size_bits; /* intersection sizes that trigger arcs */
  int directed; /* should the graph be directed? */
{@+@<Vanilla local variables@>@;
  @<Normalize the simplex parameters@>;
  @<Create a graph with one vertex for each subset@>;
  @<Name the subsets and create the arcs or edges@>;
  if (gb_alloc_trouble) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* rats, we ran out of memory somewhere back there */
  }
  return new_graph;
}

@ @<Create a graph with one vertex for each subset@>=
@<Determine the number of feasible $(x_0,\ldots,x_d)$, and allocate the graph@>;
sprintf(new_graph->id,"subsets(%u,%d,%d,%d,%d,%d,0x%lx,%d)",
    n,n0,n1,n2,n3,n4,size_bits,directed?1:0);
strcpy(new_graph->format,"ZZZIIIZZZZZZZZ"); /* hash table will not be used */

@ We generate the vertices with exactly the logic used in |simplex|.

@<Name the subsets and create the arcs or edges@>=
v=new_graph->vertices;
yy[d+1]=0;@+sig[0]=n;
for (k=d;k>=0;k--) yy[k]=yy[k+1]+nn[k];
if (yy[0]>=n) {
  k=0;@+xx[0]=(yy[1]>=n? 0: n-yy[1]);
  while (1) {
    @<Complete the partial solution $(x_0,\ldots,x_k)$@>;
    @<Assign a symbolic name for $(x_0,\ldots,x_d)$ to vertex~|v|@>;
    @<Create arcs or edges from previous subsets to~|v|@>;
    v++;
    @<Advance to the next partial solution $(x_0,\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 */

@ The only difference is that we generate the arcs or edges by brute
force, examining each pair of vertices to see if they are adjacent or not.

The code here is character-set dependent: It assumes that `\..' and null
have a character code less than `\.0', as in ASCII. It also assumes
that characters occupy exactly eight bits.
@^system dependencies@>

@d UL_BITS 8*sizeof(unsigned long) /* the number of bits in |size_bits| */

@<Create arcs or edges from previous subsets to~|v|@>=
{@+register Vertex *u;
  for (u=new_graph->vertices;u<=v;u++) {@+register char *p=u->name;
    int ss=0; /* the number of elements common to |u| and |v| */
    for (j=0;j<=d;j++,p++) {
      for (s=(*p++)-'0';*p>='0';p++) s=10*s+*p-'0'; /* |sscanf(p,"%d",&s)| */
@^character-set dependencies@>
      if (xx[j]<s) ss+=xx[j];
      else ss+=s;
    }
    if ((size_bits&(((unsigned long)1)<<ss))&& ss<UL_BITS) {
      if (directed) gb_new_arc(u,v,1);
      else gb_new_edge(u,v,1);
    }
  }
}

@* Permutation graphs. The subroutine call
`|perms(n0,n1,n2,n3,n4,max_inv,directed)|'
creates a graph whose vertices represent the permutations of a
multiset, having at most |max_inv| inversions. Two permutations are adjacent
in the graph if one is obtained from the other by interchanging two
adjacent elements. Each arc has length~1.

For example, the multiset $\{0,0,1,2\}$ has twelve permutations:
$$\vcenter{\halign{#&&\quad#\cr
0012,&0021,&0102,&0120,&0201,&0210,\cr
1002,&1020,&1200,&2001,&2010,&2100.\cr}}$$
The first of these, 0012, has two neighbors, 0021 and 0102.

The number of inversions is the number of pairs of elements $xy$ such
that $x>y$ and $x$ precedes $y$ from left to right, counting
multiplicity.  For example, 2010 has four inversions, corresponding to
$xy\in\{20,21,20,10\}$. When two permutations are adjacent, one of
them has exactly one more inversion than the other. It is not
difficult to verify that the number of inversions of a permutation is
equal to the distance in the graph from that permutation to the
lexicographically first permutation.

Parameters |n0| through |n4| specify the composition of the multiset,
just as in the |subsets| routine.
Roughly speaking, there are |n0| elements equal to~0, |n1| elements
equal to~1, and so on. The multiset $\{0,0,1,2,3,3\}$ would, for example,
be represented by |(n0,n1,n2,n3,n4)=(2,1,1,2,0)|.

Of course, we sometimes want to have multisets with more than five distinct
elements; when there are $d+1$ distinct elements, the multiset should have
$n_k$ elements equal to~$k$ and $n=n_0+n_1+\cdots+n_d$ elements in all.
The value of $d$ can be specified by making |n0=-d| (in which case
each multiplicity $n_k$ is taken to be~1); or by making |n0>0| and |n1=-d|
(in which case each multiplicity $n_k$ is taken to be equal to~|n0|); or
|n0>0|, |n1>0|, |n2=-d| (in which case the multiplicities are alternately
$(|n0|,|n1|,|n0|,|n1|,|n0|,\ldots\,)$); or |n0>0|, |n1>0|, |n2>0|, |n3=-d|,
(in which case the multiplicities are the first~|d+1| elements of the
periodic sequence $(|n0|,|n1|,|n2|,|n0|,|n1|,\ldots\,)$); or all
but |n4| are positive, while |n4=-d| (in which case the multiplicities again
are periodic).

An example like |(n0,n1,n2,n3,n4)=(1,2,3,4,-8)| is about as tricky
as you can get. It specifies the multiset $\{0,1,1,2,2,2,3,3,3,3,4,5,5,
6,6,6,7,7,7,7,8\}$.

If any of the multiplicity parameters is negative or zero, the
remaining multiplicities are ignored. For example, if |n2<=0|, the
subroutine does not look at |n3| or~|n4|.

You probably don't want to try |perms(n0,0,0,0,0,max_inv,directed)|
when |n0>0|, because a multiset with |n0| identical elements has only
one permutation.

The special case when you want all $n!$ permutations of an $n$-element set
can be obtained by calling |all_perms(n,directed)|.

@(gb_basic.h@>=
#define all_perms(n,directed) @[perms(1-n,0,0,0,0,0,directed)@]

@ If |max_inv=0|, all permutations will be considered, regardless of
the number of inversions. In that case the total number of vertices in
the graph will be the multinomial coefficient $${n\choose
n_0,n_1,\ldots,n_d}\,,\qquad n=n_0+n_1+\cdots+n_d.$$ The maximum
number of inversions in general is the number of inversions of the
lexicographically last permutation, namely ${n\choose2}-{n_0\choose2}-
{n_1\choose2}-\cdots-{n_d\choose2}=\sum_{0\le j<k\le d}n_jn_k$.

If |directed| is nonzero, the graph will contain only arcs that are directed
from permutations to their neighbors having exactly one more inversion.

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

@<Basic subroutines@>=
Graph *perms(n0,n1,n2,n3,n4,max_inv,directed)
  int n0,n1,n2,n3,n4; /* composition of the multiset */
  unsigned long max_inv; /* maximum number of inversions */
  int directed; /* should the graph be directed? */
{@+@<Vanilla local variables@>@;
  register int n; /* total number of elements in multiset */
  @<Normalize the permutation parameters@>;
  @<Create a graph with one vertex for each permutation@>;
  @<Name the permutations and create the arcs or edges@>;
  if (gb_alloc_trouble) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* shucks, we ran out of memory somewhere back there */
  }
  return new_graph;
}

@ @<Normalize the permutation parameters@>=
if (n0==0) {@+n0=1;@+n1=0;@+} /* convert the empty set into $\{0\}$ */
else if (n0<0) {@+n1=n0;@+n0=1;@+}
n=BUF_SIZE; /* this allows us to borrow code from |simplex|, already written */
@<Normalize the simplex parameters@>;
@<Determine |n| and the maximum possible number of inversions@>;

@ Here we want to set |max_inv| to the maximum possible number of
inversions, if it is zero or if it exceeds that number.

@<Determine |n| and the maximum possible number of inversions@>=
{@+register long ss; /* max inversions known to be possible */
  for (k=0,s=ss=0;k<=d;ss+=s*nn[k],s+=nn[k],k++)
    if (nn[k]>=BUF_SIZE) panic(bad_specs);
      /* too many elements in the multiset */
  if (s>=BUF_SIZE) panic(bad_specs+1); /* too many elements in the multiset */
  n=s;
  if (max_inv==0 || max_inv>ss) max_inv=ss;  
}

@ To determine the number of vertices, we sum the first |max_inv+1|
coefficients of a power series in which the coefficient of~$z^j$
is the number of permutations having $j$ inversions. It is known
[{\sl Sorting and Searching}, exercise 5.1.2--16] that this power series
is the ``$z$-multinomial coefficient''
$${n\choose n_0,\ldots,n_d}_z={n!_z\over n_0!_z\ldots n_d!_z}\,,
\qquad\hbox{where}\qquad m!_z=\prod_{k=1}^m{1-z^k\over 1-z}\,.$$

@<Create a graph with one vertex for each permutation@>=
{@+long nverts; /* the number of vertices */
  register long *coef=gb_alloc_type(max_inv+1,@[long@],working_storage);
  if (gb_alloc_trouble) panic(no_room+1); /* can't allocate |coef| array */
  coef[0]=1;
  for (j=1,s=nn[0];j<=d;s+=nn[j],j++)
    @<Multiply the power series coefficients by
       $\prod_{1\le k\le n_j}(1-z^{s+k})/(1-z^k)$@>;
  for (k=1,nverts=1;k<=max_inv;k++) {
    nverts+=coef[k];
    if (nverts>1000000000) panic(very_bad_specs); /* way too big */
  }
  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,"perms(%d,%d,%d,%d,%d,%lu,%d)",
    n0,n1,n2,n3,n4,max_inv,directed?1:0);
  strcpy(new_graph->format,"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 int 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 int *xtab,*ytab,*ztab; /* permutations and their inversions */
  int 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_alloc_type(3*n+3,@[int@],working_storage);
if (gb_alloc_trouble) { /* 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 |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 int *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,1);
    else gb_new_edge(u,v,1);
    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 having 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(n,0,0,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 n; /* the number being partitioned */
  unsigned max_parts; /* maximum number of parts */
  unsigned max_size; /* maximum size of each part */ 
  int 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_alloc_trouble) {
    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_alloc_type(n+1,@[long@],working_storage);
  if (gb_alloc_trouble) 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(%u,%u,%u,%d)",
    n,max_parts,max_size,directed?1:0);
  strcpy(new_graph->format,"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 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$, given $(x_1,\ldots,x_{j-1})$.

@<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,"+%d",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]) {@+int 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,"+%d",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,1);
  else gb_new_edge(v,u,1);
}

@* 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 distance 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 3 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,$$
if we use infix notation and name the leaves $(a,b,c,d)$ from left to right.
Here each tree is related to its two neighbors by associativity, and the
first and last trees are also related in the same way.

If |max_height=0|, it is changed to |n|, which means there is no
restriction on the height of a leaf. In this case the graph will have
exactly ${2n+1\choose n}/ (2n+1)$ vertices; furthermore, each vertex
will have exactly $n-1$ neighbors, because a rotation will be possible
just above every internal node except the root.  The graph in this
graph can also be interpreted geometrically: The vertices are in one
to one correspondence with the triangulations of a regular
$(n+2)$-gon; two triangulations are adjacent if and only if one is obtained
from the other by replacing the pair of adjacent triangles $ABC,DCB$
by the pair $ADC,BDA$.

@(gb_basic.h@>=
#define all_trees(n,directed) @[binary(n,0,directed)@]


@ The program for |binary| is very similar in structure to the program
for |parts| already considered. But the details are more exciting.

@<Basic subroutines@>=
Graph *binary(n,max_height,directed)
  unsigned n; /* the number of internal nodes */
  unsigned max_height; /* maximum height of a leaf */
  int directed; /* should the graph be directed? */
{@+@<Vanilla local variables@>@;
  if (2*n+2>BUF_SIZE) panic(bad_specs); /* |n| is too huge for us */
  if (max_height==0 || max_height>n) max_height=n;
  if (max_height>30) panic(very_bad_specs); /* more than a billion vertices */
  @<Create a graph with one vertex for each binary tree@>;
  @<Name the trees and create the arcs or edges@>;
  if (gb_alloc_trouble) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* uff da, we ran out of memory somewhere back there */
  }
  return new_graph;
}

@ The number of vertices is the coefficient of $z^n$
in the power series $G_h$, where $h=|max_height|$ and the recurrence
$$G_0=1,\qquad G_{h+1}=1+z G_h^2$$
defines $G_h$.

The coefficients of $G_5$ are $\le55308$, but the
coefficients of $G_6$ are much larger; they exceed one billion when
$28\le n\le49$, and they exceed one million when $17\le n\le 56$.
In order to avoid overflow during this calculation, we use a
special method when $h\ge6$ and $n\ge20$: In such cases, graphs
of reasonable size arise only if $n\ge 2^h-7$, and we look at the
coefficient of $z^{-(2^h-1-n)}$ in $R_h=G_h/z^{2^h-1}$, which is a
power series in $z^{-1}$ defined by the recurrence
$$R_0=1,\qquad R_{h+1}=R_h^2+z^{1-2^{h+1}}.$$

@<Create a graph with one vertex for each binary tree@>=
{@+long nverts; /* the number of vertices */
  if (n>=20 && max_height>=6) @<Compute |nverts| using the $R$ series@>@;
  else {
    nn[0]=nn[1]=1;
    for (k=2;k<=n;k++) nn[k]=0;
    for (j=2;j<=max_height;j++)
      for (k=n-1;k;k--) {
        for (s=0,i=k;i>=0;i--) s+=nn[i]*nn[k-i]; /* overflow impossible */
        nn[k+1]=s;
      }
    nverts=nn[n];
  }
  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,"binary(%u,%u,%d)",
    n,max_height,directed?1:0);
  strcpy(new_graph->format,"VVZZZZZZZZZZZZ"); /* hash table will be used */
}

@ The smallest nontrivial graph that is unilaterally disallowed by
this procedure on the grounds of size limitations occurs when |max_height=6|
and |n=20|; it has 14,162,220 vertices.

@<Compute |nverts| using the $R$ series@>=
{@+register float ss;
  d=(1<<max_height)-1-n;
  if (d>8) panic(bad_specs+1); /* too many vertices */
  if (d<0) nverts=0;
  else {
    nn[0]=nn[1]=1;
    for (k=2;k<=d;k++) nn[k]=0;
    for (j=2;j<=max_height;j++) {
      for (k=d;k;k--) {
        for (ss=0.0,i=k;i>=0;i--) ss+=((float)nn[i])*((float)nn[k-i]);
        if (ss>MAX_NNN) panic(very_bad_specs+1); /* way too big */
        for (s=0,i=k;i>=0;i--) s+=nn[i]*nn[k-i]; /* overflow impossible */
        nn[k]=s;
      }
      i=(1<<j)-1;
      if (i<=d) nn[i]++; /* add $z^{1-2^j}$ */
    }
    nverts=nn[d];
  }
}

@ We generate the trees in lexicographic order of their Polish prefix
notation, encoded in binary notation as $x_0x_1\ldots x_{2n}$, using `1'
for an internal node and `0' for a leaf. For example, the five
trees when $n=3$ are
$$1010100,\quad 1011000,\quad 1100100,\quad 1101000,\quad 1110000,$$
in lexicographic order. The algorithm for lexicographic generation maintains
three auxiliary arrays $l_j$, $y_j$, and $\sigma_j$, where
$$\sigma_j\;=\;n-j+\sum_{i=0}^{j-1}x_i\;=\;-1+\sum_{i=j}^{2n}(1-x_i)$$
is one less than the number of 0's (leaves) in $(x_j,\ldots,x_{2n})$.
The values of $l_j$ and $y_j$ are harder
to describe formally; $l_j$ is $2^{h-l}$ when $h=|max_height|$ and when
$x_j$ represents a node at level~$l$ of the tree, based on the values
of $(x_0,\ldots,x_{j-1})$. The value of $y_j$ is a binary encoding of
tree levels in which an internal node has not yet received a right child;
$y_j$ is also the maximum number of future leaves that can be produced by
previously specified internal nodes, without exceeding the maximum height.
The number of 1-bits in $y_j$ is the minimum number of future leaves,
based on previous specifications.

Therefore if $\sigma_j>y_j$, $x_j$ is forced to be~1. If $l_j=1$,
$x_j$ is forced to be~0. If the number of 1-bits of $y_j$ is equal
to $\sigma_j$, $x_j$ is forced to be~0. Otherwise $x_j$ can be
either 0 or~1, and it will be possible to complete the partial
solution $x_0\ldots x_j$ to a full Polish prefix code $x_0\ldots x_{2n}$.

For example, here are the arrays for one of the binary trees
that is generated when $n=h=3$:
$$\vcenter{\halign{$\hfil#$\quad=&&\quad#\cr
j       &0&1&2&3&4&5&6\cr
l_j     &8&4&2&2&1&1&4\cr
y_j     &0&4&6&4&5&4&0\cr
\sigma_j&3&3&3&2&2&1&0\cr
x_j     &1&1&0&1&0&0&0\cr}}$$
If $x_j=1$ and $j<2n$, we have $l_{j+1}=l_j/2$, $y_{j+1}=y_j+l_{j+1}$,
and $\sigma_{j+1}=\sigma_j$. If $x_j=0$ and $j<2n$, we have $l_{j+1}=
2^t$, $y_{j+1}=y_j-2^t$, and $\sigma_{j+1}=\sigma_j-1$, where $2^t$ is the
least power of~2 in the binary representation of~$y_j$. It is not difficult to
prove by induction that $\sigma_j<y_j+l_j$, assuming that $n<2^h$.

@<Name the trees and create the arcs or edges@>=
{@+register long *xtab,*ytab,*ltab,*stab;
  @<Initialize |xtab|, |ytab|, |ltab|, and |stab|; also set |d=2n|@>;
  v=new_graph->vertices;
  if (ltab[0]>n) {
    k=0;@+xtab[0]=n?1:0;
    while (1) {
      @<Complete the partial tree $x_0\ldots x_k$@>;
      @<Assign a Polish prefix code name to vertex~|v|@>;
      @<Create arcs or edges from |v| to previous trees@>;
      v++;
      @<Advance to the next partial tree $x_0\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 */
gb_free(working_storage);

@ @<Initialize |xtab|, |ytab|, |ltab|, and |stab|...@>=
xtab=gb_alloc_type(8*n+4,@[int@],working_storage);
if (gb_alloc_trouble) { /* no room for |xtab| */
  gb_recycle(new_graph);@+panic(no_room+2);@+}
d=n+n;
ytab=xtab+(d+1);
ltab=ytab+(d+1);
stab=ltab+(d+1);
ltab[0]=1<<max_height;
stab[0]=n; /* |ytab[0]=0| */

@ @<Complete the partial tree...@>=
for (j=k+1;j<=d;j++) {
  if (xtab[j-1]) {
    ltab[j]=ltab[j-1]>>1;
    ytab[j]=ytab[j-1]+ltab[j];
    stab[j]=stab[j-1];
  } else {
    ytab[j]=ytab[j-1]&(ytab[j-1]-1); /* remove least significant 1-bit */
    ltab[j]=ytab[j-1]-ytab[j];
    stab[j]=stab[j-1]-1;
  }
  if (stab[j]<=ytab[j]) xtab[j]=0;
  else xtab[j]=1;  /* this is the lexicographically smallest completion */
}

@ As in previous routines, 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 tree...@>=
for (k=d-1;;k--) {
  if (k<=0) goto last; /* this happens only when |n<=1| */
  if (xtab[k]) break; /* find rightmost 1 */
}
for (k--;;k--) {
  if (xtab[k]==0 && ltab[k]>1) break;
  if (k==0) goto last;
}
xtab[k]++;

@ In the |name| field, we encode internal nodes of the binary tree by
`\..' and leaves by `\.x'. Thus the five trees shown above in binary
code will be named
$$\.{.x.x.xx},\quad \.{.x..xxx},\quad \.{..xx.xx},\quad \.{..x.xxx},\quad
\.{...xxxx},$$
respectively.

@<Assign a Polish prefix...@>=
{@+register char *p=buffer; /* string pointer */
  for (k=0;k<=d;k++,p++) *p=(xtab[k]? '.': 'x');
  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 trees in lexicographic order of their
Polish prefix notation, it is relatively easy to find all pairs of trees that
are adjacent via one application of the associative law: We simply
replace a substring of the form $\..\..\alpha\beta$ by
$\..\alpha\..\beta$, when $\alpha$ and $\beta$ are Polish prefix
strings. The result comes earlier in lexicographic order, so it will
be an existing vertex unless it violates the |max_height| restriction.

@<Create arcs or edges from |v| to previous trees@>=
for (j=0;j<d;j++)
  if (xtab[j]==1 && xtab[j+1]==1) {
    for (i=j+1,s=0;s>=0;s+=(xtab[i+1]<<1)-1,i++) xtab[i]=xtab[i+1];
    xtab[i]=1;
    {@+register char *p=buffer; /* string pointer */
      register Vertex *u;
      for (k=0;k<=d;k++,p++) *p=(xtab[k]? '.': 'x');
      u=hash_out(buffer);
      if (u) {
        if (directed) gb_new_arc(v,u,1);
        else gb_new_edge(v,u,1);
      }
    }
    for (i--;i>j;i--) xtab[i+1]=xtab[i]; /* restore |xtab| */
    xtab[i+1]=1;
  }

@* Complementing and copying. We have seen how to create a wide
variety of basic graphs with the |board|, |simplex|, |subsets|,
|perms|, |parts|, and |binary| procedures. The remaining routines
of |gb_basic| are somewhat different. They transform existing
graphs into new ones, thereby presenting us with an almost
mind-boggling array of further possibilities.

The first of these transformations is perhaps the simplest: It
complements a given graph, i.e., makes vertices adjacent if and only if
they were previously non-adjacent. More precisely, the subroutine call
`|complement(g,copy,self,directed)|' returns a graph with the
same vertices as |g|, but with complemented arcs.
If |self!=0|, the new graph will have a self-loop from a vertex |v| to itself
when the original graph did not; if |self=0|, the new graph will
have no self-loops. If |directed!=0|, the new graph will have
an arc from |u| to |v| when the original graph did not; if |directed=0|,
the new graph will be undirected, and it will have an edge between |u|
and~|v| when the original graph did not. In the latter case, the original
graph should also be undirected (i.e., its arcs should come in pairs,
as described in the |gb_new_edge| routine of |gb_graph|).

If |copy!=0|, a double complement will actually be done. This means that
the new graph will essentially be a copy of the old, except that
duplicate arcs (and possibly self-loops) will be removed. Information
that may have been in the utility fields is not copied, and arc lengths
are all set to~1.

One possibly useful feature of the graphs returned by |complement| is
worth noting: The vertices adjacent to~|v|, namely the list
$$\hbox{|v->arcs->tip|,\quad |v->arcs->next->tip|,\quad
 |v->arcs->next->next->tip|,\quad \dots\thinspace,}$$
will be in strictly decreasing order (except in the case of an
undirected self-loop, when |v| itself will appear twice in succession).

@ @<Basic subroutines@>=
Graph *complement(g,copy,self,directed)
  Graph *g; /* graph to be complemented */
  int copy; /* should we double-complement? */
  int self; /* should we produce self-loops? */
  int directed; /* should the graph be directed? */
{@+@<Vanilla local variables@>@;
  register int n;
  register Vertex *u;
  register unsigned long delta; /* difference in memory addresses */
  if (g==NULL) panic(missing_operand); /* where's |g|? */
  @<Set up a graph with the vertices of |g|@>;
  sprintf(buffer,",%d,%d,%d)",copy?1:0,self?1:0,directed?1:0);
  make_compound_id(new_graph,"complement(",g,buffer);
  @<Insert complementary arcs or edges@>;
  if (gb_alloc_trouble) {
    gb_recycle(new_graph);
    panic(alloc_fault);
     /* worse luck, we ran out of memory somewhere back there */
  }
  return new_graph;
}

@ In several of the following routines it is efficient to circumvent
\Cee's normal rules for pointer arithmetic, and to use the
fact that the vertices of a graph being copied are a constant distance away
in memory from the vertices of its clone.

@d vert_offset(v,delta) ((Vertex*)(((unsigned long)v)+delta))
@^pointer hacks@>

@<Set up a graph with the vertices of |g|@>=
n=g->n;
new_graph=gb_new_graph(n);
if (new_graph==NULL)
  panic(no_room); /* out of memory before we're even started */
delta=((unsigned long)(new_graph->vertices))-((unsigned long)(g->vertices));
for (u=new_graph->vertices,v=g->vertices;v<g->vertices+n;u++,v++)
  u->name=gb_save_string(v->name);

@ A temporary utility field in the new graph is used to remember which
vertices are adjacent to a given vertex in the old one. We stamp the |tmp|
field of~|v| with a pointer to~|u| when there's an arc from |u| to~|v|.

@d tmp u.v /* utility field |u| for temporary use as a vertex pointer */

@<Insert comp...@>=
for (v=g->vertices;v<g->vertices+n;v++) {@+register Vertex *vv;
  u=vert_offset(v,delta);
           /* vertex in |new_graph| corresponding to |v| in |g| */
  {@+register Arc *a;
    for (a=v->arcs;a;a=a->next) vert_offset(a->tip,delta)->tmp=u;
  }
  if (directed) {
    for (vv=new_graph->vertices;vv<new_graph->vertices+n;vv++)
      if ((vv->tmp==u && copy) || (vv->tmp!=u && !copy))
        if (vv!=u || self) gb_new_arc(u,vv,1);
  } else {
    for (vv=(self?u:u+1);vv<new_graph->vertices+n;vv++)
      if ((vv->tmp==u && copy) || (vv->tmp!=u && !copy))
        gb_new_edge(u,vv,1);
  }
}
for (v=new_graph->vertices;v<new_graph->vertices+n;v++) v->tmp=NULL;

@* Graph union and intersection. Another simple way to get new graphs
from old ones is to take the union or intersection of their sets of arcs. The
subroutine call `|gunion(g,gg,multi,directed)|' produces a graph
with the vertices and arcs of |g| together with the
arcs of another graph~|gg|. The subroutine call `|intersection(g,gg,multi,
directed)|' produces a graph with the vertices of |g| but with only the
arcs that appear in both |g| and |gg|. In both cases we assume
that |gg| has the same vertices as |g|, in the sense that vertices
in the same relative position from the beginning of the vertex array
are considered identical. If the actual number of vertices in |gg| exceeds
the number in |g|, the extra vertices and all arcs touching them in~|gg| are
suppressed.

The input graphs are assumed to be undirected, unless the |directed|
parameter is nonzero. Peculiar results may occur if you mix directed
and undirected graphs, but the subroutines will not ``crash''
when they are asked to produce undirected output from directed input.

If |multi| is nonzero, the new graph may have multiple edges: Suppose
there are $k_1$ arcs from $u$ to $v$ in |g|, and $k_2$ in |gg|. Then
there will be $k_1+k_2$ in the union and $\min(k_1,k_2)$ in the
intersection when |multi!=0|, but at most
one in the union or intersection when |multi=0|.

The lengths of arcs are copied to the union graph when |multi!=0|;
the minimum length of multiple arcs is retained in the union when |multi=0|.

The lengths of arcs in the intersection graph are a bit trickier.
If multiple arcs occur in |g|, their minimum length, |l|, is computed. Then
we compute the maximum of |l| and the lengths of corresponding arcs
in |gg|. If |multi=0|, only the minimum of those maxima will survive.

@ @<Basic subroutines@>=
Graph *gunion(g,gg,multi,directed)
  Graph *g,*gg; /* graphs to be united */
  int multi; /* should we reproduce multiple arcs? */
  int directed; /* should the graph be directed? */
{@+@<Vanilla local variables@>@;
  register int n;
  register Vertex *u;
  register unsigned long delta,ddelta; /* differences in memory addresses */
  if (g==NULL || gg==NULL) panic(missing_operand);
    /* where are |g| and |gg|? */
  @<Set up a graph with the vertices of |g|@>;
  sprintf(buffer,",%d,%d)",multi?1:0,directed?1:0);
  make_double_compound_id(new_graph,"gunion(",g,",",gg,buffer);
  ddelta=((unsigned long)(new_graph->vertices))-((unsigned long)(gg->vertices));
  @<Insert arcs or edges present in either |g| or |gg|@>;
  if (gb_alloc_trouble) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* uh oh, we ran out of memory somewhere back there */
  }
  return new_graph;
}

@ @<Insert arcs or edges present in either |g| or |gg|@>=
for (v=g->vertices;v<g->vertices+n;v++) {@+register Arc *a;
  register Vertex *vv=vert_offset(v,delta);
      /* vertex in |new_graph| corresponding to |v| in |g| */
  register Vertex *vvv=vert_offset(vv,-ddelta);
      /* vertex in |gg| corresponding to |v| in  |g| */
  for (a=v->arcs;a;a=a->next) {
    u=vert_offset(a->tip,delta);
    @<Insert a union arc or edge from |vv| to |u|, if appropriate@>;
  }
  if (vvv<gg->vertices+gg->n) for (a=vvv->arcs;a;a=a->next) {
    u=vert_offset(a->tip,ddelta);
    if (u<new_graph->vertices+n)
      @<Insert a union arc or edge from |vv| to |u|, if appropriate@>;
  }
}
for (v=new_graph->vertices;v<new_graph->vertices+n;v++)
  v->tmp=NULL,v->tlen=NULL;

@ We use the |tmp| trick of |complement| to remember which arcs have
already been recorded from |u|, and we extend it so that we can maintain
minimum lengths. Namely, |uu->tmp| will equal |u| if and only
if we have already seen an arc from |u| to |uu|; and if so, |uu->tlen|
will be one such arc. In the undirected case, |uu->tlen| will point to the
first arc of an edge pair that touches~|u|.

The only thing slightly nontrivial here is the way we keep undirected
edges grouped in pairs. We generate a new edge from |vv| to |u| only
if |vv<=u|, and if equality holds we advance~|a| so that we don't
see the self-loop in both directions. Similar logic will be repeated 
in many of the programs below.

@d tlen z.a /* utility field |z| regarded as a pointer to an arc */

@<Insert a union arc or edge from |vv| to |u|, if appropriate@>=
{@+register Arc *b;
  if (directed) {
    if (multi || u->tmp!=vv) gb_new_arc(vv,u,a->len);
    else {
      b=u->tlen;
      if (a->len<b->len) b->len=a->len;
    }
    u->tmp=vv; /* remember that we've seen this */
    u->tlen=vv->arcs;
  } else if (u>=vv) {
    if (multi || u->tmp!=vv) gb_new_edge(vv,u,a->len);
    else {
      b=u->tlen;
      if (a->len<b->len) b->len=(b+1)->len=a->len;
    }
    u->tmp=vv;
    u->tlen=vv->arcs;
    if (u==vv && a->next==a+1) a++; /* bypass second half of self-loop */
  }
}

@ @<Basic subroutines@>=
Graph *intersection(g,gg,multi,directed)
  Graph *g,*gg; /* graphs to be intersected */
  int multi; /* should we reproduce multiple arcs? */
  int directed; /* should the graph be directed? */
{@+@<Vanilla local variables@>@;
  register int n;
  register Vertex *u;
  register unsigned long delta,ddelta; /* differences in memory addresses */
  if (g==NULL || gg==NULL) panic(no_room+1); /* where are |g| and |gg|? */
  @<Set up a graph with the vertices of |g|@>;
  sprintf(buffer,",%d,%d)",multi?1:0,directed?1:0);
  make_double_compound_id(new_graph,"intersection(",g,",",gg,buffer);
  ddelta=((unsigned long)(new_graph->vertices))-((unsigned long)(gg->vertices));
  @<Insert arcs or edges present in both |g| and |gg|@>;
  if (gb_alloc_trouble) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* whoops, we ran out of memory somewhere back there */
  }
  return new_graph;
}

@ Two more temporary utility fields are needed here.

@d mult v.i /* utility field |v|, counts multiplicity of arcs */
@d minlen w.i /* utility field |w|, records the smallest length */

@<Insert arcs or edges present in both |g| and |gg|@>=
for (v=g->vertices;v<g->vertices+n;v++) {@+register Arc *a;
  register Vertex *vv=vert_offset(v,delta);
      /* vertex in |new_graph| corresponding to |v| in |g| */
  register Vertex *vvv=vert_offset(vv,-ddelta);
      /* vertex in |gg| corresponding to |v| in  |g| */
  if (vvv>=gg->vertices+gg->n) continue;
  @<Take note of all arcs from |v|@>;
  for (a=vvv->arcs;a;a=a->next) {
    u=vert_offset(a->tip,ddelta);
    if (u>=new_graph->vertices+n) continue;
    if (u->tmp==vv) {@+int l=u->minlen;
      if (a->len>l) l=a->len; /* maximum */
      if (u->mult<0) @<Update minimum of multiple maxima@>@;
      else @<Generate a new arc or edge for the intersection,
                and reduce the multiplicity@>;
    }
  }
}
@<Clear out the temporary utility fields@>;

@ @<Generate a new arc or edge for the intersection...@>=
{
  if (directed) gb_new_arc(vv,u,l);
  else {
    if (vv<=u) gb_new_edge(vv,u,l);
    if (vv==u && a->next==a+1) a++; /* skip second half of self-loop */
  }
  if (!multi) {
    u->tlen=vv->arcs;
    u->mult=-1;
  } else if (u->mult==0) u->tmp=NULL;
  else u->mult--;
}

@ We get here if and only |multi=0| and |gg|~has more than one arc from |vv|
to~|u| and |g|~has at least one arc from |vv| to~|u|.

@<Update minimum of multiple maxima@>=
{@+register Arc *b=u->tlen; /* previous arc or edge from |vv| to |u| */
  if (l<b->len) {
    b->len=l;
    if (!directed) (b+1)->len=l;
  }
}

@ @<Take note of all arcs from |v|@>=
for (a=v->arcs;a;a=a->next) {
  u=vert_offset(a->tip,delta);
  if (u->tmp==vv) {
    u->mult++;
    if (a->len<u->minlen) u->minlen=a->len;
  } else u->tmp=vv, u->mult=0, u->minlen=a->len;
  if (u==vv && !directed && a->next==a+1) a++;
           /* skip second half of self-loop */
}

@ @<Clear out the temporary utility fields@>=
for (v=new_graph->vertices;v<new_graph->vertices+n;v++) {
  v->tmp=NULL;
  v->tlen=NULL;
  v->mult=0;
  v->minlen=0;
}

@* Line graphs. The next operation in |gb_basic|'s repertoire constructs
the so-called line graph of a given graph~$g$. The subroutine that does
this is invoked by calling `|lines(g,directed)|'.

If |directed=0|, the line graph has one vertex for each edge of~|g|;
two vertices are adjacent if and only if the corresponding edges
have a common vertex.

If |directed!=0|, the line graph has one vertex for each arc of~|g|;
there is an arc from vertex |u| to vertex |v| if and only if the
arc corresponding to~|u| ends at the vertex that begins the arc
corresponding to~|v|.

All arcs of the line graph will have length~1.

Utility fields |u.v| and |v.v| of each vertex in the line graph will point to
the vertices of |g| that define the corresponding arc or edge, and |w.a| will
point to the arc from |u.v| to |v.v| in~|g|. In the undirected case we will
have |u.v<=v.v|.

@<Basic subroutines@>=
Graph *lines(g,directed)
  Graph *g; /* graph whose lines will become vertices */
  int directed; /* should the graph be directed? */
{@+@<Vanilla local variables@>@;
  register int m; /* the number of lines */
  register Vertex *u;
  if (g==NULL) panic(no_room+1); /* where is |g|? */
  @<Set up a graph whose vertices are the lines of |g|@>;
  if (directed) @<Insert arcs of a directed line graph@>@;
  else @<Insert edges of an undirected line graph@>;
  @<Restore |g| to its pristine original condition@>;
  if (gb_alloc_trouble) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* (sigh) we ran out of memory somewhere back there */
  }
  return new_graph;
near_panic:@<Recover from potential disaster due to bad data@>;
}

@ We want to add a data structure to |g| so that the line graph can be
built efficiently. But we also want to preserve |g| so that it
exhibits no traces of occupation when |lines| has finished its
work.  To do this, we will move utility field~|v->z| temporarily into
a utility field~|u->z| of the line graph, where |u| is the first
vertex having |u->u.v==v|, whenever such a |u| exists. Then we'll
set |v->map=u|. We will then be able to find |u| when |v|
is given, and we'll be able to cover our tracks later.

In the undirected case further structure is needed. We will temporarily
change the |tip| field in the second arc of each edge pair so that
it points to the line-graph vertex that points to the first arc of the pair.

The |format| field of the graph does not indicate the fact that utility
fields |u.v|, |v.v|, and |w.a| of each vertex will be set, because those
utility fields are pointers from the new graph to the original graph.
The |save_graph| procedure does not deal with pointers between graphs.

@d map z.v /* the |z| field treated as a vertex pointer */

@<Restore |g| to its pristine original condition@>=
for (u=new_graph->vertices,v=NULL;u<new_graph->vertices+m;u++) {
  if (u->u.v!=v) {
    v=u->u.v; /* original vertex of |g| */
    v->map=u->map; /* restore original value of |v->z| */
    u->map=NULL;
  }
  if (!directed) ((u->w.a)+1)->tip=v;
}

@ Special care must be taken to avoid chaos when the user is trying to
construct the undirected line graph of a directed graph.  Otherwise we
might trash the memory, or leave the original graph in a garbled state
with pointers leading into supposedly free space.

@<Set up a graph whose vertices are the lines of |g|@>=
m=(directed? g->m: (g->m)/2);
new_graph=gb_new_graph(m);
if (new_graph==NULL)
  panic(no_room); /* out of memory before we're even started */
make_compound_id(new_graph,"lines(",g,directed? ",1)": ",0)");
u=new_graph->vertices;
for (v=g->vertices+g->n-1;v>=g->vertices;v--) {@+register Arc *a;
  register int mapped=0; /* has |v->map| been set? */
  for (a=v->arcs;a;a=a->next) {@+register Vertex *vv=a->tip;
    if (!directed) {
      if (vv<v) continue;
      if (vv>=g->vertices+g->n) goto near_panic;
            /* original graph not undirected */
    }
    @<Make |u| a vertex representing the arc |a| from |v| to |vv|@>;
    if (!mapped) {
      u->map=v->map;  /* |z.v=map| incorporates all bits of utility field |z|,
                           whatever its type */
      v->map=u;
      mapped=1;
    }
    u++;
  }
}
if (u!=new_graph->vertices+m) goto near_panic;

@ @<Recover...@>=
m=u-new_graph->vertices;
@<Restore |g| to its pristine...@>;
gb_recycle(new_graph);
panic(invalid_operand);
 /* |g| did not obey the conventions for an undirected graph */

@ The vertex names in the line graph are pairs of original vertex names,
separated by `\.{--}' when undirected, `\.{->}' when directed. If either
of the original names is horrendously long, the villainous Procrustes
chops it off arbitrarily so that it fills at most half of the name buffer.

@<Make |u| a vertex representing the arc |a| from |v| to |vv|@>=
u->u.v=v;
u->v.v=vv;
u->w.a=a;
if (!directed) {
  if (u>=new_graph->vertices+m || (a+1)->tip!=v) goto near_panic;
  if (v==vv && a->next==a+1) a++; /* skip second half of self-loop */
  else (a+1)->tip=u;
}
sprintf(buffer,"%.*s-%c%.*s",(BUF_SIZE-3)/2,v->name,@|
   directed? '>': '-',BUF_SIZE/2-1,vv->name);
u->name=gb_save_string(buffer);

@ @<Insert arcs of a directed line graph@>=
for (u=new_graph->vertices;u<new_graph->vertices+m;u++) {
  v=u->v.v;
  if (v->arcs) { /* |v->map| has been set up */
    v=v->map;
    do@+{gb_new_arc(u,v,1);
    v++;
    }@+while (v->u.v==u->v.v);
  }
}

@ An undirected line graph will contain no self-loops. It contains
multiple edges only if the original graph did; in that case, there
are two edges joining a line to each of its parallel mates, because
each mate hits both of its endpoints.

The details of this section are worthy of careful study.  We use the
fact that the first vertices of the lines occur in nonincreasing order.

@<Insert edges of an undirected line graph@>=
for (u=new_graph->vertices;u<new_graph->vertices+m;u++) {@+register Vertex *vv;
  register Arc *a;@+register int mapped=0;
  v=u->u.v; /* we look first for prior lines that touch the first vertex */
  for (vv=v->map;vv<u;vv++) gb_new_edge(u,vv,1);
  v=u->v.v; /* then we look for prior lines that touch the other one */
  for (a=v->arcs;a;a=a->next) {
    vv=a->tip;
    if (vv<u && vv>=new_graph->vertices) gb_new_edge(u,vv,1);
    else if (vv>=v && vv<g->vertices+g->n) mapped=1;
  }
  if (mapped && v>u->u.v)
    for (vv=v->map;vv->u.v==v;vv++) gb_new_edge(u,vv,1);
}

@* Graph products. Three ways have traditionally been used to define the
product of two graphs. In all three cases the vertices of the product graph
are ordered pairs $(v,v')$, where $v$ and $v'$ are vertices of the original
graphs; the difference occurs in the definition of arcs. Suppose $g$ has
$m$ arcs and $n$ vertices, while $g'$ has $m'$ arcs and $n'$ vertices. The
{\it cartesian product\/} of $g$ and~$g'$ has $mn'+m'n$ arcs, namely from
$(u,u')$ to $(v,u')$ whenever there's an arc from $u$ to $v$ in~$g$, and from
$(u,u')$ to $(u,v')$ whenever there's an arc from $u'$ to $v'$ in~$g'$.
The {\it direct product\/} has $mm'$ arcs, namely from $(u,u')$ to
$(v,v')$ in the same circumstances. The {\it strong product\/}
has both the arcs of the cartesian product and the direct product.

Notice that an undirected graph with $m$ edges has $2m$ arcs. Thus the
number of edges in the direct product of two undirected graphs is
twice the product of the number of edges in the individual graphs.
A self-loop in~$g$ will combine with an edge in~$g'$ to make
two parallel edges in the direct product.

The subroutine call `|product(g,gg,type,directed)|' produces the product
graph of one of these three types, where |type=0| for cartesian product,
|type=1| for direct product, and |type=2| for strong product.
The length of an arc in the cartesian product is copied from the length
of the original arc that it replicates; the length of an arc in the direct
product is the minimum of the two arc lengths that induce it. If |directed=0|,
the product graph will be an undirected graph, with its edges consisting
of consecutive arc pairs according to the standard GraphBase conventions,
and the input graphs should adhere to  the same conventions.

@(gb_basic.h@>=
#define cartesian 0
#define direct 1
#define strong 2

@ @<Basic subroutines@>=
Graph *product(g,gg,type,directed)
  Graph *g,*gg; /* graphs to be multiplied */
  int type; /* |cartesian|, |direct|, or |strong| */
  int directed; /* should the graph be directed? */
{@+@<Vanilla local variables@>@;
  register Vertex *u,*vv;
  register long n; /* the number of vertices in the product graph */
  if (g==NULL || gg==NULL) panic(no_room+1); /* where are |g| and |gg|? */
  @<Set up a graph with ordered pairs of vertices@>;
  if ((type&1)==0) @<Insert arcs or edges for cartesian product@>;
  if (type) @<Insert arcs or edges for direct product@>;
  if (gb_alloc_trouble) {
    gb_recycle(new_graph);
    panic(alloc_fault);
      /* @@?`$*$\#!\&, we ran out of memory somewhere back there */
  }
  return new_graph;
}

@ We must be constantly on guard against running out of memory, especially
when multiplying information.

The vertex names in the product are pairs of original vertex names, separated
by a comma.

@<Set up a graph with ordered pairs of vertices@>=
{@+float test_product=((float)(g->n))*((float)(gg->n));
  if (test_product>MAX_NNN) panic(very_bad_specs); /* way too many vertices */
}
n=(g->n)*(gg->n);
new_graph=gb_new_graph(n);
if (new_graph==NULL)
  panic(no_room); /* out of memory before we're even started */
for (u=new_graph->vertices,v=g->vertices,vv=gg->vertices;@|
    u<new_graph->vertices+n;u++) {
  sprintf(buffer,"%.*s,%.*s",BUF_SIZE/2-1,v->name,(BUF_SIZE-1)/2,vv->name);
  u->name=gb_save_string(buffer);
  if (++vv==gg->vertices+gg->n) vv=gg->vertices,v++; /* ``carry'' */
}
sprintf(buffer,",%d,%d)",(type?2:0)-(type&1),directed?1:0);
make_double_compound_id(new_graph,"product(",g,",",gg,buffer);

@ @<Insert arcs or edges for cartesian product@>=
{@+register Vertex *uu,*uuu;
  register Arc *a;
  register unsigned long delta; /* difference in memory addresses */
  delta=((unsigned long)(new_graph->vertices))-((unsigned long)(gg->vertices));
  for (u=gg->vertices;u<gg->vertices+gg->n;u++)
    for (a=u->arcs;a;a=a->next) {
      v=a->tip;
      if (!directed) {
        if (u>v) continue;
        if (u==v && a->next==a+1) a++; /* skip second half of self-loop */
      }
      for (uu=vert_offset(u,delta),vv=vert_offset(v,delta);@|
             uu<new_graph->vertices+n;uu+=gg->n,vv+=gg->n)
        if (directed) gb_new_arc(uu,vv,a->len);
        else gb_new_edge(uu,vv,a->len);      
    }
  @<Insert arcs or edges for first component of cartesian product@>;
}

@ @<Insert arcs or edges for first component...@>=
for (u=g->vertices,uu=new_graph->vertices;uu<new_graph->vertices+n;
          u++,uu+=gg->n)
  for (a=u->arcs;a;a=a->next) {
    v=a->tip;
    if (!directed) {
      if (u>v) continue;
      if (u==v && a->next==a+1) a++; /* skip second half of self-loop */
    }
    vv=new_graph->vertices+((gg->n)*(v-g->vertices));
    for (uuu=uu;uuu<uu+gg->n;uuu++,vv++)
      if (directed) gb_new_arc(uuu,vv,a->len);
      else gb_new_edge(uuu,vv,a->len);
  }

@ @<Insert arcs or edges for direct product@>=
{@+Vertex *uu;@+Arc *a;
  unsigned long delta0=
   ((unsigned long)(new_graph->vertices))-((unsigned long)(gg->vertices));
  unsigned long del=(gg->n)*sizeof(Vertex);
  register unsigned long delta,ddelta;
  for (uu=g->vertices,delta=delta0;uu<g->vertices+g->n;uu++,delta+=del)
    for (a=uu->arcs;a;a=a->next) {
      vv=a->tip;
      if (!directed) {
        if (uu>vv) continue;
        if (uu==vv && a->next==a+1) a++; /* skip second half of self-loop */
        ddelta=delta0+del*(vv-g->vertices);
        for (u=gg->vertices;u<gg->vertices+gg->n;u++) {@+register Arc *aa;
          for (aa=u->arcs;aa;aa=aa->next) {@+long length=a->len;
            if (length>aa->len) length=aa->len;
            v=aa->tip;
            if (directed)
              gb_new_arc(vert_offset(u,delta),vert_offset(v,ddelta),length);
            else gb_new_edge(vert_offset(u,delta),vert_offset(v,ddelta),length);
          }
        }
      }
    }
}

@* Induced graphs. Another important way to transform a graph is to
remove, identify, or split some of its vertices. All of these
operations are performed by the |induced| routine, which users can
invoke by calling `|induced(g,description,self,multi,directed)|'.

Each vertex |v| of |g| should first be assigned an ``induction code'' in
its field |v->ind|, which is actually utility field~|z|. The
induction code is 0~if |v| is to be eliminated; it is 1~if |v| is to be
retained; it is |k>1| if |v| is to be split into $k$ nonadjacent vertices
having the same neighbors as~|v| did; and it is |k<0| if |v| is to be
identified with all other vertices having the same value of~|k|.

For example, suppose |g| is a circuit with vertices $\{0,1,\ldots,9\}$,
where |j| is adjacent to~|k| if and only if $k=(j\pm1)\bmod10$.
If we set
$$\vcenter{\halign{\hbox{\hfil#\hfil}\cr
|0->ind=0|,\quad |1->ind=5->ind=9->ind=-1|,\quad |2->ind=3->ind=-2|,\cr
|4->ind=6->ind=8->ind=1|,\quad and |7->ind=3|,\cr}}$$
the induced graph will have vertices $\{-1,-2,4,6,7,7',7'',8\}$.
The vertices adjacent to 6, say, will be $-1$ (formerly~5), 7, $7'$,
and~$7''$. The vertices adjacent to $-1$ will be those formerly
adjacent to 1, 5, or~9, namely $-2$ (formerly~2), 4, 6, and~8. The
vertices adjacent to $-2$ will be those formerly adjacent to 2 or~3,
namely $-1$ (formerly~1), $-2$ (formerly~3), $-2$ (formerly~2), and~4.
Duplicate edges will be discarded if |multi==0|, and self-loops will
be discarded if |self==0|.

The total number of vertices in the induced graph will be the sum
of the positive |ind| fields plus the absolute value of the most
negative |ind| field. This rule implies, for example, that if at least
one vertex has |ind=-5|, the induced graph will always have a vertex $-4$,
even though no |ind| field has been set to~$-4$.

The |description| parameter is a string that will appear as part of
the name of the induced graph; if |description=0|, this string will
be empty. In the latter case, users are encouraged to assign a suitable
name to the |id| field of the induced graph, characterizing the method
by which the |ind| codes were set.

If the |directed| parameter is zero, the input graph will be assumed to
be undirected, and the output graph will be undirected.

When |multi=0|, the length of an arc that represents multiple arcs
will be the minimum of the multiple arc lengths.

@d ind z.i

@(gb_basic.h@>=
#define ind @[z.i /* utility field |z| when used to induce a graph */@]

@ Here's a simple example: To get a complete bipartite graph with
parts of sizes |n1| and |n2|, we can start with a trivial two-point
graph and split its vertices into |n1| and |n2| parts.

@<Applications...@>=
Graph *complete_bipartite(n1,n2,directed)
  unsigned n1; /* size of first part */
  unsigned n2; /* size of second part */
  int directed; /* should all arcs go from first part to second? */
{@+Graph *new_graph=board(2,0,0,0,1,0,directed);
  if (new_graph) {
    new_graph->vertices->ind=n1;
    (new_graph->vertices+1)->ind=n2;
    new_graph=induced(new_graph,0,0,0,directed);
    if (new_graph) {
      sprintf(new_graph->id,"complete_bipartite(%u,%u,%d)",@|
                              n1,n2,directed?1:0);
      mark_bipartite(new_graph,n1);
    }
  }
  return new_graph;
}

@ The |induced| routine also provides a special feature not mentioned
above: If the |ind| field of any vertex |v| is |IND_GRAPH| or greater
(where |IND_GRAPH| is a large constant, much larger than the number
of vertices that would fit in computer memory), then utility field |v->subst|
should point to a graph. A copy of the vertices of
that graph will then be substituted for |v| in the induced graph.

This feature extends the ordinary case when |v->ind>0|, which essentially
substitutes an empty graph for~|v|.

If substitution is being used to replace all of $g$'s vertices
by disjoint copies of some other graph~$g'$,
the induced graph will be somewhat similar to
a product graph. But it will not be the same as any of the three
types of output produced by |product|, because the relation between
$g$ and $g'$ is not symmetrical. Assuming that no self-loops are
present, and that graphs $(g,g')$ have respectively $(m,m')$ arcs and
$(n,n')$ vertices, the resulting of substituting $g'$ for all
vertices of~$g$ has $m'n+mn'^2$ arcs.


@d IND_GRAPH 1000000000 /* when |ind| is a billion or more, */
@d subst y.g /* we'll look at the |subst| field */

@(gb_basic.h@>=
#define IND_GRAPH 1000000000
#define subst @[y.g@]

@ For example, we can use the |IND_GRAPH| feature to create a
``wheel'' of $n$ vertices arranged cyclically, all connected to one or
more center points. In the directed case, the arcs will run from the
center(s) to a cycle; in the undirected case, the edges will join the
center(s) to a circuit.

@<Applications...@>=
Graph *wheel(n,n1,directed)
  unsigned n; /* size of the rim */
  unsigned n1; /* number of center points */
  int directed; /* should all arcs go from center to rim and around? */
{@+Graph *new_graph=board(2,0,0,0,1,0,directed); /* trivial 2-vertex graph */
  if (new_graph) {
    new_graph->vertices->ind=n1;
    (new_graph->vertices+1)->ind=IND_GRAPH;
    (new_graph->vertices+1)->subst=board(n,0,0,0,1,1,directed);
                             /* cycle or circuit */
    new_graph=induced(new_graph,0,0,0,directed);
    if (new_graph) {
      sprintf(new_graph->id,"wheel(%u,%u,%d)",@|
                              n,n1,directed?1:0);
    }
  }
  return new_graph;
}

@ @(gb_basic.h@>=
extern Graph *complete_bipartite();
extern Graph *wheel(); /* standard applications of |induced| */

@ @<Basic subroutines@>=
Graph *induced(g,description,self,multi,directed)
  Graph *g; /* graph marked for induction in its |ind| fields */
  char *description; /* string to be mentioned in |new_graph->id| */
  int self; /* should self-loops be permitted? */
  int multi; /* should multiple arcs be permitted? */
  int directed; /* should the graph be directed? */
{@+@<Vanilla local variables@>@;
  register Vertex *u;
  register long n=0; /* total number of vertices in induced graph */
  register long nn=0; /* number of negative vertices in induced graph */
  if (g==NULL) panic(missing_operand); /* where is |g|? */
  @<Set up a graph with the induced vertices@>;
  @<Insert arcs or edges for induced vertices@>;
  @<Restore |g| to its original state@>;
  if (gb_alloc_trouble) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* aargh, we ran out of memory somewhere back there */
  }
  return new_graph;
}

@ @<Set up a graph with the induced vertices@>=
@<Determine |n| and |nn|@>;
new_graph=gb_new_graph(n);
if (new_graph==NULL)
  panic(no_room); /* out of memory before we're even started */
@<Assign names to the new vertices, and create a map from |g| to |new_graph|@>;
sprintf(buffer,",%s,%d,%d,%d)",@|description?description:null_string,@|
                    self?1:0,multi?1:0,directed?1:0);
make_compound_id(new_graph,"induced(",g,buffer);

@ @<Determine |n| and |nn|@>=
for (v=g->vertices;v<g->vertices+g->n;v++)
  if (v->ind>0) {
    if (n>IND_GRAPH) panic(very_bad_specs); /* way too big */
    if (v->ind>=IND_GRAPH) {
      if (v->subst==NULL) panic(missing_operand+1);
        /* substitute graph is missing */
      n+=v->subst->n;
    } else n+=v->ind;
  } else if (v->ind<-nn) nn=-(v->ind);
if (n>IND_GRAPH || nn>IND_GRAPH) panic(very_bad_specs+1); /* gigantic */
n+=nn;

@ The negative vertices get the negative number as their name. Split vertices
get names with an optional prime appended, if the |ind| field is 2;
otherwise split vertex names are obtained by appending a colon and an index
number between |0| and~|ind-1|. The name of a vertex within a
graph |v->subst| is composed of the name of |v| followed by a
colon and the name within that graph.

We store the original |ind| field in the |mult| field of the first
corresponding vertex in the new graph, and change |ind| to point to
that vertex. That convention will make it easy
to determine the location of each vertex's clone or clones.
Of course, if the original |ind| field is zero, we leave it zero (|NULL|),
because it has no corresponding vertex in the new graph.

@<Assign names to the new vertices, and create a map from |g| to |new_graph|@>=
for (k=1,u=new_graph->vertices;k<=nn;k++,u++) {
  u->mult=-k;
  sprintf(buffer,"%d",-k);
  u->name=gb_save_string(buffer);
}  
for (v=g->vertices;v<g->vertices+g->n;v++)
  if ((k=v->ind)<0) v->map=(new_graph->vertices)-(k+1);
  else if (k>0) {
    u->mult=k;
    v->map=u;
    if (k<=2) {
      u->name=gb_save_string(v->name);
      u++;
      if (k==2) {
        sprintf(buffer,"%s'",v->name);
        u->name=gb_save_string(buffer);
        u++;
      }
    } else if (k>=IND_GRAPH) @<Make names and arcs for a substituted graph@>@;
    else for (j=0;j<k;j++,u++) {
      sprintf(buffer,"%.*s:%d",BUF_SIZE-12,v->name,j);
      u->name=gb_save_string(buffer);
    }
  }

@ @<Restore |g| to its original state@>=
for (v=g->vertices;v<g->vertices+g->n;v++)
  if (v->map) v->ind=v->map->mult;
for (v=new_graph->vertices;v<new_graph->vertices+n;v++)
  v->u.i=v->v.i=v->z.i=0; /* clear |tmp|, |mult|, |tlen| */

@ The heart of the procedure to construct an induced graph is, of course,
the part where we map the arcs of |g| into arcs of |new_graph|.

Notice that if |v| has a self-loop
in the original graph and if |v| is being split into several vertices,
it will produce arcs between different clones of itself, but it will not
produce self-loops unless |self!=0|. In an undirected graph, a loop
from a vertex to itself will not produce multiple edges among its clones,
even if |multi!=0|.

More precisely, if |v| has |k| clones |u| through |u+k-1|, an original
directed arc from |v| to~|v| will generate all $k^2$ possible arcs between
them, except that the |k| self-loops will be eliminated when
|self==0|.  An original undirected edge from |v| to~|v| will generate
$k\choose2$ edges between distinct clones, together with |k|
undirected self-loops if |self!=0|.

@<Insert arcs or edges for induced vertices@>=
for (v=g->vertices;v<g->vertices+g->n;v++) {
  u=v->map;
  if (u) {@+register Arc *a;@+register Vertex *uu,*vv;
    k=u->mult;
    if (k<0) k=1; /* |k| is the number of clones of |v| */
    else if (k>=IND_GRAPH) k=v->subst->n;
    for (;k;k--,u++) {
      if (!multi)
        @<Take note of existing edges that touch |u|@>;
      for (a=v->arcs;a;a=a->next) {
        vv=a->tip;
        uu=vv->map;
        if (uu==NULL) continue;
        j=uu->mult;
        if (j<0) j=1; /* |j| is the number of clones of |vv| */
        else if (j>=IND_GRAPH) j=vv->subst->n;
        if (!directed) {
          if (vv<v) continue;
          if (vv==v) {
            if (a->next==a+1) a++; /* skip second half of self-loop */
            j=k,uu=u; /* also skip duplicate edges generated by self-loop */
          }
        }
        @<Insert arcs or edges from vertex |u| to vertices
              |uu| through |uu+j-1|@>;
      }
    }
  }
}

@ Again we use the |tmp| and |tlen| trick of |gunion| to handle
multiple arcs. (This trick explains why the code in the previous
section tries to generate as many arcs as possible from a single
vertex |u|, before changing~|u|.)

@<Take note of existing edges that touch |u|@>=
for (a=u->arcs;a;a=a->next) {
  a->tip->tmp=u;
  if (directed || a->tip>u || a->next==a+1) a->tip->tlen=a;
  else a->tip->tlen=a+1;
}

@ @<Insert arcs or edges from vertex |u| to vertices |uu|...@>=
for (;j;j--,uu++) {
  if (u==uu && !self) continue;
  if (uu->tmp==u && !multi)
    @<Update the minimum arc length from |u| to |uu|, then |continue|@>;
  if (directed) gb_new_arc(u,uu,a->len);
  else gb_new_edge(u,uu,a->len);
  uu->tmp=u;
  uu->tlen=((directed || u<=uu)? u->arcs: uu->arcs);
}

@ @<Update the minimum arc length from |u| to |uu|, then |continue|@>=
{@+register Arc *b=uu->tlen; /* existing arc or edge from |u| to |uu| */
  if (a->len<b->len) {
    b->len=a->len; /* remember the minimum length */
    if (!directed) (b+1)->len=a->len;
  }
  continue;
}

@ We have now accumulated enough experience to finish off the one
remaining piece of program with ease.

@<Make names and arcs for a sub...@>=
{@+register Graph *gg=v->subst;
  register Vertex *vv=gg->vertices;
  register Arc *a;
  unsigned long delta=((unsigned long)u)-((unsigned long)vv);
  for (j=0;j<v->subst->n;j++,u++,vv++) {
    sprintf(buffer,"%.*s:%.*s",BUF_SIZE/2-1,v->name,(BUF_SIZE-1)/2,vv->name);
    u->name=gb_save_string(buffer);
    for (a=vv->arcs;a;a=a->next) {@+register Vertex *vvv=a->tip;
      Vertex *uu=vert_offset(vvv,delta);
      if (vvv==vv && !self) continue;
      if (uu->tmp==u && !multi) @<Update the minimum arc length...@>;
      if (!directed) {
        if (vvv<vv) continue;
        if (vvv==vv && a->next==a+1) a++; /* skip second half of self-loop */
        gb_new_edge(u,uu,a->len);
      } else gb_new_arc(u,uu,a->len);
      uu->tmp=u;
      uu->tlen=((directed || u<=uu)? u->arcs: uu->arcs);
    }
  }
}

@* Index. As usual, we close with an index that
shows where the identifiers of \\{gb\_basic} are defined and used.