计算几何模板

自用计算几何模版

退役了,整理一下自己几何的板子菜队不需要几何选手
或许你可以在我的hexo博客里找到一些证明

一些写在前面的代码

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 50000 + 10;
const double pi = acos(-1.0);
const double inf = DBL_MAX;
const double eps = 1e-12;

判断浮点数的正负

#define zhengfu(x) (((x) > eps) - ((x) < -eps))

• 正数返回1
• 负数返回-1
• \(0\)返回\(0\)

基础类型的表示

写在前面:

• 直线使用两点式
  • 如果是\(ax+by+c=0\)形式,方向向量为\((b,-a)\)
  • 可根据方向向量写出两点
  • 两点也阔以表示一个线段
• 点和向量共用一个\(struct\)
• \(poe\)表示极角,\(poe=atan2(y,x)\),在\(atn2()\)函数中
• 圆用一个点表示圆心,以及一个半径

点的定义

typedef struct point vec;
struct point { //整数点的基本数据结构
    int x, y;
    point(int _x = 0, int _y = 0)
        : x(_x), y(_y) {
    }
    double len() //模长
    {
        return sqrt(1.0 * x * x + y * y);
    }
    vec chuizhi() {
        return vec(-y, x);
    }
    int operator*(const point &i_T) const //点积
    {
        return x * i_T.x + y * i_T.y;
    }
    int operator^(const point &i_T) const //叉积
    {
        return x * i_T.y - y * i_T.x;
    }
    point operator*(int u) const {
        return point(x * u, y * u);
    }
    bool operator==(const point &i_T) const {
        return x == i_T.x && y == i_T.y;
    }
    point operator/(int u) const {
        return point(x / u, y / u);
    }
    point operator+(const point &i_T) {
        return point(x + i_T.x, y + i_T.y);
    }
    point operator-(const point &i_T) {
        return point(x - i_T.x, y - i_T.y);
    }
    friend bool operator<(point a, point b) {
        return a.y == b.y ? a.x < b.x : a.y < b.y;
    }
    friend ostream &operator<<(ostream &out, point &a) {
        //cout << a.x << ' ' << a.y;
        printf("%d %d", a.x, a.y);
        return out;
    }
    friend istream &operator>>(istream &in, point &a) {
        scanf("%d%d", &a.x, &a.y);
        return in;
    }
};
typedef struct point vec;
struct point { //点的基本数据结构
    double x, y;
    double poe;
    point(double _x = 0, double _y = 0)
        : x(_x), y(_y) {
    }
    double len() //模长
    {
        return sqrt(x * x + y * y);
    }
    vec chuizhi() {
        return vec(-y, x);
    }
    double operator*(const point &i_T) const //点积
    {
        return x * i_T.x + y * i_T.y;
    }
    double operator^(const point &i_T) const //叉积
    {
        return x * i_T.y - y * i_T.x;
    }
    point operator*(double u) const {
        return point(x * u, y * u);
    }
    bool operator==(const point &i_T) const {
        return fabs(x - i_T.x) < eps && fabs(y - i_T.y) < eps;
    }
    point operator/(double u) const {
        return point(x / u, y / u);
    }
    point operator+(const point &i_T) {
        return point(x + i_T.x, y + i_T.y);
    }
    point operator-(const point &i_T) {
        return point(x - i_T.x, y - i_T.y);
    }
    friend bool operator<(point a, point b) {
        return fabs(a.y - b.y) < eps ? a.x < b.x : a.y < b.y;
    }
    void atn2() {
        poe = atan2(y, x);
    }
    friend ostream &operator<<(ostream &out, point &a) {
        //cout << a.x << ' ' << a.y;
        printf("%.8f %.8f", a.x, a.y);
        return out;
    }
    friend istream &operator>>(istream &in, point &a) {
        scanf("%lf%lf", &a.x, &a.y);
        return in;
    }
};

圆的定义

struct circle {
    point o;
    double r;
    circle(point _o = point(), double _r = 0.0)
        : r(_r), o(_o) {
    }
    point Point(double t) //圆上任意一点,t表示与x的夹角
    {
        return point(o.x + r * cos(t), o.y + r * sin(t));
    }
    friend istream &operator>>(istream &in, circle &a) {
        in >> a.o >> a.r;
        return in;
    }
    friend ostream &operator<<(ostream &out, circle &a) {
        out << a.o << ' ';
        printf("%.8f", a.r);
        return out;
    }
};

直线的定义

这里\(poe\)代表着直线/线段的倾角

typedef struct Line Ray;	//射线
typedef struct Line Segment;	//线段Segment
struct Line {                //直线
    point a, b;
    double poe;
    Line(point _a = point(), point _b = point())
        : a(_a), b(_b) {
    }
    vec Vec() {
        return (b - a);
    }
    double len() {
        return (b - a).len();
    }
    //与x轴的交点
    double x_len() {
        return b.x - (a.x - b.x) * b.y / (a.y - b.y);
    }
    // 点pi是否在线段ab内
    bool in(point pi) {
        if (b < a)
            return ((b == pi || b < pi) && (a == pi || pi < a));
        return ((a == pi || a < pi) && (b == pi || pi < b));
    }
    bool operator==(Line i_t) {
        return (fabs(Vec() ^ i_t.Vec()) < eps && fabs(x_len() - i_t.x_len()) < eps);
    }
    bool operator<(Line i_t) {
        return fabs(poe - i_t.poe) < eps ? (Vec() ^ (i_t.b - a)) > 0 : poe < i_t.poe;
    }
    friend istream &operator>>(istream &in, Line &a) {
        cin >> a.a >> a.b;
        return in;
    }
    friend ostream &operator<<(ostream &out, Line &a) {
        out << a.a << ' ' << a.b;
        return out;
    }
    void atn2() {
        poe = atan2(b.y - a.y, b.x - a.x);
    }
};

函数

从各处搜集到的方法,以及一些自己写的函数。

杂七杂八，写法证明的话，看心情补。（能用就行）

浮点型GCD

从\(codeforces\)的一场div2上抄的

double gcd(double x, double y) //浮点型gcd
{
    while (fabs(x) > eps && fabs(y) > eps) {
        if (x > y)
            x -= floor(x / y) * y;
        else
            y -= floor(y / x) * x;
    }
    return x + y;
}

浮点型快读

不记得从哪个大佬的板子里找到的,不过一般用不到.

inline double read() //浮点型快读
{
    double num;
    char in;
    double Dec = 0.1;
    bool IsN = false, IsD = false;
    in = getchar();
    while (in != '-' && in != '.' && (in < '0' || in > '9'))
        in = getchar();
    if (in == '-') {
        IsN = true;
        num = 0;
    } else if (in == '.') {
        IsD = true;
        num = 0;
    } else
        num = in - '0';
    if (!IsD) {
        while (in = getchar(), in >= '0' && in <= '9') {
            num *= 10;
            num += in - '0';
        }
    }
    if (in != '.') {
        if (IsN)
            num = -num;
        return num;
    } else {
        while (in = getchar(), in >= '0' && in <= '9') {
            num += Dec * (in - '0');
            Dec *= 0.1;
        }
    }
    if (IsN)
        num = -num;
    return num;
}

比较两个数的大小关系

//这个函数其实完全可以用上面的zhengfu(x-y)代替
int bijiao(double x, double y) {
    if (fabs(x - y) < eps)
        return 0;
    if (x > y)
        return 1;
    return -1;
}
//返回x>=y
int dayu_dengyu(double x, double y) {
    if (fabs(x - y) < eps || x > y)
        return 1;
    return 0;
}

求线段的中垂线

原理很简单,求两点的中点a

然后求个线段的向量v,再求垂直向量n,

第二个点b=a+n

Line zhongchuixian(Segment l) //求线段中垂线
{
    point a = (l.a + l.b) / 2.0;
    vec tp = l.b - l.a;
    return Line(a, a + point(-tp.y, tp.x));
}

点a关于点b的对称点

point zhongxin_duichen(point a, point b) //a为点，b为对称中心
{
    return point(2 * b.x - a.x, 2 * b.y - a.y);
}

求三点形成的三角形外心/三点共圆的圆心

或者直接套两个中垂线函数,然后求交点

point sanjiaoxing_waixin(point a, point b, point c) //三角形外心
{
    double A = a.x * a.x - b.x * b.x + a.y * a.y - b.y * b.y, B = (b.x - a.x) * 2, C = (b.y - a.y) * 2;
    double D = c.x * c.x - b.x * b.x + c.y * c.y - b.y * b.y, E = (b.x - c.x) * 2, F = (b.y - c.y) * 2;
    return point((D * C - A * F) / (B * F - E * C), (D * B - A * E) / (C * E - F * B));
}

求\(\overrightarrow{ab}\)和\(\overrightarrow{ac}\)的叉积

double xuanzhuan(vec a, vec b, vec c) //求三点叉积，正为(b-a)->(c-a)左旋
{
    return (b - a) ^ (c - a);
}

通过叉积进行极角排序

bool cmp(vec a, vec b) //极角排序
{
    vec c(0, 0);
    if ((a - c) ^ (b - c) == 0)
        return a.x < b.x;
    return ((a - c) ^ (b - c)) > 0;
}

求夹角大小

double jiajiao(point a, point b, point c) //a为角的顶点
{
    vec A = b - c, B = a - c, C = b - a;
    return acos((C * C + B * B - A * A) / (B.len() * C.len() * 2));
}
double jiajiao(vec u, vec v) //向量夹角
{
    return acos((u * v) / (u.len() * v.len()));
}

求向量长度

直接用\(struct\)里的\(.len()\)方法也是一样的,只是在某些时候不如这么写方便罢了

double changdu(vec a) //向量长度
{
    return sqrt(a * a);
}

求c在直线ab上的投影点

aoj上对应习题CGL_1_A: Projection

point touying(Line l, point c) //c投影在直线ab上的位置
{
    vec A = l.b - l.a;
    vec B = c - l.a;
    double La = A.len();
    double Lc = (A * B) / (La * La);
    return A * Lc + l.a;
}

将向量按照起点旋转一定角度

vec xiangliang_xuanzhuan(vec a, double pi) //将向量按照起点旋转逆时针pi
{
    double x, y;
    x = cos(pi) * a.x - sin(pi) * a.y;
    y = sin(pi) * a.x + cos(pi) * a.y;
    return vec(x, y);
}

求c关于直线ab的对称的c'

aoj上对应题目CGL_1_B: Reflection

point fanshe(Line l, point c) //求c关于直线ab的对称点c'
{
    point A = touying(l, c);
    return (A - c) * 2.0 + c;
}

判断两线段是否有交点

bool xianduan_xiangjiao(Segment l1, Segment l2) //两线段是否有交点
{
    point a = l1.a, b = l1.b, c = l2.a, d = l2.b;
    if (!(min(a.x, b.x) <= max(c.x, d.x) && min(c.y, d.y) <= max(a.y, b.y) && min(c.x, d.x) <= max(a.x, b.x) && min(a.y, b.y) <= max(c.y, d.y)))
        return 0;
    double u, v, w, z;
    u = (c.x - a.x) * (b.y - a.y) - (b.x - a.x) * (c.y - a.y);
    v = (d.x - a.x) * (b.y - a.y) - (b.x - a.x) * (d.y - a.y);
    w = (a.x - c.x) * (d.y - c.y) - (d.x - c.x) * (a.y - c.y);
    z = (b.x - c.x) * (d.y - c.y) - (d.x - c.x) * (b.y - c.y);
    return (u * v <= eps && w * z <= eps);
}

判断向量平行

bool pingxing(vec a, vec b) //向量平行
{
    // return bijiao(a.x * b.y, a.y * b.x) == 0;
    return zhengfu(a ^ b) == 0;
}

求两个线段的交点

