#include <stdio.h>
/* #include <stdlib.h> */

/*@
  terminates \false;
  assigns \nothing;
  ensures \false;
*/

void exit(int status);


#define true 1
#define false 0

static char pathdata[] = "style=\"fill:none;stroke:#000000;stroke-width:1px;stroke-opacity:1\"";


/* %define DEBUG 1 */

struct point {
  int x, y;
};

struct edge {
/* all fields are indices, used instead of pointers. 
   p1 is the source point and p2 is the target point of the edge
   hopefully, t1 should be on the counter-clockwise side of 
   the edge, and t2 should be on the other side.
   Edges on the convex hull should have a valid triangle as t1 field,
   the t2 field should be -1 (the field t2 is thus an optional field. */
  int p1, p2, t1, t2;
};

/* We hope to maintain the invariant that the counter-clockwise
   ordering of the points is respected by the edge e1 (resp. e2, e3)
   when the orientation o1 (resp. o2, o3) is true, and by the reverse
   of edge e1 when the orientation o2 (res. o2, o3) is false. */
   
struct triangle {
  int e1, e2, e3, o1, o2, o3;
};

#define array_size 1000
static struct triangle t_map[array_size];
static struct edge e_map[array_size];
static struct point p_map[array_size];
static int h_map[array_size];
static int rotate_tmp[array_size];
static int work_stack[array_size];

static int t_bound = 0;
static int e_bound = 0;
static int p_bound = 0;
static int h_bound = 0;
static int w_bound = 0;
static int outer_edge = 0;

void dump_array(int a[], int bound, char text[]) {
  int i = 0;
  printf("%% %s\n", text);
  while(i < bound) {
    printf("%% %d\n", a[i]);
    i ++;
  }
}
void dump_state () {
  int j = 0;
  printf("%% p_bound %d\n", p_bound);
  printf("%% dumping points\n");
  while (j < p_bound) {
    printf("%% %d %d\n", p_map[j].x, p_map[j].y);
    j++;
  }
  printf("%% dumping edges\n");
  j=0;
  while (j < e_bound) {
    printf("%% %d %d %d %d\n", e_map[j].p1, e_map[j].p2, e_map[j].t1, e_map[j].t2);
    j++;
  }
  printf("%% dumping triangles\n");
  j=0;
  while (j < t_bound) {
    printf("%% edges       %d %d %d\n%% orientation %d %d %d\n",
           t_map[j].e1, t_map[j].e2, t_map[j].e3,
           t_map[j].o1, t_map[j].o2, t_map[j].o3);
    j++;
  }
  dump_array(h_map, h_bound, "dumping convex hull");
  dump_array(work_stack, w_bound, "dumpin working stack");
}

void add_triangle(int i);

void make_triangles(int i);

int find_red_edge(int i);

/* First section: reading in data */

int read_point (int i) {
  int the_x, the_y;
  int success;
  success = scanf("%d %d", &the_x, &the_y);
  if (2 == success) {
    p_map[i].x = the_x;
    p_map[i].y = the_y;
    return true;
  } else {
    return false;
  }
}

void read_points () {
  int escape = false;
  int i;
  for(p_bound = 0; true != escape; ) {
    if(true == read_point(p_bound)) {
      i = p_bound;
      p_bound++;
      make_triangles(i);
    } else {
      escape = true;
    };
  }
}

static int det_debug = 0;

int det3(int a, int b, int c, int d, int e, int f, int g, int h, int i) {
  int value = a * e * i + d * h * c + g * b * f -
          a * f * h - b * d * i - c * e * g;
  if (0 == value) {
    det_debug = 0;
    printf("%% computed a zero determinant\n%d %d %d\n%d %d %d\n%d %d %d\n",
	   a, b, c, d, e, f, g, h, i);
  }
  return value;
}

void add_work(int the_edge) {
  work_stack[w_bound] = the_edge;
  w_bound = w_bound + 1;
}

void add_red_edges_work(int buffer_hull_edge) {
  while (buffer_hull_edge < h_bound) {
    add_work(h_map[buffer_hull_edge]);
    buffer_hull_edge = buffer_hull_edge + 1;
  }
}

