unit geompack;

interface

uses
  arrtypes, Math, SysUtils;


function diaedg(x0: Real; y0: Real; x1: Real; y1: Real; x2: Real; y2: Real; x3: Real; y3: Real): Integer;

function dtris2(point_num: Integer; var point_xy: TRealArray; var tri_num: Integer; var tri_vert: TIntegerArray; var tri_nabe: TIntegerArray): Integer;

function i4_modp(i: Integer; j: Integer): Integer;

function i4_sign(i: Integer): Integer;

function i4_wrap(ival: Integer; ilo: Integer; ihi: Integer): Integer;

procedure i4vec_heap_d(n: Integer; var a: TIntegerArray);

function i4vec_indicator(n: Integer): TIntegerArray;

procedure i4vec_sort_heap_a(n: Integer; var a: TIntegerArray);

procedure i4vec_sorted_unique(n: Integer; var a: TIntegerArray; var nuniq: Integer);

function lrline(xu: Real; yu: Real; xv1: Real; yv1: Real; xv2: Real; yv2: Real; dv: Real): Integer;

function perm_check(n: Integer; var p: TIntegerArray): Boolean;

procedure perm_inv(n: Integer; var p: TIntegerArray);

function points_delaunay_naive_2d(n: Integer; var p: TRealArray; var ntri: Integer): TIntegerArray;

function r8_epsilon(): Real;

procedure r82vec_part_quick_a(from: Integer; n: Integer; var a: TRealArray; var l: Integer; var r: Integer);

procedure r82vec_permute(n: Integer; var a: TRealArray; var p: TIntegerArray);

function r82vec_sort_heap_index_a(n: Integer; var a: TRealArray): TIntegerArray;

procedure r82vec_sort_quick_a(n: Integer; var a: TRealArray);

procedure r8mat_uniform(m: Integer; n: Integer; b: Real; c: Real; var seed: Integer; var r: TRealArray);

function r8vec_eq(from1: Integer; from2: Integer; n: Integer; var a1: TRealArray; var a2: TRealArray): Boolean;

function r8vec_gt(from1: Integer; from2: Integer; n: Integer; var a1: TRealArray; var a2: TRealArray): Boolean;

function r8vec_lt(from1: Integer; from2: Integer; n: Integer; var a1: TRealArray; var a2: TRealArray): Boolean;

procedure r8vec_swap(from1: Integer; from2: Integer; n: Integer; var a1: TRealArray; var a2: TRealArray);

function swapec(i: Integer; var top: Integer; var btri: Integer; var bedg: Integer; point_num: Integer; var point_xy: TRealArray; tri_num: Integer; var tri_vert: TIntegerArray; var tri_nabe: TIntegerArray; var stack: TIntegerArray): Integer;

function triangle_circumcenter_2d(var t: TRealArray): TRealArray;

procedure vbedg(x: Real; y: Real; point_num: Integer; var point_xy: TRealArray; tri_num: Integer; var tri_vert: TIntegerArray; var tri_nabe: TIntegerArray; var ltri: Integer; var ledg: Integer; var rtri: Integer; var redg: Integer);



implementation


//****************************************************************************80

function diaedg(x0: Real; y0: Real; x1: Real; y1: Real; x2: Real; y2: Real; x3: Real; y3: Real): Integer;

//****************************************************************************80
//
//  Purpose:
//
//    DIAEDG chooses a diagonal edge.
//
//  Discussion:
//
//    The routine determines whether 0--2 or 1--3 is the diagonal edge
//    that should be chosen, based on the circumcircle criterion, where
//(X0,Y0),(X1,Y1),(X2,Y2),(X3,Y3) are the vertices of a simple
//    quadrilateral in counterclockwise order.
//
//  Modified:
//
//    28 August 2003
//
//  Author:
//
//    Barry Joe,
//    Department of Computing Science,
//    University of Alberta,
//    Edmonton, Alberta, Canada  T6G 2H1
//
//  Reference:
//
//    Barry Joe,
//    GEOMPACK - a software package for the generation of meshes
//    using geometric algorithms,
//    Advances in Engineering Software,
//    Volume 13, pages 325-331, 1991.
//
//  Parameters:
//
//    Input, Real X0, Y0, X1, Y1, X2, Y2, X3, Y3, the coordinates of the
//    vertices of a quadrilateral, given in counter clockwise order.
//
//    Output, Integer DIAEDG, chooses a diagonal:
//    +1, if diagonal edge 02 is chosen;
//    -1, if diagonal edge 13 is chosen;
//     0, if the four vertices are cocircular.
//
var
  ca: Real;
  cb: Real;
  dx10: Real;
  dx12: Real;
  dx30: Real;
  dx32: Real;
  dy10: Real;
  dy12: Real;
  dy30: Real;
  dy32: Real;
  s: Real;
  tol: Real;
  tola: Real;
  tolb: Real;
  value: Integer;
begin
  tol:= 100.0 * r8_epsilon();

  dx10:= x1 - x0;
  dy10:= y1 - y0;
  dx12:= x1 - x2;
  dy12:= y1 - y2;
  dx30:= x3 - x0;
  dy30:= y3 - y0;
  dx32:= x3 - x2;
  dy32:= y3 - y2;

  tola:= tol * Max(Abs(dx10), Max(Abs(dy10), Max(Abs(dx30), Abs(dy30))));

  tolb:= tol * Max(Abs(dx12), Max(Abs(dy12), Max(Abs(dx32), Abs(dy32))));

  ca:= dx10 * dx30 + dy10 * dy30;
  cb:= dx12 * dx32 + dy12 * dy32;

  if ((tola < ca) and(tolb < cb)) then begin
    value:= -1;
  end
  else if ((ca < -tola) and(cb < -tolb)) then begin
    value:= 1;
  end
  else begin
    tola:= Max(tola, tolb);
    s:= (dx10 * dy30 - dx30 * dy10) * cb +(dx32 * dy12 - dx12 * dy32) * ca;

    if (tola < s) then begin
      value:= -1;
    end
    else if (s < -tola) then begin
      value:= 1;
    end
    else begin
      value:= 0;
    end;
  end;

  Result:= value;
end;
//****************************************************************************80

function dtris2(point_num: Integer; var point_xy: TRealArray; var tri_num: Integer; var tri_vert: TIntegerArray; var tri_nabe: TIntegerArray): Integer;

//****************************************************************************80
//
//  Purpose:
//
//    DTRIS2 constructs a Delaunay triangulation of 2D vertices.
//
//  Discussion:
//
//    The routine constructs the Delaunay triangulation of a set of 2D vertices
//    using an incremental approach and diagonal edge swaps.  Vertices are
//    first sorted in lexicographically increasing(X,Y) order, and
//    then are inserted one at a time from outside the convex hull.
//
//  Modified:
//
//    15 January 2004
//
//  Author:
//
//    Barry Joe,
//    Department of Computing Science,
//    University of Alberta,
//    Edmonton, Alberta, Canada  T6G 2H1
//
//  Reference:
//
//    Barry Joe,
//    GEOMPACK - a software package for the generation of meshes
//    using geometric algorithms,
//    Advances in Engineering Software,
//    Volume 13, pages 325-331, 1991.
//
//  Parameters:
//
//    Input, Integer POINT_NUM, the number of vertices.
//
//    Input/output, Real POINT_XY[POINT_NUM*2], the coordinates of the vertices.
//    On output, the vertices have been sorted into dictionary order.
//
//    Output, Integer *TRI_NUM, the number of triangles in the triangulation;
//    TRI_NUM is equal to 2*POINT_NUM - NB - 2, where NB is the number
//    of boundary vertices.
//
//    Output, Integer TRI_VERT[TRI_NUM*3], the nodes that make up each triangle.
//    The elements are indices of POINT_XY.  The vertices of the triangles are
//    in counter clockwise order.
//
//    Output, Integer TRI_NABE[TRI_NUM*3], the triangle neighbor list.
//    Positive elements are indices of TIL; negative elements are used for links
//    of a counter clockwise linked list of boundary edges; LINK:= -(3*I + J-1)
//    where I, J:= triangle, edge index; TRI_NABE[I,J] refers to
//    the neighbor along edge from vertex J to J+1(mod 3).
//
//    Output, Integer DTRIS2, is 0 for no error.
var
  cmax: Real;
  e: Integer;
  error: Integer;
  i: Integer;
  indx: TIntegerArray;
  j: Integer;
  k: Integer;
  l: Integer;
  ledg: Integer;
  lr: Integer;
  ltri: Integer;
  m: Integer;
  m1: Integer;
  m2: Integer;
  n: Integer;
  redg: Integer;
  rtri: Integer;
  stack: TIntegerArray;
  t: Integer;
  tol: Real;
  top: Integer;