point xianduan_jiaodian(Segment l1, Segment l2) //两线段交点
{
    double tmpLeft, tmpRight, x = inf, y = inf;
    if (xianduan_xiangjiao(l1, l2)) {
        tmpLeft = (l2.b.x - l2.a.x) * (l1.a.y - l1.b.y) - (l1.b.x - l1.a.x) * (l2.a.y - l2.b.y);
        tmpRight = (l1.a.y - l2.a.y) * (l1.b.x - l1.a.x) * (l2.b.x - l2.a.x) + l2.a.x * (l2.b.y - l2.a.y) * (l1.b.x - l1.a.x) - l1.a.x * (l1.b.y - l1.a.y) * (l2.b.x - l2.a.x);
        x = tmpRight / tmpLeft;
        tmpLeft = (l1.a.x - l1.b.x) * (l2.b.y - l2.a.y) - (l1.b.y - l1.a.y) * (l2.a.x - l2.b.x);
        tmpRight = l1.b.y * (l1.a.x - l1.b.x) * (l2.b.y - l2.a.y) + (l2.b.x - l1.b.x) * (l2.b.y - l2.a.y) * (l1.a.y - l1.b.y) - l2.b.y * (l2.a.x - l2.b.x) * (l1.b.y - l1.a.y);
        y = tmpRight / tmpLeft;
    }
    return point(x, y);
}

求三点形成的三角形面积

double mianji(point a, point b, point c) //三点面积
{
    return fabs((b - a) ^ (c - a)) * 0.5;
}

点到线段的距离

double dian_dao_xianduan(Segment l, point c) //点到线段的距离
{
    double L = changdu(l.b - l.a);
    double r = (l.b - l.a) * (c - l.a) / (L * L);
    if (zhengfu(r) == -1) {
        return changdu(c - l.a);
    } else if (dayu_dengyu(r, 1)) {
        return changdu(c - l.b);
    } else {
        double L = r * changdu(l.b - l.a), l2 = changdu(c - l.a);
        return sqrt(l2 * l2 - L * L);
    }
}

两线段距离

double xianduanjuli(Segment l1, Segment l2) //两线段距离
{
    if (xianduan_xiangjiao(l1, l2))
        return 0.0;
    double minn = inf;
    double l = dian_dao_xianduan(l1, l2.a);
    minn = dayu_dengyu(minn, l) ? l : minn;
    l = dian_dao_xianduan(l1, l2.b);
    minn = dayu_dengyu(minn, l) ? l : minn;
    l = dian_dao_xianduan(l2, l1.a);
    minn = dayu_dengyu(minn, l) ? l : minn;
    l = dian_dao_xianduan(l2, l1.b);
    minn = dayu_dengyu(minn, l) ? l : minn;
    return minn;
}

多边形重心

point duobianxing_zhongxin(point po[], int num) //多边形重心
{
    point c;
    double ans = 0;
    po[0] = po[num];
    for (int i = 0; i < num; i++) {
        double tmp = (po[i] ^ po[i + 1]);
        c = c + (po[i] + po[i + 1]) * tmp / 3.0;
        ans += tmp;
    }
    c = c / ans;
    return c;
}

求多边形面积

double duobianxingmianji(vector<point> po) //多边形面积
{
    double ans = 0;
    po.push_back(po[0]);
    for (int i = 0; i < po.size(); i++)
        ans += (po[i] ^ po[i + 1]);
    ans = fabs(ans) * 0.5;
    return ans;
}
double duobianxingmianji(point po[], int num) //多边形面积
{
    double ans = 0;
    po[0] = po[num];
    for (int i = 0; i < num; i++)
        ans += (po[i] ^ po[i + 1]);
    ans = fabs(ans) * 0.5;
    return ans;
}

多边形周长

double duobianxing_zhouchang(point po[], int n) {
    po[0] = po[n];
    double sum = 0;
    for (int i = 0; i < n; i++)
        sum += (po[i] - po[i + 1]).len();
    return sum;
}

判断点是否在多边形内

bool duobianxingnei(point t, point po[], int n) //点在多边形内
{
    int t1, t2, sum, i;
    double f;
    po[0] = po[n];
    for (i = 0; i <= n; i++)
        po[i] = po[i] - t;                                               // 坐标平移
    t1 = po[0].x >= 0 ? (po[0].y >= 0 ? 0 : 3) : (po[0].y >= 0 ? 1 : 2); // 计算象限

    for (sum = 0, i = 1; i <= n; i++) {
        if (fabs(po[i].x) < eps && fabs(po[i].y) < eps)
            break;
        // 被测点为多边形顶点
        f = (po[i] ^ po[i - 1]);

        // 计算叉积
        if (fabs(f) < eps && po[i - 1].x * po[i].x <= 0 && po[i - 1].y * po[i].y <= 0)
            break;                                                           // 点在边上
        t2 = po[i].x >= 0 ? (po[i].y >= 0 ? 0 : 3) : (po[i].y >= 0 ? 1 : 2); // 计算象限
        if (t2 == (t1 + 1) % 4)
            sum += 1;
        // 情况1
        else if (t2 == (t1 + 3) % 4)
            sum -= 1;
        // 情况2
        else if (t2 == (t1 + 2) % 4) // 情况3
            if (f > 0)
                sum += 2;
            else
                sum -= 2;

        t1 = t2;
    }
    bool tf = 0;
    if (i <= n || sum)
        tf = 1;
    for (i = 0; i <= n; i++)
        po[i] = po[i] + t; // 恢复坐标
    return tf;
}

判断点是否在凸包内,包含边界

二分查找

bool dian_in_tubao(vec a, vec p[], int tail) //点在凸包内,包含边界
{
    if ((xuanzhuan(p[1], p[2], a)) < 0 || xuanzhuan(p[1], p[tail], a) > 0) //if ((xuanzhuan(p[1], p[2], a)) <= eps || xuanzhuan(p[1], p[tail], a) >= -eps)
        return false;
    int left = 2, right = tail;
    while (right - left != 1) {
        int mid = (right + left) >> 1;
        if (xuanzhuan(p[1], p[mid], a) >= eps)
            left = mid;
        else
            right = mid;
    }
    return (zhengfu(xuanzhuan(p[left], p[right], a)) > 0); //return (zhengfu(xuanzhuan(p[left], p[right], a)) > eps);
}

点在线段内

bool dian_zai_xianshang(Segment l1, vec c) //点在线段上
{
    int x1 = l1.a.x, x2 = l1.b.x, y1 = l1.a.y, y2 = l1.b.y;
    if (x1 > x2)
        swap(x1, x2);
    if (y1 > y2)
        swap(y1, y2);
    if (xuanzhuan(l1.a, l1.b, c) == 0 && x1 <= c.x && c.x <= x2 && y1 <= c.y && c.y <= y2)
        return 1;
    return 0;
}

点在多边形边上

bool zaibianshang(point &t, point po[], int n) //点在多边形边上否
{
    for (int i = 1; i <= n; i++)
        if (dian_zai_xianshang(Line(po[i], po[(i + 1) % n]), t))
            return 1;
    return 0;
}

向量共线关系

int xiangliang_gongxian(vec a, vec b) { //-1反向，1同向，0不共线
    if (zhengfu(a ^ b) == 0)
        return zhengfu(a * b);
    return 0;
}

求两直线的交点

point zhixian_zhixian_jiaodian(Line l1, Line l2) //两直线交点
{
    double t = ((l1.a.x - l2.a.x) * (l2.a.y - l2.b.y) - (l1.a.y - l2.a.y) * (l2.a.x - l2.b.x)) / ((l1.a.x - l1.b.x) * (l2.a.y - l2.b.y) - (l1.a.y - l1.b.y) * (l2.a.x - l2.b.x));
    return l1.a + (l1.b - l1.a) * t;
}

判断直线和线段是否有交点

int zhixian_xianduan_xiangjiao(Line l1, Segment l2) //直线与线段是否有交点
{
    double c1 = chaji(l1.b - l1.a, l2.a - l1.a), d1 = chaji(l1.b - l1.a, l2.b - l1.a);
    if (zhengfu(c1) == 0 && zhengfu(d1) == 0) { //重合
        return -1;
    } else if (zhengfu(c1 * d1) <= 0) //有交点
        return 1;
    return 0; //平行
}

求直线和线段的交点

point zhixian_xianduan_jiaodian(Line l1, Segment l2) //直线与线段的交点
{
    point ret;
    ret.x = ((l1.b.x - l1.a.x) * (l2.a.y * (l2.b.x - l2.a.x) + l2.a.x * (l2.a.y - l2.b.y)) - (l2.b.x - l2.a.x) * (l1.a.y * (l1.b.x - l1.a.x) + l1.a.x * (l1.a.y - l1.b.y))) / ((l1.b.y - l1.a.y) * (l2.b.x - l2.a.x) - (l2.b.y - l2.a.y) * (l1.b.x - l1.a.x));
    if (l1.a.x != l1.b.x) {
        ret.y = ((l1.b.y - l1.a.y) * ret.x + l1.a.y * (l1.b.x - l1.a.x) + l1.a.x * (l1.a.y - l1.b.y)) / (l1.b.x - l1.a.x);
    } else {
        ret.y = ((l2.b.y - l2.a.y) * ret.x + l2.a.y * (l2.b.x - l2.a.x) + l2.a.x * (l2.a.y - l2.b.y)) / (l2.b.x - l2.a.x);
    }
    return ret;
}

点到直线距离

double zhixian_dian_juli(Line l, point a) //直线到点的距离
{
    return fabs((a - l.a) ^ (l.b - l.a)) / (l.b - l.a).len();
}

线段和点p最近的端点

point xianduan_duandian_dian(point p, Segment l) //线段距离点p最近的端点
{
    point t = p;
    t.x += l.a.y - l.b.y;
    t.y += l.b.x - l.a.x;
    if (((l.a - p) ^ (t - p)) * ((l.b - p) ^ (t - p)) > eps)
        return (p - l.a).len() < (p - l.b).len() ? l.a : l.b;
    return zhixian_zhixian_jiaodian(Line(p, t), l);
}

射线和线段是否有交点

如果有交点,就直接套线段和直线交点求出来即可

bool shexian_xianduan_jiaodian(Ray l1, Segment l2, point &ret) {
    if (zhixian_xianduan_xiangjiao(l1, l2) != 1)
        return false;
    ret = zhixian_xianduan_jiaodian(l1, l2);
    vec A = l1.b - l1.a;
    vec B = ret - l1.a;
    return (A * B > 0);
}

直线与圆交点以及个数

int zhixian_yuan_jiaodian(circle c, Line l, point &p1, point &p2) { //直线与圆交点以及个数
    if (zhixian_dian_juli(l, c.o) > c.r)
        return 0;
    point pi = c.o;
    pi.x += l.a.y - l.b.y;
    pi.y += l.b.x - l.a.x;
    pi = zhixian_zhixian_jiaodian(Line(pi, c.o), l);
    double t = sqrt(c.r * c.r - (pi - c.o).len() * (pi - c.o).len()) / (l.a - l.b).len();
    p1 = pi + (l.b - l.a) * t;
    p2 = pi - (l.b - l.a) * t;
    return 1 + !(p1 == p2);
}
pair<point, point> zhixian_yuan_jiaodian(circle c, Line l) //直线与圆相交的交点
{
    point p = c.o;
    double t;
    p.x += l.a.y - l.b.y;
    p.y += l.b.x - l.a.x;
    p = zhixian_zhixian_jiaodian(Line(p, c.o), l);
    t = sqrt(c.r * c.r - (p - c.o).len() * (p - c.o).len()) / (l.a - l.b).len();
    return make_pair<point, point>(p + (l.b - l.a) * t, p - (l.b - l.a) * t);
}

线段和圆是否相交

int xianduan_yuan_xiangjiao(Segment l, circle a) //圆和线段是否相交
{
    if (bijiao(a.r, zhixian_dian_juli(l, a.o)) == 1) {
        vec d = touying(l, a.o);
        if (zhengfu(((d - a.o) ^ (l.a - a.o)) * ((d - a.o) ^ (l.b - a.o))) <= 0)
            return 1;
        if (bijiao(a.r, (l.b - a.o).len()) || bijiao(a.r, (l.a - a.o).len()))
            return 1;
    }
    return 0;
}

线段与圆交点及个数

int xianduan_yuan_jiaodian(circle a, Segment l, point &p1, point &p2) { //线段与圆交点以及个数
    int tmp = zhixian_yuan_jiaodian(a, l, p1, p2);
    if (tmp == 0)
        return 0;
    if (l.b < l.a)
        swap(l.a, l.b);
    return l.in(p1) + l.in(p2) - (p1 == p2);
}