int flip_cond(int p1, int p2, int p3, int p4) {
  int det, x1, y1, x2, y2, x3, y3, x4, y4;

  x1 = p_map[p1].x;
  y1 = p_map[p1].y;
  x2 = p_map[p2].x;
  y2 = p_map[p2].y;
  x3 = p_map[p3].x;
  y3 = p_map[p3].y;
  x4 = p_map[p4].x;
  y4 = p_map[p4].y;

  det = det3(x1, y1, x1 * x1 + y1 * y1, x2, y2, x2 * x2 + y2 * y2, x3, y3, x3 * x3 + y3 * y3) +
    det3(x1, y1, x1 * x1 + y1 * y1, x3, y3, x3 * x3 + y3 * y3, x4, y4, x4 * x4 + y4 * y4) -
   det3(x1, y1, x1 * x1 + y1 * y1, x2, y2, x2 * x2 + y2 * y2, x4, y4, x4 * x4 + y4 * y4) -
    det3(x2, y2, x2 * x2 + y2 * y2, x3, y3, x3 * x3 + y3 * y3, x4, y4, x4 * x4 + y4 * y4);
  if(det > 0) {
    return true;
  } else if (det < 0) {
    return false;
  } else {
    printf("%% %s\n%d %d %d %d %d %d %d %d\n",
           "computed a degenerate incircle condition",
	   x1, y1, x2, y2, x3, y3, x4, y4);
    exit(-1);
  }
}


void make_delaunay() {
  int the_edge, p1, p2, p3, p4, t1, t2, edgeA, edgeB, edgeC, edgeD,
    oA, oB, oC, oD;

  while (0 < w_bound) {
    printf("%% in make_delaunay working stack size %d\n", w_bound);
    dump_state();

    the_edge = work_stack[w_bound-1];

    p1 = e_map[the_edge].p2;
    p3 = e_map[the_edge].p1;

    t1 = e_map[the_edge].t1;
    t2 = e_map[the_edge].t2;
    if (-1 == t2) {
      w_bound = w_bound - 1;
    } else {
      if(the_edge == t_map[t1].e1) {
	oA = t_map[t1].o2;
	edgeA = t_map[t1].e2;
	oB = t_map[t1].o3;
	edgeB = t_map[t1].e3;
      } else if(the_edge == t_map[t1].e2) {
	oA = t_map[t1].o3;
	edgeA = t_map[t1].e3;
	oB = t_map[t1].o1;
	edgeB = t_map[t1].e1;
      } else {
	oA = t_map[t1].o1;
	edgeA = t_map[t1].e1;
	oB = t_map[t1].o2;
	edgeB = t_map[t1].e2;
      }
      if(true == oA) {
	p2 = e_map[edgeA].p2;
      } else {
	p2 = e_map[edgeA].p1;
      }
      if(the_edge == t_map[t2].e1) {
	oC = t_map[t2].o2;
	edgeC = t_map[t2].e2;
	oD = t_map[t2].o3;
	edgeD = t_map[t2].e3;
      } else if (the_edge == t_map[t2].e2) {
	oC = t_map[t2].o3;
	edgeC = t_map[t2].e3;
	oD = t_map[t2].o1;
	edgeD = t_map[t2].e1;
      } else {
	oC = t_map[t2].o1;
	edgeC = t_map[t2].e1;
	oD = t_map[t2].o2;
	edgeD = t_map[t2].e2;
      }
      if(true == oC) {
	p4 = e_map[edgeC].p2;
      } else {
	p4 = e_map[edgeC].p1;
      }
      if(true == flip_cond(p1, p2, p3, p4)) {
	e_map[the_edge].p1 = p4;
	e_map[the_edge].p2 = p2;
	if(true == oC) {
	  e_map[edgeC].t1 = t1;
	} else {
	  e_map[edgeC].t2 = t1;
	}
	if(true == oA) {
	  e_map[edgeA].t1 = t2;
	} else {
	  e_map[edgeA].t2 = t2;
	}
	t_map[t1].e1 = the_edge;
	t_map[t1].o1 = true;
	t_map[t1].e2 = edgeB;
	t_map[t1].o2 = oB;
	t_map[t1].e3 = edgeC;
	t_map[t1].o3 = oC;
	t_map[t2].e1 = the_edge;
	t_map[t2].o1 = false;
	t_map[t2].e2 = edgeD;
	t_map[t2].o2 = oD;
	t_map[t2].e3 = edgeA;
	t_map[t2].o3 = oA;
	w_bound = w_bound - 1;
	add_work(edgeA);
	add_work(edgeB);
	add_work(edgeC);
	add_work(edgeD);
      } else {
	w_bound = w_bound - 1;
      }
    }
  }
}