begin
  SetLength(stack, point_num);

  tol:= 100.0 * r8_epsilon();
//
//  Sort the vertices by increasing(x,y).
//
  indx:= r82vec_sort_heap_index_a(point_num, point_xy);

  r82vec_permute(point_num, point_xy, indx);
//
//  Make sure that the data points are "reasonably" distinct.
//
  m1:= 1;

  for i:= 2 to point_num do begin
    m:= m1;
    m1:= i;

    k:= -1;

    for j:= 0 to 1 do begin
      cmax:= Max(Abs(point_xy[2*(m-1)+j]), Abs(point_xy[2*(m1-1)+j]));

      if (tol *(cmax + 1.0) < Abs(point_xy[2*(m-1)+j] - point_xy[2*(m1-1)+j])) then begin
        k:= j;
        break;
      end;
    end;

    if (k = -1) then begin
      raise Exception.CreateFmt('DTRIS2 - Fatal error #%d! Fails for point number I=%d. M=%d, M1=%d, X,Y(M)=(%f,%f), X,Y(M1)=(%f,%f).',[224,i,m,m1,point_xy[2*(m-1)+0],point_xy[2*(m-1)+1],point_xy[2*(m1-1)+0],point_xy[2*(m1-1)+1]]);
     end;
   end;
//
//  Starting from points M1 and M2, search for a third point M that
//  makes a "healthy" triangle(M1,M2,M)
//
  m1:= 1;
  m2:= 2;
  j:= 3;
  m:= j;
  lr:= -1;

  while(true) do begin
    if (point_num < j) then begin
      raise Exception.CreateFmt('DTRIS2 - Fatal error #%d!',[225]);
    end;

    m:= j;

    lr:= lrline(point_xy[2*(m-1)+0], point_xy[2*(m-1)+1], point_xy[2*(m1-1)+0], point_xy[2*(m1-1)+1], point_xy[2*(m2-1)+0], point_xy[2*(m2-1)+1], 0.0);

    if (lr <> 0) then begin
      break;
    end;

    inc(j);
  end;
//
//  Set up the triangle information for(M1,M2,M), and for any other
//  triangles you created because points were collinear with M1, M2.
//
  tri_num:= j - 2;

  if (lr = -1) then begin
    tri_vert[3*0+0]:= m1;
    tri_vert[3*0+1]:= m2;
    tri_vert[3*0+2]:= m;
    tri_nabe[3*0+2]:= -3;

    for i:= 2 to tri_num do begin
      m1:= m2;
      m2:= i+1;
      tri_vert[3*(i-1)+0]:= m1;
      tri_vert[3*(i-1)+1]:= m2;
      tri_vert[3*(i-1)+2]:= m;
      tri_nabe[3*(i-1)+0]:= -3 * i;
      tri_nabe[3*(i-1)+1]:= i;
      tri_nabe[3*(i-1)+2]:= i - 1;
    end;

    tri_nabe[3*(tri_num-1)+0]:= -3 *(tri_num) - 1;
    tri_nabe[3*(tri_num-1)+1]:= -5;
    ledg:= 2;
    ltri:= tri_num;
  end
  else begin
    tri_vert[3*0+0]:= m2;
    tri_vert[3*0+1]:= m1;
    tri_vert[3*0+2]:= m;
    tri_nabe[3*0+0]:= -4;

    for i:= 2 to tri_num do begin
      m1:= m2;
      m2:= i+1;
      tri_vert[3*(i-1)+0]:= m2;
      tri_vert[3*(i-1)+1]:= m1;
      tri_vert[3*(i-1)+2]:= m;
      tri_nabe[3*(i-2)+2]:= i;
      tri_nabe[3*(i-1)+0]:= -3 * i - 3;
      tri_nabe[3*(i-1)+1]:= i - 1;
    end;

    tri_nabe[3*(tri_num-1)+2]:= -3 *(tri_num);
    tri_nabe[3*0+1]:= -3 *(tri_num) - 2;
    ledg:= 2;
    ltri:= 1;
  end;
//
//  Insert the vertices one at a time from outside the convex hull,
//  determine visible boundary edges, and apply diagonal edge swaps until
//  Delaunay triangulation of vertices(so far) is obtained.
//
  top:= 0;

  for i:= j+1 to point_num do begin
    m:= i;
    m1:= tri_vert[3*(ltri-1)+ledg-1];

    if (ledg <= 2) then begin
      m2:= tri_vert[3*(ltri-1)+ledg];
    end
    else begin
      m2:= tri_vert[3*(ltri-1)+0];
    end;

    lr:= lrline(point_xy[2*(m-1)+0], point_xy[2*(m-1)+1], point_xy[2*(m1-1)+0], point_xy[2*(m1-1)+1], point_xy[2*(m2-1)+0], point_xy[2*(m2-1)+1], 0.0);

    if (0 < lr) then begin
      rtri:= ltri;
      redg:= ledg;
      ltri:= 0;
    end
    else begin
      l:= -tri_nabe[3*(ltri-1)+ledg-1];
      rtri:= l div 3;
      redg:= (l mod 3) + 1;
    end;

    vbedg(point_xy[2*(m-1)+0], point_xy[2*(m-1)+1], point_num, point_xy, tri_num, tri_vert, tri_nabe, ltri, ledg, rtri, redg);

    n:= tri_num + 1;
    l:= -tri_nabe[3*(ltri-1)+ledg-1];

    while (true) do begin
      t:= l div 3;
      e:= (l mod 3) + 1;
      l:= -tri_nabe[3*(t-1)+e-1];
      m2:= tri_vert[3*(t-1)+e-1];

      if (e <= 2) then begin
        m1:= tri_vert[3*(t-1)+e];
      end
      else begin
        m1:= tri_vert[3*(t-1)+0];
      end;

      tri_num:= tri_num + 1;

      tri_nabe[3*(t-1)+e-1]:= tri_num;
      tri_vert[3*(tri_num-1)+0]:= m1;
      tri_vert[3*(tri_num-1)+1]:= m2;
      tri_vert[3*(tri_num-1)+2]:= m;
      tri_nabe[3*(tri_num-1)+0]:= t;
      tri_nabe[3*(tri_num-1)+1]:= tri_num - 1;
      tri_nabe[3*(tri_num-1)+2]:= tri_num + 1;
      inc(top);

      if (point_num < top) then begin
        raise Exception.CreateFmt('DTRIS2 - Fatal error #%d! Stack overflow.',[8]);
       end;

      stack[top-1]:= tri_num;

      if ((t = rtri) and(e = redg)) then begin
        break;
      end;
    end;

    tri_nabe[3*(ltri-1)+ledg-1]:= -3 * n - 1;
    tri_nabe[3*(n-1)+1]:= -3 *(tri_num) - 2;
    tri_nabe[3*(tri_num-1)+2]:= -l;
    ltri:= n;
    ledg:= 2;

    error:= swapec(m, top, ltri, ledg, point_num, point_xy, tri_num, tri_vert, tri_nabe, stack);

    if (error <> 0) then begin
      raise Exception.CreateFmt('DTRIS2 - Fatal error #%d! Error return from SWAPEC',[error]);
    end;
  end;
//
//  Now account for the sorting that we did.
//
  for i:= 0 to 2 do begin
    for j:= 0 to tri_num - 1 do begin
      tri_vert[i+j*3]:= indx [ tri_vert[i+j*3] - 1 ];
    end;
  end;

  perm_inv(point_num, indx);

  r82vec_permute(point_num, point_xy, indx);

  Result:= 0;
 end;
//****************************************************************************80

function i4_modp(i: Integer; j: Integer): Integer;