圆心与ab两点形成的夹角

double yuanxin_liangdian_jiajiao(point a, point b, circle c) //圆心与两点形成的夹角
{
    point _a = touying(Line(c.o, b), a);
    double l = (a - _a).len();
    return asin(l / c.r);
}

求圆和圆的交点

pair<point, point> yuan_yuan_jiaodian(circle a, circle b) //求圆和圆的交点
{
    double l = changdu(a.o - b.o);
    double x = acos((a.r * a.r + l * l - b.r * b.r) / (2.0 * a.r * l));
    double t = atan2((b.o - a.o).y, (b.o - a.o).x);
    return make_pair(a.o + vec(cos(t - x) * a.r, sin(t - x) * a.r), a.o + vec(cos(x + t) * a.r, sin(t + x) * a.r));
}

过一点做圆的切线，求切点

pair<point, point> guodian_yuan_qiedian(circle a, point p) //过一点做圆的切线求切点
{
    double l = changdu(a.o - p);
    double t = asin(a.r / l);
    double lb = l * cos(t);
    vec x = a.o - p;
    x = x / l * lb;
    return make_pair(p + xiangliang_xuanzhuan(x, t), p + xiangliang_xuanzhuan(x, -t));
}

求圆和圆相交面积

double yuan_yuan_xiangjiao(circle a, circle b) { //圆与圆相交面积
    double l = (a.o - b.o).len();
    double ang1 = 2.0 * acos((a.r * a.r + l * l - b.r * b.r) / (2.0 * a.r * l));
    double ang2 = 2.0 * acos((b.r * b.r + l * l - a.r * a.r) / (2.0 * b.r * l));
    double s1 = 0.5 * a.r * a.r * sin(ang1), s2 = 0.5 * b.r * b.r * sin(ang2);
    double h1 = 0.5 * ang1 * a.r * a.r, h2 = 0.5 * ang2 * b.r * b.r;
    return h1 + h2 - s1 - s2;
}

求两圆的外公切线以及个数

int yuanyuan_waigongqiexian(point p, circle a, circle b, Line *l) //求两圆的外公切线以及个数
{
    int cnt = 0;
    point _o = b.o - a.o;
    _o.atn2();
    double base = _o.poe;
    double d = (a.o - b.o).len(), ang = acos((a.r - b.r) / d);
    l[cnt] = Line(a.Point(base - ang), b.Point(base - ang));
    if (zhengfu((a.o - l[cnt].a) ^ (l[cnt].b - l[cnt].a)) == zhengfu((p - l[cnt].a) ^ (l[cnt].b - l[cnt].a))) //p和圆心在切线的同侧
        cnt++;
    l[cnt] = Line(a.Point(base + ang), b.Point(base + ang));
    if (zhengfu((a.o - l[cnt].a) ^ (l[cnt].b - l[cnt].a)) == zhengfu((p - l[cnt].a) ^ (l[cnt].b - l[cnt].a)))
        cnt++;
    return cnt;
}

求两圆公切线以及个数

int yuanyuan_gongqiexian(circle a, circle b, point *u, point *v) //求圆和圆公切线以及切线个数
{
    int cnt = 0;
    if (a.r < b.r) {
        swap(a, b);
        swap(u, v);
    }
    double l = changdu(a.o - b.o);
    double rdiff = a.r - b.r;
    double rsum = a.r + b.r;
    if (zhengfu(l - rdiff) < 0)
        return 0;
    double base = atan2((b.o - a.o).y, (b.o - a.o).x);
    if (zhengfu(l) == 0)
        return -1;
    if (zhengfu(l - rdiff) == 0) {
        u[cnt] = v[cnt] = a.Point(base);
        cnt++;
        return 1;
    }
    double ang = acos((a.r - b.r) / l);
    u[cnt] = a.Point(base + ang);
    v[cnt] = b.Point(base + ang);
    cnt++;
    u[cnt] = a.Point(base - ang);
    v[cnt] = b.Point(base - ang);
    cnt++;
    if (zhengfu(l - rsum) == 0) {
        u[cnt] = v[cnt] = a.Point(base);
        cnt++;
    } else if (zhengfu(l - rsum) > 0) {
        double ang = acos((a.r + b.r) / l);
        u[cnt] = a.Point(base + ang);
        v[cnt] = b.Point(pi + base + ang);
        cnt++;
        u[cnt] = a.Point(base - ang);
        v[cnt] = b.Point(pi + base - ang);
        cnt++;
    }
    return cnt;
}

旋转卡壳法求凸包直径

double tubao_zhijing(int tail) //求出凸包直径，卡壳
{
    double re = 0;
    if (tail == 2) //仅有两个点
        return changdu(q[2] - q[1]);
    q[0] = q[tail];                       //把最后的点放到最前面
    for (int i = 0, j = 2; i < tail; i++) //枚举边
    {
        while (xuanzhuan(q[i], q[i + 1], q[j]) < xuanzhuan(q[i], q[i + 1], q[j + 1]))
            j = (j + 1) % tail;
        re = max(re, max(changdu(q[j] - q[i]), changdu(q[j] - q[i + 1])));
    }
    return re;
}

使用Andrew求凸包

void Andrew(point po[], int n, point qo[], int &tail) //求凸包,<0则包含凸包边上的点,<=0则只输出拐点
{
    sort(po + 1, po + 1 + n);
    tail = 1;
    qo[1] = po[1];
    for (int i = 2; i <= n; i++) {
        while (tail > 1 && xuanzhuan(qo[tail - 1], qo[tail], po[i]) <= 0)
            tail--;
        qo[++tail] = po[i];
    }
    int basic = tail;
    for (int i = n - 1; i >= 1; i--) {
        while (tail > basic && xuanzhuan(qo[tail - 1], qo[tail], po[i]) <= 0)
            tail--;
        qo[++tail] = po[i];
    }
}

分治求平面最近点对

bool cmp(vec a, vec b) {
    return a.y < b.y;
}
double solve(int l, int r) //平面最近点对
{
    if (l == r)
        return inf;
    int mid = (l + r) >> 1;

    double ans = solve(l, mid);
    ans = min(ans, solve(mid + 1, r));

    int tot = 0;
    for (int i = l; i <= r; i++)
        if (fabs(p[mid].x - p[i].x) <= ans)
            temp[tot++] = p[i];

    sort(temp, temp + tot, cmp);

    for (int i = 0; i < tot; i++)
        for (int j = i + 1; j < tot; j++) {
            if (temp[j].y - temp[i].y > ans)
                break;
            ans = min(ans, changdu(temp[j] - temp[i]));
        }
    return ans;
}

多边形和圆相交的面积

将交点与圆心相连,形成三角形,求面积

对三角形进行分类讨论

double distp(Line l) //长度的平方
{
    return (l.a.x - l.b.x) * (l.a.x - l.b.x) + (l.a.y - l.b.y) * (l.a.y - l.b.y);
}
double yuanxin_dian_sanjiao(Line l, circle c) //求圆心与两点所成三角形有向面积
{
    double sign = 1.0;
    l.a = l.a - c.o;
    l.b = l.b - c.o;
    c.o = point(0.0, 0.0);
    if (fabs((l.a - c.o) ^ (l.b - c.o)) < eps)
        return 0.0;
    if (distp(Line(l.a, c.o)) > distp(Line(l.b, c.o))) {
        swap(l.a, l.b);
        sign = -1.0;
    }
    point p1, p2;
    pair<point, point> q = zhixian_yuan_jiaodian(c, l); //oa和ob与圆的交点
    p1 = q.first, p2 = q.second;
    if (distp(Line(l.a, c.o)) < c.r * c.r + eps) {   //a在圆内
        if (distp(Line(l.b, c.o)) < c.r * c.r + eps) //b也在圆内,返回叉积/2
            return ((l.a - c.o) ^ (l.b - c.o)) / 2.0 * sign;
        if ((p1 - l.b).len() > (p2 - l.b).len())
            swap(p1, p2);
        double ret1 = fabs((l.a - c.o) ^ (p1 - c.o));
        double ret2 = acos((p1.x * l.b.x + p1.y * l.b.y) / p1.len() / l.b.len()) * c.r * c.r;
        double ret = (ret1 + ret2) / 2.0;
        if (((l.a - c.o) ^ (l.b - c.o)) < eps && sign > 0.0 || ((l.a - c.o) ^ (l.b - c.o)) > eps && sign < 0.0)
            ret = -ret;
        return ret;
    }
    point ins = xianduan_duandian_dian(c.o, l);
    if (distp(Line(c.o, ins)) > c.r * c.r - eps) {
        double ret = acos((l.a.x * l.b.x + l.a.y * l.b.y) / l.a.len() / l.b.len()) * c.r * c.r / 2.0;
        if (((l.a - c.o) ^ (l.b - c.o)) < eps && sign > 0.0 || ((l.a - c.o) ^ (l.b - c.o)) > eps && sign < 0.0)
            ret = -ret;
        return ret;
    }
    double cm = c.r / ((c.o - l.a).len() - c.r);
    point m = point((c.o.x + cm * l.a.x) / (1 + cm), (c.o.y + cm * l.a.y) / (1 + cm));
    double cn = c.r / ((c.o - l.b).len() - c.r);
    point n = point((c.o.x + cn * l.b.x) / (1 + cn), (c.o.y + cn * l.b.y) / (1 + cn));
    double ret1 = acos((m.x * n.x + m.y * n.y) / m.len() / n.len()) * c.r * c.r;
    double ret2 = acos((p1.x * p2.x + p1.y * p2.y) / p1.len() / p2.len()) * c.r * c.r - fabs((p1 - c.o) ^ (p2 - c.o));
    double ret = (ret1 - ret2) / 2.0;
    if (((l.a - c.o) ^ (l.b - c.o)) < eps && sign > 0.0 || ((l.a - c.o) ^ (l.b - c.o)) > eps && sign < 0.0)
        ret = -ret;
    return ret;
}
double duobianxing_yuan_xiangjiao(circle c, point p[], int n) //多边形与圆相交面积
{
    double sum = 0;
    p[0] = p[n];
    for (int i = 0; i < n; i++)
        sum += yuanxin_dian_sanjiao(Line(p[i], p[i + 1]), c);
    return fabs(sum);
}

最小圆覆盖

circle zuixiaoyuan_fugai(point p[], int n) //最小圆覆盖
{
    random_shuffle(p + 1, p + 1 + n);
    circle c(p[1], 0.0);
    for (int i = 2; i <= n; i++)
        if (zhengfu(changdu(c.o - p[i]) - c.r) > 0) {
            c = circle(p[i], 0.0);
            for (int j = 1; j < i; j++)
                if (zhengfu(changdu(c.o - p[j]) - c.r) > 0) {
                    c.o = (p[i] + p[j]) / 2;
                    c.r = changdu(c.o - p[i]);
                    for (int k = 1; k < j; k++)
                        if (zhengfu(changdu(c.o - p[k]) - c.r) > 0) {
                            c.o = triangle(p[i], p[j], p[k]).waixin();
                            c.r = changdu(c.o - p[i]);
                        }
                }
        }
    return c;
}

辛普森积分法

double fun(double x) //原积分函数
{
    return x;
}
double simpson(double l, double r) //辛普森法则
{
    double mid = (l + r) / 2.0;
    return (fun(l) + fun(r) + 4.0 * fun(mid)) * (r - l) / 6.0;
}
double asr(double l, double r, double eps, double ans) {
    double mid = (l + r) / 2;
    double l_ = simpson(l, mid), r_ = simpson(mid, r);
    if (fabs(l_ + r_ - ans) <= 15 * eps)
        return l_ + r_ + (l_ + r_ - ans) / 15;
    return asr(l, mid, eps / 2, l_) + asr(mid, r, eps / 2, r_);
}
inline double _asr(double l, double r, double eps) {
    return asr(l, r, eps, simpson(l, r));
}

半平面交