void add_point_in_triangle(int the_triangle, int the_point) {
  int p1, p2, p3;
  int edge1, edge2, edge3, edge4, edge5, edge6;
  int o1, o2, o3, t2, t3;

  o1 = t_map[the_triangle].o1;
  o2 = t_map[the_triangle].o2;
  o3 = t_map[the_triangle].o3;

  edge4 = t_map[the_triangle].e1;
  edge5 = t_map[the_triangle].e2;
  edge6 = t_map[the_triangle].e3;

  p1 = (true == o1 ? e_map[edge4].p1 : e_map[edge4].p2);
  p2 = (true == o2 ? e_map[edge5].p1 : e_map[edge5].p2);
  p3 = (true == o3 ? e_map[edge6].p1 : e_map[edge6].p2);

  e_map[e_bound].p1 = the_point;
  e_map[e_bound].p2 = p1;
  edge1 = e_bound;
  e_bound = e_bound + 1;
  e_map[e_bound].p1 = the_point;
  e_map[e_bound].p2 = p2;
  edge2 = e_bound;
  e_bound = e_bound + 1;
  e_map[e_bound].p1 = the_point;
  e_map[e_bound].p2 = p3;
  edge3 = e_bound;
  e_bound = e_bound + 1;

  t_map[the_triangle].e1 = edge1;
  t_map[the_triangle].e2 = edge4;
  t_map[the_triangle].e3 = edge2;
  t_map[the_triangle].o1 = true;
  t_map[the_triangle].o2 = o1;
  t_map[the_triangle].o3 = false;

  t_map[t_bound].e1 = edge2;
  t_map[t_bound].e2 = edge5;
  t_map[t_bound].e3 = edge3;
  t_map[t_bound].o1 = true;
  t_map[t_bound].o2 = o2;
  t_map[t_bound].o3 = false;
  t2 = t_bound;
  t_bound = t_bound + 1;

  t_map[t_bound].e1 = edge3;
  t_map[t_bound].e2 = edge6;
  t_map[t_bound].e3 = edge1;
  t_map[t_bound].o1 = true;
  t_map[t_bound].o2 = o3;
  t_map[t_bound].o3 = false;
  t3 = t_bound;
  t_bound = t_bound + 1;

  e_map[edge1].t1 = the_triangle;
  e_map[edge1].t2 = t3;
  e_map[edge2].t1 = t2;
  e_map[edge2].t2 = the_triangle;
  e_map[edge3].t1 = t3;
  e_map[edge3].t2 = t2;

  if(true == o1) {
    e_map[edge4].t1 = the_triangle;
  } else {
    e_map[edge4].t2 = the_triangle;
  }
  if(true == o2) {
    e_map[edge5].t1 = t2;
  } else {
    e_map[edge5].t2 = t2;
  }
  if(true == o3) {
    e_map[edge6].t1 = t3;
  } else {
    e_map[edge6].t2 = t3;
  }
  add_work(edge1);
  add_work(edge2);
  add_work(edge3);
  add_work(edge4);
  add_work(edge5);
  add_work(edge6);
}

