二维计算几何

#define double long double

const double eps = 1e-12;  //可根据输出要求调整spj
const double inf = 1e9;
const double pi = M_PI;
const int N = 510;

double random_double(int l, int r) {
  int tmp = rand() % (r - l + 1);
  double ans = (tmp + l) * 0.01;
  return ans;
}  //随机生成[l/100,r/100]范围内的实数

double abss(double tmp) { return tmp < 0 ? -tmp : tmp; }  //取绝对值

double Helen(double a, double b, double c) {
  double p = (a + b + c) * 0.5;
  return sqrt(p * (p - a) * (p - b) * (p - c));
}  //海伦公式

double cos(double a, double b, double c) {
  return (a * a + b * b - c * c) / (2 * a * b);
}  //余弦公式

struct point {
  double x, y;

  point() { x = 0, y = 0; }

  point(double p, double q) {
    x = p;
    y = q;
  }

  void rotate_self(double seta) {
    x = x * cos(seta) + y * sin(seta);
    y = y * cos(seta) - x * sin(seta);
  }  //(向量)顺时针自旋seta

  void rotate_in_centernode(point p, double seta) {
    point q(x - p.x, y - p.y);
    q.rotate_self(seta);
    x = q.x + p.x;
    y = q.y + p.y;
  }  //(点)以p为旋转中心顺时针自旋seta

  double length() const { return sqrt(x * x + y * y); }

  point operator-(const point &p) const { return {x - p.x, y - p.y}; }

  point operator+(const point &p) const { return {x + p.x, y + p.y}; }

} orz(0, 0);

point operator*(const double &p, const point q) {
  return {p * q.x, p * q.y};
}  //数乘

double dot_product(point a, point b) { return a.x * b.x + a.y * b.y; }  //点积

double cross_product(point a, point b) { return a.x * b.y - a.y * b.x; }  //叉积

double distance(point p, point q) {
  return sqrt((p.x - q.x) * (p.x - q.x) + (p.y - q.y) * (p.y - q.y));
}  //两点之间距离

struct line {
  double a, b, c;
  //射线模块
  double l, r;
  point direction;
  double polar_angle;

  line() {
    a = 0, b = 0;
    c = 0;
    l = -inf, r = inf;
    polar_angle = 0;
  }

  //表示自变量范围为[l,r]
  line(double aa, double bb, double cc, double ll = -inf, double rr = inf) {
    if (ll != -inf || rr != inf) l = ll, r = rr;  //射线模块开启
    a = aa, b = bb, c = cc;
    polar_angle = 0;
  }

  line(point p, point q, double ll = -inf, double rr = inf) {
    if (ll != -inf || rr != inf) l = ll, r = rr;  //射线模块开启
    direction = q - p;
    polar_angle = atan2(direction.y, direction.x);
    if (abss(p.y - q.y) < eps) {
      a = 0.0, b = 1.0, c = -p.y;
      return;
    }
    // a = random_double(50, 500);//0.50~5.00随机出现一个实数
    a = 1;
    b = a * (q.x - p.x) / (p.y - q.y);
    c = -a * p.x - b * p.y;
  }  //两点确定一条直线

  point get_node(double x) const {
    if (abss(b) <= eps) {  //直线垂直于x轴
      point ans(-c / a, random_double(-200, 200));
      return ans;
    } else {
      point ans(x, -(c + a * x) / b);
      return ans;
    }
  }  //给出一个横坐标x，求直线上对应点
  //若直线垂直于x轴，随机给出一个纵坐标

  bool on_right(point p) const {
    return cross_product(direction, p - get_node(1)) <= eps;
  }  //判断一个点是否在这条直线对应的向量右边

  void rotate_in_centernode(point q, double seta) {
    if (abss(q.x * a + q.y * b + c) <= eps) {  // q是直线上的点
      double a_x = q.x + random_double(300, 1000);
      point tmp = get_node(a_x);
      tmp.rotate_in_centernode(q, seta);
      line(tmp, q);
      return;
    } else {  // q不是直线上的点
      double a_x = random_double(-1000, -100);
      double b_x = random_double(100, 1000);
      point tmp = get_node(a_x);
      point nmp = get_node(b_x);
      tmp.rotate_in_centernode(q, seta);
      nmp.rotate_in_centernode(q, seta);
      line(tmp, nmp);
    }
  }                 //某直线以某点为转轴顺时针旋转seta度
                    //考虑中心点在直线上和不在直线上
} x_axis(0, 1, 0);  // ax+by+c=0