bool check(Line a, Line b, Line c) //判断a,b交点是否在c的右边
{
    point o = zhixian_zhixian_jiaodian(a, b);
    return zhengfu(c.Vec() ^ (o - c.a)) < 0;
}
double banpingmianjiao_mianji(Line qu[], int head, int tail) {
    double sum = 0;
    vector<point> q;
    if (tail - head + 1 >= 3) {
        q.push_back(zhixian_zhixian_jiaodian(qu[head], qu[tail]));
        for (int i = head; i < tail; i++)
            q.push_back(zhixian_zhixian_jiaodian(qu[i], qu[i + 1]));
        q.erase(unique(q.begin(), q.end()), q.end());
        sum = duobianxingmianji(q);
    }
    return sum;
}
void banpingmianjiao(int &head, int &tail, Line L[], Line qu[], int t) //求半平面交,q为双端队列
{
    sort(L + 1, L + t + 1);
    head = 1, tail = 0;
    int cnt = 0;
    for (int i = 1; i < t; i++) {
        if (zhengfu(L[i].poe - L[i + 1].poe) == 0)
            continue;
        L[++cnt] = L[i]; //因为排过序，即使极角相同，后面的也比前面的优
    }
    L[++cnt] = L[t];
    for (int i = 1; i <= cnt; i++) {
        while (head < tail && check(qu[tail - 1], qu[tail], L[i]))
            tail--;
        while (head < tail && check(qu[head + 1], qu[head], L[i]))
            head++;
        qu[++tail] = L[i];
    }
    while (head < tail && check(qu[tail - 1], qu[tail], qu[head]))
        tail--;
    while (head < tail && check(qu[head + 1], qu[head], qu[tail]))
        head++;
}

闵可夫斯基和

/**
 * 闵可夫斯基和
 * pl1+pl2=s
 * tail1为pl1的大小
 * tail2为pl2的大小
 * cnt为所得s凸包的大小
 * s[1]==s[cnt]
 * */
void Minkefusiji(point s[], int &cnt, point pl1[], int tail1, point pl2[], int tail2) {
    s[cnt = 1] = pl1[1] + pl2[1];
    pl1[tail1 + 1] = pl1[1];
    pl2[tail2 + 1] = pl2[1];
    for (int i = 1; i <= tail1; i++)
        pl1[i] = pl1[i + 1] - pl1[i];
    for (int i = 1; i <= tail2; i++)
        pl2[i] = pl2[i + 1] - pl2[i];
    int a1 = 1, a2 = 1;
    while (a1 <= tail1 && a2 <= tail2) {
        ++cnt;
        if ((pl1[a1] ^ pl2[a2]) >= 0)
            s[cnt] = s[cnt - 1] + pl1[a1++];
        else
            s[cnt] = s[cnt - 1] + pl2[a2++];
    }
    while (a1 <= tail1)
        ++cnt, s[cnt] = s[cnt - 1] + pl1[a1++];
    while (a2 <= tail2)
        ++cnt, s[cnt] = s[cnt - 1] + pl2[a2++];
}

皮克定理

/**
 * 皮克定理相关函数
 * 线段上整数点
 * 多边形内整数点个数
 * 多边形面积为:S=(2*in+on-2)*0.5
 * */
int gcd(int x, int y) {
    return y == 0 ? x : gcd(y, x % y);
}
int xianduan_shang_zhengdian(Segment l1) //线段内整数个数为,最大公约数
{
    point pl(fabs(l1.b.x - l1.a.x), fabs(l1.b.y - l1.a.y));
    return gcd((int)pl.x, (int)pl.y);
}
//多边形面积直接叉积算
int duobianxing_nei_zhengshu(int s, int on) //内部整数点个数将公式变形即可
{
    return max(s + 2 - on, 0) / 2;
}

欧拉公式

/** 
 * 欧拉公式
 * V-E+F=2
 * V:顶点
 * E:边
 * F:面
 * 假设nn个点都在凸包上，那么V=n，每个面有三条边，每条边被算了两次，即2E=3F
 * 通过上面的公式可以得到F=2n−4,E=3n−6。
 * 增量构造法的复杂度是面数×点数，所以是O(n^2)级别
 * */

反演

/**
 * 反演变换
 * 1.不经过O的直线反演后成为一个经过O的圆
 * 2.经过O的圆，反演后成为不经过O的一条直线
 * 3.不经过O的圆，反演后成为另一个不过O的圆，且两个圆关于O位似
 * 4.过O的直线反演后不变
 * 5.反演不改变相切性
 * */
circle buguo_O_zhixian_fanyan(Line l, point o) //不过中心的直线,反演后为过中心的圆
{
    double r1 = zhixian_dian_juli(l, o);
    double r2 = 1.0 / r1;
    double r = r2 / 2.0;
    vec l1 = touying(l, o) - o;
    point _o = l1 * r / l1.len() + o;
    return circle(_o, r);
}
Line guo_O_zhixian_fanyan(Line l, point o) //过反演中心的直线,是其本身
{
    return l;
}
Line guo_O_yuan_fanyan(circle c, point o) //过反演中心的圆,为一个不过中心的直线
{
    double r2 = c.r * 2.0;
    double r1 = 1.0 / r2;
    vec l1 = (c.o - o) / c.r * r1;
    point _o = o + l1;
    return Line(_o, _o + l1.chuizhi());
}
point buguo_O_dian_fanyan(point c, point o) //不过反演中心的点
{
    vec co = c - o;
    double l1 = co.len();
    double l2 = 1.0 / l1;
    return o + co * (l2 / l1);
}
circle buguo_O_yuan_fanyan(circle c, point o) //不过反演中心的圆,为一个形似的圆
{
    vec co = c.o - o;
    double l1 = co.len();
    double r = 0.5 * (1.0 / (l1 - c.r) - 1.0 / (l1 + c.r));
    double l2 = 1.0 / (l1 + c.r) + r;
    point c2 = o + co * (l2 / l1);
    return circle(c2, r);
}

三分

/**
 * 三分
 * 凹函数
 * judge()目标函数
 * */
double find(double l, double r) {
    while (l + eps < r) {
        double lm = (r - l) / 3.0 + l, rm = r - (r - l) / 3.0;
        if (judge(lm) > judge(rm))
            l = lm;
        else
            r = rm;
    }
    return (l + r) / 2;
}

综合

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 50000 + 10;
const double pi = acos(-1.0);
const double inf = DBL_MAX;
const double eps = 1e-12;
#define zhengfu(x) (((x) > eps) - ((x) < -eps))

/**
 * 直线使用两点式
 * 如果是ax+by+c=0形式,方向向量为(b,-a)
 * 可根据方向向量写出两点
 * 两点也阔以表示一个线段
 * 
 * 点和向量共用一个struct
 * 
 * poe表示极角,poe=atan2(y,x),在ata2函数中
 * 
 * 圆用一个点表示圆心,以及一个半径
 * */
