gb_plane.w

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

@ @<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$