point get_intersection(line p, line q) {
  double xx, yy;
  xx = p.c * q.b - q.c * p.b;
  yy = p.c * q.a - q.c * p.a;
  double tmp = p.a * q.b - q.a * p.b;
  xx /= (-tmp), yy /= tmp;
  point ans(xx, yy);
  return ans;
}  //两条直线交点

line get_vertical(point p, line q) {
  double aa, bb, cc;
  aa = -q.b;
  bb = q.a;
  cc = q.b * p.x - q.a * p.y;
  line z(aa, bb, cc);
  return z;
}  //过一点垂直于已知直线

point get_pedal(point p, line q) {
  return get_intersection(get_vertical(p, q), q);
}  //求点到直线的垂足

line get_parallel(point p, line q) {
  line ans(p.x, p.y, -q.a * p.x - q.b * p.y);
  return ans;
}  //过直线外一点与已知直线平行

line get_vertical_parallel(point p, point q) {
  point mid((p.x + q.x) * 0.5, (p.y + q.y) * 0.5);  //中点
  line pq(p, q);
  return get_vertical(mid, pq);
}  //求平面上两点的中垂线

struct circle {
  point center;
  double r;

  circle(point p, double q) {
    center = p;
    r = q;
  }

  circle(point p, point q, double k) {
    point c = p + (k / (k + 1)) * (q - p);
    point d;
    if (k > 1)
      d = p + (k / (k - 1)) * (q - p);
    else
      d = q + (1 / (1 - k)) * (p - q);
    center = c;
    r = distance(c, d);
  }  //到p的距离是到q的距离的k倍

  circle(point p, point q) {
    center = point((p.x + q.x) * 0.5, (p.y + q.y) * 0.5);
    r = distance(p, q) * 0.5;
  }  //直径

  circle(point p, point q, point s) {
    /*double delta = cross_product(p - q, s - q);
    if (abss(delta) < eps) {
        double ps = distance(p, s);
        double pq = distance(p, q);
        double qs = distance(q, s);
        if (ps >= pq - eps && ps >= qs - eps)circle(p, s);
        else if (pq >= ps - eps && pq >= qs - eps)circle(p, q);
        else circle(q, s);
        return;
    }*/
    line a = get_vertical_parallel(p, q);
    line b = get_vertical_parallel(q, s);
    point c = get_intersection(a, b);
    center = c;
    r = distance(c, p);
  }  //三点共圆

  bool not_in_this_circle(point p) const {
    return distance(p, center) > r + eps;
  }  //不在圆内
};

int rela_between_cri(circle p, circle q) {
  if (p.r < q.r) swap(p, q);
  double dis = distance(p.center, q.center);
  if (abss(dis - p.r - q.r) < eps)
    return -1;  //外切
  else if (abss(dis - p.r + q.r) < eps)
    return 1;  //内切
  else if (dis > p.r + q.r)
    return -2;  //相离
  else if (dis < p.r - q.r)
    return 2;  //内含
  else
    return 0;  //相交
}  //任意两圆之间的关系

int rela_between_line_and_cri(line p, circle q) {
  double dis = distance(q.center, get_pedal(q.center, p));  //点到直线距离
  if (abss(dis - q.r) < eps)
    return 0;  //相切
  else if (dis > q.r)
    return -1;  //相离
  else
    return 1;  //相交
}  //直线和圆之间的关系

point get_tangent_point(line p, circle q) {
  if (!rela_between_line_and_cri(p, q)) return get_pedal(q.center, p);
  return orz;
}  //直线和圆的切点

point get_tangent_point(circle p, circle q) {
  if (p.r < q.r) swap(p, q);
  int x = rela_between_cri(p, q);
  if (abss(x) == 1)
    if (x == 1)
      return (p.r / (p.r + q.r)) * (q.center - p.center) + p.center;  //外切
    else
      return (p.r / (p.r - q.r)) * (q.center - p.center) + p.center;  //内切
  else
    return orz;
}  //圆和圆的切点

line get_tangent_line(circle p, circle q) {
  if (abss(rela_between_cri(p, q)) == 1) {  //确认相切
    line qq(p.center, q.center);
    point pp = get_tangent_point(p, q);
    return get_vertical(pp, qq);
  } else
    return x_axis;
}  //圆和圆的切线