typedef struct point vec;
struct point { //整数点的基本数据结构
    int x, y;
    point(int _x = 0, int _y = 0)
        : x(_x), y(_y) {
    }
    double len() //模长
    {
        return sqrt(1.0 * x * x + y * y);
    }
    vec chuizhi() {
        return vec(-y, x);
    }
    int operator*(const point &i_T) const //点积
    {
        return x * i_T.x + y * i_T.y;
    }
    int operator^(const point &i_T) const //叉积
    {
        return x * i_T.y - y * i_T.x;
    }
    point operator*(int u) const {
        return point(x * u, y * u);
    }
    bool operator==(const point &i_T) const {
        return x == i_T.x && y == i_T.y;
    }
    point operator/(int u) const {
        return point(x / u, y / u);
    }
    point operator+(const point &i_T) {
        return point(x + i_T.x, y + i_T.y);
    }
    point operator-(const point &i_T) {
        return point(x - i_T.x, y - i_T.y);
    }
    friend bool operator<(point a, point b) {
        return a.y == b.y ? a.x < b.x : a.y < b.y;
    }
    friend ostream &operator<<(ostream &out, point &a) {
        //cout << a.x << ' ' << a.y;
        printf("%d %d", a.x, a.y);
        return out;
    }
    friend istream &operator>>(istream &in, point &a) {
        scanf("%d%d", &a.x, &a.y);
        return in;
    }
};
typedef struct point vec;
struct point { //点的基本数据结构
    double x, y;
    double poe;
    point(double _x = 0, double _y = 0)
        : x(_x), y(_y) {
    }
    double len() //模长
    {
        return sqrt(x * x + y * y);
    }
    vec chuizhi() {
        return vec(-y, x);
    }
    double operator*(const point &i_T) const //点积
    {
        return x * i_T.x + y * i_T.y;
    }
    double operator^(const point &i_T) const //叉积
    {
        return x * i_T.y - y * i_T.x;
    }
    point operator*(double u) const {
        return point(x * u, y * u);
    }
    bool operator==(const point &i_T) const {
        return fabs(x - i_T.x) < eps && fabs(y - i_T.y) < eps;
    }
    point operator/(double u) const {
        return point(x / u, y / u);
    }
    point operator+(const point &i_T) {
        return point(x + i_T.x, y + i_T.y);
    }
    point operator-(const point &i_T) {
        return point(x - i_T.x, y - i_T.y);
    }
    friend bool operator<(point a, point b) {
        return fabs(a.y - b.y) < eps ? a.x < b.x : a.y < b.y;
    }
    void atn2() {
        poe = atan2(y, x);
    }
    friend ostream &operator<<(ostream &out, point &a) {
        //cout << a.x << ' ' << a.y;
        printf("%.8f %.8f", a.x, a.y);
        return out;
    }
    friend istream &operator>>(istream &in, point &a) {
        scanf("%lf%lf", &a.x, &a.y);
        return in;
    }
};
struct circle {
    point o;
    double r;
    circle(point _o = point(), double _r = 0.0)
        : r(_r), o(_o) {
    }
    point Point(double t) //圆上任意一点
    {
        return point(o.x + r * cos(t), o.y + r * sin(t));
    }
    friend istream &operator>>(istream &in, circle &a) {
        in >> a.o >> a.r;
        return in;
    }
    friend ostream &operator<<(ostream &out, circle &a) {
        out << a.o << ' ';
        printf("%.8f", a.r);
        return out;
    }
};
typedef struct Line Ray;
typedef struct Line Segment; //线段Segment
struct Line {                //直线
    point a, b;
    double poe;
    Line(point _a = point(), point _b = point())
        : a(_a), b(_b) {
    }
    vec Vec() {
        return (b - a);
    }
    double len() {
        return (b - a).len();
    }
    double x_len() {
        return b.x - (a.x - b.x) * b.y / (a.y - b.y);
    }
    bool in(point pi) {
        if (b < a)
            return ((b == pi || b < pi) && (a == pi || pi < a));
        return ((a == pi || a < pi) && (b == pi || pi < b));
    }
    bool operator==(Line i_t) {
        return (fabs(Vec() ^ i_t.Vec()) < eps && fabs(x_len() - i_t.x_len()) < eps);
    }
    bool operator<(Line i_t) {
        return fabs(poe - i_t.poe) < eps ? (Vec() ^ (i_t.b - a)) > 0 : poe < i_t.poe;
    }
    friend istream &operator>>(istream &in, Line &a) {
        cin >> a.a >> a.b;
        return in;
    }
    friend ostream &operator<<(ostream &out, Line &a) {
        out << a.a << ' ' << a.b;
        return out;
    }
    void atn2() {
        poe = atan2(b.y - a.y, b.x - a.x);
    }
};
double gcd(double x, double y) //浮点型gcd
{
    while (fabs(x) > eps && fabs(y) > eps) {
        if (x > y)
            x -= floor(x / y) * y;
        else
            y -= floor(y / x) * x;
    }
    return x + y;
}
inline double read() //浮点型快读
{
    double num;
    char in;
    double Dec = 0.1;
    bool IsN = false, IsD = false;
    in = getchar();
    while (in != '-' && in != '.' && (in < '0' || in > '9'))
        in = getchar();
    if (in == '-') {
        IsN = true;
        num = 0;
    } else if (in == '.') {
        IsD = true;
        num = 0;
    } else
        num = in - '0';
    if (!IsD) {
        while (in = getchar(), in >= '0' && in <= '9') {
            num *= 10;
            num += in - '0';
        }
    }
    if (in != '.') {
        if (IsN)
            num = -num;
        return num;
    } else {
        while (in = getchar(), in >= '0' && in <= '9') {
            num += Dec * (in - '0');
            Dec *= 0.1;
        }
    }
    if (IsN)
        num = -num;
    return num;
}
int bijiao(double x, double y) {
    if (fabs(x - y) < eps)
        return 0;
    if (x > y)
        return 1;
    return -1;
}
int dayu_dengyu(double x, double y) {
    if (fabs(x - y) < eps || x > y)
        return 1;
    return 0;
}
Line zhongchuixian(Segment l) //求线段中垂线
{
    point a = (l.a + l.b) / 2.0;
    vec tp = l.b - l.a;
    return Line(a, a + point(-tp.y, tp.x));
}
point zhongxin_duichen(point a, point b) //a为点，b为对称中心
{
    return point(2 * b.x - a.x, 2 * b.y - a.y);
}
point sanjiaoxing_waixin(point a, point b, point c) //三角形外心
{
    double A = a.x * a.x - b.x * b.x + a.y * a.y - b.y * b.y, B = (b.x - a.x) * 2, C = (b.y - a.y) * 2;
    double D = c.x * c.x - b.x * b.x + c.y * c.y - b.y * b.y, E = (b.x - c.x) * 2, F = (b.y - c.y) * 2;
    return point((D * C - A * F) / (B * F - E * C), (D * B - A * E) / (C * E - F * B));
}
double xuanzhuan(vec a, vec b, vec c) //求三点叉积，正为(b-a)->(c-a)左旋
{
    return (b - a) ^ (c - a);
}
bool cmp(vec a, vec b) //极角排序
{
    vec c(0, 0);
    if ((a - c) ^ (b - c) == 0)
        return a.x < b.x;
    return ((a - c) ^ (b - c)) > 0;
}
double jiajiao(point a, point b, point c) //a为角的顶点
{
    vec A = b - c, B = a - c, C = b - a;
    return acos((C * C + B * B - A * A) / (B.len() * C.len() * 2));
}
double jiajiao(vec u, vec v) //向量夹角
{
    return acos((u * v) / (u.len() * v.len()));
}
double changdu(vec a) //向量长度
{
    return sqrt(a * a);
}
point touying(Line l, point c) //c投影在直线ab上的位置
{
    vec A = l.b - l.a;
    vec B = c - l.a;
    double La = A.len();
    double Lc = (A * B) / (La * La);
    return A * Lc + l.a;
}
vec xiangliang_xuanzhuan(vec a, double pi) //将向量按照起点旋转逆时针pi,
{
    double x, y;
    x = cos(pi) * a.x - sin(pi) * a.y;
    y = sin(pi) * a.x + cos(pi) * a.y;
    return vec(x, y);
}
point fanshe(Line l, point c) //求c关于直线ab的对称点c'
{
    point A = touying(l, c);
    return (A - c) * 2.0 + c;
}
bool xianduan_xiangjiao(Segment l1, Segment l2) //两线段是否有交点
{
    point a = l1.a, b = l1.b, c = l2.a, d = l2.b;
    if (!(min(a.x, b.x) <= max(c.x, d.x) && min(c.y, d.y) <= max(a.y, b.y) && min(c.x, d.x) <= max(a.x, b.x) && min(a.y, b.y) <= max(c.y, d.y)))
        return 0;
    double u, v, w, z;
    u = (c.x - a.x) * (b.y - a.y) - (b.x - a.x) * (c.y - a.y);
    v = (d.x - a.x) * (b.y - a.y) - (b.x - a.x) * (d.y - a.y);
    w = (a.x - c.x) * (d.y - c.y) - (d.x - c.x) * (a.y - c.y);
    z = (b.x - c.x) * (d.y - c.y) - (d.x - c.x) * (b.y - c.y);
    return (u * v <= eps && w * z <= eps);
}
bool pingxing(vec a, vec b) //向量平行
{
    // return bijiao(a.x * b.y, a.y * b.x) == 0;
    return zhengfu(a ^ b) == 0;
}
point xianduan_jiaodian(Segment l1, Segment l2) //两线段交点
{
    double tmpLeft, tmpRight, x = inf, y = inf;
    if (xianduan_xiangjiao(l1, l2)) {
        tmpLeft = (l2.b.x - l2.a.x) * (l1.a.y - l1.b.y) - (l1.b.x - l1.a.x) * (l2.a.y - l2.b.y);
        tmpRight = (l1.a.y - l2.a.y) * (l1.b.x - l1.a.x) * (l2.b.x - l2.a.x) + l2.a.x * (l2.b.y - l2.a.y) * (l1.b.x - l1.a.x) - l1.a.x * (l1.b.y - l1.a.y) * (l2.b.x - l2.a.x);

        x = tmpRight / tmpLeft;

        tmpLeft = (l1.a.x - l1.b.x) * (l2.b.y - l2.a.y) - (l1.b.y - l1.a.y) * (l2.a.x - l2.b.x);
        tmpRight = l1.b.y * (l1.a.x - l1.b.x) * (l2.b.y - l2.a.y) + (l2.b.x - l1.b.x) * (l2.b.y - l2.a.y) * (l1.a.y - l1.b.y) - l2.b.y * (l2.a.x - l2.b.x) * (l1.b.y - l1.a.y);
        y = tmpRight / tmpLeft;
    }
    return point(x, y);
}
double mianji(point a, point b, point c) //三点面积
{
    return fabs((b - a) ^ (c - a)) * 0.5;
}
double dian_dao_xianduan(Segment l, point c) //点到线段的距离
{
    double L = changdu(l.b - l.a);
    double r = (l.b - l.a) * (c - l.a) / (L * L);

    if (zhengfu(r) == -1) {
        return changdu(c - l.a);
    } else if (dayu_dengyu(r, 1)) {
        return changdu(c - l.b);
    } else {
        double L = r * changdu(l.b - l.a), l2 = changdu(c - l.a);
        return sqrt(l2 * l2 - L * L);
    }
}
double xianduanjuli(Segment l1, Segment l2) //两线段距离
{
    if (xianduan_xiangjiao(l1, l2))
        return 0.0;
    double minn = inf;
    double l = dian_dao_xianduan(l1, l2.a);
    minn = dayu_dengyu(minn, l) ? l : minn;
    l = dian_dao_xianduan(l1, l2.b);
    minn = dayu_dengyu(minn, l) ? l : minn;
    l = dian_dao_xianduan(l2, l1.a);
    minn = dayu_dengyu(minn, l) ? l : minn;
    l = dian_dao_xianduan(l2, l1.b);
    minn = dayu_dengyu(minn, l) ? l : minn;
    return minn;
}
point duobianxing_zhongxin(point po[], int num) //多边形重心
{
    point c;
    double ans = 0;
    po[0] = po[num];
    for (int i = 0; i < num; i++) {
        double tmp = (po[i] ^ po[i + 1]);
        c = c + (po[i] + po[i + 1]) * tmp / 3.0;
        ans += tmp;
    }
    c = c / ans;
    return c;
}
double duobianxingmianji(vector<point> po) //多边形面积
{
    double ans = 0;
    po.push_back(po[0]);
    for (int i = 0; i < po.size(); i++)
        ans += (po[i] ^ po[i + 1]);
    ans = fabs(ans) * 0.5;
    return ans;
}
double duobianxingmianji(point po[], int num) //多边形面积
{
    double ans = 0;
    po[0] = po[num];
    for (int i = 0; i < num; i++)
        ans += (po[i] ^ po[i + 1]);
    ans = fabs(ans) * 0.5;
    return ans;
}
double duobianxing_zhouchang(point po[], int n) {
    po[0] = po[n];
    double sum = 0;
    for (int i = 0; i < n; i++)
        sum += (po[i] - po[i + 1]).len();
    return sum;
}
bool duobianxingnei(point t, point po[], int n) //点在多边形内
{
    int t1, t2, sum, i;
    double f;
    po[0] = po[n];
    for (i = 0; i <= n; i++)
        po[i] = po[i] - t;                                               // 坐标平移
    t1 = po[0].x >= 0 ? (po[0].y >= 0 ? 0 : 3) : (po[0].y >= 0 ? 1 : 2); // 计算象限

    for (sum = 0, i = 1; i <= n; i++) {
        if (fabs(po[i].x) < eps && fabs(po[i].y) < eps)
            break;
        // 被测点为多边形顶点
        f = (po[i] ^ po[i - 1]);

        // 计算叉积
        if (fabs(f) < eps && po[i - 1].x * po[i].x <= 0 && po[i - 1].y * po[i].y <= 0)
            break;                                                           // 点在边上
        t2 = po[i].x >= 0 ? (po[i].y >= 0 ? 0 : 3) : (po[i].y >= 0 ? 1 : 2); // 计算象限
        if (t2 == (t1 + 1) % 4)
            sum += 1;
        // 情况1
        else if (t2 == (t1 + 3) % 4)
            sum -= 1;
        // 情况2
        else if (t2 == (t1 + 2) % 4) // 情况3
            if (f > 0)
                sum += 2;
            else
                sum -= 2;

        t1 = t2;
    }
    bool tf = 0;
    if (i <= n || sum)
        tf = 1;
    for (i = 0; i <= n; i++)
        po[i] = po[i] + t; // 恢复坐标
    return tf;
}
bool dian_in_tubao(vec a, vec p[], int tail) //点在凸包内,包含边界
{
    if ((xuanzhuan(p[1], p[2], a)) < 0 || xuanzhuan(p[1], p[tail], a) > 0) //if ((xuanzhuan(p[1], p[2], a)) <= eps || xuanzhuan(p[1], p[tail], a) >= -eps)
        return false;
    int left = 2, right = tail;
    while (right - left != 1) {
        int mid = (right + left) >> 1;
        if (xuanzhuan(p[1], p[mid], a) >= eps)
            left = mid;
        else
            right = mid;
    }
    return (zhengfu(xuanzhuan(p[left], p[right], a)) > 0); //return (zhengfu(xuanzhuan(p[left], p[right], a)) > eps);
}

bool dian_zai_xianshang(Segment l1, vec c) //点在线段上
{
    int x1 = l1.a.x, x2 = l1.b.x, y1 = l1.a.y, y2 = l1.b.y;
    if (x1 > x2)
        swap(x1, x2);
    if (y1 > y2)
        swap(y1, y2);
    if (xuanzhuan(l1.a, l1.b, c) == 0 && x1 <= c.x && c.x <= x2 && y1 <= c.y && c.y <= y2)
        return 1;
    return 0;
}
bool zaibianshang(point &t, point po[], int n) //点在多边形边上否
{
    for (int i = 1; i <= n; i++)
        if (dian_zai_xianshang(Line(po[i], po[(i + 1) % n]), t))
            return 1;
    return 0;
}
int xiangliang_gongxian(vec a, vec b) { //-1反向，1同向，0不共线
    if (zhengfu((a ^ b)) == 0)
        return zhengfu(a * b);
    return 0;
}
bool shexian_xianduan_jiaodian(Ray l1, Segment l2, point &ret) {
    if (zhixian_xianduan_xiangjiao(l1, l2) != 1)
        return false;
    ret = zhixian_xianduan_jiaodian(l1, l2);
    vec A = l1.b - l1.a;
    vec B = ret - l1.a;
    return (A * B > 0);
}

vec q[MAXN];
double tubao_zhijing(int tail) //求出凸包直径，卡壳
{
    double re = 0;
    if (tail == 2) //仅有两个点
        return changdu(q[2] - q[1]);
    q[0] = q[tail];                       //把最后的点放到最前面
    for (int i = 0, j = 2; i < tail; i++) //枚举边
    {
        while (xuanzhuan(q[i], q[i + 1], q[j]) < xuanzhuan(q[i], q[i + 1], q[j + 1]))
            j = (j + 1) % tail;
        re = max(re, max(changdu(q[j] - q[i]), changdu(q[j] - q[i + 1])));
    }
    return re;
}