int find_triangle(int i) {
  int the_triangle, p1, x1, y1, p2, x2, y2, p3, x3, y3;
  int d1, d2, d3;
  int found = false;
  
  if(0 == t_bound) {
    printf("looking for a triangle in an empty triangle array\n");
    exit(-1);
  } else {
    the_triangle = t_bound - 1;
  }
  while(false == found) {
    if(true == t_map[the_triangle].o1) {
      p1 = e_map[t_map[the_triangle].e1].p1;
      } else {
      p1 = e_map[t_map[the_triangle].e1].p2;
    }
    if(true == t_map[the_triangle].o2) {
      p2 = e_map[t_map[the_triangle].e2].p1;
      } else {
      p2 = e_map[t_map[the_triangle].e2].p2;
    }
    if(true == t_map[the_triangle].o3) {
      p3 = e_map[t_map[the_triangle].e3].p1;
      } else {
      p3 = e_map[t_map[the_triangle].e3].p2;
    }
    x1 = p_map[p1].x; y1 = p_map[p1].y;
    x2 = p_map[p2].x; y2 = p_map[p2].y;
    x3 = p_map[p3].x; y3 = p_map[p3].y;
    d1 = det3(1, x1, y1, 1, x2, y2, 1,
		p_map[i].x, p_map[i].y);
    d2 = det3(1, x2, y2, 1, x3, y3, 1,
		p_map[i].x, p_map[i].y);
    d3 = det3(1, x3, y3, 1, x1, y1, 1,
		p_map[i].x, p_map[i].y);
    if(0 == d1) {
      exit(-1);
    } else if (d1 > 0) {
      if (0 == d2) {
	exit(-1);
      } else if (d2 > 0) {
	if (0 == d3) {
	  exit(-1);
	} else if (d3 > 0) {
	  return the_triangle;
	} else {
	  if(true == t_map[the_triangle].o3) {
	    the_triangle = e_map[t_map[the_triangle].e3].t2;
	  } else {
	    the_triangle = e_map[t_map[the_triangle].e3].t1;
	  }
	}
      } else {
	if(true == t_map[the_triangle].o2) {
	  the_triangle = e_map[t_map[the_triangle].e2].t2;
	} else {
	  the_triangle = e_map[t_map[the_triangle].e2].t1;
	}
      }
    } else {
      if(true == t_map[the_triangle].o1) {
	the_triangle = e_map[t_map[the_triangle].e1].t2;
      } else {
	the_triangle = e_map[t_map[the_triangle].e1].t1;
      }
    }
  }
}