line get_angle_bisector(line p, line q) {
  point dp = (1 / p.direction.length()) * p.direction;
  point dq = (1 / q.direction.length()) * q.direction;
  point st = get_intersection(p, q);
  return {st, st + 0.5 * (dp + dq)};
}  //两条直线的角平分线

double area_of_intersected_cricle(circle p, circle q) {
  if (p.r < q.r) swap(p, q);
  int fl = rela_between_cri(p, q);
  if (!fl) {
    double dis = distance(p.center, q.center);
    double s = Helen(p.r, q.r, dis);
    double alpha = acos(cos(p.r, dis, q.r));
    double beta = acos(cos(q.r, dis, p.r));
    double ans = alpha * p.r * p.r + beta * q.r * q.r - 2 * s;
    return ans;
  } else if (fl < 0)
    return 0;
  else
    return q.r * q.r * pi;
}  //两圆的相交面积

bool cmpyud(point p, point q) {
  return abss(p.y - q.y) <= eps ? p.x < q.x : p.y < q.y;
}  //按照纵坐标从上到下排序
bool cmpydu(point p, point q) {
  return abss(p.y - q.y) <= eps ? p.x > q.x : p.y > q.y;
}  //按照纵坐标从下到上排序
bool cmpxlr(point p, point q) {
  return abss(p.x - q.x) <= eps ? p.y < q.y : p.x < q.x;
}  //横坐标从左到右排序
bool cmpxrl(point p, point q) {
  return abss(p.x - q.x) <= eps ? p.y > q.y : p.x > q.x;
}  //横坐标从右到左排序
bool is_parallel(line p, line q) {
  return abss(p.a * q.b - p.b * q.a) < eps;
}  //两条直线平行
bool cmp_polar_angle(line p, line q) {
  return p.polar_angle < q.polar_angle;
}  //按直线的极角排序

void convex_to_lines(point an[], int n, line bn[], int &num) {
  for (int i = 1; i <= n; i++) bn[++num] = line(an[i], an[i == n ? 1 : i + 1]);
}  //凸多边形点转线

circle get_minimum_circle_coverage(point an[], int num) {
  for (int i = num; i >= 1; i--) swap(an[i], an[rand() % i + 1]);  //随机化
  circle ans(an[1], 0);  //初始化答案圆
  for (int i = 2; i <= num; i++)
    if (ans.not_in_this_circle(an[i])) {
      ans = circle(an[i], 0);
      for (int j = 1; j < i; j++)
        if (ans.not_in_this_circle(an[j])) {
          ans = circle(an[i], an[j]);
          for (int k = 1; k < j; k++)
            if (ans.not_in_this_circle(an[k]))
              ans = circle(an[i], an[j], an[k]);
        }
    }
  return ans;
}  //最小圆覆盖

circle get_inside_cutter(point A, point B, point C) {
  line p = get_angle_bisector(line(A, B), line(A, C));
  line q = get_angle_bisector(line(B, C), line(B, A));
  point c = get_intersection(p, q);
  return {c, distance(c, get_pedal(c, p))};
}  //三角形内切圆

void get_half_convex(point an[], int num, point bn[], int &m) {
  stack<point> s;
  s.push(an[1]);
  s.push(an[2]);  //塞入前两个点
  double p;
  for (int i = 3; i <= num; i++) {
    while (!s.empty()) {
      point tmp = s.top();
      s.pop();
      if (s.empty()) {
        s.push(tmp);
        break;
      }  //如果只有一个点,新点入栈
      p = (an[i].x - s.top().x) * (tmp.y - s.top().y) -
          (tmp.x - s.top().x) * (an[i].y - s.top().y);
      if (p < 0) {
        s.push(tmp);
        break;
      }  //如果新点和top形成的极角小于tmp与top形成的极角，说明新点无法代替top
    }
    s.push(an[i]);
  }
  //现在栈里的点就是凸壳
  while (!s.empty()) {
    bn[++m] = s.top();
    s.pop();
  }
  m--;
}
//求半个凸壳

bool if_intersect(line p, line q) {
  //先判断定义域是否相交
  double ll = max(p.l, q.l);
  double rr = min(p.r, q.r);
  if (rr < ll) return false;  //定义域无交集则不可能相交
  double lp = p.get_node(ll).y, lq = q.get_node(ll).y;
  double rp = p.get_node(rr).y, rq = q.get_node(rr).y;
  if ((lp - lq) * (rp - rq) >= 0)
    return false;
  else
    return true;
}  // p与q是否有交点