void Andrew(point po[], int n, point qo[], int &tail) //求凸包,<0则包含凸包边上的点,<=0则只输出拐点
{
    sort(po + 1, po + 1 + n);
    tail = 1;
    qo[1] = po[1];
    for (int i = 2; i <= n; i++) {
        while (tail > 1 && xuanzhuan(qo[tail - 1], qo[tail], po[i]) <= 0)
            tail--;
        qo[++tail] = po[i];
    }
    int basic = tail;
    for (int i = n - 1; i >= 1; i--) {
        while (tail > basic && xuanzhuan(qo[tail - 1], qo[tail], po[i]) <= 0)
            tail--;
        qo[++tail] = po[i];
    }
}
int zhixian_xianduan_xiangjiao(Line l1, Segment l2) //直线与线段是否有交点
{
    double c1 = chaji(l1.b - l1.a, l2.a - l1.a), d1 = chaji(l1.b - l1.a, l2.b - l1.a);
    if (zhengfu(c1) == 0 && zhengfu(d1) == 0) { //重合
        return -1;
    } else if (zhengfu(c1 * d1) <= 0) //有交点
        return 1;
    return 0; //平行
}
bool cmp(vec a, vec b) {
    return a.y < b.y;
}
double solve(int l, int r) //平面最近点对
{
    if (l == r)
        return inf;
    int mid = (l + r) >> 1;

    double ans = solve(l, mid);
    ans = min(ans, solve(mid + 1, r));

    int tot = 0;
    for (int i = l; i <= r; i++)
        if (fabs(p[mid].x - p[i].x) <= ans)
            temp[tot++] = p[i];

    sort(temp, temp + tot, cmp);

    for (int i = 0; i < tot; i++)
        for (int j = i + 1; j < tot; j++) {
            if (temp[j].y - temp[i].y > ans)
                break;
            ans = min(ans, changdu(temp[j] - temp[i]));
        }
    return ans;
}
point zhixian_xianduan_jiaodian(Line l1, Segment l2) //直线与线段的交点
{
    point ret;
    ret.x = ((l1.b.x - l1.a.x) * (l2.a.y * (l2.b.x - l2.a.x) + l2.a.x * (l2.a.y - l2.b.y)) - (l2.b.x - l2.a.x) * (l1.a.y * (l1.b.x - l1.a.x) + l1.a.x * (l1.a.y - l1.b.y))) / ((l1.b.y - l1.a.y) * (l2.b.x - l2.a.x) - (l2.b.y - l2.a.y) * (l1.b.x - l1.a.x));
    if (l1.a.x != l1.b.x) {
        ret.y = ((l1.b.y - l1.a.y) * ret.x + l1.a.y * (l1.b.x - l1.a.x) + l1.a.x * (l1.a.y - l1.b.y)) / (l1.b.x - l1.a.x);
    } else {
        ret.y = ((l2.b.y - l2.a.y) * ret.x + l2.a.y * (l2.b.x - l2.a.x) + l2.a.x * (l2.a.y - l2.b.y)) / (l2.b.x - l2.a.x);
    }
    return ret;
}
int yuan_yuan_xiangjiao(circle a, circle b) //询问圆和圆的切线个数
{
    double l = changdu(b.o - a.o);
    if (bijiao(l, a.r + b.r) == 1)
        return 4;
    else if (bijiao(l, a.r + b.r) == 0)
        return 3;
    else {
        if (bijiao(l + min(a.r, b.r), max(a.r, b.r)) == 1)
            return 2;
        else if (bijiao(l + min(a.r, b.r), max(a.r, b.r)) == 0)
            return 1;
        else
            return 0;
    }
}
double zhixian_dian_juli(Line l, point a) //直线到点的距离
{
    return fabs((a - l.a) ^ (l.b - l.a)) / (l.b - l.a).len();
}
int zhixian_yuan_jiaodian(circle c, Line l, point &p1, point &p2) { //直线与圆交点以及个数
    if (zhixian_dian_juli(l, c.o) > c.r)
        return 0;
    point pi = c.o;
    pi.x += l.a.y - l.b.y;
    pi.y += l.b.x - l.a.x;
    pi = zhixian_zhixian_jiaodian(Line(pi, c.o), l);
    double t = sqrt(c.r * c.r - (pi - c.o).len() * (pi - c.o).len()) / (l.a - l.b).len();
    p1 = pi + (l.b - l.a) * t;
    p2 = pi - (l.b - l.a) * t;
    return 1 + !(p1 == p2);
}
pair<point, point> zhixian_yuan_jiaodian(circle c, Line l) //直线与圆相交的交点
{
    point p = c.o;
    double t;
    p.x += l.a.y - l.b.y;
    p.y += l.b.x - l.a.x;
    p = zhixian_zhixian_jiaodian(Line(p, c.o), l);
    t = sqrt(c.r * c.r - (p - c.o).len() * (p - c.o).len()) / (l.a - l.b).len();
    return make_pair<point, point>(p + (l.b - l.a) * t, p - (l.b - l.a) * t);
}

int xianduan_yuan_xiangjiao(Segment l, circle a) //圆和线段是否相交
{
    if (bijiao(a.r, zhixian_dian_juli(l, a.o)) == 1) {
        vec d = touying(l, a.o);
        if (zhengfu(((d - a.o) ^ (l.a - a.o)) * ((d - a.o) ^ (l.b - a.o))) <= 0)
            return 1;
        if (bijiao(a.r, (l.b - a.o).len()) || bijiao(a.r, (l.a - a.o).len()))
            return 1;
    }
    return 0;
}
int xianduan_yuan_jiaodian(circle a, Segment l, point &p1, point &p2) { //线段与圆交点以及个数
    int tmp = zhixian_yuan_jiaodian(a, l, p1, p2);
    if (tmp == 0)
        return 0;
    if (l.b < l.a)
        swap(l.a, l.b);
    return l.in(p1) + l.in(p2) - (p1 == p2);
}
double yuanxin_liangdian_jiajiao(point a, point b, circle c) //圆心与两点形成的夹角
{
    point _a = touying(Line(c.o, b), a);
    double l = (a - _a).len();
    return asin(l / c.r);
}
pair<point, point> yuan_yuan_jiaodian(circle a, circle b) //求圆和圆的交点
{
    double l = changdu(a.o - b.o);
    double x = acos((a.r * a.r + l * l - b.r * b.r) / (2.0 * a.r * l));
    double t = atan2((b.o - a.o).y, (b.o - a.o).x);
    return make_pair(a.o + vec(cos(t - x) * a.r, sin(t - x) * a.r), a.o + vec(cos(x + t) * a.r, sin(t + x) * a.r));
}
pair<point, point> guodian_yuan_qiedian(circle a, point p) //过一点做圆的切线求切点
{
    double l = changdu(a.o - p);
    double t = asin(a.r / l);
    double lb = l * cos(t);
    vec x = a.o - p;
    x = x / l * lb;
    return make_pair(p + xiangliang_xuanzhuan(x, t), p + xiangliang_xuanzhuan(x, -t));
}
double yuan_yuan_xiangjiao(circle a, circle b) { //圆与圆相交面积
    double l = (a.o - b.o).len();
    double ang1 = 2.0 * acos((a.r * a.r + l * l - b.r * b.r) / (2.0 * a.r * l));
    double ang2 = 2.0 * acos((b.r * b.r + l * l - a.r * a.r) / (2.0 * b.r * l));
    double s1 = 0.5 * a.r * a.r * sin(ang1), s2 = 0.5 * b.r * b.r * sin(ang2);
    double h1 = 0.5 * ang1 * a.r * a.r, h2 = 0.5 * ang2 * b.r * b.r;
    return h1 + h2 - s1 - s2;
}
int yuanyuan_waigongqiexian(point p, circle a, circle b, Line *l) //求两圆的外公切线以及个数
{
    int cnt = 0;
    point _o = b.o - a.o;
    _o.atn2();
    double base = _o.poe;
    double d = (a.o - b.o).len(), ang = acos((a.r - b.r) / d);
    l[cnt] = Line(a.Point(base - ang), b.Point(base - ang));
    if (zhengfu((a.o - l[cnt].a) ^ (l[cnt].b - l[cnt].a)) == zhengfu((p - l[cnt].a) ^ (l[cnt].b - l[cnt].a))) //p和圆心在切线的同侧
        cnt++;
    l[cnt] = Line(a.Point(base + ang), b.Point(base + ang));
    if (zhengfu((a.o - l[cnt].a) ^ (l[cnt].b - l[cnt].a)) == zhengfu((p - l[cnt].a) ^ (l[cnt].b - l[cnt].a)))
        cnt++;
    return cnt;
}
int yuanyuan_gongqiexian(circle a, circle b, point *u, point *v) //求圆和圆公切线以及切线个数
{
    int cnt = 0;
    if (a.r < b.r) {
        swap(a, b);
        swap(u, v);
    }
    double l = changdu(a.o - b.o);
    double rdiff = a.r - b.r;
    double rsum = a.r + b.r;
    if (zhengfu(l - rdiff) < 0)
        return 0;
    double base = atan2((b.o - a.o).y, (b.o - a.o).x);
    if (zhengfu(l) == 0)
        return -1;
    if (zhengfu(l - rdiff) == 0) {
        u[cnt] = v[cnt] = a.Point(base);
        cnt++;
        return 1;
    }
    double ang = acos((a.r - b.r) / l);
    u[cnt] = a.Point(base + ang);
    v[cnt] = b.Point(base + ang);
    cnt++;
    u[cnt] = a.Point(base - ang);
    v[cnt] = b.Point(base - ang);
    cnt++;
    if (zhengfu(l - rsum) == 0) {
        u[cnt] = v[cnt] = a.Point(base);
        cnt++;
    } else if (zhengfu(l - rsum) > 0) {
        double ang = acos((a.r + b.r) / l);
        u[cnt] = a.Point(base + ang);
        v[cnt] = b.Point(pi + base + ang);
        cnt++;
        u[cnt] = a.Point(base - ang);
        v[cnt] = b.Point(pi + base - ang);
        cnt++;
    }
    return cnt;
}

point zhixian_zhixian_jiaodian(Line l1, Line l2) //两直线交点
{
    double t = ((l1.a.x - l2.a.x) * (l2.a.y - l2.b.y) - (l1.a.y - l2.a.y) * (l2.a.x - l2.b.x)) / ((l1.a.x - l1.b.x) * (l2.a.y - l2.b.y) - (l1.a.y - l1.b.y) * (l2.a.x - l2.b.x));
    return l1.a + (l1.b - l1.a) * t;
}