//****************************************************************************80
//
//  Purpose:
//
//    I4_MODP returns the nonnegative remainder of I4 division.
//
//  Formula:
//
//    If
//      NREM:= I4_MODP(I, J)
//      NMULT:= (I - NREM) / J
//    then
//      I:= J * NMULT + NREM
//    where NREM is always nonnegative.
//
//  Comments:
//
//    The MOD function computes a result with the same sign as the
//    quantity being divided.  Thus, suppose you had an angle A,
//    and you wanted to ensure that it was between 0 and 360.
//    Then mod(A,360) would do, if A was positive, but if A
//    was negative, your result would be between -360 and 0.
//
//    On the other hand, I4_MODP(A,360) is between 0 and 360, always.
//
//  Examples:
//
//        I         J     MOD  I4_MODP   I4_MODP Factorization
//
//      107        50       7       7    107:=  2 *  50 + 7
//      107       -50       7       7    107:= -2 * -50 + 7
//     -107        50      -7      43   -107:= -3 *  50 + 43
//     -107       -50      -7      43   -107:=  3 * -50 + 43
//
//  Modified:
//
//    26 May 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer I, the number to be divided.
//
//    Input, Integer J, the number that divides I.
//
//    Output, Integer I4_MODP, the nonnegative remainder when I is
//    divided by J.
//
var
  value: Integer;
 begin

  if (j = 0) then begin
    raise EDivByZero.Create('I4_MODP - Fatal error! I4_MODP(I, J) called with J = 0');
  end;

  value:= i mod j;

  if (value < 0) then begin
    value:= value + Abs(j);
  end;

  Result:= value;
end;
//****************************************************************************80

function i4_sign(i: Integer): Integer;

//****************************************************************************80
//
//  Purpose:
//
//    I4_SIGN returns the sign of an I4.
//
//  Discussion:
//
//    The sign of 0 and all positive integers is taken to be +1.
//    The sign of all negative integers is -1.
//
//  Modified:
//
//    06 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer I, the integer whose sign is desired.
//
//    Output, Integer I4_SIGN, the sign of I.
begin
  if (i < 0) then begin
    Result:= -1;
  end
  else begin
    Result:= 1;
  end;
end;
//****************************************************************************80*

function i4_wrap(ival: Integer; ilo: Integer; ihi: Integer): Integer;

//****************************************************************************80*
//
//  Purpose:
//
//    I4_WRAP forces an I4 to lie between given limits by wrapping.
//
//  Example:
//
//    ILO:= 4, IHI:= 8
//
//    I  I4_WRAP
//
//    -2     8
//    -1     4
//     0     5
//     1     6
//     2     7
//     3     8
//     4     4
//     5     5
//     6     6
//     7     7
//     8     8
//     9     4
//    10     5
//    11     6
//    12     7
//    13     8
//    14     4
//
//  Modified:
//
//    19 August 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer IVAL, an integer value.
//
//    Input, Integer ILO, IHI, the desired bounds for the integer value.
//
//    Output, Integer I4_WRAP, a "wrapped" version of IVAL.
//
var
  jhi: Integer;
  jlo: Integer;
  value: Integer;
  wide: Integer;
begin
  jlo:= Min(ilo, ihi);
  jhi:= Max(ilo, ihi);

  wide:= jhi + 1 - jlo;

  if (wide = 1) then begin
    value:= jlo;
  end
  else begin
    value:= jlo + i4_modp(ival - jlo, wide);
  end;

  Result:= value;
end;

//****************************************************************************80

procedure i4vec_heap_d(n: Integer; var a: TIntegerArray);

//****************************************************************************80
//
//  Purpose:
//
//    I4VEC_HEAP_D reorders an I4VEC into a descending heap.
//
//  Definition:
//
//    A heap is an array A with the property that, for every index J,
//    A[J] >= A[2*J+1] and A[J] >= A[2*J+2],(as long as the indices
//    2*J+1 and 2*J+2 are legal).
//
//  Diagram:
//
//                  A(0)
//                /      \
//            A(1)         A(2)
//          /     \        /  \
//      A(3)       A(4)  A(5) A(6)
//      /  \       /   \
//    A(7) A(8)  A(9) A(10)
//
//  Reference:
//
//    Albert Nijenhuis, Herbert Wilf,
//    Combinatorial Algorithms,
//    Academic Press, 1978, second edition,
//    ISBN 0-12-519260-6.
//
//  Modified:
//
//    30 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the size of the input array.
//
//    Input/output, Integer A[N].
//    On input, an unsorted array.
//    On output, the array has been reordered into a heap.
//
var
  i: Integer;
  ifree: Integer;
  key: Integer;
  m: Integer;
begin
//
//  Only nodes(N/2)-1 down to 0 can be "parent" nodes.
//
  for i:= (n div 2)-1 downto 0 do begin
//
//  Copy the value out of the parent node.
//  Position IFREE is now "open".
//
    key:= a[i];
    ifree:= i;

    while(true) do begin
//
//  Positions 2*IFREE + 1 and 2*IFREE + 2 are the descendants of position
//  IFREE.(One or both may not exist because they equal or exceed N.)
//
      m:= 2 * ifree + 1;
//
//  Does the first position exist?
//
      if (n <= m) then begin
        break;
      end
      else begin
//
//  Does the second position exist?
//
        if (m + 1 < n) then begin
//
//  If both positions exist, take the larger of the two values,
//  and update M if necessary.
//
          if (a[m] < a[m+1]) then begin
            inc(m);
          end;
        end;
//
//  If the large descendant is larger than KEY, move it up,
//  and update IFREE, the location of the free position, and
//  consider the descendants of THIS position.
//
        if (key < a[m]) then begin
          a[ifree]:= a[m];
          ifree:= m;
        end
        else begin
          break;
        end;
      end;
    end;
//
//  When you have stopped shifting items up, return the item you
//  pulled out back to the heap.
//
    a[ifree]:= key;
  end;
end;
//****************************************************************************80

function i4vec_indicator(n: Integer): TIntegerArray;

//****************************************************************************80
//
//  Purpose:
//
//    I4VEC_INDICATOR sets an I4VEC to the indicator vector.
//
//  Modified:
//
//    13 January 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the number of elements of A.
//
//    Output, Integer I4VEC_INDICATOR(N), the initialized array.
//
var
  a: TIntegerArray;
  i: Integer;
begin
  SetLength(a,n);

  for i:= 0 to n-1 do begin
    a[i]:= i + 1;
  end;

  Result:= a;
end;
//****************************************************************************80

procedure i4vec_sort_heap_a(n: Integer; var a: TIntegerArray);

//****************************************************************************80
//
//  Purpose:
//
//    I4VEC_SORT_HEAP_A ascending sorts an I4VEC using heap sort.
//
//  Reference:
//
//    Albert Nijenhuis, Herbert Wilf,
//    Combinatorial Algorithms,
//    Academic Press, 1978, second edition,
//    ISBN 0-12-519260-6.
//
//  Modified:
//
//    30 April 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the number of entries in the array.
//
//    Input/output, Integer A[N].
//    On input, the array to be sorted;
//    On output, the array has been sorted.
//
var
  n1: Integer;
  temp: Integer;
begin

  if (n <= 1) then begin
    exit;
  end;
//
//  1: Put A into descending heap form.
//
  i4vec_heap_d(n, a);
//
//  2: Sort A.
//
//  The largest object in the heap is in A[0].
//  Move it to position A[N-1].
//
  temp:= a[0];
  a[0]:= a[n-1];
  a[n-1]:= temp;
//
//  Consider the diminished heap of size N1.
//
  for n1:= n-1 downto 2 do begin
//
//  Restore the heap structure of the initial N1 entries of A.
//
    i4vec_heap_d(n1, a);
//
//  Take the largest object from A[0] and move it to A[N1-1].
//
    temp:= a[0];
    a[0]:= a[n1-1];
    a[n1-1]:= temp;

  end;
end;
//****************************************************************************80

procedure i4vec_sorted_unique(n: Integer; var a: TIntegerArray; var nuniq: Integer);

//****************************************************************************80
//
//  Purpose:
//
//    I4VEC_SORTED_UNIQUE finds unique elements in a sorted I4VEC.
//
//  Modified:
//
//    02 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the number of elements in A.
//
//    Input/output, Integer A[N].  On input, the sorted
//    integer array.  On output, the unique elements in A.
//
//    Output, Integer *NUNIQ, the number of unique elements in A.
//
var
  i: Integer;