int count_intersect_point(line p, point an[], int num) {
  int ans = 0;
  for (int i = 1; i < num; i++) {
    double ll = min(an[i].x, an[i + 1].x), rr = max(an[i].x, an[i + 1].x);
    ans += (int)(if_intersect(p, line(an[i], an[i + 1], ll, rr)));
  }
  double ll = min(an[num].x, an[1].x), rr = max(an[num].x, an[1].x);
  ans += (int)(if_intersect(p, line(an[num], an[1], ll, rr)));
  return ans;
}
// p与多边形an的交点数(考虑射线）
// p为直线一定有偶数个交点

bool in_the_polygon(point p, point an[], int num) {
  point tmp;  //平面上再随机生成一个点
  tmp.x = random_double(-10000, 10000);
  tmp.y = random_double(-10000, 10000);
  line L(p, tmp, p.x);  //生成一条直线（开启射线模块）
  if (count_intersect_point(L, an, num) & 1)
    return true;
  else
    return false;
}  //判断点p在多边形an内

void get_convex(point an[], int num, point bn[], int &m)  //求凸包,bn用来存凸包
{
  sort(an + 1, an + num + 1, cmpxlr);
  get_half_convex(an, num, bn, m);  //按横坐标从左到右得下凸壳
  sort(an + 1, an + num + 1, cmpxrl);
  get_half_convex(an, num, bn, m);  //按横坐标从右到左得上凸壳
}

double get_area_of(point an[], int num) {
  double area = 0.0;
  for (int i = 2; i <= num; i++) area += cross_product(an[i], an[i - 1]);
  area += cross_product(an[1], an[num]);
  area *= 0.5;
  return abss(area);
}  //求凸多边形的面积

double get_length_of(point an[], int num) {
  double ans = 0.0;
  for (int i = 1; i < num; i++) ans += distance(an[i], an[i + 1]);
  ans += distance(an[num], an[1]);
  return ans;
}  //求凸多边形周长

double rotate_jam(point an[], int num) {
  int t = 1;
  double ans = 0.0;
  int nxt[100010];
  for (int i = 1; i <= num; i++) nxt[i] = i + 1;
  nxt[num] = 1;
  for (int i = 1; i <= num; i++) {
    while (true) {
      double dis1 = distance(an[i], an[nxt[t]]);
      double dis2 = distance(an[i], an[nxt[nxt[t]]]);
      double dis = distance(an[i], an[t]);
      if (max(dis1, dis2) > dis) {
        if (dis1 > dis2)
          t = nxt[t];
        else
          t = nxt[nxt[t]];
      } else
        break;
    }
    ans = max(ans, distance(an[i], an[t]));
  }
  return ans;
}  //旋转卡壳求最远点对

int hed = 1, tal = 0;
line que[N];  //求半平面交的时候用到的队列

double half_plane_intersection(line an[], int n, point its[]) {
  sort(an + 1, an + n + 1, cmp_polar_angle);
  hed = 1, tal = 0;
  que[++tal] = an[1];
  for (int i = 2; i <= n; i++) {
    while (hed < tal && an[i].on_right(its[tal - 1])) tal--;
    while (hed < tal && an[i].on_right(its[hed])) hed++;
    que[++tal] = an[i];
    if (que[tal].polar_angle - que[tal - 1].polar_angle <= eps) {
      tal--;
      if (!que[tal].on_right(que[tal + 1].get_node(0))) que[tal] = que[tal + 1];
    }
    if (hed < tal) its[tal - 1] = get_intersection(que[tal], que[tal - 1]);
  }
  while (hed < tal && que[hed].on_right(its[tal - 1])) --tal;
  if (tal - hed <= 1) return 0;
  its[tal] = get_intersection(que[hed], que[tal]);
  double area = 0.0;
  for (int i = hed + 1; i <= tal; i++)
    area += cross_product(its[i], its[i - 1]);
  area += cross_product(its[hed], its[tal]);
  area *= 0.5;
  if (area < 0) area = -area;
  return area;
}
// Half-plane intersection
//逆时针给出n个凸多边形的顶点坐标，求它们交的面积