File: gb_plane.w

package info (click to toggle)
sgb 1:20030623-3
  • links: PTS
  • area: non-free
  • in suites: sarge
  • size: 1,868 kB
  • ctags: 28
  • sloc: makefile: 197; sh: 15
file content (994 lines) | stat: -rw-r--r-- 39,922 bytes parent folder | download | duplicates (4)
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
% This file is part of the Stanford GraphBase (c) Stanford University 1993
@i boilerplate.w %<< legal stuff: PLEASE READ IT BEFORE MAKING ANY CHANGES!
@i gb_types.w

\def\title{GB\_\,PLANE}

\prerequisite{GB\_\,MILES}
@* Introduction. This GraphBase module contains the |plane| subroutine,
which constructs undirected planar graphs from vertices located randomly
in a rectangle,
as well as the |plane_miles| routine, which constructs planar graphs
based on the mileage and coordinate data in \.{miles.dat}. Both
routines use a general-purpose |delaunay| subroutine,
which computes the Delaunay triangulation of a given set of points.

@d plane_miles p_miles /* abbreviation for Procrustean external linkage */

@(gb_plane.h@>=
#define plane_miles p_miles
extern Graph *plane();
extern Graph *plane_miles();
extern void delaunay();

@ The subroutine call |plane(n,x_range,y_range,extend,prob,seed)| constructs
a planar graph whose vertices have integer coordinates
uniformly distributed in the rectangle
$$\{\,(x,y)\;\mid\;0\le x<|x_range|, \;0\le y<|y_range|\,\}\,.$$
The values of |x_range| and |y_range| must be at most $2^{14}=16384$; the
latter value is the default, which is substituted if |x_range| or |y_range|
is given as zero. If |extend==0|, the graph will have |n| vertices; otherwise
it will have |n+1| vertices, where the |(n+1)|st is assigned the coordinates
$(-1,-1)$ and may be regarded as a point at~$\infty$.
Some of the |n|~finite vertices might have identical coordinates, particularly
if the point density |n/(x_range*y_range)| is not very small.

The subroutine works by first constructing the Delaunay triangulation
of the points, then discarding
each edge of the resulting graph with probability |prob/65536|. Thus,
for example, if |prob| is zero the full Delaunay triangulation will be
returned; if |prob==32768|, about half of the Delaunay edges will remain.
Each finite edge is assigned a length equal to the Euclidean distance between
points, multiplied by $2^{10}$ and
rounded to the nearest integer. If |extend!=0|, the
Delaunay triangulation will also contain edges between $\infty$ and
all points of the convex hull; such edges, if not discarded, are
assigned length $2^{28}$, otherwise known as |INFTY|.

If |extend!=0| and |prob==0|, the graph will have $n+1$ vertices and
$3(n-1)$ edges; this is the maximum number of edges that a planar graph
on $n+1$ vertices can have. In such a case the average degree of a vertex will
be $6(n-1)/(n+1)$, slightly less than~6; hence, if |prob==32768|,
the average degree of a vertex will usually be near~3.

As with all other GraphBase routines that rely on random numbers,
different values of |seed| will produce different graphs, in a
machine-independent fashion that is reproducible on many different
computers. Any |seed| value between 0 and $2^{31}-1$ is permissible.

@d INFTY 0x10000000L /* ``infinite'' length */

@(gb_plane.h@>=
#define INFTY @t\quad@> 0x10000000L

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

@d panic(c) @+{@+panic_code=c;@+gb_trouble_code=0;@+return NULL;@+}

@ Here is the overall shape of the \CEE/ file \.{gb\_plane.c}\kern.2em:

@p
#include "gb_flip.h"
 /* we will use the {\sc GB\_\,FLIP} routines for random numbers */
#include "gb_graph.h" /* we will use the {\sc GB\_\,GRAPH} data structures */
#include "gb_miles.h" /* and we might use {\sc GB\_\,MILES} for mileage data */
#include "gb_io.h"
 /* and {\sc GB\_\,MILES} uses {\sc GB\_\,IO}, which has |str_buf| */
@h@#
@<Type declarations@>@;
@<Global variables@>@;
@<Subroutines for arithmetic@>@;
@<Other subroutines@>@;
@<The |delaunay| routine@>@;
@<The |plane| routine@>@;
@<The |plane_miles| routine@>@;

@ @<The |plane| routine@>=
Graph *plane(n,x_range,y_range,extend,prob,seed)
  unsigned long n; /* number of vertices desired */
  unsigned long x_range,y_range; /* upper bounds on rectangular coordinates */
  unsigned long extend; /* should a point at infinity be included? */
  unsigned long prob; /* probability of rejecting a Delaunay edge */
  long seed; /* random number seed */
{@+Graph *new_graph; /* the graph constructed by |plane| */
  register Vertex *v; /* the current vertex of interest */
  register long k; /* the canonical all-purpose index */
  gb_init_rand(seed);
  if (x_range>16384 || y_range>16384) panic(bad_specs); /* range too large */
  if (n<2) panic(very_bad_specs); /* don't make |n| so small, you fool */
  if (x_range==0) x_range=16384; /* default */
  if (y_range==0) y_range=16384; /* default */
  @<Set up a graph with |n| uniformly distributed vertices@>;
  @<Compute the Delaunay triangulation and
    run through the Delaunay edges; reject them with probability
    |prob/65536|, otherwise append them with their Euclidean length@>;
  if (gb_trouble_code) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* oops, we ran out of memory somewhere back there */
  }
  if (extend) new_graph->n++; /* make the ``infinite'' vertex legitimate */
  return new_graph;
}

@ The coordinates are placed into utility fields |x_coord| and |y_coord|.
A random ID number is also stored in utility field~|z_coord|; this number is
used by the |delaunay| subroutine to break ties when points are equal or
collinear or cocircular. No two vertices have the same ID number.
(The header file \.{gb\_miles.h} defines |x_coord|, |y_coord|, and
|index_no| to be |x.I|, |y.I|, and |z.I| respectively.)

@d z_coord z.I

@<Set up a graph with |n| uniform...@>=
if (extend) extra_n++; /* allocate one more vertex than usual */
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,"plane(%lu,%lu,%lu,%lu,%lu,%ld)",
  n,x_range,y_range,extend,prob,seed);