begin
  nuniq:= 0;

  if (n <= 0) then begin
    exit;
  end;

  nuniq:= 1;

  for i:= 1 to n-1 do begin
    if (a[i] <> a[nuniq]) then begin
      nuniq:= nuniq + 1;
      a[nuniq]:= a[i];
    end;
  end;
end;
//****************************************************************************80

function lrline(xu: Real; yu: Real; xv1: Real; yv1: Real; xv2: Real; yv2: Real; dv: Real): Integer;

//****************************************************************************80
//
//  Purpose:
//
//    LRLINE determines where a point lies in relation to a directed line.
//
//  Discussion:
//
//    LRLINE determines whether a point is to the left of, right of,
//    or on a directed line parallel to a line through given points.
//
//  Modified:
//
//    28 August 2003
//
//  Author:
//
//    Barry Joe,
//    Department of Computing Science,
//    University of Alberta,
//    Edmonton, Alberta, Canada  T6G 2H1
//
//  Reference:
//
//    Barry Joe,
//    GEOMPACK - a software package for the generation of meshes
//    using geometric algorithms,
//    Advances in Engineering Software,
//    Volume 13, pages 325-331, 1991.
//
//  Parameters:
//
//    Input, Real XU, YU, XV1, YV1, XV2, YV2, are vertex coordinates; the
//    directed line is parallel to and at signed distance DV to the left of
//    the directed line from(XV1,YV1) to(XV2,YV2);(XU,YU) is the vertex for
//    which the position relative to the directed line is to be determined.
//
//    Input, Real DV, the signed distance, positive for left.
//
//    Output, Integer LRLINE, is +1, 0, or -1 depending on whether(XU,YU) is
//    to the right of, on, or left of the directed line.  LRLINE is 0 if
//    the line degenerates to a point.
//
var
  dx: Real;
  dxu: Real;
  dy: Real;
  dyu: Real;
  t: Real;
  tolabs: Real;
  value: Integer;
const
  tol: Real = 0.0000001;
begin

  dx:= xv2 - xv1;
  dy:= yv2 - yv1;
  dxu:= xu - xv1;
  dyu:= yu - yv1;

  tolabs:= tol * Max(Abs(dx), Max(Abs(dy), Max(Abs(dxu), Max(Abs(dyu), Abs(dv)))));

  t:= dy * dxu - dx * dyu + dv * sqrt(dx * dx + dy * dy);

  if (tolabs < t) then begin
    value:= 1;
  end
  else if (-tolabs <= t) then begin
    value:= 0;
  end
  else if (t < -tolabs) then begin
    value:= -1;
  end
  else begin
    raise Exception.Create('lrline - Fatal error! Point position could not be calculated.');
  end;

  Result:= value;
end;
//****************************************************************************80

function perm_check(n: Integer; var p: TIntegerArray): Boolean;

//****************************************************************************80
//
//  Purpose:
//
//    PERM_CHECK checks that a vector represents a permutation.
//
//  Discussion:
//
//    The routine verifies that each of the integers from 1
//    to N occurs among the N entries of the permutation.
//
//  Modified:
//
//    13 January 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the number of entries.
//
//    Input, Integer P[N], the array to check.
//
//    Output, bool PERM_CHECK, is TRUE if the permutation is OK.
//
var
  found: Boolean;
  i: Integer;
  seek: Integer;
begin
  for seek:= 1 to n do begin
    found:= false;

    for i:= 0 to n-1 do begin
      if (p[i] = seek) then begin
        found:= true;
        break;
      end;
    end;

    if (not found) then begin
      Result:= false;
      exit;
    end;
  end;

  Result:= true;
 end;
//****************************************************************************80

procedure perm_inv(n: Integer; var p: TIntegerArray);

//****************************************************************************80
//
//  Purpose:
//
//    PERM_INV inverts a permutation "in place".
//
//  Modified:
//
//    13 January 2004
//
//  Parameters:
//
//    Input, Integer N, the number of objects being permuted.
//
//    Input/output, Integer P[N], the permutation, in standard index form.
//    On output, P describes the inverse permutation
//
var
  i: Integer;
  i0: Integer;
  i1: Integer;
  i2: Integer;
  iss: Integer;
begin
  if (n <= 0) then begin
    raise Exception.CreateFmt('PERM_INV - Fatal error! Input value of N=%d.',[n]);
  end;

  if (not perm_check(n, p)) then begin
    raise Exception.Create('PERM_INV - Fatal error! The input array does not represent a proper permutation.');
  end;

  for i:= 1 to n do begin
    i1:= p[i-1];

    while(i < i1) do begin
      i2:= p[i1-1];
      p[i1-1]:= -i2;
      i1:= i2;
    end;

    iss:= - i4_sign(p[i-1]);
    p[i-1]:= i4_sign(iss) * Abs(p[i-1]);
  end;

  for i:= 1 to n do begin
    i1:= -p[i-1];

    if (0 <= i1) then begin
      i0:= i;

      while(true) do begin
        i2:= p[i1-1];
        p[i1-1]:= i0;

        if (i2 < 0) then begin
          break;
        end;

        i0:= i1;
        i1:= i2;
      end;
    end;
  end;
end;
//****************************************************************************80

function points_delaunay_naive_2d(n: Integer; var p: TRealArray; var ntri: Integer): TIntegerArray;

//****************************************************************************80
//
//  Purpose:
//
//    POINTS_DELAUNAY_NAIVE_2D computes the Delaunay triangulation in 2D.
//
//  Discussion:
//
//    A naive and inefficient(but extremely simple) method is used.
//
//    This routine is only suitable as a demonstration code for small
//    problems.  Its running time is of order N^4.  Much faster algorithms
//    are available.
//
//    Given a set of nodes in the plane, a triangulation is a set of
//    triples of distinct nodes, forming triangles, so that every
//    point with the convex hull of the set of  nodes is either one
//    of the nodes, or lies on an edge of one or more triangles,
//    or lies within exactly one triangle.
//
//  Modified:
//
//    05 February 2005
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Joseph O'Rourke,
//    Computational Geometry,
//    Cambridge University Press,
//    Second Edition, 1998, page 187.
//
//  Parameters:
//
//    Input, Integer N, the number of nodes.  N must be at least 3.
//
//    Input, Real P[2*N], the coordinates of the nodes.
//
//    Output, Integer *NTRI, the number of triangles.
//
//    Output, Integer POINTS_DELAUNAY_NAIVE_2D[3*NTRI], the indices of the
//    nodes making each triangle.
//
var
  count: Integer;
  flag: Integer;
  i: Integer;
  j: Integer;
  k: Integer;
  m: Integer;
  pass: Integer;
  tri: TIntegerArray;
  xn: Real;
  yn: Real;
  zn: Real;
  z: TRealArray;
begin
  count:= 0;

  SetLength(z,n);

  for i:= 0 to n-1 do begin
    z[i]:= p[0+i*2] * p[0+i*2] + p[1+i*2] * p[1+i*2];
  end;
//
//  First pass counts triangles,
//  Second pass allocates triangles and sets them.
//
  for pass:= 1 to 2 do begin
    if (pass = 2) then  begin
      SetLength(tri, 3*count);
    end;
    count:= 0;
//
//  For each triple(I,J,K):
//
    for i:= 0 to n - 3 do begin
      for j:= i+1 to n-1 do begin
        for k:= i+1 to n-1 do begin
          if (j <> k) then begin
            xn:= (p[1+j*2] - p[1+i*2]) *(z[k] - z[i]) -(p[1+k*2] - p[1+i*2]) *(z[j] - z[i]);
            yn:= (p[0+k*2] - p[0+i*2]) *(z[j] - z[i]) -(p[0+j*2] - p[0+i*2]) *(z[k] - z[i]);
            zn:= (p[0+j*2] - p[0+i*2]) *(p[1+k*2] - p[1+i*2]) -(p[0+k*2] - p[0+i*2]) *(p[1+j*2] - p[1+i*2]);

            flag:= Integer(zn < 0);

            if (flag<>0) then begin
              for m:= 0 to n-1 do begin
                flag:= Integer(Boolean(flag) and((p[0+m*2] - p[0+i*2]) * xn +(p[1+m*2] - p[1+i*2]) * yn +(z[m] - z[i]) * zn <= 0));
              end
            end;

            if (flag<>0) then begin
              if (pass = 2) then begin
                tri[0+count*3]:= i;
                tri[1+count*3]:= j;
                tri[2+count*3]:= k;
              end;
              count:= count + 1;
            end;
          end;
        end;
      end;
    end;
  end;

  ntri:= count;

  Result:= tri;