void make_triangles(int i) {
  int the_det, red_edge, buffer_h_bound, the_triangle, orientation;
  int buffer_hull_edge, new_edge1, new_edge0;

  int j = 0;

  if (i < 2) {
    return;
  } else if(2 == i) {
    the_det = det3 (1, p_map[0].x, p_map[0].y,
                    1, p_map[1].x, p_map[1].y,
                    1, p_map[2].x, p_map[2].y);
    if(0 < the_det ) {
      t_map[0].e1 = 0;
      t_map[0].e2 = 1;
      t_map[0].e3 = 2;
      e_map[0].p1 = 0;
      e_map[0].p2 = 1;
      e_map[0].t1 = 0;
      e_map[0].t2 = -1;
      e_map[1].p1 = 1;
      e_map[1].p2 = 2;
      e_map[1].t1 = 0;
      e_map[1].t2 = -1;
      e_map[2].p1 = 2;
      e_map[2].p2 = 0;
      e_map[2].t1 = 0;
      e_map[2].t2 = -1;
      t_map[0].o1 = true;
      t_map[0].o2 = true;
      t_map[0].o3 = true;
      h_map[0] = 0;
      h_map[1] = 1;
      h_map[2] = 2;
      h_bound = 3;
      t_bound = 1;
      e_bound = 3;
    } else if (the_det < 0) {
      t_map[0].e1 = 0;
      t_map[0].e2 = 1;
      t_map[0].e3 = 2;
      e_map[0].p1 = 0;
      e_map[0].p2 = 2;
      e_map[0].t1 = 0;
      e_map[0].t2 = -1;
      e_map[1].p1 = 2;
      e_map[1].p2 = 1;
      e_map[1].t1 = 0;
      e_map[1].t2 = -1;
      e_map[2].p1 = 1;
      e_map[2].p2 = 0;
      e_map[2].t1 = 0;
      e_map[2].t2 = -1;
      t_map[0].o1 = true;
      t_map[0].o2 = true;
      t_map[0].o3 = true;
      h_map[0] = 0;
      h_map[1] = 1;
      h_map[2] = 2;
      t_bound = 1;
      e_bound = 3;
      h_bound = 3;
    } else {
      exit (-1);
    }
  } else {
    printf("%% inside make_triangles\n");
    dump_state();
   red_edge = find_red_edge(i);
   dump_array(h_map, h_bound, "dumping hull");
   printf("%% found red_edge %d\n", red_edge);
   if (-1 == red_edge) {
     the_triangle = find_triangle(i);
     add_point_in_triangle(the_triangle, i);
     make_delaunay();
   } else {
     red_edge = rotate_hull(red_edge, i);
     dump_array(h_map, h_bound, "hull after rotation");
     printf("%% red edge after rotation %d\n", red_edge);
     /* At this point, all edges in h_map from 0 to */
     /* (red_edge-1) are blue */
     /* all edges from red_edge to h_bound-1 are red */
     buffer_hull_edge = red_edge;
     /* initializing a first edge from the first purple point */
     /* to the new point. */
     new_edge1 = e_bound;
     e_bound = e_bound + 1;
     e_map[new_edge1].p1 = e_map[h_map[red_edge]].p1;
     e_map[new_edge1].p2 = i;
     printf("%% making new edge from point %d to point %d\n",
            e_map[new_edge1].p1, e_map[new_edge1].p2);
     /* This edge is created in such a way that the next
        created triangle will be on the good side. */
     orientation = true;
     new_edge0 = new_edge1;
     e_map[new_edge1].t2 = -1;
     /* for each red edge, we add a new triangle, where one
        of the edges is the red edge, another edge is the
        last created edge, and the third one is created at
        this moment. */
     buffer_h_bound = h_bound;
     while(red_edge < buffer_h_bound) {
       t_map[t_bound].e1 = h_map[red_edge];
       t_map[t_bound].e2 = new_edge1;
       t_map[t_bound].o1 = false;
       t_map[t_bound].o2 = orientation;
       e_map[h_map[red_edge]].t2 = t_bound;
       if (true == orientation) {
	 e_map[new_edge1].t1 = t_bound;
       } else {
	 e_map[new_edge1].t2 = t_bound;
       }
       new_edge1 = e_bound;
       add_work(new_edge1);
       e_map[new_edge1].p1 = i;
       e_map[new_edge1].p2 = e_map[h_map[red_edge]].p2;
       e_bound = e_bound + 1;
       /* the new edge is created in such a way that the current
          triangle is on the good side. */
       t_map[t_bound].e3 = new_edge1;
       e_map[new_edge1].t1 = t_bound;
       t_map[t_bound].o3 = true;
       /* But the next triangle (if any) will be on the wrong side */
       orientation = false;
       e_map[new_edge1].t2 = -1;
       t_bound = t_bound + 1;
       red_edge = red_edge + 1;
     }
     /* All the red edges must be added to the working list */
     add_red_edges_work(buffer_hull_edge);
     /* Here we maintain the convex hull */
     h_map[buffer_hull_edge] = new_edge0;
     buffer_hull_edge = buffer_hull_edge + 1;
     h_map[buffer_hull_edge] = new_edge1;
     h_bound = buffer_hull_edge + 1;
     make_delaunay ();
   }
  }
}

void rotate_buffer(int buffer[], int h_bound, int index) {
  int i = 0;
  while(i + index < h_bound) {
    rotate_tmp[i] = buffer[i + index];
    i = i + 1;
  }
  i = 0;
  while(i < index) {
    rotate_tmp[i + h_bound - index] = buffer[i];
    i = i + 1;
  }
  i = 0;
  while(i < h_bound) {
    buffer[i] = rotate_tmp[i];
    i = i + 1;
  }
}