point xianduan_duandian_dian(point p, Segment l) //线段距离点p最近的端点
{
    point t = p;
    t.x += l.a.y - l.b.y;
    t.y += l.b.x - l.a.x;
    if (((l.a - p) ^ (t - p)) * ((l.b - p) ^ (t - p)) > eps)
        return (p - l.a).len() < (p - l.b).len() ? l.a : l.b;
    return zhixian_zhixian_jiaodian(Line(p, t), l);
}
double distp(Line l) //长度的平方
{
    return (l.a.x - l.b.x) * (l.a.x - l.b.x) + (l.a.y - l.b.y) * (l.a.y - l.b.y);
}
double yuanxin_dian_sanjiao(Line l, circle c) //求圆心与两点所成三角形有向面积
{
    double sign = 1.0;
    l.a = l.a - c.o;
    l.b = l.b - c.o;
    c.o = point(0.0, 0.0);
    if (fabs((l.a - c.o) ^ (l.b - c.o)) < eps)
        return 0.0;
    if (distp(Line(l.a, c.o)) > distp(Line(l.b, c.o))) {
        swap(l.a, l.b);
        sign = -1.0;
    }
    point p1, p2;
    pair<point, point> q = zhixian_yuan_jiaodian(c, l); //oa和ob与圆的交点
    p1 = q.first, p2 = q.second;
    if (distp(Line(l.a, c.o)) < c.r * c.r + eps) {   //a在圆内
        if (distp(Line(l.b, c.o)) < c.r * c.r + eps) //b也在圆内,返回叉积/2
            return ((l.a - c.o) ^ (l.b - c.o)) / 2.0 * sign;
        if ((p1 - l.b).len() > (p2 - l.b).len())
            swap(p1, p2);
        double ret1 = fabs((l.a - c.o) ^ (p1 - c.o));
        double ret2 = acos((p1.x * l.b.x + p1.y * l.b.y) / p1.len() / l.b.len()) * c.r * c.r;
        double ret = (ret1 + ret2) / 2.0;
        if (((l.a - c.o) ^ (l.b - c.o)) < eps && sign > 0.0 || ((l.a - c.o) ^ (l.b - c.o)) > eps && sign < 0.0)
            ret = -ret;
        return ret;
    }
    point ins = xianduan_duandian_dian(c.o, l);
    if (distp(Line(c.o, ins)) > c.r * c.r - eps) {
        double ret = acos((l.a.x * l.b.x + l.a.y * l.b.y) / l.a.len() / l.b.len()) * c.r * c.r / 2.0;
        if (((l.a - c.o) ^ (l.b - c.o)) < eps && sign > 0.0 || ((l.a - c.o) ^ (l.b - c.o)) > eps && sign < 0.0)
            ret = -ret;
        return ret;
    }
    double cm = c.r / ((c.o - l.a).len() - c.r);
    point m = point((c.o.x + cm * l.a.x) / (1 + cm), (c.o.y + cm * l.a.y) / (1 + cm));
    double cn = c.r / ((c.o - l.b).len() - c.r);
    point n = point((c.o.x + cn * l.b.x) / (1 + cn), (c.o.y + cn * l.b.y) / (1 + cn));
    double ret1 = acos((m.x * n.x + m.y * n.y) / m.len() / n.len()) * c.r * c.r;
    double ret2 = acos((p1.x * p2.x + p1.y * p2.y) / p1.len() / p2.len()) * c.r * c.r - fabs((p1 - c.o) ^ (p2 - c.o));
    double ret = (ret1 - ret2) / 2.0;
    if (((l.a - c.o) ^ (l.b - c.o)) < eps && sign > 0.0 || ((l.a - c.o) ^ (l.b - c.o)) > eps && sign < 0.0)
        ret = -ret;
    return ret;
}
double duobianxing_yuan_xiangjiao(circle c, point p[], int n) //多边形与圆相交面积
{
    double sum = 0;
    p[0] = p[n];
    for (int i = 0; i < n; i++)
        sum += yuanxin_dian_sanjiao(Line(p[i], p[i + 1]), c);
    return fabs(sum);
}
circle zuixiaoyuan_fugai(point p[], int n) //最小圆覆盖
{
    random_shuffle(p + 1, p + 1 + n);
    circle c(p[1], 0.0);
    for (int i = 2; i <= n; i++)
        if (zhengfu(changdu(c.o - p[i]) - c.r) > 0) {
            c = circle(p[i], 0.0);
            for (int j = 1; j < i; j++)
                if (zhengfu(changdu(c.o - p[j]) - c.r) > 0) {
                    c.o = (p[i] + p[j]) / 2;
                    c.r = changdu(c.o - p[i]);
                    for (int k = 1; k < j; k++)
                        if (zhengfu(changdu(c.o - p[k]) - c.r) > 0) {
                            c.o = triangle(p[i], p[j], p[k]).waixin();
                            c.r = changdu(c.o - p[i]);
                        }
                }
        }
    return c;
}

/**
 * 以下
 * 皆为
 * 辛普森函数
 * */
double fun(double x) //原积分函数
{
    return x;
}
double simpson(double l, double r) //辛普森法则
{
    double mid = (l + r) / 2.0;
    return (fun(l) + fun(r) + 4.0 * fun(mid)) * (r - l) / 6.0;
}
double asr(double l, double r, double eps, double ans) {
    double mid = (l + r) / 2;
    double l_ = simpson(l, mid), r_ = simpson(mid, r);
    if (fabs(l_ + r_ - ans) <= 15 * eps)
        return l_ + r_ + (l_ + r_ - ans) / 15;
    return asr(l, mid, eps / 2, l_) + asr(mid, r, eps / 2, r_);
}
inline double _asr(double l, double r, double eps) {
    return asr(l, r, eps, simpson(l, r));
}
/**
 * 以下
 * 皆为半平面交
 * */
bool check(Line a, Line b, Line c) //判断a,b交点是否在c的右边
{
    point o = zhixian_zhixian_jiaodian(a, b);
    return zhengfu(c.Vec() ^ (o - c.a)) < 0;
}
double banpingmianjiao_mianji(Line qu[], int head, int tail) {
    double sum = 0;
    vector<point> q;
    if (tail - head + 1 >= 3) {
        q.push_back(zhixian_zhixian_jiaodian(qu[head], qu[tail]));
        for (int i = head; i < tail; i++)
            q.push_back(zhixian_zhixian_jiaodian(qu[i], qu[i + 1]));
        q.erase(unique(q.begin(), q.end()), q.end());
        sum = duobianxingmianji(q);
    }
    return sum;
}
void banpingmianjiao(int &head, int &tail, Line L[], Line qu[], int t) //求半平面交,q为双端队列
{
    sort(L + 1, L + t + 1);
    head = 1, tail = 0;
    int cnt = 0;
    for (int i = 1; i < t; i++) {
        if (zhengfu(L[i].poe - L[i + 1].poe) == 0)
            continue;
        L[++cnt] = L[i]; //因为排过序，即使极角相同，后面的也比前面的优
    }
    L[++cnt] = L[t];

    for (int i = 1; i <= cnt; i++) {
        while (head < tail && check(qu[tail - 1], qu[tail], L[i]))
            tail--;
        while (head < tail && check(qu[head + 1], qu[head], L[i]))
            head++;
        qu[++tail] = L[i];
    }

    while (head < tail && check(qu[tail - 1], qu[tail], qu[head]))
        tail--;
    while (head < tail && check(qu[head + 1], qu[head], qu[tail]))
        head++;
}
/**
 * 闵可夫斯基和
 * pl1+pl2=s
 * tail1为pl1的大小
 * tail2为pl2的大小
 * cnt为所得s凸包的大小
 * s[1]==s[cnt]
 * */
void Minkefusiji(point s[], int &cnt, point pl1[], int tail1, point pl2[], int tail2) {
    s[cnt = 1] = pl1[1] + pl2[1];
    pl1[tail1 + 1] = pl1[1];
    pl2[tail2 + 1] = pl2[1];
    for (int i = 1; i <= tail1; i++)
        pl1[i] = pl1[i + 1] - pl1[i];
    for (int i = 1; i <= tail2; i++)
        pl2[i] = pl2[i + 1] - pl2[i];
    int a1 = 1, a2 = 1;
    while (a1 <= tail1 && a2 <= tail2) {
        ++cnt;
        if ((pl1[a1] ^ pl2[a2]) >= 0)
            s[cnt] = s[cnt - 1] + pl1[a1++];
        else
            s[cnt] = s[cnt - 1] + pl2[a2++];
    }
    while (a1 <= tail1)
        ++cnt, s[cnt] = s[cnt - 1] + pl1[a1++];
    while (a2 <= tail2)
        ++cnt, s[cnt] = s[cnt - 1] + pl2[a2++];
}
/**
 * 皮克定理相关函数
 * 线段上整数点
 * 多边形内整数点个数
 * 多边形面积为:S=(2*in+on-2)*0.5
 * */
int gcd(int x, int y) {
    return y == 0 ? x : gcd(y, x % y);
}
int xianduan_shang_zhengdian(Segment l1) //线段内整数个数为,最大公约数
{
    point pl(fabs(l1.b.x - l1.a.x), fabs(l1.b.y - l1.a.y));
    return gcd((int)pl.x, (int)pl.y);
}
//多边形面积直接叉积算
int duobianxing_nei_zhengshu(int s, int on) //内部整数点个数将公式变形即可
{
    return max(s + 2 - on, 0) / 2;
}
/** 
 * 欧拉公式
 * V-E+F=2
 * V:顶点
 * E:边
 * F:面
 * 假设nn个点都在凸包上，那么V=n，每个面有三条边，每条边被算了两次，即2E=3F
 * 通过上面的公式可以得到F=2n−4,E=3n−6。
 * 增量构造法的复杂度是面数×点数，所以是O(n^2)级别
 * */

/**
 * 反演变换
 * 1.不经过O的直线反演后成为一个经过O的圆
 * 2.经过O的圆，反演后成为不经过O的一条直线
 * 3.不经过O的圆，反演后成为另一个不过O的圆，且两个圆关于O位似
 * 4.过O的直线反演后不变
 * 5.反演不改变相切性
 * */
circle buguo_O_zhixian_fanyan(Line l, point o) //不过中心的直线,反演后为过中心的圆
{
    double r1 = zhixian_dian_juli(l, o);
    double r2 = 1.0 / r1;
    double r = r2 / 2.0;
    vec l1 = touying(l, o) - o;
    point _o = l1 * r / l1.len() + o;
    return circle(_o, r);
}
Line guo_O_zhixian_fanyan(Line l, point o) //过反演中心的直线,是其本身
{
    return l;
}
Line guo_O_yuan_fanyan(circle c, point o) //过反演中心的圆,为一个不过中心的直线
{
    double r2 = c.r * 2.0;
    double r1 = 1.0 / r2;
    vec l1 = (c.o - o) / c.r * r1;
    point _o = o + l1;
    return Line(_o, _o + l1.chuizhi());
}
point buguo_O_dian_fanyan(point c, point o) //不过反演中心的点
{
    vec co = c - o;
    double l1 = co.len();
    double l2 = 1.0 / l1;
    return o + co * (l2 / l1);
}
circle buguo_O_yuan_fanyan(circle c, point o) //不过反演中心的圆,为一个形似的圆
{
    vec co = c.o - o;
    double l1 = co.len();
    double r = 0.5 * (1.0 / (l1 - c.r) - 1.0 / (l1 + c.r));
    double l2 = 1.0 / (l1 + c.r) + r;
    point c2 = o + co * (l2 / l1);
    return circle(c2, r);
}
/**
 * 三分
 * 凹函数
 * judge()目标函数
 * */
double find(double l, double r) {
    while (l + eps < r) {
        double lm = (r - l) / 3.0 + l, rm = r - (r - l) / 3.0;
        if (judge(lm) > judge(rm))
            l = lm;
        else
            r = rm;
    }
    return (l + r) / 2;
}

3D板子

#include <bits/stdc++.h>
using namespace std;

const int MAXN = 50000 + 10;
const double pi = acos(-1.0);
const double inf = 1e100;
const double eps = 1e-12;

#define zero(x) (((x) > 0 ? (x) : -(x)) < eps)
typedef struct point3d vec3d;
struct point3d {
    double x, y, z;
    point3d(double _x = 0.0, double _y = 0.0, double _z = 0.0)
        : x(_x)
        , y(_y)
        , z(_z)
    {
    }
    double len() //模长
    {
        return sqrt(x * x + y * y + z * z);
    }
    double operator*(const vec3d& i_T) const //点积
    {
        return x * i_T.x + y * i_T.y + z * i_T.z;
    }
    vec3d operator*(double u) const
    {
        return vec3d(x * u, y * u, z * u);
    }
    bool operator==(const vec3d& i_T) const
    {
        return fabs(x - i_T.x) < eps && fabs(y - i_T.y) < eps && fabs(z - i_T.z) < eps;
    }
    vec3d operator/(double u) const
    {
        return vec3d(x / u, y / u, z / u);
    }
    vec3d operator^(const vec3d& i_T) //叉积
    {
        return vec3d(y * i_T.z - z * i_T.y, z * i_T.x - x * i_T.z, x * i_T.y - y * i_T.x);
    }
    vec3d operator+(const vec3d& i_T)
    {
        return vec3d(x + i_T.x, y + i_T.y, z + i_T.z);
    }
    vec3d operator-(const point3d& i_T)
    {
        return vec3d(x - i_T.x, y - i_T.y, z - i_T.z);
    }
    friend ostream& operator<<(ostream& out, point3d& a)
    {
        //cout << a.x << ' ' << a.y;
        printf("%.8f %.8f %.8f", a.x, a.y, a.z);
        return out;
    }
    friend istream& operator>>(istream& in, point3d& a)
    {
        scanf("%lf%lf%lf", &a.x, &a.y, &a.z);
        return in;
    }
}p[MAXN];
typedef struct line3d Segment3d;
struct line3d {
    point3d a, b;
    line3d(point3d _a = point3d(), point3d _b = point3d())
        : a(_a)
        , b(_b)
    {
    }
};
struct plane { //3点确定一个平面
    point3d a, b, c; //具体三点
    plane(point3d _a = 0.0, point3d _b = 0.0, point3d _c = 0.0)
        : a(_a)
        , b(_b)
        , c(_c)
    {
    }
    vec3d faxiangliang() //求该平面的法向量
    {
        return ((b - a) ^ (c - a));
    }
};
struct mian {
    int v[3];
    point3d u[3];
    vec3d mj()
    {
        return ((p[v[1]] - p[v[0]]) ^ (p[v[2]] - p[v[0]]));
    }
    bool kejian(point3d as)
    {
        return (as - p[v[0]]) * mj() > 0 ? 1 : 0;
    }
};
int zhengfu(double x)
{
    if (fabs(x) < eps)
        return 0;
    return x > 0 ? 1 : -1;
}
//判三点共线
bool sandian_gongxian(point3d p1, point3d p2, point3d p3)
{
    return zhengfu(((p1 - p2) ^ (p2 - p3)).len()) == 0;
}