end;
//****************************************************************************80

function r8_epsilon(): Real;

//****************************************************************************80
//
//  Purpose:
//
//    R8_EPSILON returns the round off unit for Real precision arithmetic.
//
//  Discussion:
//
//    R8_EPSILON is a number R which is a power of 2 with the property that,
//    to the precision of the computer's arithmetic,
//      1 < 1 + R
//    but
//      1:= (1 + R / 2)
//
//  Modified:
//
//    06 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, Real R8_EPSILON, the floating point round-off unit.
//
var
  r: Real;
begin

  r:= 1.0;

  while(1.0 < 1.0 + r) do begin
    r:= r / 2.0;
  end;

  Result:= (2.0 * r);
end;
//****************************************************************************80

procedure r82vec_part_quick_a(from: Integer; n: Integer; var a: TRealArray; var l: Integer; var r: Integer);

//****************************************************************************80
//
//  Purpose:
//
//    R82VEC_PART_QUICK_A reorders an R2 vector as part of a quick sort.
//
//  Discussion:
//
//    The routine reorders the entries of A.  Using A(1:2,1) as a
//    key, all entries of A that are less than or equal to the key will
//    precede the key, which precedes all entries that are greater than the key.
//
//  Example:
//
//    Input:
//
//      N:= 8
//
//      A:= ((2,4),(8,8),(6,2),(0,2),(10,6),(10,0),(0,6),(4,8))
//
//    Output:
//
//      L:= 2, R:= 4
//
//      A:= ((0,2),(0,6),(2,4),(8,8),(6,2),(10,6),(10,0),(4,8))
//             -----------          ----------------------------------
//                LEFT          KEY           RIGHT
//
//  Modified:
//
//    01 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the number of entries of A.
//
//    Input/output, Real A[N*2].  On input, the array to be checked.
//    On output, A has been reordered as described above.
//
//    Output, Integer *L, *R, the indices of A that define the three segments.
//    Let KEY:= the input value of A(1:2,1).  Then
//    I <= L                 A(1:2,I) < KEY;
//         L < I < R         A(1:2,I):= KEY;
//                 R <= I    A(1:2,I) > KEY.
//
var
  i: Integer;
  j: Integer;
  key: TRealArray;
  ll: Integer;
  m: Integer;
  rr: Integer;
begin
  SetLength(key,2);

  if (n < 1) then begin
    raise Exception.Create('R82VEC_PART_QUICK_A - Fatal error! N < 1.');
  end;

  if (n = 1) then begin
    l:= 0;
    r:= 2;

    exit;
  end;

  key[0]:= a[from+(2*0+0)];
  key[1]:= a[from+(2*0+1)];
  m:= 1;
//
//  The elements of unknown size have indices between L+1 and R-1.
//
  ll:= 1;
  rr:= n + 1;

  for i:= 2 to n do begin
    if (r8vec_gt(from+(2*ll), 0, 2, a, key)) then begin
      rr:= rr - 1;
      r8vec_swap(from+(2*(rr-1)), from+(2*ll), 2, a, a);
    end
    else if (r8vec_eq(from+(2*ll), 0, 2, a, key)) then begin
      m:= m + 1;
      r8vec_swap(from+(2*(m-1)), from+(2*ll), 2, a, a);
      ll:= ll + 1;
    end
    else if (r8vec_lt(from+2*ll, 0, 2, a, key)) then begin
      ll:= ll + 1;
    end;
  end;
//
//  Now shift small elements to the left, and KEY elements to center.
//
  for i:= 0 to ll - m -1 do begin
    for j:= 0 to 1 do begin
      a[from+(2*i+j)]:= a[from+(2*(i+m)+j)];
    end;
  end;

  ll:= ll - m;

  for i:= ll to ll+m -1 do begin
    for j:= 0 to 1 do begin
      a[from+(2*i+j)]:= key[j];
    end;
  end;

  l:= ll;
  r:= rr;
end;
//****************************************************************************80*

procedure r82vec_permute(n: Integer; var a: TRealArray; var p: TIntegerArray);

//****************************************************************************80*
//
//  Purpose:
//
//    R82VEC_PERMUTE permutes an R2 vector in place.
//
//  Discussion:
//
//    This routine permutes an TRealArray "objects", but the same
//    logic can be used to permute an array of objects of any arithmetic
//    type, or an array of objects of any complexity.  The only temporary
//    storage required is enough to store a single object.  The number
//    of data movements made is N + the number of cycles of order 2 or more,
//    which is never more than N + N/2.
//
//  Example:
//
//    Input:
//
//      N:= 5
//      P:= (2,    4,    5,    1,    3)
//      A:= (1.0,  2.0,  3.0,  4.0,  5.0)
//(11.0, 22.0, 33.0, 44.0, 55.0)
//
//    Output:
//
//      A:= (2.0,  4.0,  5.0,  1.0,  3.0)
//(22.0, 44.0, 55.0, 11.0, 33.0).
//
//  Modified:
//
//    19 February 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the number of objects.
//
//    Input/output, Real A[2*N], the array to be permuted.
//
//    Input, Integer P[N], the permutation.  P(I):= J means
//    that the I-th element of the output array should be the J-th
//    element of the input array.  P must be a legal permutation
//    of the integers from 1 to N, otherwise the algorithm will
//    fail catastrophically.
//
var
  a_temp: array [0..1] of Real;
  i: Integer;
  iget: Integer;
  iput: Integer;
  istart: Integer;
begin
  if (not perm_check(n, p)) then begin
    raise Exception.Create('R82VEC_PERMUTE - Fatal error! The input array does not represent a proper permutation.');
  end;
//
//  Search for the next element of the permutation that has not been used.
//
  for istart:= 1 to n do begin
    if (p[istart-1] < 0) then begin
      continue;
    end
    else if (p[istart-1] = istart) then begin
      p[istart-1]:= -p[istart-1];
      continue;
    end
    else begin
      a_temp[0]:= a[0+(istart-1)*2];
      a_temp[1]:= a[1+(istart-1)*2];
      iget:= istart;
//
//  Copy the new value into the vacated entry.
//
      while(true) do begin
        iput:= iget;
        iget:= p[iget-1];

        p[iput-1]:= -p[iput-1];

        if ((iget < 1) or(n < iget)) then begin
          raise Exception.Create('R82VEC_PERMUTE - Fatal error!');
        end;

        if (iget = istart) then begin
          a[0+(iput-1)*2]:= a_temp[0];
          a[1+(iput-1)*2]:= a_temp[1];
          break;
        end;

        a[0+(iput-1)*2]:= a[0+(iget-1)*2];
        a[1+(iput-1)*2]:= a[1+(iget-1)*2];
      end;
    end;
  end;
//
//  Restore the signs of the entries.
//
  for i:= 0 to n-1 do begin
    p[i]:= -p[i];
  end;
end;
//****************************************************************************80

function r82vec_sort_heap_index_a(n: Integer; var a: TRealArray): TIntegerArray;

//****************************************************************************80
//
//  Purpose:
//
//    R82VEC_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an R2 vector.
//
//  Discussion:
//
//    The sorting is not actually carried out.  Rather an index array is
//    created which defines the sorting.  This array may be used to sort
//    or index the array, or to sort or index related arrays keyed on the
//    original array.
//
//    Once the index array is computed, the sorting can be carried out
//    "implicitly:
//
//      A(1:2,INDX(I)), I:= 1 to N is sorted,
//
//    or explicitly, by the call
//
//      call R82VEC_PERMUTE(N, A, INDX)
//
//    after which A(1:2,I), I:= 1 to N is sorted.
//
//  Modified:
//
//    13 January 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the number of entries in the array.
//
//    Input, Real A[2*N], an array to be index-sorted.
//
//    Output, Integer R82VEC_SORT_HEAP_INDEX_A[N], the sort index.  The
//    I-th element of the sorted array is A(0:1,R82VEC_SORT_HEAP_INDEX_A(I-1)).
//
var
  aval: array [0..1] of Real;
  i: Integer;
  indx: TIntegerArray;
  indxt: Integer;
  ir: Integer;
  j: Integer;
  l: Integer;