strcpy(new_graph->util_types,"ZZZIIIZZZZZZZZ");
for (k=0,v=new_graph->vertices; k<n; k++,v++) {
  v->x_coord=gb_unif_rand(x_range);
  v->y_coord=gb_unif_rand(y_range);
  v->z_coord=((long)(gb_next_rand()/n))*n+k;
  sprintf(str_buf,"%ld",k);@+v->name=gb_save_string(str_buf);
}
if (extend) {
  v->name=gb_save_string("INF");
  v->x_coord=v->y_coord=v->z_coord=-1;
  extra_n--;
}

@ @(gb_plane.h@>=
#define x_coord @t\quad@> x.I
#define y_coord @t\quad@> y.I
#define z_coord @t\quad@> z.I

@* Delaunay triangulation. The Delaunay triangulation of a set of
vertices in the plane consists of all line segments $uv$ such that
there exists a circle passing through $u$ and~$v$ containing no other
vertices. Equivalently, $uv$ is a Delaunay edge if and only if the
Voronoi regions for $u$ and $v$ are adjacent; the Voronoi region of a
vertex~$u$ is the polygon with the property that all points inside it
are closer to $u$ than to any other vertex. In this sense, we can say
that Delaunay edges connect vertices with their ``neighbors.''
@^Delaunay [Delone], Boris Nikolaevich@>

The definitions in the previous paragraph assume that no two vertices are
equal, that no three vertices lie on a straight line, and that no four vertices
lie on a circle. If those nondegeneracy conditions aren't satisfied, we can
perturb the points very slightly so that the assumptions do hold.

Another way to characterize the Delaunay triangulation is to consider
what happens when we map a given set of
points onto the unit sphere via stereographic projection: Point $(x,y)$ is
mapped to
$$\bigl(2x/(r^2+1),2y/(r^2+1),(r^2-1)/(r^2+1)\bigr)\,,$$ where $r^2=x^2+y^2$.
If we now extend the configuration by adding $(0,0,1)$,
which is the limiting point on the sphere when $r$ approaches infinity,
the Delaunay edges of the original points
turn out to be edges of the polytope defined by the mapped
points. This polytope, which is the 3-dimensional convex hull of $n+1$ points
on the sphere, also has edges from $(0,0,1)$ to the mapped points
that correspond to the 2-dimensional convex hull of the original points. Under
our assumption of nondegeneracy, the faces of this polytope are all
triangles; hence its edges are said to form a triangulation.

A self-contained presentation of all the relevant theory, together with
an exposition and proof of correctness of the algorithm below, can be found
in the author's monograph {\sl Axioms and Hulls}, Lecture Notes in
Computer Science {\bf606} (Springer-Verlag, 1992).
@^Axioms and Hulls@>
For further references, see Franz Aurenhammer, {\sl ACM Computing Surveys\/
@^Aurenhammer, Franz@>
\bf23} (1991), 345--405.

@ The |delaunay| procedure, which finds the Delaunay triangulation of
a given set of vertices, is the key ingredient in \\{gb\_plane}'s
algorithms for generating planar graphs. The given vertices should
appear in a GraphBase graph~|g| whose edges, if any, are ignored by
|delaunay|. The coordinates of each vertex appear in utility fields
|x_coord| and~|y_coord|, which must be nonnegative and less than
$2^{14}=16384$. The utility field~|z_coord| must contain a unique ID
number, distinct for every vertex, so that the algorithm can break
ties in cases of degeneracy. (Note: These assumptions about the input
data are the responsibility of the calling procedure; |delaunay| does not
double-check them. If they are violated, catastrophic failure is possible.)

Instead of returning the Delaunay triangulation as a graph, |delaunay|
communicates its answer implicitly by performing the procedure call
|f(u,v)| on every pair of vertices |u| and~|v| joined by a Delaunay edge.
Here |f|~is a procedure supplied as a parameter; |u| and~|v| are either
pointers to vertices or |NULL| (i.e., \.{NULL}), where |NULL| denotes the
vertex ``$\infty$.'' As remarked above, edges run between $\infty$ and all
vertices on the convex hull of the given points. The graph of all edges,
including the infinite edges, is planar.

For example, if the vertex at infinity is being ignored, the user can
declare
$$\vcenter{\halign{#\hfil\cr
|void ins_finite(u,v)|\cr
\qquad|Vertex *u,*v;|\cr
|{@+if (u&&v)@+gb_new_edge(u,v,1L);@+}|\cr}}$$
Then the procedure call |delaunay(g,ins_finite)| will add all the finite
Delaunay edges to the current graph~|g|, giving them all length~1.

If |delaunay| is unable to allocate enough storage to do its work, it
will set |gb_trouble_code| nonzero and there will be no edges in
the triangulation.

@<The |delaunay| routine@>=
void delaunay(g,f)
  Graph *g; /* vertices in the plane */
  void @[@] (*f)(); /* procedure that absorbs the triangulated edges */
{@+@<Local variables for |delaunay|@>;@#
  @<Find the Delaunay triangulation of |g|, or return with |gb_trouble_code|
    nonzero if out of memory@>;
  @<Call |f(u,v)| for each Delaunay edge \\{uv}@>;
  gb_free(working_storage);
}

@ The procedure passed to |delaunay| will communicate with |plane| via
global variables called |gprob| and |inf_vertex|.

@<Glob...@>=
static unsigned long gprob; /* copy of the |prob| parameter */
static Vertex *inf_vertex; /* pointer to the vertex $\infty$, or |NULL| */

@ @<Compute the Delaunay triangulation and
    run through the Delaunay edges; reject them with probability
    |prob/65536|, otherwise append them with their Euclidean length@>=
gprob=prob;
if (extend) inf_vertex=new_graph->vertices+n;
else inf_vertex=NULL;
delaunay(new_graph,new_euclid_edge);

@ @<Other...@>=
static void new_euclid_edge(u,v)
  Vertex *u,*v;
{@+register long dx,dy;
  if ((gb_next_rand()>>15)>=gprob) {
    if (u) {
      if (v) {
        dx=u->x_coord-v->x_coord;
        dy=u->y_coord-v->y_coord;
        gb_new_edge(u,v,int_sqrt(dx*dx+dy*dy));
      }@+else if (inf_vertex) gb_new_edge(u,inf_vertex,INFTY);
    }@+else if (inf_vertex) gb_new_edge(inf_vertex,v,INFTY);
  }
}

@* Arithmetic. Before we lunge into the world of geometric algorithms,
let's build up some confidence by polishing off some subroutines that
will be needed to ensure correct results. We assume that |long| integers
are less than $2^{31}$.

First is a routine to calculate $s=\lfloor2^{10}\sqrt x+{1\over2}\rfloor$,
the nearest integer to $2^{10}$ times the square root of a given nonnegative
integer~|x|. If |x>0|, this
is the unique integer such that $2^{20}x-s\le s^2<2^{20}x+s$.

The following routine appears to work by magic, but the mystery goes
away when one considers the invariant relations
$$ m=\lfloor 2^{2k-21}\rfloor,\qquad
   0<y=\lfloor 2^{20-2k}x\rfloor-s^2+s\le q=2s.$$
(Exception: We might actually have $y=0$ for a short time when |q=2|.)

@<Subroutines for arith...@>=
static long int_sqrt(x)
  long x;
{@+register long y, m, q=2; long k;
  if (x<=0) return 0;
  for (k=25,m=0x20000000;x<m;k--,m>>=2) ; /* find the range */
  if (x>=m+m) y=1;
  else y=0;
  do @<Decrease |k| by 1, maintaining the invariant relations
       between |x|, |y|, |m|, and |q|@>@;
  while (k);
  return q>>1;
}

@ @<Decrease |k| by 1, maintaining the invariant relations...@>=
{
  if (x&m) y+=y+1;
  else y+=y;
  m>>=1;
  if (x&m) y+=y-q+1;
  else y+=y-q;
  q+=q;
  if (y>q)
    y-=q,q+=2;
  else if (y<=0)
    q-=2,y+=q;
  m>>=1;
  k--;
}

@ We are going to need multiple-precision arithmetic in order to
calculate certain geometric predicates properly, but it turns out
that we do not need to implement general-purpose subroutines for bignums.
It suffices to have a single special-purpose routine called
$|sign_test|(|x1|,|x2|,|x3|,\allowbreak|y1|,|y2|,|y3|)$, which
computes a single-precision integer having the same sign as the dot product
$$\hbox{|x1*y1+x2*y2+x3*y3|}$$
when we have $-2^{29}<|x1|,|x2|,|x3|<2^{29}$ and $0\le|y1|,|y2|,|y3|<2^{29}$.

@<Subroutines for arith...@>=
static long sign_test(x1,x2,x3,y1,y2,y3)
  long x1,x2,x3,y1,y2,y3;
{@+long s1,s2,s3; /* signs of individual terms */
  long a,b,c; /* components of a redundant representation of the dot product */
  register long t; /* temporary register for swapping */
  @<Determine the signs of the terms@>;
  @<If the answer is obvious, return it without further ado; otherwise,
    arrange things so that |x3*y3| has the opposite sign to |x1*y1+x2*y2|@>;
  @<Compute a redundant representation of |x1*y1+x2*y2+x3*y3|@>;
  @<Return the sign of the redundant representation@>;
}

@ @<Determine the signs of the terms@>=
if (x1==0 || y1==0) s1=0;
else {
  if (x1>0) s1=1;
  else x1=-x1,s1=-1;
}
if (x2==0 || y2==0) s2=0;
else {
  if (x2>0) s2=1;
  else x2=-x2,s2=-1;
}
if (x3==0 || y3==0) s3=0;
else {
  if (x3>0) s3=1;
  else x3=-x3,s3=-1;
}

@ The answer is obvious unless one of the terms is positive and one
of the terms is negative.

@<If the answer is obvious, return it without further ado; otherwise,
    arrange things so that |x3*y3| has the opposite sign to |x1*y1+x2*y2|@>=
if ((s1>=0 && s2>=0 && s3>=0) || (s1<=0 && s2<=0 && s3<=0))
  return (s1+s2+s3);
if (s3==0 || s3==s1) {
  t=s3;@+s3=s2;@+s2=t;
  t=x3;@+x3=x2;@+x2=t;
  t=y3;@+y3=y2;@+y2=t;
}@+else if (s3==s2) {
  t=s3;@+s3=s1;@+s1=t;
  t=x3;@+x3=x1;@+x1=t;
  t=y3;@+y3=y1;@+y1=t;
}

@ We make use of a redundant representation $2^{28}a+2^{14}b+c$, which
can be computed by brute force. (Everything is understood to be multiplied
by |-s3|.)

@<Compute a redundant...@>=
{@+register long lx,rx,ly,ry;
  lx=x1/0x4000;@+rx=x1%0x4000; /* split off the least significant 14 bits */
  ly=y1/0x4000;@+ry=y1%0x4000;
  a=lx*ly;@+b=lx*ry+ly*rx;@+c=rx*ry;
  lx=x2/0x4000;@+rx=x2%0x4000;
  ly=y2/0x4000;@+ry=y2%0x4000;
  a+=lx*ly;@+b+=lx*ry+ly*rx;@+c+=rx*ry;
  lx=x3/0x4000;@+rx=x3%0x4000;
  ly=y3/0x4000;@+ry=y3%0x4000;
  a-=lx*ly;@+b-=lx*ry+ly*rx;@+c-=rx*ry;
}

@ Here we use the fact that $\vert c\vert<2^{29}$.

@<Return the sign...@>=
if (a==0) goto ez;
if (a<0)
  a=-a,b=-b,c=-c,s3=-s3;
while (c<0) {
  a--;@+c+=0x10000000;
  if (a==0) goto ez;
}
if (b>=0) return -s3; /* the answer is clear when |a>0 && b>=0 && c>=0| */
b=-b;
a-=b/0x4000;
if (a>0) return -s3;
if (a<=-2) return s3;
return -s3*((a*0x4000-b%0x4000)*0x4000+c);
ez:@+ if (b>=0x8000) return -s3;
if (b<=-0x8000) return s3;
return -s3*(b*0x4000+c);

@*Determinants. The |delaunay| routine bases all of its decisions on
two geometric predicates, which depend on whether certain determinants
are positive or negative.

The first predicate, |ccw(u,v,w)|, is true if and only if the three points
$(u,v,w)$ have a counterclockwise orientation. This means that if we draw the
unique circle through those points, and if we travel along that circle
in the counterclockwise direction starting at~|u|, we will encounter
|v| before~|w|.

It turns out that |ccw(u,v,w)| holds if and only if the determinant
$$\left\vert\matrix{x_u&y_u&1\cr x_v&y_v&1\cr x_w&y_w&1\cr}
 \right\vert=\left\vert\matrix{x_u-x_w&y_u-y_w\cr x_v-x_w&y_v-y_w\cr}
 \right\vert$$
is positive. The evaluation must be exact; if the answer is zero, a special
tie-breaking rule must be used because the three points were collinear.
The tie-breaking rule is tricky (and necessarily so, according to the
theory in {\sl Axioms and Hulls\/}).

Integer evaluation of that determinant will not cause |long| integer
overflow, because we have assumed that all |x| and |y| coordinates lie
between 0 and~$2^{14}-1$, inclusive. In fact, we could go up to
$2^{15}-1$ without risking overflow; but the limitation to 14 bits will
be helpful when we consider a more complicated determinant below.

@<Other...@>=
static long ccw(u,v,w)
  Vertex *u,*v,*w;
{@+register long wx=w->x_coord, wy=w->y_coord; /* $x_w$, $y_w$ */
  register long det=(u->x_coord-wx)*(v->y_coord-wy)
                     -(u->y_coord-wy)*(v->x_coord-wx);
  Vertex *t;
  if (det==0) {
    det=1;
    if (u->z_coord>v->z_coord) {
           t=u;@+u=v;@+v=t;@+det=-det;
    }
    if (v->z_coord>w->z_coord) {
           t=v;@+v=w;@+w=t;@+det=-det;
    }
    if (u->z_coord>v->z_coord) {
           t=u;@+u=v;@+v=t;@+det=-det;
    }
    if (u->x_coord>v->x_coord || (u->x_coord==v->x_coord &&@|
        (u->y_coord>v->y_coord || (u->y_coord==v->y_coord &&@|
         (w->x_coord>u->x_coord ||
          (w->x_coord==u->x_coord && w->y_coord>=u->y_coord))))))
           det=-det;
  }
  return (det>0);
}

@ The other geometric predicate, |incircle(t,u,v,w)|, is true if and only
if point |t| lies outside the circle passing through |u|, |v|, and~|w|,
assuming that |ccw(u,v,w)| holds. This predicate makes us work harder, because
it is equivalent to the sign of a $4\times4$ determinant that requires
twice as much precision:
$$\left\vert\matrix{x_t&y_t&x_t^2+y_t^2&1\cr
                    x_u&y_u&x_u^2+y_u^2&1\cr
                    x_v&y_v&x_v^2+y_v^2&1\cr
                    x_w&y_w&x_w^2+y_w^2&1\cr}\right\vert=
\left\vert\matrix{x_t-x_w&y_t-y_w&(x_t-x_w)^2+(y_t-y_w)^2\cr
                  x_u-x_w&y_u-y_w&(x_u-x_w)^2+(y_u-y_w)^2\cr
                  x_v-x_w&y_v-y_w&(x_v-x_w)^2+(y_v-y_w)^2\cr}
 \right\vert\,.$$
The sign can, however, be deduced by the |sign_test| subroutine we had
the foresight to provide earlier.

@<Other...@>=
static long incircle(t,u,v,w)
  Vertex *t,*u,*v,*w;
{@+register long wx=w->x_coord, wy=w->y_coord; /* $x_w$, $y_w$ */
  long tx=t->x_coord-wx, ty=t->y_coord-wy; /* $x_t-x_w$, $y_t-y_w$ */
  long ux=u->x_coord-wx, uy=u->y_coord-wy; /* $x_u-x_w$, $y_u-y_w$ */
  long vx=v->x_coord-wx, vy=v->y_coord-wy; /* $x_v-x_w$, $y_v-y_w$ */
  register long det=sign_test(tx*uy-ty*ux,ux*vy-uy*vx,vx*ty-vy*tx,@|
                            vx*vx+vy*vy,tx*tx+ty*ty,ux*ux+uy*uy);
  Vertex *s;
  if (det==0) {
    @<Sort |(t,u,v,w)| by ID number@>;
    @<Remove incircle degeneracy@>;
  }
  return (det>0);
}

@ @<Sort...@>=
det=1;
if (t->z_coord>u->z_coord) {
   s=t;@+t=u;@+u=s;@+det=-det;
}
if (v->z_coord>w->z_coord) {
   s=v;@+v=w;@+w=s;@+det=-det;
}
if (t->z_coord>v->z_coord) {
   s=t;@+t=v;@+v=s;@+det=-det;
}
if (u->z_coord>w->z_coord) {
   s=u;@+u=w;@+w=s;@+det=-det;
}
if (u->z_coord>v->z_coord) {
   s=u;@+u=v;@+v=s;@+det=-det;
}

@ By slightly perturbing the points, we can always make them nondegenerate,
although the details are complicated. A sequence of 12 steps, involving
up to four auxiliary functions
$$\openup3\jot
\eqalign{|ff|(t,u,v,w)&=\left\vert
  \matrix{x_t-x_v&(x_t-x_w)^2+(y_t-y_w)^2-(x_v-x_w)^2-(y_v-y_w)^2\cr
          x_u-x_v&(x_u-x_w)^2+(y_u-y_w)^2-(x_v-x_w)^2-(y_v-y_w)^2\cr}
     \right\vert\,,\cr
|gg|(t,u,v,w)&=\left\vert
  \matrix{y_t-y_v&(x_t-x_w)^2+(y_t-y_w)^2-(x_v-x_w)^2-(y_v-y_w)^2\cr
          y_u-y_v&(x_u-x_w)^2+(y_u-y_w)^2-(x_v-x_w)^2-(y_v-y_w)^2\cr}
     \right\vert\,,\cr
|hh|(t,u,v,w)&=(x_u-x_t)(y_v-y_w)\,,\cr
|jj|(t,u,v,w)&=(x_u-x_v)^2+(y_u-y_w)^2-(x_t-x_v)^2-(y_t-y_w)^2\,,\cr}
$$
does the trick, as explained in {\sl Axioms and Hulls}.

@<Remove incircle degeneracy@>=
{@+long dd;
  if ((dd=ff(t,u,v,w))<0 || (dd==0 &&@|
       ((dd=gg(t,u,v,w))<0 || (dd==0 &&@|
         ((dd=ff(u,t,w,v))<0 || (dd==0 &&@|
           ((dd=gg(u,t,w,v))<0 || (dd==0 &&@|
             ((dd=ff(v,w,t,u))<0 || (dd==0 &&@|
               ((dd=gg(v,w,t,u))<0 || (dd==0 &&@|
                 ((dd=hh(t,u,v,w))<0 || (dd==0 &&@|
                   ((dd=jj(t,u,v,w))<0 || (dd==0 &&@|
                     ((dd=hh(v,t,u,w))<0 || (dd==0 &&@|
                       ((dd=jj(v,t,u,w))<0 || (dd==0 &&
                         jj(t,w,u,v)<0))))))))))))))))))))
    det=-det;
}

@ @<Subroutines for arith...@>=
static long ff(t,u,v,w)
  Vertex *t,*u,*v,*w;
{@+register long wx=w->x_coord, wy=w->y_coord; /* $x_w$, $y_w$ */
  long tx=t->x_coord-wx, ty=t->y_coord-wy; /* $x_t-x_w$, $y_t-y_w$ */
  long ux=u->x_coord-wx, uy=u->y_coord-wy; /* $x_u-x_w$, $y_u-y_w$ */
  long vx=v->x_coord-wx, vy=v->y_coord-wy; /* $x_v-x_w$, $y_v-y_w$ */
  return sign_test(ux-tx,vx-ux,tx-vx,vx*vx+vy*vy,tx*tx+ty*ty,ux*ux+uy*uy);
}
static long gg(t,u,v,w)
  Vertex *t,*u,*v,*w;
{@+register long wx=w->x_coord, wy=w->y_coord; /* $x_w$, $y_w$ */
  long tx=t->x_coord-wx, ty=t->y_coord-wy; /* $x_t-x_w$, $y_t-y_w$ */
  long ux=u->x_coord-wx, uy=u->y_coord-wy; /* $x_u-x_w$, $y_u-y_w$ */
  long vx=v->x_coord-wx, vy=v->y_coord-wy; /* $x_v-x_w$, $y_v-y_w$ */
  return sign_test(uy-ty,vy-uy,ty-vy,vx*vx+vy*vy,tx*tx+ty*ty,ux*ux+uy*uy);
}
static long hh(t,u,v,w)
  Vertex *t,*u,*v,*w;
{
  return (u->x_coord-t->x_coord)*(v->y_coord-w->y_coord);
}
static long jj(t,u,v,w)
  Vertex *t,*u,*v,*w;
{@+register long vx=v->x_coord, wy=w->y_coord;
  return (u->x_coord-vx)*(u->x_coord-vx)+(u->y_coord-wy)*(u->y_coord-wy)@|
        -(t->x_coord-vx)*(t->x_coord-vx)-(t->y_coord-wy)*(t->y_coord-wy);
}

@* Delaunay data structures. Now we have the primitive predicates
we need, and we can get on with the geometric aspects of |delaunay|.
As mentioned above, each vertex is represented by two coordinates and an
ID number, stored in the utility fields |x_coord|, |y_coord|, and~|z_coord|.

Each edge of the current triangulation is represented by two arcs
pointing in opposite directions; the two arcs are called {\sl mates}. Each
arc conceptually has a triangle on its left and a mate on its right.

An \&{arc} record differs from an |Arc|; it has three fields:
\smallskip
|vert| is the vertex this arc leads to, or |NULL| if that vertex is $\infty$;
\smallskip
|next| is the next arc having the same triangle at the left;
\smallskip
|inst| is the branch node that points to the triangle at the left, as
explained below.

\smallskip\noindent
If |p| points to an arc, then |p->next->next->next==p|, because a triangle
is bounded by three arcs. We also have |p->next->inst==p->inst|, for
all arcs~|p|.

@<Type...@>=
typedef struct a_struct {
  Vertex *vert; /* |v|, if this arc goes from |u| to |v| */
  struct a_struct *next; /* the arc from |v| that shares
         a triangle with this one */
  struct n_struct *inst; /* instruction to change
          when the triangle is modified */
} arc;

@ Storage is allocated in such a way that, if |p| and |q| point respectively
to an arc and its mate, then |p+q=&arc_block[0]+&arc_block[m-1]|, where |m| is
the total number of arc records allocated in the |arc_block| array. This
convention saves us one pointer field in each arc.

When setting |q| to the mate of |p|, we need to do the calculation
cautiously using an auxiliary register, because the constant
|&arc_block[0]+&arc_block[m-1]| might be too large to evaluate without
integer overflow on some systems.
@^pointer hacks@>

@d mate(a,b) { /* given |a|, set |b| to its mate */
  reg=max_arc-(siz_t)a;
  b=(arc*)(reg+min_arc);
}

@<Local variables for |delaunay|@>=
register siz_t reg; /* used while computing mates */
siz_t min_arc,max_arc; /* |&arc_block[0]|, |&arc_block[m-1]| */
arc *next_arc; /* the first arc record that hasn't yet been used */

@ @<Initialize the array of arcs@>=
next_arc=gb_typed_alloc(6*g->n-6,arc,working_storage);
if (next_arc==NULL) return; /* |gb_trouble_code| is nonzero */
min_arc=(siz_t)next_arc;
max_arc=(siz_t)(next_arc+(6*g->n-7));

@ @<Call |f(u,v)| for each Delaunay edge...@>=
a=(arc *)min_arc;
b=(arc *)max_arc;
for (; a<next_arc; a++,b--)
  (*f)(a->vert,b->vert);

@ The last and probably most crucial component of the data structure
is the collection of {\sl branch nodes}, which will be linked together
into a binary tree.  Given a new vertex |w|, we will ascertain what
triangle it belongs to by starting at the root of this tree and
executing a sequence of instructions, each of which has the form `if
|w| lies to the right of the straight line from |u| to~|v| then go to
$\alpha$ else go to~$\beta$', where $\alpha$ and~$\beta$ are nodes
that continue the search. This process continues until we reach a
terminal node, which says `congratulations, you're done, |w|~is in
triangle such-and-such'. The terminal node points to one of the three
arcs bounding that triangle. If a vertex of the triangle is~$\infty$,
the terminal node points to the arc whose |vert| pointer is~|NULL|.

@<Type...@>=
typedef struct n_struct {
  Vertex *u; /* first vertex, or |NULL| if this is a terminal node */
  Vertex *v; /* second vertex, or pointer to the triangle
                corresponding to a terminal node */
  struct n_struct *l; /* go here if |w| lies to the left of $uv$ */
  struct n_struct *r; /* go here if |w| lies to the right of $uv$ */
} node;

@ The search tree just described is actually a dag (a directed acyclic
graph), because it has overlapping subtrees. As the algorithm proceeds,
the dag gets bigger and bigger, since the number of triangles keeps
growing. Instructions are never deleted; we just extend the dag by
substituting new branches for nodes that once were terminal.

The expected number of nodes in this dag is $O(n)$ when there are $n$~vertices,
if we input the vertices in random order. But it can be as high as order~$n^2$
in the worst case. So our program will allocate blocks of nodes dynamically
instead of assuming a maximum size.

@d nodes_per_block 127 /* on most computers we want it $\equiv 15$ (mod 16) */
@d new_node(x)
  if (next_node==max_node) {
    x=gb_typed_alloc(nodes_per_block,node,working_storage);
    if (x==NULL) {
      gb_free(working_storage); /* release |delaunay|'s auxiliary memory */
      return; /* |gb_trouble_code| is nonzero */
    }
    next_node=x+1; max_node=x+nodes_per_block;
  }@+else x=next_node++;
@#
@d terminal_node(x,p) {@+new_node(x); /* allocate a new node */
   x->v=(Vertex*)(p); /* make it point to a given arc from the triangle */
}  /* note that |x->u==NULL|, representing a terminal node */

@<Local variables for |delaunay|@>=
node *next_node; /* the first yet-unused node slot
   in the current block of nodes */
node *max_node; /* address of nonexistent node following the current
   block of nodes */
node root_node; /* start here to locate a vertex in its triangle */
Area working_storage; /* where |delaunay| builds its triangulation */

@ The algorithm begins with a trivial triangulation that contains
only the first two vertices, together with two ``triangles'' extending
to infinity at their left and right.

@<Initialize the data structures@>=
next_node=max_node=NULL;
init_area(working_storage);
@<Initialize the array of arcs@>;
u=g->vertices;
v=u+1;
@<Make two ``triangles'' for |u|, |v|, and $\infty$@>;

@ We'll need a bunch of local variables to do elementary operations on
data structures.

@<Local variables for |delaunay|@>=
Vertex *p, *q, *r, *s, *t, *tp, *tpp, *u, *v;
arc *a,*aa,*b,*c,*d, *e;
node *x,*y,*yp,*ypp;

@ @<Make two ``triangles'' for |u|, |v|, and $\infty$@>=
root_node.u=u; root_node.v=v;
a=next_arc;
terminal_node(x,a+1);
root_node.l=x;
a->vert=v;@+a->next=a+1;@+a->inst=x;
(a+1)->next=a+2;@+(a+1)->inst=x;
 /* |(a+1)->vert=NULL|, representing $\infty$ */
(a+2)->vert=u;@+(a+2)->next=a;@+(a+2)->inst=x;
mate(a,b);
terminal_node(x,b-2);
root_node.r=x;
b->vert=u;@+b->next=b-2;@+b->inst=x;
(b-2)->next=b-1;@+(b-2)->inst=x;
 /* |(b-2)->vert=NULL|, representing $\infty$ */
(b-1)->vert=v;@+(b-1)->next=b;@+(b-1)->inst=x;
next_arc+=3;

@*Delaunay updating.
The main loop of the algorithm updates the data structure incrementally
by adding one new vertex at a time. The new vertex will always be connected
by an edge (i.e., by two arcs) to each of the vertices of the triangle that
previously enclosed it. It might also deserve to be connected to other
nearby vertices.

@<Find the Delaunay triangulation...@>=
if (g->n<2) return; /* no edges unless there are at least 2 vertices */
@<Initialize the data structures@>;
for (p=g->vertices+2;p<g->vertices+g->n;p++) {
  @<Find an arc |a| on the boundary of the triangle containing |p|@>;
  @<Divide the triangle left of |a| into three triangles surrounding |p|@>;
  @<Explore the triangles surrounding |p|, ``flipping'' their neighbors
    until all triangles that should touch |p| are found@>;
}

@ We have set up the branch nodes so that they solve the triangle location
problem.

@<Find an arc |a| on the boundary of the triangle containing |p|@>=
x=&root_node;
do@+{
  if (ccw(x->u,x->v,p))
    x = x->l;
  else x = x->r;
}@+while (x->u);
a = (arc*) x->v; /* terminal node points to the arc we want */

@ Subdividing a triangle is an easy exercise in data structure manipulation,
except that we must do something special when one of the vertices is
infinite. Let's look carefully at what needs to be done.

Suppose the triangle containing |p| has the vertices |q|, |r|, and |s|
in counterclockwise order. Let |x| be the terminal node that points to
the triangle~$\Delta qrs$. We want to change |x| so that we will be
able to locate a future point of $\Delta qrs$ within either $\Delta pqr$,
$\Delta prs$, or $\Delta psq$.

If |q|, |r|, and |s| are finite, we will change |x| and add five new nodes
as follows:
$$\vbox{\halign{\hfil#:\enspace&#\hfil\cr
$x$&if left of $rp$, go to $x''$, else go to $x'$;\cr
$x'$&if left of $sp$, go to $y$, else go to $y'$;\cr
$x''$&if left of $qp$, go to $y'$, else go to $y''$;\cr
$y$&you're in $\Delta prs$;\cr
$y'$&you're in $\Delta psq$;\cr
$y''$&you're in $\Delta pqr$.\cr}}$$

But if, say, $q=\infty$, such instructions make no sense,
because there are lines in all directions that run from $\infty$ to any point.
In such a case we use ``wedges'' instead of triangles, as explained below.

At the beginning of the following code, we have |x==a->inst|.

@<Divide the triangle left of |a| into three triangles surrounding |p|@>=
b=a->next;@+c=b->next;
q=a->vert;@+r=b->vert;@+s=c->vert;
@<Create new terminal nodes |y|, |yp|, |ypp|, and new arcs pointing to them@>;
if (q==NULL) @<Compile instructions to update convex hull@>
else {@+register node *xp;
  x->u=r;@+x->v=p;
  new_node(xp);
  xp->u=q;@+xp->v=p;@+xp->l=yp;@+xp->r=ypp; /* instruction $x''$ above */
  x->l=xp;
  new_node(xp);
  xp->u=s;@+xp->v=p;@+xp->l=y;@+xp->r=yp; /* instruction $x'$ above */
  x->r=xp;
}

@ The only subtle point here is that |q=a->vert| might be |NULL|. A terminal
node must point to the proper arc of an infinite triangle.

@<Create new terminal nodes |y|, |yp|, |ypp|, and new arcs pointing to them@>=
terminal_node(yp,a);@+terminal_node(ypp,next_arc);@+terminal_node(y,c);
c->inst=y;@+a->inst=yp;@+b->inst=ypp;
mate(next_arc,e);
a->next=e;@+b->next=e-1;@+c->next=e-2;
next_arc->vert=q;@+next_arc->next=b;@+next_arc->inst=ypp;
(next_arc+1)->vert=r;@+(next_arc+1)->next=c;@+(next_arc+1)->inst=y;
(next_arc+2)->vert=s;@+(next_arc+2)->next=a;@+(next_arc+2)->inst=yp;
e->vert=(e-1)->vert=(e-2)->vert=p;
e->next=next_arc+2;@+(e-1)->next=next_arc;@+(e-2)->next=next_arc+1;
e->inst=yp;@+(e-1)->inst=ypp;@+(e-2)->inst=y;
next_arc += 3;

@ Outside of the current convex hull, we have ``wedges'' instead of
triangles. Wedges are exterior angles whose points lie outside an
edge $rs$ of the convex hull, but not outside the next edge on the other
side of point |r|. When a new point lies in such a wedge, we have to
see if it also lies outside the edges $st$, $tu$, etc., in the
clockwise direction, in which case the convex hull loses points
$s$, $t$, etc., and we must update the new wedges accordingly.

This was the hardest part of the program to prove correct; a complete
proof can be found in {\sl Axioms and Hulls}.

@<Compile...@>=
{@+register node *xp;
  x->u=r;@+x->v=p;@+x->l=ypp;
  new_node(xp);
  xp->u=s;@+xp->v=p;@+xp->l=y;@+xp->r=yp;
  x->r=xp;
  mate(a,aa);@+d=aa->next;@+t=d->vert;
  while (t!=r && (ccw(p,s,t))) {@+register node *xpp;
    terminal_node(xpp,d);
    xp->r=d->inst;
    xp=d->inst;
    xp->u=t;@+xp->v=p;@+xp->l=xpp;@+xp->r=yp;
    flip(a,aa,d,s,NULL,t,p,xpp,yp);
    a=aa->next;@+mate(a,aa);@+d=aa->next;
    s=t;@+t=d->vert;
    yp->v=(Vertex*)a;
  }
  terminal_node(xp,d->next);
  x=d->inst;@+x->u=s;@+x->v=p;@+x->l=xp;@+x->r=yp;
  d->inst=xp;@+d->next->inst=xp;@+d->next->next->inst=xp;
  r=s; /* this value of |r| shortens the exploration step that follows */
}

@ The updating process finishes by walking around the triangles
that surround |p|, making sure that none of them are adjacent to
triangles containing |p| in their circumcircle. (Such triangles are
no longer in the Delaunay triangulation, by definition.)

@<Explore...@>=
while(1) {
  mate(c,d);@+e=d->next;
  t=d->vert;@+tp=c->vert;@+tpp=e->vert;
  if (tpp && incircle(tpp,tp,t,p)) { /* triangle $tt''t'$ no longer Delaunay */
    register node *xp, *xpp;
    terminal_node(xp,e);
    terminal_node(xpp,d);
    x=c->inst;@+x->u=tpp;@+x->v=p;@+x->l=xp;@+x->r=xpp;
    x=d->inst;@+x->u=tpp;@+x->v=p;@+x->l=xp;@+x->r=xpp;
    flip(c,d,e,t,tp,tpp,p,xp,xpp);
    c=e;
  }
  else if (tp==r) break;
  else {
    mate(c->next,aa);
    c=aa->next;
  }
}

@ Here |d| is the mate of |c|, |e=d->next|, |t=d->vert|, |tp=c->vert|,
and |tpp=e->vert|. The triangles $\Delta tt'p$ and $\Delta t'tt''$ to the
left and right of arc~|c| are being replaced in the current triangulation
by $\Delta ptt''$ and $\Delta t''t'p$, corresponding to terminal nodes
|xp| and |xpp|. (The values of |t| and |tp| are not actually used, so
some optimization is possible.)

@<Other...@>=
static void flip(c,d,e,t,tp,tpp,p,xp,xpp)
  arc *c,*d,*e;
  Vertex *t,*tp,*tpp,*p;
  node *xp,*xpp;
{@+register arc *ep=e->next, *cp=c->next, *cpp=cp->next;
  e->next=c;@+c->next=cpp;@+cpp->next=e;
  e->inst=c->inst=cpp->inst=xp;
  c->vert=p;
  d->next=ep;@+ep->next=cp;@+cp->next=d;
  d->inst=ep->inst=cp->inst=xpp;
  d->vert=tpp;
}

@*Use of mileage data. The |delaunay| routine is now complete, and the
only missing piece of code is the promised routine that generates
planar graphs based on data from the real world.

The subroutine call {\advance\thinmuskip 0mu plus 2mu
|plane_miles(n,north_weight,west_weight,pop_weight, extend,prob,seed)|}
will construct a planar graph with min$(128,n)$ vertices, where the
vertices are exactly the same as the cities produced by the subroutine call
|miles(n,north_weight,west_weight, pop_weight,0,0,seed)|. (As
explained in module {\sc GB\_\,MILES}, the weight parameters |north_weight|,
|west_weight|, and |pop_weight| are used to rank the cities by
location and/or population.)  The edges of the new graph are obtained
by first constructing the Delaunay triangulation of those cities,
based on a simple projection onto the plane using their latitude and
longitude, then discarding each Delaunay edge with probability
|prob/65536|. The length of each surviving edge is the same as the
mileage between cities that would appear in the complete graph
produced by |miles|.

If |extend!=0|, an additional vertex representing $\infty$ is also
included. The Delaunay triangulation includes edges of length |INFTY|
connecting this vertex with all cities on the convex hull; these edges,
like the others, are subject to being discarded with probability |prob/65536|.
(See the description of |plane| for further comments about using
|prob| to control the sparseness of the graph.)

The weight parameters must satisfy
$$ \vert|north_weight|\vert\le100{,}000,\quad
   \vert|west_weight|\vert\le100{,}000,\quad
   \vert|pop_weight|\vert\le100.$$
Vertices of the graph will appear in order of decreasing weight.
The |seed| parameter defines the pseudo-random numbers used wherever
a ``random'' choice between equal-weight vertices needs to be made,
or when deciding whether to discard a Delaunay edge.

@<The |plane_miles| routine@>=
Graph *plane_miles(n,north_weight,west_weight,pop_weight,extend,prob,seed)
  unsigned long n; /* number of vertices desired */
  long north_weight; /* coefficient of latitude in the weight function */
  long west_weight; /* coefficient of longitude in the weight function */
  long pop_weight; /* coefficient of population in the weight function */
  unsigned long extend; /* should a point at infinity be included? */
  unsigned long prob; /* probability of rejecting a Delaunay edge */
  long seed; /* random number seed */
{@+Graph *new_graph; /* the graph constructed by |plane_miles| */
  @<Use |miles| to set up the vertices of a graph@>;
  @<Compute the Delaunay triangulation and
    run through the Delaunay edges; reject them with probability
    |prob/65536|, otherwise append them with the road length in miles@>;
  if (gb_trouble_code) {
    gb_recycle(new_graph);
    panic(alloc_fault); /* oops, we ran out of memory somewhere back there */
  }
  gb_free(new_graph->aux_data); /* recycle special memory used by |miles| */
  if (extend) new_graph->n++; /* make the ``infinite'' vertex legitimate */
  return new_graph;
}

@ By setting the |max_distance| parameter to~1, we cause |miles|
to produce a graph having the desired vertices but no edges.
The vertices of this graph will have appropriate coordinate fields
|x_coord|, |y_coord|, and~|z_coord|.

@<Use |miles|...@>=
if (extend) extra_n++; /* allocate one more vertex than usual */
if (n==0 || n>MAX_N) n=MAX_N; /* compute true number of vertices */
new_graph=miles(n,north_weight,west_weight,pop_weight,1L,0L,seed);
if (new_graph==NULL) return NULL; /* |panic_code| has been set by |miles| */
sprintf(new_graph->id,"plane_miles(%lu,%ld,%ld,%ld,%lu,%lu,%ld)",
  n,north_weight,west_weight,pop_weight,extend,prob,seed);
if (extend) extra_n--; /* restore |extra_n| to its previous value */

@ @<Compute the Delaunay triangulation and
    run through the Delaunay edges; reject them with probability
    |prob/65536|, otherwise append them with the road length in miles@>=
gprob=prob;
if (extend) {
  inf_vertex=new_graph->vertices+new_graph->n;
  inf_vertex->name=gb_save_string("INF");
  inf_vertex->x_coord=inf_vertex->y_coord=inf_vertex->z_coord= -1;
}@+else inf_vertex=NULL;
delaunay(new_graph,new_mile_edge);

@ The mileages will all have been negated by |miles|, so we make them
positive again.

@<Other...@>=
static void new_mile_edge(u,v)
  Vertex *u,*v;
{
  if ((gb_next_rand()>>15)>=gprob) {
    if (u) {
      if (v) gb_new_edge(u,v,-miles_distance(u,v));
      else if (inf_vertex) gb_new_edge(u,inf_vertex,INFTY);
    }@+else if (inf_vertex) gb_new_edge(inf_vertex,v,INFTY);
  }
}

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