//判四点共面
int sidian_gongmian(point3d a, point3d b, point3d c, point3d d)
{
    return zhengfu((plane(a, b, c).faxiangliang() * (d - a))) == 0;
}

//判点是否在线段上,包括端点和共线
int dian_zai_xianduan(point3d p, Segment3d l)
{
    return zero(((p - l.a) ^ (p - l.b)).len()) && (l.a.x - p.x) * (l.b.x - p.x) < eps && (l.a.y - p.y) * (l.b.y - p.y) < eps && (l.a.z - p.z) * (l.b.z - p.z) < eps;
}

//判点是否在线段上,不包括端点
int dian_zai_xianduan(point3d p, Segment3d l)
{
    return dian_zai_xianduan(p, l) && (!zero(p.x - l.a.x) || !zero(p.y - l.a.y) || !zero(p.z - l.a.z)) && (!zero(p.x - l.b.x) || !zero(p.y - l.b.y) || !zero(p.z - l.b.z));
}

//判点是否在空间三角形上,包括边界,三点共线无意义
int dian_zai_sanjiaoxing(point3d p, plane s)
{
    return zero(((s.a - s.b) ^ (s.a - s.c)).len() - ((p - s.a) ^ (p - s.b)).len() - ((p - s.b) ^ (p - s.c)).len() - ((p - s.c) ^ (p - s.a)).len());
}

//判点是否在空间三角形上,不包括边界,三点共线无意义
int dian_zai_sanjiaoxingnei(point3d p, plane s)
{
    return dian_zai_sanjiaoxing(p, s) && ((p - s.a) ^ (p - s.b)).len() > eps && ((p - s.b) ^ (p - s.c)).len() > eps && ((p - s.c) ^ (p - s.a)).len() > eps;
}

//判两点在线段同侧,点在线段上返回0,不共面无意义
int liangdian_pingmian_tongce(point3d p1, point3d p2, Segment3d l)
{
    return (((l.a - l.b) ^ (p1 - l.b)) * ((l.a - l.b) ^ (p2 - l.b))) > eps;
}

//判两点在线段异侧,点在线段上返回0,不共面无意义
int liangdian_xianduan_yice(point3d p1, point3d p2, Segment3d l)
{
    return (((l.a - l.b) ^ (p1 - l.b)) * ((l.a - l.b) ^ (p2 - l.b))) < -eps;
}

//判两点在平面同侧,点在平面上返回0
int liangdian_pingmian_tongce(point3d p1, point3d p2, plane s)
{
    vec3d w = s.faxiangliang();
    return (w * (p1 - s.a)) * (w * (p2 - s.a)) > eps;
}

//判两点在平面异侧,点在平面上返回0
int liangdian_xianduan_yice(point3d p1, point3d p2, plane s)
{
    vec3d w = s.faxiangliang();
    return (w * (p1 - s.a)) * (w * (p2 - s.a)) < -eps;
}

//判两直线平行
int zhixian_pingxing(line3d u, line3d v)
{
    return ((u.a - u.b) ^ (v.a - v.b)).len() < eps;
}

//判两平面平行
int pingmian_pingxing(plane u, plane v)
{
    return (u.faxiangliang() ^ v.faxiangliang()).len() < eps;
}

//判直线与平面平行
int zhixian_pingmian_pingxing(line3d l, plane s)
{
    return zero(((l.a - l.b) * s.faxiangliang()));
}

//判两直线垂直
int zhixian_chuizhi(line3d u, line3d v)
{
    return zero(((u.a - u.b) * (v.a - v.b)));
}

//判两平面垂直
int pingmian_chuizhi(plane u, plane v)
{
    return zero((u.faxiangliang() * v.faxiangliang()));
}

//判直线与平面平行
int zhixian_pingmian_chuizhi(line3d l, plane s)
{
    return ((l.a - l.b) ^ s.faxiangliang()).len() < eps;
}

//判两线段相交,包括端点和部分重合
int xianduan_xiangjiao(Segment3d u, Segment3d v)
{
    if (!sidian_gongmian(u.a, u.b, v.a, v.b))
        return 0;
    if (!sandian_gongxian(u.a, u.b, v.a) || !sandian_gongxian(u.a, u.b, v.b))
        return !liangdian_pingmian_tongce(u.a, u.b, v) && !liangdian_pingmian_tongce(v.a, v.b, u);
    return dian_zai_xianduan(u.a, v) || dian_zai_xianduan(u.b, v) || dian_zai_xianduan(v.a, u) || dian_zai_xianduan(v.b, u);
}

//判两线段相交,不包括端点和部分重合
int xianduan_mofanxiangjiao(Segment3d u, Segment3d v)
{
    return sidian_gongmian(u.a, u.b, v.a, v.b) && liangdian_xianduan_yice(u.a, u.b, v) && liangdian_xianduan_yice(v.a, v.b, u);
}

//判线段与空间三角形相交,包括交于边界和(部分)包含
int xianduan_sanjiaoxing_xiangjiao(Segment3d l, plane s)
{
    return !liangdian_pingmian_tongce(l.a, l.b, s) && !liangdian_pingmian_tongce(s.a, s.b, plane(l.a, l.b, s.c)) && !liangdian_pingmian_tongce(s.b, s.c, plane(l.a, l.b, s.a)) && !liangdian_pingmian_tongce(s.c, s.a, plane(l.a, l.b, s.b));
}

//判线段与空间三角形相交,不包括交于边界和(部分)包含
int xianduan_sanjiaoxing_mofanxiangjiao(Segment3d l, plane s)
{
    return liangdian_xianduan_yice(l.a, l.b, s) && liangdian_xianduan_yice(s.a, s.b, plane(l.a, l.b, s.c)) && liangdian_xianduan_yice(s.b, s.c, plane(l.a, l.b, s.a)) && liangdian_xianduan_yice(s.c, s.a, plane(l.a, l.b, s.b));
}

//计算两直线交点,注意事先判断直线是否共面和平行!
//线段交点请另外判线段相交(同时还是要判断是否平行!)
point3d zhixian_jiaodian(line3d u, line3d v)
{
    point3d ret = u.a;
    double t = ((u.a.x - v.a.x) * (v.a.y - v.b.y) - (u.a.y - v.a.y) * (v.a.x - v.b.x)) / ((u.a.x - u.b.x) * (v.a.y - v.b.y) - (u.a.y - u.b.y) * (v.a.x - v.b.x));
    ret = ret + (u.b - u.a) * t;
    return ret;
}

//计算直线与平面交点,注意事先判断是否平行,并保证三点不共线!
//线段和空间三角形交点请另外判断
point3d zhixian_pingmian_jiaodian(line3d l, plane s)
{
    point3d ret = s.faxiangliang();
    double t = (ret.x * (s.a.x - l.a.x) + ret.y * (s.a.y - l.a.y) + ret.z * (s.a.z - l.a.z)) / (ret.x * (l.b.x - l.a.x) + ret.y * (l.b.y - l.a.y) + ret.z * (l.b.z - l.a.z));
    ret = l.a + (l.b - l.a) * t;
    return ret;
}

//计算两平面交线,注意事先判断是否平行,并保证三点不共线!
line3d pingmian_jiaoxian(plane u, plane v)
{
    line3d ret;
    ret.a = zhixian_pingmian_pingxing(line3d(v.a, v.b), u) ? zhixian_pingmian_jiaodian(line3d(v.b, v.c), u) : zhixian_pingmian_jiaodian(line3d(v.a, v.b), u);
    ret.b = zhixian_pingmian_pingxing(line3d(v.c, v.a), u) ? zhixian_pingmian_jiaodian(line3d(v.b, v.c), u) : zhixian_pingmian_jiaodian(line3d(v.c, v.a), u);
    return ret;
}

//点到直线距离
double dian_zhixian_juli(point3d p, line3d l)
{
    return ((p - l.a) ^ (l.b - l.a)).len() / (l.a - l.b).len();
}

//点到平面距离
double dina_pingmian_juli(point3d p, plane s)
{
    vec3d w = s.faxiangliang();
    return fabs((w * (p - s.a))) / w.len();
}

//直线到直线距离
double zhixian_zhixian_juli(line3d u, line3d v)
{
    point3d n = ((u.a - u.b) ^ (v.a - v.b));
    return fabs(((u.a - v.a) * n)) / n.len();
}

//两直线夹角cos值
double zhixian_cos(line3d u, line3d v)
{
    return ((u.a - u.b) * (v.a - v.b)) / (u.a - u.b).len() / (v.a - v.b).len();
}

//两平面夹角cos值
double pingmian_cos(plane u, plane v)
{
    vec3d uw = u.faxiangliang(), vw = v.faxiangliang();
    return (uw * vw) / (uw.len() * vw.len());
}

//直线平面夹角sin值
double zhixian_pingmian_sin(line3d l, plane s)
{
    vec3d w = s.faxiangliang();
    return ((l.a - l.b) * w) / ((l.a - l.b).len() * w.len());
}

/**
 * 以下为三维凸包
 * 在求解前需对每一个点进行raodong
 * 点集的名称为p
 * */
struct mian { //空间中3点确定的三角形
    int v[3]; //用于三维凸包时的point3d中3点的位置
    vec3d mianji() //用于求三维凸包的表面积
    {
        return ((p[v[0]] - p[v[1]]) ^ (p[v[2]] - p[v[1]]));
    }
    vec3d faxiangliang() //求该三角形所在平面的法向量
    {
        return ((v[1] - v[0]) ^ (v[2] - v[0]));
    }
    bool kejian(point3d as) //用于求三维凸包的函数,判断面是否可见
    {
        return (as - p[v[0]]) * mianji() > 0 ? 1 : 0;
    }
};
void raodong(point3d& a) //蜜汁扰动,使得三点不共线
{
    a = a + ((double)rand() / (double)RAND_MAX - 0.5) * eps * 10;
}
bool vist[MAXN][MAXN];
void sanwei_tubao(point3d* pl, int n, vector<mian>& m) //求解三维凸包
{
    m.push_back(mian{ 1, 2, 3 });
    m.push_back(mian{ 3, 2, 1 });
    for (int i = 4; i <= n; i++) {
        vector<mian> ml;
        for (int j = 0; j < m.size(); j++) {
            bool r = m[j].kejian(p[i]);
            if (!r)
                ml.push_back(m[j]);
            for (int k = 0; k < 3; k++)
                vist[m[j].v[k]][m[j].v[(k + 1) % 3]] = r;
        }
        for (int j = 0; j < m.size(); j++)
            for (int k = 0; k < 3; k++) {
                int a = m[j].v[k], b = m[j].v[(k + 1) % 3];
                if (vist[a][b] != vist[b][a] && vist[a][b])
                    ml.push_back(mian{ a, b, i });
            }
        m = ml;
    }
}
double sanwei_tubao_biaominaji(vector<mian> m) //三维凸包表面积
{
    double ans = 0;
    for (int i = 0; i < m.size(); i++)
        ans += m[i].mianji().len() / 2;
    return ans;
}

/*
模拟退火求最小覆盖球
*/
double solve(int n) {
    int start_T = 100000;
    double step = start_T, ans = inf;
    point3d z;
    int s = 0;
    while (step > eps) {
        double l = (p[s] - z).len();
        for (int i = 0; i < n; i++) {
            double len = (p[i] - z).len();
            if (l < len)
                s = i,l=len;
        }
        ans = min(ans, l);
        z = z + (p[s] - z) / start_T * step;
        step *= 0.97;
    }
    return ans;
}