begin
  if (n < 1) then begin
    Result:= nil;
  end;

  if (n = 1) then begin
    SetLength(indx,1);
    indx[0]:= 1;

    Result:= indx;
    exit;
  end;

  indx:= i4vec_indicator(n);

  l:= (n div 2) + 1;
  ir:= n;

  while(true) do begin
    if (1 < l) then begin
      l:= l - 1;
      indxt:= indx[l-1];
      aval[0]:= a[0+(indxt-1)*2];
      aval[1]:= a[1+(indxt-1)*2];
    end
    else begin
      indxt:= indx[ir-1];
      aval[0]:= a[0+(indxt-1)*2];
      aval[1]:= a[1+(indxt-1)*2];
      indx[ir-1]:= indx[0];
      ir:= ir - 1;

      if (ir = 1) then begin
        indx[0]:= indxt;
        break;
      end;
    end;

    i:= l;
    j:= l + l;

    while(j <= ir) do begin
      if (j < ir) then begin
        if ((a[0+(indx[j-1]-1)*2] <  a[0+(indx[j]-1)*2])
        or ((a[0+(indx[j-1]-1)*2] =  a[0+(indx[j]-1)*2]) and (a[1+(indx[j-1]-1)*2] <  a[1+(indx[j]-1)*2]))) then begin
          j:= j + 1;
        end
      end;

      if ((aval[0] <  a[0+(indx[j-1]-1)*2])
      or ((aval[0] =  a[0+(indx[j-1]-1)*2]) and (aval[1] <  a[1+(indx[j-1]-1)*2]))) then begin
        indx[i-1]:= indx[j-1];
        i:= j;
        j:= j + j;
      end
      else begin
        j:= ir + 1;
      end;
    end;

    indx[i-1]:= indxt;
   end;

  Result:= indx;
end;
//****************************************************************************80

procedure r82vec_sort_quick_a(n: Integer; var a: TRealArray);

//****************************************************************************80
//
//  Purpose:
//
//    R82VEC_SORT_QUICK_A ascending sorts an R2 vector using quick sort.
//
//  Discussion:
//
//    The data structure is a set of N pairs of real numbers.
//    These values are stored in a one dimensional array, by pairs.
//
//  Modified:
//
//    01 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the number of entries in the array.
//
//    Input/output, Real A[N*2].
//    On input, the array to be sorted.
//    On output, the array has been sorted.
//
var
  base: Integer;
  l_segment: Integer;
  level: Integer;
  n_segment: Integer;
  rsave: array of Integer;
  r_segment: Integer;
const
  LEVEL_MAX: Integer = 25;
begin
  SetLength(rsave, LEVEL_MAX);

  if (n < 1) then begin
    raise Exception.Create('R82VEC_SORT_QUICK_A - Fatal error! N < 1.');
  end;

  if (n = 1) then begin
    exit;
  end;

  level:= 1;
  rsave[level-1]:= n + 1;
  base:= 1;
  n_segment:= n;

  while(0 < n_segment) do begin
//
//  Partition the segment.
//
    r82vec_part_quick_a(2*(base-1)+0, n_segment, a, l_segment, r_segment);
//
//  If the left segment has more than one element, we need to partition it.
//
    if (1 < l_segment) then begin
      if (LEVEL_MAX < level) then begin
        raise Exception.Create(Format('R82VEC_SORT_QUICK_A - Fatal error! Exceeding recursion maximum of %d.', [LEVEL_MAX]));
      end;

      level:= level + 1;
      n_segment:= l_segment;
      rsave[level-1]:= r_segment + base - 1;
    end
//
//  The left segment and the middle segment are sorted.
//  Must the right segment be partitioned?
//
    else if (r_segment < n_segment) then begin
      n_segment:= n_segment + 1 - r_segment;
      base:= base + r_segment - 1;
    end
//
//  Otherwise, we back up a level if there is an earlier one.
//
    else begin
      while (true) do begin
        if (level <= 1) then begin
          n_segment:= 0;
          break;
        end;

        base:= rsave[level-1];
        n_segment:= rsave[level-2] - rsave[level-1];
        level:= level - 1;

        if (0 < n_segment) then begin
          break;
        end;
      end;
    end;
  end
end;
//****************************************************************************80

procedure r8mat_uniform(m: Integer; n: Integer; b: Real; c: Real; var seed: Integer; var r: TRealArray);

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_UNIFORM fills a Real precision array with scaled pseudorandom values.
//
//  Discussion:
//
//    This routine implements the recursion
//
//      seed:= 16807 * seed mod(2**31 - 1)
//      unif:= seed /(2**31 - 1)
//
//    The integer arithmetic never requires more than 32 bits,
//    including a sign bit.
//
//  Modified:
//
//    30 January 2005
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Paul Bratley, Bennett Fox, Linus Schrage,
//    A Guide to Simulation,
//    Springer Verlag, pages 201-202, 1983.
//
//    Bennett Fox,
//    Algorithm 647:
//    Implementation and Relative Efficiency of Quasirandom
//    Sequence Generators,
//    ACM Transactions on Mathematical Software,
//    Volume 12, Number 4, pages 362-376, 1986.
//
//    Peter Lewis, Allen Goodman, James Miller,
//    A Pseudo-Random Number Generator for the System/360,
//    IBM Systems Journal,
//    Volume 8, pages 136-143, 1969.
//
//  Parameters:
//
//    Input, Integer M, N, the number of rows and columns.
//
//    Input, Real B, C, the limits of the pseudorandom values.
//
//    Input/output, Integer *SEED, the "seed" value.  Normally, this
//    value should not be 0, otherwise the output value of SEED
//    will still be 0, and R8_UNIFORM will be 0.  On output, SEED has
//    been updated.
//
//    Output, Real R8MAT_UNIFORM[M*N], a matrix of pseudorandom values.
//
var
  i: Integer;
  j: Integer;
  k: Integer;
begin
  for j:= 0 to n-1 do begin
    for i:= 0 to m-1 do begin
      k:= seed div 127773;

      seed:= 16807 *(seed - k * 127773) - k * 2836;

      if (seed < 0) then begin
        seed:= seed + 2147483647;
      end;
//
//  Although SEED can be represented exactly as a 32 bit integer,
//  it generally cannot be represented exactly as a 32 bit real number!
//
      r[i+j*m]:= b +(c - b) * seed * 4.656612875E-10;
    end;
  end;
end;
//****************************************************************************80

function r8vec_eq(from1: Integer; from2: Integer; n: Integer; var a1: TRealArray; var a2: TRealArray): Boolean;

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_EQ is true if every pair of entries in two vectors is equal.
//
//  Modified:
//
//    28 August 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the number of entries in the vectors.
//
//    Input, Real A1[N], A2[N], two vectors to compare.
//
//    Output, bool R8VEC_EQ.
//    R8VEC_EQ is TRUE if every pair of elements A1(I) and A2(I) are equal,
//    and FALSE otherwise.
//
var
  i: Integer;
  retval: Boolean;
begin
  retval:= true;

  for i:= 0 to n-1 do begin
    if (a1[from1+i] <> a2[from2+i]) then begin
      retval:= false;
      break;
    end;
  end;

  Result:= retval;
end;
//****************************************************************************80

function r8vec_gt(from1: Integer; from2: Integer; n: Integer; var a1: TRealArray; var a2: TRealArray): Boolean;

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_GT =(A1 > A2) for real vectors.
//
//  Discussion:
//
//    The comparison is lexicographic.
//
//    A1 > A2  <=>                              A1(1) > A2(1) or
//(A1(1)     = A2(1)     and A1(2) > A2(2)) or
//                 ...
//(A1(1:N-1) = A2(1:N-1) and A1(N) > A2(N)
//
//  Modified:
//
//    28 August 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the dimension of the vectors.
//
//    Input, Real A1[N], A2[N], the vectors to be compared.
//
//    Output, bool R8VEC_GT, is TRUE if and only if A1 > A2.
//
var
  i: Integer;
  retval: Boolean;
begin
  retval:= false;

  for i:= 0 to n-1 do begin
    if (a2[from2+i] < a1[from1+i]) then begin
      retval:= true;
      break;
    end
    else if (a1[from1+i] < a2[from2+i]) then begin
      retval:= false;
      break;
    end;
  end;

  Result:= retval;