int is_red_edge(struct edge the_edge, int the_point) {
  /* edge color is computed with respect to the point
     at p_map[the_point].  The edge is on the hull, so its
     orientation is taken in the normal direction. */

  int d = det3(1, p_map[the_edge.p1].x, p_map[the_edge.p1].y,
	       1, p_map[the_edge.p2].x, p_map[the_edge.p2].y,
	       1, p_map[the_point].x, p_map[the_point].y);
  return (0 < d ? false : true);
}

int find_red_edge(int the_point) {
  /* We just want to circulate the hull and find
     the index in h_map that designates the first red edge. */
  int result = -1;
  int j = 0;

  while (j < h_bound) {
    if(true == is_red_edge(e_map[h_map[j]], the_point)) {
      return j;
    }
    j++;
  }
  return (-1);
}


int rotate_hull(int red_edge, int the_point) {
  /* assumes that red_edge is the first red_edge in the
     h_map.  Thus red_edge may be 0, and there
     may be more red edges at the end of the map. */

  /* The function must rotate the hull so that all
     blue edges appear before the red edges, the
     returned value must be the index of the first
     red edge after rotation.  This value can be computed
     by dichotomy, but it is probably more efficient
     to just look from from end. */
  int first = red_edge;
  int count = 0;
  int second = first + 1;

  while (true == is_red_edge(e_map[h_map[second]], the_point)
         && second < h_bound) {
    second = second + 1;
  }
  if (second == h_bound) {
    first = red_edge;
  } else {
    rotate_buffer(h_map, h_bound, second);
    first = red_edge + h_bound - second;
  };
  printf("%% looking for first red edge %d\n", first);
  while (true == is_red_edge(e_map[h_map[first]], the_point)) {
    printf("%% looking for first red edge %d\n", first);
    first = first - 1;
  }
  return (first + 1);
}

void display_map_ps() {
  int i;
  for(i = 0; i < e_bound; i++) {
    printf("%d %d moveto %d %d lineto\n",
	   p_map[e_map[i].p1].x,
	   p_map[e_map[i].p1].y,
	   p_map[e_map[i].p2].x,
	   p_map[e_map[i].p2].y);
  }
  printf("stroke\nshowpage\n");
}
int main() {
  printf("%%!PS\n");
  /*
  read_darts();
#ifdef DEBUG
  printf("%% %d darts\n", bound);
#endif
  mkdelaunay();
  add_point_triangle(120, 160, 5);
  printf("%% %d darts after adding a point\n", bound);
  mkdelaunay();
  add_point_out(20, 210);
  dump_map();
  mkdelaunay();
#ifdef DEBUG
  check_map();
#endif
  */
  read_points();
  display_map_ps();
  return 0;
}

/* Coq accessory code:

  Require Import ZArith.
  Open Scope Z_scope.

  Definition d3 (a b c
                 d e f
                 g h i : Z) : Z :=
   a * e * i + d * h * c + g * b * f -
   a * h * f - d * b * i - g * e * c.

Definition ccw_d (p1 p2 p3 : Z * Z) : Z :=
    let (x1, y1) := p1 in
    let (x2, y2) := p2 in
    let (x3, y3) := p3 in
    d3 1 x1 y1 1 x2 y2 1 x3 y3.

Definition incircle (p1 p2 p3 p4 : Z * Z) : Z :=
    let (x1, y1) := p1 in
    let (x2, y2) := p2 in
    let (x3, y3) := p3 in
    let (x4, y4) := p4 in
  d3 x2 y2 (x2^2 + y2^2) x3 y3 (x3^2 + y3^2) x4 y4 (x4^2 + y4^2) -
  d3 x1 y1 (x1^2 + y1^2) x3 y3 (x3^2 + y3^2) x4 y4 (x4^2 + y4^2) +
  d3 x1 y1 (x1^2 + y1^2) x2 y2 (x2^2 + y2^2) x4 y4 (x4^2 + y4^2) -
  d3 x1 y1 (x1^2 + y1^2) x2 y2 (x2^2 + y2^2) x3 y3 (x3^2 + y3^2). 
     
      
 */