end;
//****************************************************************************80

function r8vec_lt(from1: Integer; from2: Integer; n: Integer; var a1: TRealArray; var a2: TRealArray): Boolean;

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_LT =(A1 < A2) for real vectors.
//
//  Discussion:
//
//    The comparison is lexicographic.
//
//    A1 < A2  <=>                              A1(1) < A2(1) or
//(A1(1)     = A2(1)     and A1(2) < A2(2)) or
//                 ...
//(A1(1:N-1) = A2(1:N-1) and A1(N) < A2(N)
//
//  Modified:
//
//    28 August 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the dimension of the vectors.
//
//    Input, Real A1[N], A2[N], the vectors to be compared.
//
//    Output, bool R8VEC_LT, is TRUE if and only if A1 < A2.
//
var
  i: Integer;
  retval: Boolean;
begin
  retval:= false;

  for i:= 0 to n-1 do begin
    if (a1[from1+i] < a2[from2+i]) then begin
      retval:= true;
      break;
    end
    else if (a2[from2+i] < a1[from1+i]) then begin
      retval:= false;
      break;
    end;
  end;

  Result:= retval;
end;
//****************************************************************************80

procedure r8vec_swap(from1: Integer; from2: Integer; n: Integer; var a1: TRealArray; var a2: TRealArray);

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_SWAP swaps the entries of two R8VEC's.
//
//  Modified:
//
//    28 August 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Integer N, the number of entries in the arrays.
//
//    Input/output, Real A1[N], A2[N], the vectors to swap.
//
var
  i: Integer;
  temp: Real;
begin
  for i:= from1 to from1+n-1 do begin
    temp:= a1[i];
    a1[i]:= a2[i-from1+from2];
    a2[i-from1+from2]:= temp;
  end
end;
//****************************************************************************80

function swapec(i: Integer; var top: Integer; var btri: Integer; var bedg: Integer; point_num: Integer; var point_xy: TRealArray; tri_num: Integer; var tri_vert: TIntegerArray; var tri_nabe: TIntegerArray; var stack: TIntegerArray): Integer;

//****************************************************************************80
//
//  Purpose:
//
//    SWAPEC swaps diagonal edges until all triangles are Delaunay.
//
//  Discussion:
//
//    The routine swaps diagonal edges in a 2D triangulation, based on
//    the empty circumcircle criterion, until all triangles are Delaunay,
//    given that I is the index of the new vertex added to the triangulation.
//
//  Modified:
//
//    03 September 2003
//
//  Author:
//
//    Barry Joe,
//    Department of Computing Science,
//    University of Alberta,
//    Edmonton, Alberta, Canada  T6G 2H1
//
//  Reference:
//
//    Barry Joe,
//    GEOMPACK - a software package for the generation of meshes
//    using geometric algorithms,
//    Advances in Engineering Software,
//    Volume 13, pages 325-331, 1991.
//
//  Parameters:
//
//    Input, Integer I, the index of the new vertex.
//
//    Input/output, Integer *TOP, the index of the top of the stack.
//    On output, TOP is zero.
//
//    Input/output, Integer *BTRI, *BEDG; on input, if positive, are the
//    triangle and edge indices of a boundary edge whose updated indices
//    must be recorded.  On output, these may be updated because of swaps.
//
//    Input, Integer POINT_NUM, the number of points.
//
//    Input, Real POINT_XY[POINT_NUM*2], the coordinates of the points.
//
//    Input, Integer TRI_NUM, the number of triangles.
//
//    Input/output, Integer TRI_VERT[TRI_NUM*3], the triangle incidence list.
//    May be updated on output because of swaps.
//
//    Input/output, Integer TRI_NABE[TRI_NUM*3], the triangle neighbor list;
//    negative values are used for links of the counter-clockwise linked
//    list of boundary edges; May be updated on output because of swaps.
//
//      LINK:= -(3*I + J-1) where I, J:= triangle, edge index.
//
//    Workspace, Integer STACK[MAXST]; on input, entries 1 through TOP
//    contain the indices of initial triangles(involving vertex I)
//    put in stack; the edges opposite I should be in interior; entries
//    TOP+1 through MAXST are used as a stack.
//
//    Output, Integer SWAPEC, is set to 8 for abnormal return.
//
var
  a: Integer;
  b: Integer;
  c: Integer;
  e: Integer;
  ee: Integer;
  em1: Integer;
  ep1: Integer;
  f: Integer;
  fm1: Integer;
  fp1: Integer;
  l: Integer;
  r: Integer;
  s: Integer;
  swap: Integer;
  t: Integer;
  tt: Integer;
  u: Integer;
  x: Real;
  y: Real;

begin
//
//  Determine whether triangles in stack are Delaunay, and swap
//  diagonal edge of convex quadrilateral if not.
//
  x:= point_xy[2*(i-1)+0];
  y:= point_xy[2*(i-1)+1];

  while(true) do begin
    if (top <= 0)  then begin
      break;
    end;

    t:= stack[(top)-1];
    top:= top - 1;

    if (tri_vert[3*(t-1)+0] = i) then begin
      e:= 2;
      b:= tri_vert[3*(t-1)+2];
    end
    else if (tri_vert[3*(t-1)+1] = i) then begin
      e:= 3;
      b:= tri_vert[3*(t-1)+0];
    end
    else begin
      e:= 1;
      b:= tri_vert[3*(t-1)+1];
    end;

    a:= tri_vert[3*(t-1)+e-1];
    u:= tri_nabe[3*(t-1)+e-1];

    if (tri_nabe[3*(u-1)+0] = t) then begin
      f:= 1;
      c:= tri_vert[3*(u-1)+2];
    end
    else if (tri_nabe[3*(u-1)+1] = t) then begin
      f:= 2;
      c:= tri_vert[3*(u-1)+0];
    end
    else begin
      f:= 3;
      c:= tri_vert[3*(u-1)+1];
    end;

    swap:= diaedg(x, y, point_xy[2*(a-1)+0], point_xy[2*(a-1)+1], point_xy[2*(c-1)+0], point_xy[2*(c-1)+1], point_xy[2*(b-1)+0], point_xy[2*(b-1)+1]);

    if (swap = 1) then begin
      em1:= i4_wrap(e - 1, 1, 3);
      ep1:= i4_wrap(e + 1, 1, 3);
      fm1:= i4_wrap(f - 1, 1, 3);
      fp1:= i4_wrap(f + 1, 1, 3);

      tri_vert[3*(t-1)+ep1-1]:= c;
      tri_vert[3*(u-1)+fp1-1]:= i;
      r:= tri_nabe[3*(t-1)+ep1-1];
      s:= tri_nabe[3*(u-1)+fp1-1];
      tri_nabe[3*(t-1)+ep1-1]:= u;
      tri_nabe[3*(u-1)+fp1-1]:= t;
      tri_nabe[3*(t-1)+e-1]:= s;
      tri_nabe[3*(u-1)+f-1]:= r;

      if (0 < tri_nabe[3*(u-1)+fm1-1]) then begin
        top:= top + 1;
        stack[top-1]:= u;
      end;

      if (0 < s) then begin
        if (tri_nabe[3*(s-1)+0] = u) then begin
          tri_nabe[3*(s-1)+0]:= t;
        end
        else if (tri_nabe[3*(s-1)+1] = u) then begin
          tri_nabe[3*(s-1)+1]:= t;
        end
        else begin
          tri_nabe[3*(s-1)+2]:= t;
        end;

        top:= top + 1;

        if (point_num < top) then begin
          Result:= 8;
          exit;
        end;

        stack[top-1]:= t;
      end
      else begin
        if ((u = btri) and(fp1 = bedg)) then begin
          btri:= t;
          bedg:= e;
        end;

        l:= -(3 * t + e - 1);
        tt:= t;
        ee:= em1;

        while(0 < tri_nabe[3*(tt-1)+ee-1]) do begin
          tt:= tri_nabe[3*(tt-1)+ee-1];

          if (tri_vert[3*(tt-1)+0] = a) then begin
            ee:= 3;
          end
          else if (tri_vert[3*(tt-1)+1] = a) then begin
            ee:= 1;
          end
          else begin
            ee:= 2;
          end;
        end;

        tri_nabe[3*(tt-1)+ee-1]:= l;
      end;

      if (0 < r) then begin
        if (tri_nabe[3*(r-1)+0] = t) then begin
          tri_nabe[3*(r-1)+0]:= u;
        end
        else if (tri_nabe[3*(r-1)+1] = t) then begin
          tri_nabe[3*(r-1)+1]:= u;
        end
        else begin
          tri_nabe[3*(r-1)+2]:= u;
        end
      end
      else begin
        if ((t = btri) and(ep1 = bedg)) then begin
          btri:= u;
          bedg:= f;
        end;

        l:= -(3 * u + f - 1);
        tt:= u;
        ee:= fm1;

        while(0 < tri_nabe[3*(tt-1)+ee-1]) do begin
          tt:= tri_nabe[3*(tt-1)+ee-1];

          if (tri_vert[3*(tt-1)+0] = b) then begin
            ee:= 3;
          end
          else if (tri_vert[3*(tt-1)+1] = b) then begin
            ee:= 1;
          end
          else begin
            ee:= 2;
          end;
        end;

        tri_nabe[3*(tt-1)+ee-1]:= l;
      end;
    end;
  end;

  Result:= 0;
end;
//****************************************************************************80

function triangle_circumcenter_2d(var t: TRealArray): TRealArray;

//****************************************************************************80
//
//  Purpose:
//
//    TRIANGLE_CIRCUMCENTER_2D computes the circumcenter of a triangle in 2D.
//
//  Discussion:
//
//    The circumcenter of a triangle is the center of the circumcircle, the
//    circle that passes through the three vertices of the triangle.
//
//    The circumcircle contains the triangle, but it is not necessarily the
//    smallest triangle to do so.
//
//    If all angles of the triangle are no greater than 90 degrees, then
//    the center of the circumscribed circle will lie inside the triangle.
//    Otherwise, the center will lie outside the circle.
//
//    The circumcenter is the intersection of the perpendicular bisectors
//    of the sides of the triangle.
//
//    In geometry, the circumcenter of a triangle is often symbolized by "O".
//
//  Modified:
//
//    09 February 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, Real T[2*3], the triangle vertices.
//
//    Output, Real *X, *Y, the coordinates of the circumcenter of the triangle.
//
const
  DIM_NUM: Integer = 2;
var
  asq: Real;
  bot: Real;
  center: TRealArray;
  csq: Real;
  top1: Real;
  top2: Real;
begin
  SetLength(center, DIM_NUM);

  asq:= (t[0+1*2] - t[0+0*2]) *(t[0+1*2] - t[0+0*2]) +(t[1+1*2] - t[1+0*2]) *(t[1+1*2] - t[1+0*2]);

  csq:= (t[0+2*2] - t[0+0*2]) *(t[0+2*2] - t[0+0*2]) +(t[1+2*2] - t[1+0*2]) *(t[1+2*2] - t[1+0*2]);

  top1:= (t[1+1*2] - t[1+0*2]) * csq -(t[1+2*2] - t[1+0*2]) * asq;
  top2:= (t[0+1*2] - t[0+0*2]) * csq -(t[0+2*2] - t[0+0*2]) * asq;

  bot:= (t[1+1*2] - t[1+0*2]) *(t[0+2*2] - t[0+0*2]) -(t[1+2*2] - t[1+0*2]) *(t[0+1*2] - t[0+0*2]);

  center[0]:= t[0+0*2] + 0.5 * top1 / bot;
  center[1]:= t[1+0*2] + 0.5 * top2 / bot;

  Result:= center;
end;
//****************************************************************************80

procedure vbedg(x: Real; y: Real; point_num: Integer; var point_xy: TRealArray; tri_num: Integer; var tri_vert: TIntegerArray; var tri_nabe: TIntegerArray; var ltri: Integer; var ledg: Integer; var rtri: Integer; var redg: Integer);

//****************************************************************************80
//
//  Purpose:
//
//    VBEDG determines which boundary edges are visible to a point.
//
//  Discussion:
//
//    The point(X,Y) is assumed to be outside the convex hull of the
//    region covered by the 2D triangulation.
//
//  Author:
//
//    Barry Joe,
//    Department of Computing Science,
//    University of Alberta,
//    Edmonton, Alberta, Canada  T6G 2H1
//
//  Reference:
//
//    Barry Joe,
//    GEOMPACK - a software package for the generation of meshes
//    using geometric algorithms,
//    Advances in Engineering Software,
//    Volume 13, pages 325-331, 1991.
//
//  Modified:
//
//    02 September 2003
//
//  Parameters:
//
//    Input, Real X, Y, the coordinates of a point outside the convex hull
//    of the current triangulation.
//
//    Input, Integer POINT_NUM, the number of points.
//
//    Input, Real POINT_XY[POINT_NUM*2], the coordinates of the vertices.
//
//    Input, Integer TRI_NUM, the number of triangles.
//
//    Input, Integer TRI_VERT[TRI_NUM*3], the triangle incidence list.
//
//    Input, Integer TRI_NABE[TRI_NUM*3], the triangle neighbor list; negative
//    values are used for links of a counter clockwise linked list of boundary
//    edges;
//      LINK:= -(3*I + J-1) where I, J:= triangle, edge index.
//
//    Input/output, Integer *LTRI, *LEDG.  If LTRI <> 0 then these values are
//    assumed to be already computed and are not changed, else they are updated.
//    On output, LTRI is the index of boundary triangle to the left of the
//    leftmost boundary triangle visible from(X,Y), and LEDG is the boundary
//    edge of triangle LTRI to the left of the leftmost boundary edge visible
//    from(X,Y).  1 <= LEDG <= 3.
//
//    Input/output, Integer *RTRI.  On input, the index of the boundary triangle
//    to begin the search at.  On output, the index of the rightmost boundary
//    triangle visible from(X,Y).
//
//    Input/output, Integer *REDG, the edge of triangle RTRI that is visible
//    from(X,Y).  1 <= REDG <= 3.
//
var
  a: Integer;
  ax: Real;
  ay: Real;
  b: Integer;
  bx: Real;
  by: Real;
  done: Boolean;
  e: Integer;
  l: Integer;
  lr: Integer;
  t: Integer;
begin
//
//  Find the rightmost visible boundary edge using links, then possibly
//  leftmost visible boundary edge using triangle neighbor information.
//
  if (ltri = 0) then begin
    done:= false;
    ltri:= rtri;
    ledg:= redg;
  end
  else begin
    done:= true;
  end;

  while(true) do begin
    l:= -tri_nabe[3*((rtri)-1)+(redg)-1];
    t:= l div 3;
    e:= 1 + l mod 3;
    a:= tri_vert[3*(t-1)+e-1];

    if (e <= 2) then begin
      b:= tri_vert[3*(t-1)+e];
    end
    else begin
      b:= tri_vert[3*(t-1)+0];
    end;

    ax:= point_xy[2*(a-1)+0];
    ay:= point_xy[2*(a-1)+1];

    bx:= point_xy[2*(b-1)+0];
    by:= point_xy[2*(b-1)+1];

    lr:= lrline(x, y, ax, ay, bx, by, 0.0);

    if (lr <= 0) then begin
      break;
    end;

    rtri:= t;
    redg:= e;
  end;

  if (done) then begin
    exit;
  end;

  t:= ltri;
  e:= ledg;

  while(true) do begin
    b:= tri_vert[3*(t-1)+e-1];
    e:= i4_wrap(e-1, 1, 3);

    while(0 < tri_nabe[3*(t-1)+e-1]) do begin
      t:= tri_nabe[3*(t-1)+e-1];

      if (tri_vert[3*(t-1)+0] = b) then begin
        e:= 3;
      end
      else if (tri_vert[3*(t-1)+1] = b) then begin
        e:= 1;
      end
      else begin
        e:= 2;
      end;
    end;

    a:= tri_vert[3*(t-1)+e-1];
    ax:= point_xy[2*(a-1)+0];
    ay:= point_xy[2*(a-1)+1];

    bx:= point_xy[2*(b-1)+0];
    by:= point_xy[2*(b-1)+1];

    lr:= lrline(x, y, ax, ay, bx, by, 0.0);

    if (lr <= 0) then begin
      break;
    end;
  end;

  ltri:= t;
  ledg:= e;
end;

end.