//reference link：https://www.acwing.com/solution/content/61667/
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;

typedef pair<double, double> PDD;
typedef double db;
const double PI = acos(-1), eps = 1e-8;
const int N = 60;
db R;
int n;
PDD q[N];

PDD operator - (PDD a, PDD b) {return {a.x - b.x, a.y - b.y};}
PDD operator + (PDD a, PDD b) {return {a.x + b.x, a.y + b.y};}
PDD operator * (PDD a, db b) {return {a.x * b, a.y * b};}
PDD operator / (PDD a, db b) {return {a.x / b, a.y / b};}
db operator & (PDD a, PDD b) {return a.x * b.x + a.y * b.y;}
db operator * (PDD a, PDD b) {return a.x * b.y - a.y * b.x;}

int dcmp(db x, db y) {if(fabs(x - y) < eps) return 0; if(x < y) return -1; return 1;}
int sign(db x) {if(fabs(x) < eps) return 0; if(x < 0) return -1; return 1;}
db get_len(PDD a) {return sqrt(a & a);}
PDD norm(PDD a) {return a / get_len(a);}
PDD rotate(PDD a, db b) {return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};}
PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w) {PDD u = p - q; db t = w * u / (v * w); return p + v * t;}
db get_dis(PDD a, PDD b) {db dx = a.x - b.x, dy = a.y - b.y; return sqrt(dx * dx + dy * dy);}
db get_sector_area(PDD a, PDD b) {db c = fabs(acos((a & b) / get_len(a) / get_len(b))); if(sign(a * b) < 0) c = -c;return c * R * R / 2;}
// double get_sector_area(PDD a, PDD b)
// {
//     auto angle = acos((a & b) / get_len(a) / get_len(b));
//     if (sign(a * b) < 0) angle = -angle;
//     return R * R * angle / 2;
// }
bool on_segment(PDD a, PDD b, PDD c) {return !sign((c - a) * (c - b)) && sign((c - a) & (c - b)) <= 0;}//c是否在线段ab上。
db get_line_circle_intersection(PDD a, PDD b, PDD &pa, PDD &pb)//找到圆心到线段的距离
{
    PDD o = {0, 0};
    auto e = get_line_intersection(a, b - a, o, rotate(b - a, PI / 2));
    db mid = get_dis(e, o);
    if(!on_segment(a, b, e)) mid = min(get_dis(o, a), get_dis(o, b));
    if(dcmp(R, mid) <= 0) return mid;
    db len = sqrt(R * R - get_dis(o, e) * get_dis(o, e));
    pa = e + norm(a - b) * len;
    pb = e + norm(b - a) * len;
    return mid;
}
db get_circle_angle_area(PDD a, PDD b)
{
    PDD o = {0, 0};
    db d_a = get_dis(o, a), d_b = get_dis(o, b);
    if(dcmp(R, d_a) >= 0 && dcmp(R, d_b) >= 0) return a * b / 2;//第一种情况：a,b都在圆内 -- 返回三角形面积
    if(!sign(a * b)) return 0;//a, b在一条直线
    PDD pa, pb;
    db mid = get_line_circle_intersection(a, b, pa, pb);//圆心到线段ab最短距离
    if(dcmp(R, mid) <= 0) return get_sector_area(a, b);//第二种情况：a，b都在圆外且有效面积为扇形。
    if(dcmp(R, d_a) >= 0) return get_sector_area(pb, b) + a * pb / 2;//第三种情况：a在里面，b在里面
    if(dcmp(R, d_b) >= 0) return get_sector_area(a, pa) + pa * b / 2;//第四种情况：b在里面，a在外面
    return get_sector_area(a, pa) + pa * pb / 2 + get_sector_area(pb, b);//第五种情况：线段与圆有交点，且a，b都在圆外
}

db work()
{
    db res = 0;
    for(int i = 0; i < n; i ++) 
        res += get_circle_angle_area(q[i], q[(i + 1) % n]);
    return fabs(res);
}

int main()
{
    while(cin >> R >> n)
    {
        for(int i = 0; i < n; i ++) cin >> q[i].x >> q[i].y;
        printf("%.2lf\n", work());
    }

    return 0;
}

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
#define pb push_back
using namespace std;

typedef pair<double, double> PDD;
typedef double db;

const int N = 110;
const db eps = 1e-8, INF = 1e6;
PDD tr[N][3], q[N];
int n;

PDD operator + (PDD a, PDD b) {return {a.x + b.x, a.y + b.y};}
PDD operator - (PDD a, PDD b) {return {a.x - b.x, a.y - b.y};}
PDD operator * (PDD a, db t) {return {a.x * t, a.y * t};}
PDD operator / (PDD a, db t) {return {a.x / t, a.y / t};}
db operator * (PDD a, PDD b) {return a.x * b.y - a.y * b.x;}
db operator & (PDD a, PDD b) {return a.x * b.x + a.y * b.y;}

int dcmp(db x, db y) {if(fabs(x - y) < eps) return 0; if(x < y) return -1; return 1;}
int sign(db x) {if(fabs(x) < eps) return 0; if(x < 0) return -1; return 1;}
bool on_segment(PDD p, PDD a, PDD b) {return sign((p - a) & (p - b)) <= 0;}//若p在ab线段外，则(p - a)与(p - b)点积>0
PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w) 
{ 
    if(!sign(v * w)) return{INF, INF};//两条线重合or平行
    PDD u = p - q; db t = (w * u) / (v * w); 
    PDD o = p + v * t;
    if(!on_segment(o, p, p + v) || !on_segment(o, q, q + w)) return{INF, INF};//交点不在两线段上
    return  o;
}

db get_strength(db a, int side)
{
    int cnt = 0;
    for(int i = 0; i < n; i ++)
    {
        auto t = tr[i];
        if(dcmp(t[0].x, a) > 0 || dcmp(t[2].x, a) < 0) continue;//三角形与直线没有交点
        //三角形一边与直线X=a重合
        if(!dcmp(t[0].x, a) && !dcmp(t[1].x, a))
        {
            if(side) q[cnt ++] = {t[0].y, t[1].y};//这里加入的纵坐标(t[0].y, t[1].y)大小不确定！
        }
        else if(!dcmp(t[1].x, a) && !dcmp(t[2].x, a))
        {
            if(!side) q[cnt ++] = {t[1].y, t[2].y};//t[1].y，t[2].y大小不确定
        }
        else
        {
            db d[3];
            int u = 0;
            for(int j = 0; j < 3; j ++)
            {
                //求三角形与直线X=a的交点 -- 2个or3个交点
                PDD o = get_line_intersection(t[j], t[(j + 1) % 3] - t[j], {a, -INF}, {0, 2 * INF});
                if(dcmp(o.x, INF)) d[u ++] = o.y;//直线与三角形有交点，且交点在三角形的边上
            }
            if(u) sort(d, d + u), q[cnt ++] = {d[0], d[u - 1]};
        }
    }

    if(!cnt) return 0;
    for(int i = 0; i < cnt; i ++) if(q[i].x > q[i].y) swap(q[i].x, q[i].y);
    sort(q, q + cnt);
    db res = 0, st = q[0].x, ed = q[0].y;
    for(int i = 1; i < cnt; i ++)
        if(q[i].x <= ed) ed = max(ed, q[i].y);
        else res += ed - st, ed = q[i].y, st = q[i].x;
    return res + ed - st;
}

db work(db a, db b)
{
    //三角要出现在a的右边or在b的左边
    return (b - a) * (get_strength(a, 1) + get_strength(b, 0)) / 2;
}

int main()
{
    cin >> n;
    vector<db> idx;
    //将顶点横坐标存入
    for(int i = 0; i < n; i ++)
    {
        for(int j = 0; j < 3; j ++) cin >> tr[i][j].x >> tr[i][j].y, idx.pb(tr[i][j].x);
        sort(tr[i], tr[i] + 3);
    }
    //将每条边的交点存入
    for(int i = 0; i < n; i ++)
        for(int j = i + 1; j < n; j ++) 
            for(int x = 0; x < 3; x ++)
                for(int y = 0; y < 3; y ++)
                {
                    PDD o = get_line_intersection(tr[i][x], tr[i][(x + 1) % 3] - tr[i][x],
                                                    tr[j][y], tr[j][(y + 1) % 3] - tr[j][y]);
                    if(dcmp(o.x, INF)) idx.pb(o.x);
                }

    sort(idx.begin(), idx.end());
    db res = 0;
    for(int i = 0; i + 1 < idx.size(); i ++)
        if(dcmp(idx[i], idx[i + 1]))
            res += work(idx[i], idx[i + 1]);
    printf("%.2lf", res);

    return 0;
}

/*
该算法瓶颈在于区间合并:
    每遍历一个区间都要遍历所有矩形求覆盖区间 -- O(n ^ 2)
    当n = 1e5时，要用线段树来维护覆盖区间 -- O(n*lgn)
*/
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath> 
#include <vector>
#define x first
#define y second
using namespace std;
typedef pair<double, double> PDD;
typedef pair<int, int> PII;
typedef long long LL;
const int N = 1010;

PII l[N], r[N], q[N];
int n;

LL get_strength(int L, int R)
{
    int cnt = 0;
    for(int i = 0; i < n; i ++)
        if(R <= l[i].x || r[i].x <= L) continue;//该矩形与区间没有交集
        else q[cnt ++] = {l[i].y, r[i].y};//有交集，将该矩形的纵坐标加入

    //区间合并
    if(!cnt) return 0;
    sort(q, q + cnt);
    LL res = 0, st = q[0].x, ed = q[0].y;
    for(int i = 0; i < cnt; i ++)
        if(q[i].x <= ed) ed = max(ed, 1ll * q[i].y);
        else res += ed - st, ed = q[i].y, st = q[i].x;

    return res + ed - st;
}

LL work(int l, int r)
{
    return 1ll * (r - l) * get_strength(l, r);
}

int main()
{
    cin >> n;
    vector<int> idx(n * 2);
    int cnt = 0;
    for(int i = 0; i < n; i ++) 
    {
        cin >> l[i].x >> l[i].y >> r[i].x >> r[i].y;
        idx[cnt ++] = l[i].x, idx[cnt ++] = r[i].x;
    }

    LL res = 0;
    sort(idx.begin(), idx.begin() + cnt);
    for(int i = 0; i + 1 < cnt; i ++ )
    {
        if(idx[i] != idx[i + 1])
            res += work(idx[i], idx[i + 1]);
    }

    cout << res << "\n";

    return 0;
}

当数据较离散时，会导致错误；（及解决方法）

原题Code：

/*
自适应辛普森积分法：
    令f(x0):直线X = x0与所有圆的长度并
    将f(X)在[-INF ~ INF]上面积分即为：所有圆的面积并


    //注：  当区间比较离散时，运气不好的话（感觉第一次的时候很有可能）l和r的x轴都没与圆有交集，
            mid也与圆没有交集，那直接就等于0<eps了,直接G了。
            这里可以打个补丁：  可以先将所有圆投影到x轴上，然后做个区间合并，
                                然后对每个区间分别做辛普森。
*/
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;
typedef pair<double, double> PDD;
typedef double db;
const int N = 1010;
const db eps = 1e-8;

int n;

struct Circle
{
    PDD o;
    db r;
}circle[N];
PDD q[2010];

int dcmp(db x, db y)
{
    if(fabs(x - y) < eps) return 0;
    if(x < y) return -1;
    return 1;
}

db f(db x)
{
    int cnt = 0;
    for(int i = 0; i < n; i ++)
    {
        db a = circle[i].o.x, b = circle[i].o.y, R = circle[i].r;
        db X = fabs(a - x);
        if(dcmp(X, R) < 0)
        {
            db Y = sqrt(R * R - X * X);
            q[cnt ++] = {b - Y, b + Y};
        }

    }
    if(!cnt) return 0;

    db res = 0;
    sort(q, q + cnt);
    db st = q[0].x, ed = q[0].y;
    for(int i = 1; i < cnt; i ++)
        if(q[i].x <= ed) ed = max(ed, q[i].y);
        else 
        {
            res += ed - st;
            st = q[i].x, ed = q[i].y;
        }
    return res + ed - st;
}

db simpson(db l, db r)
{
    db mid = (l + r) / 2;
    return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}

db asr(db l, db r, db s)
{
    db mid = (l + r) / 2;
    db left = simpson(l, mid), right = simpson(mid, r);
    if(fabs(left + right - s) < eps) return right + left;
    return asr(l, mid, left) + asr(mid, r, right);
}

int main()
{
    scanf("%d", &n);
    db l = 2000, r = -2000;
    for(int i = 0; i < n; i ++)
    {
        scanf("%lf %lf %lf", &circle[i].o.x, &circle[i].o.y, &circle[i].r);
        l = min(l, circle[i].o.x - circle[i].r), r = max(r, circle[i].o.x + circle[i].r);
    }

    printf("%.3lf\n", asr(l - 100, r + 100, simpson(l, r)));

    return 0;
}

数据较为离散时：

这个数据就会出错：

2
53 -154 2.2203502244049247
351 -459 0

解决方法：

当区间比较离散时，运气不好的话（感觉第一次的时候很有可能）l和r的x轴都没与圆有交集，
mid也与圆没有交集，那直接就等于0<eps了,直接G了。
这里可以打个补丁： 可以先将所有圆投影到x轴上，然后做个区间合并，
然后对每个区间分别做辛普森。

改进Code：（但是$eps>1e-8$会wa最后一个点，有0.001误差；$eps <= 1e-9$会T，呜呜呜）

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;
typedef pair<long double, long double> PDD;
typedef long double db;
const int N = 1010;
const db eps = 1e-8;

int n;

struct Circle
{
    PDD o;
    db r;
}circle[N];
PDD q[2010];

int dcmp(db x, db y)
{
    if(fabs(x - y) < eps) return 0;
    if(x < y) return -1;
    return 1;
}

db f(db x)
{
    int cnt = 0;
    for(int i = 0; i < n; i ++)
    {
        db a = circle[i].o.x, b = circle[i].o.y, R = circle[i].r;
        db X = fabs(a - x);
        if(dcmp(X, R) < 0)
        {
            db Y = sqrt(R * R - X * X);
            q[cnt ++] = {b - Y, b + Y};
        }

    }
    if(!cnt) return 0;

    db res = 0;
    sort(q, q + cnt);
    db st = q[0].x, ed = q[0].y;
    for(int i = 1; i < cnt; i ++)
        if(q[i].x <= ed) ed = max(ed, q[i].y);
        else 
        {
            res += ed - st;
            st = q[i].x, ed = q[i].y;
        }
    return res + ed - st;
}

db simpson(db l, db r)
{
    db mid = (l + r) / 2;
    return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}

db asr(db l, db r, db s)
{
    db mid = (l + r) / 2;
    db left = simpson(l, mid), right = simpson(mid, r);
    if(fabs(left + right - s) < eps) return right + left;
    return asr(l, mid, left) + asr(mid, r, right);
}

int main()
{
    scanf("%d", &n);
    db l = 2000, r = -2000;
    PDD p[2010];
    int cnt = 0;
    for(int i = 0; i < n; i ++)
    {
        scanf("%lf %lf %lf", &circle[i].o.x, &circle[i].o.y, &circle[i].r);
        p[cnt ++] = {circle[i].o.x - circle[i].r, circle[i].o.x + circle[i].r};
    }

    db res = 0;
    sort(p, p + cnt);
    db st = p[0].x, ed = p[0].y;
    for(int i = 1; i < cnt; i ++) 
        if(ed >= p[i].x) ed = max(ed, p[i].y);
        else 
        {
            res += asr(st, ed, simpson(st, ed));
            st = q[i].x, ed = q[i].y;
        }
    if(st != ed) res += asr(st, ed, simpson(st, ed));

    printf("%.3lf\n", res);

    return 0;
}

/*
自适应辛普森积分法：
    求f(X)在区间[l ~ r]的积分：
        1. mid = (l + r) / 2。分别求：
            S  = f(X)在[l ~ r]的积分
            S1 = f(X)在[l ~ mid]的积分
            S2 = f(X)在[mid ~ r]的积分
        2. 若|S1 + S2 - S| < eps，说明所求answer精度已经达到，直接返回S1+S2
        3. 否则递归到[l ~ mid]和[mid ~ r]进行上述步骤，知道达到精度。
*/
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;

typedef pair<double, double> PDD;
typedef double db;
const db eps = 1e-12;

db f(db x)
{
    return sin(x) / x;
}

db simpson(db l, db r)
{
    db mid = (l + r) / 2;
    return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}

db asr(db l, db r, db s)
{
    db mid = (l + r) / 2;
    db left = simpson(l, mid), right = simpson(mid, r);
    if(fabs(left + right - s) < eps) return left + right;
    return asr(l, mid, left) + asr(mid, r, right);
}

int main()
{
    db l, r;
    cin >> l >> r;
    printf("%.6lf\n", asr(l, r, simpson(l, r)));

    return 0;
}

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;

typedef double db;
typedef pair<double, double> PDD;

const int N = 5e4 + 10;
const db PI = acos(-1), eps = 1e-12, INF = 1e20;

PDD q[N], ans[N];
int n, stk[N], top;

int sign(db x) {if(fabs(x) < eps) return 0; if(x < 0) return -1; return 1;}
int dcmp(db x, db y){if(fabs(x - y) < eps) return 0; if(x < y) return -1; return 1;}
PDD operator - (PDD a, PDD b) {return {a.x - b.x, a.y - b.y};}
PDD operator + (PDD a, PDD b) {return {a.x + b.x, a.y + b.y};}
PDD operator / (PDD a, db b) {return {a.x / b, a.y / b};}
PDD operator * (PDD a, db b) {return {a.x * b, a.y * b};}
db operator & (PDD a, PDD b) {return a.x * b.x + a.y * b.y;}
db operator * (PDD a, PDD b) {return a.x * b.y - a.y * b.x;}


db get_len(PDD a) {return sqrt(a & a); }
PDD rotate(PDD a, db b) {return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};}
PDD norm(PDD a) {return a / get_len(a);}//向量a的单位向量
db project(PDD a, PDD b) {return (a & b) / get_len(b);}//向量a在b上的投影
db area(PDD a, PDD b, PDD c) {return (b - a) * (c - a);}

void andrew()
{
    sort(q, q + n);
    for(int i = 0; i < n; i ++)
    {
        while(top >= 2 && sign(area(q[stk[top - 2]], q[stk[top - 1]], q[i])) <= 0) top --;
        stk[top ++] = i;
    }

    int k = top;
    for(int i = n - 1; i >= 0; i --)
    {
        while(top > k && sign(area(q[stk[top - 2]], q[stk[top - 1]], q[i])) <= 0) top --;
        stk[top ++] = i;
    }
    // cout << top << endl;
}

db rotating_calipers()
{
    top --;
    db res_s = INF;
    for(int i = 0,a = 2, b = 1, c = 2; i < top; i ++ )
    {
        //a:到直线ed最远的点；b：在直线ed投影最大的点；c:在直线ed投影最小的点
        PDD d = q[stk[i]], e = q[stk[i + 1]];
        while(dcmp(area(d, e, q[stk[a]]), area(d, e, q[stk[a + 1]])) < 0) a = (a + 1) % top;
        while(dcmp(project((q[stk[b]] - d), e - d), project((q[stk[b + 1]] - d), e - d)) < 0) b = (b + 1) % top;
        if(!i) c = a;
        while(dcmp(project((q[stk[c]] - d), e - d), project((q[stk[c + 1]] - d), e - d)) > 0) c = (c + 1) % top;
        PDD x = q[stk[a]], y = q[stk[b]], z = q[stk[c]];
        db h = area(d, e, x) / get_len(e - d);
        db w = project(y - z, e -d);
        if(w * h < res_s)
        {
            res_s = w * h;
            ans[0] = d + norm(e - d) * project(y - d, e - d);
            ans[3] = d + norm(e - d) * project(z - d, e - d);
            PDD u = norm(rotate(e - d, -PI / 2));
            ans[1] = ans[0] + u * h;
            ans[2] = ans[3] + u * h;
        }

    }
    return res_s;
}



int main()
{
    scanf("%d", &n);
    for(int i = 0; i < n; i ++) scanf("%lf %lf", &q[i].x, &q[i].y);

    andrew();

    printf("%.5lf\n", rotating_calipers());
    // exit(0);
    int k = 0;
    for(int i = 1; i < 4; i ++)
        if(dcmp(ans[i].y, ans[k].y) < 0 || (!dcmp(ans[i].y, ans[k].y) && dcmp(ans[i].x, ans[k].x))) k = i;

    for(int i = 0; i < 4; i ++, k ++)
    {
        k %= 4;
        db u = ans[k].x, v = ans[k].y;
        if(!sign(u)) u = 0;
        if(!sign(v)) v = 0;
        printf("%.5lf %.5lf\n", u, v);
    }
    return 0;
}

旋转卡壳步骤：

（最远两点的一定在构成凸包的两点）
1. andrew求凸包
2. 逆时针求每条线对应最远的点更新答案，随着向后枚举线，其对应的点也在向后移，故可以双指针。

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define y second
#define x first
using namespace std;
typedef pair<int, int> PII;
const int N = 5e4 + 10;

PII q[N];
int stk[N], top, n;

PII operator - (PII a, PII b) {return {a.x - b.x, a.y - b.y}; }
int cross(PII a, PII b) {return a.x * b.y - a.y * b.x; }
int area(PII a, PII b, PII c) {return cross(b - a, c - a); }
int get_dis(PII a, PII b) {int dx = a.x - b.x, dy = a.y - b.y; return dx * dx + dy * dy;}

void andrew()
{
    sort(q, q + n);
    for(int i = 0; i < n; i ++)
    {
        //先找下凸壳
        while(top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0) top --;
        stk[top ++] = i;//注意这里是top++，为了totating_calipers里面%方便
    }

    int k = top;
    for(int i = n - 1; i >= 0; i --)
    {
        while(top > k && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0) top --;
        stk[top ++] = i;
    }
}

int rotating_calipers()
{
    top --;//第一个是q[0],最后一个也是q[0]
    if(top <= 2) return get_dis(q[stk[0]], q[stk[top - 1]]);//一条直线，此时top=2
    // cout << top << endl;
    int res = 0;
    for(int i = 0, j = 2; i < top; i ++)
    {
        PII u = q[stk[i]], v = q[stk[i + 1]];
        while(area(u, v, q[stk[j]]) < area(u, v, q[stk[(j + 1)]])) j = (j + 1) % top;//点j+1离直线uv距离更远
        res = max(res, max(get_dis(u, q[stk[j]]), get_dis(v, q[stk[j]])));
    }
    return res;
}

int main()
{
    scanf("%d", &n);
    for(int i = 0; i < n; i ++) scanf("%d%d", &q[i].x, &q[i].y);
    andrew();
    printf("%d\n", rotating_calipers());

    return 0;
}

求三维凸包步骤:

1. 将所有点进行微小抖动
2. 前三的点构成平面的正反面加入plane
3. 从第四个点枚举新点（增量法$O(n)$）
4. 枚举plane，能被新点“照到”的点删除，不能“照到”的点保留
5. 枚举plane中的所有边，当边ab能被“照到”，ba不能被“照到”，则该边为一条分界边，该边与新点构成的三角平面应加入plane
6. 操作进行完，plane中存的即为前i个点构成的三维凸包；i枚举到n，即为前n个点构成三维凸包。

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
const int N = 210;
const double eps = 1e-10;
typedef double db;

int n, m;
bool g[N][N];

db get_eps()
{
    return ((db)rand() / RAND_MAX - 0.5) * eps;//RAND_MAX是rand()随机的最大数；(db)rand() / RAND_MAX - 0.5 --[0 ~ 0.5]
}

struct Point
{
    db x, y, z;
    //微小抖动，防止四点共面（四点共面分割的三角形不唯一）
    void shake(){x += get_eps(), y += get_eps(), z += get_eps();}
    Point operator+ (Point t) {return {x + t.x, y + t.y, z + t.z};}
    Point operator- (Point t) {return {x - t.x, y - t.y, z - t.z};}
    db operator& (Point t) {return {x * t.x + y * t.y + z * t.z};}//点积--|X|*|Y|*cos(c)
    //叉积-- 返回向量x与向量y构成平面的法向量
    Point operator* (Point t) {return{y * t.z - z * t.y, z * t.x - x * t.z, x * t.y - y * t.x}; } 
    db len() {return sqrt(x * x + y * y + z * z);}//求模长
}q[N];

struct Plane
{
    int v[3];//逆时针存储每个三角平面上的三个点
    Point norm() {return (q[v[1]] - q[v[0]]) * (q[v[2]] - q[v[0]]); }//求平面法向量
    db area() {return norm().len() / 2; }//求面积
    bool above(Point a) {return ((a - q[v[0]]) & norm()) >= 0; } //平面一点到a的向量与法向量的点积 >= 0 --> 该点在平面里或平面上
}plane[N], bp[N];//平面和备份平面

void get_convex_3d()
{
    //前三个点构成的平面
    plane[m ++] = {0, 1, 2};//正面
    plane[m ++] = {2, 1, 0};//反面

    for(int i = 3; i < n; i ++)//从第四个点开始
    {
        int cnt = 0;
        for(int j = 0; j < m; j ++) //将不能i照到的面保留--（点i在面的下面）
        {
            bool t = plane[j].above(q[i]);
            if(!t) bp[cnt ++] = plane[j];
            for(int k = 0; k < 3; k ++) g[plane[j].v[k]][plane[j].v[(k + 1) % 3]] = t;//若面能被照到，则该面的三边也能照到；反之
        }
        for(int j = 0; j < m; j ++)//遍历每个面的三条边，边ab能被照到，ba不能被照到，则该边为分界边，该边与新点构成的面要保留
            for(int k = 0; k < 3; k ++)
            {
                int a = plane[j].v[k], b = plane[j].v[(k + 1) % 3];
                if(g[a][b] && !g[b][a]) bp[cnt ++] = {a, b, i};
            }
        m = cnt;
        for(int j = 0; j < m; j ++) plane[j] = bp[j];
    }
}

int main()
{
    scanf("%d", &n);
    for(int i = 0; i < n; i ++) scanf("%lf %lf %lf", &q[i].x, &q[i].y, &q[i].z), q[i].shake();
    get_convex_3d();

    db res = 0;
    for(int i = 0; i < m; i ++) res += plane[i].area();
    printf("%.6lf\n", res);

    return 0;
}

椭圆可以看做由圆通过将某一个轴伸缩变换得到的。

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first 
#define y second
using namespace std;
typedef pair<double, double> PDD;
typedef double db;
const int N = 5e4 + 10;
const db PI = acos(-1), eps = 1e-9;

PDD p[N];
int n;

struct Circle
{
    PDD o;
    db r;
};

PDD operator - (PDD a, PDD b) {return {a.x - b.x, a.y - b.y};}
PDD operator + (PDD a, PDD b) {return {a.x + b.x, a.y + b.y};}
PDD operator * (PDD a, db b) {return {a.x * b, a.y * b};}
PDD operator / (PDD a, db b) {return {a.x / b, a.y / b};}

int cmp(db a, db b) {if(fabs(a - b) < eps) return 0; if(a < b) return -1; return 1;}
int sign(db a) {if(fabs(a) < eps) return 0; if(a < 0) return -1; return 1;}
db cross(PDD a, PDD b){return a.x * b.y - a.y * b.x;}
PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w) {PDD u = p - q; db t = cross(w, u) / cross(v, w); return {p.x + t * v.x, p.y + t * v.y};}
PDD rotate(PDD a, db b) {return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};}
pair<PDD, PDD> get_mid_c_line(PDD a, PDD b) {return {(a + b) / 2, rotate(b - a, PI / 2)};}
db get_dis(PDD a, PDD b) {db dx = a.x - b.x, dy = a.y - b.y; return sqrt(dx * dx + dy * dy);}
Circle get_circle(PDD a, PDD b, PDD c) {auto u = get_mid_c_line(a, b), v = get_mid_c_line(a, c); PDD t = get_line_intersection(u.x, u.y, v.x, v.y); return {t, get_dis(t, a)};}


int main()
{
    scanf("%d", &n);
    for(int i = 1; i <= n; i ++) scanf("%lf %lf", &p[i].x, &p[i].y);

    db h, t;
    scanf("%lf %lf", &h, &t);

    random_shuffle(p + 1, p + 1 + n);
    for(int i = 1; i <= n; i ++) p[i] = rotate(p[i], h / 180 * PI);//这里h单位是度，要将其先变为弧度（rad）
    for(int i = 1; i <= n; i ++) p[i].x /= t;

    Circle c = {p[1], 0};
    for(int i =  1; i <= n; i ++)
        if(cmp(c.r, get_dis(c.o, p[i])) < 0)
        {
            c = {p[i], 0};
            for(int j = 1; j < i; j ++)
                if(cmp(c.r, get_dis(c.o, p[j])) < 0)
                {
                    c = {(p[i] + p[j]) / 2, get_dis(p[i] , p[j]) / 2};
                    for(int k = 1; k < j; k ++)
                        if(cmp(c.r, get_dis(c.o, p[k])) < 0)
                            c = get_circle(p[i], p[j], p[k]);
                }
        }

    printf("%.3lf\n", c.r);

    return 0;
}

压行，hh~

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;
typedef pair<double, double> PDD;
typedef double db;
const int N = 1e5 + 10;
const db PI = acos(-1), eps = 1e-8;

PDD p[N];
int n;
//定义圆。
struct Circle
{
    PDD c;//圆心
    double r;//半径
};

PDD operator - (PDD a, PDD b) {return {a.x - b.x, a.y - b.y};}
PDD operator + (PDD a, PDD b) {return {a.x + b.x, a.y + b.y};}
PDD operator * (PDD a, db b) {return {a.x * b, a.y * b};}
PDD operator / (PDD a, db b) {return {a.x / b, a.y / b};}
int cmp(db x, db y){ if(fabs(x - y) < eps) return 0; if(x < y) return -1; return 0;}
int sign(db x){ if(fabs(x) < eps) return 0; if(x < 0) return -1; return 0;}
//叉积
db cross(PDD a, PDD b){return a.x * b.y - a.y * b.x;}
//求两直线交点：参数--（线1起点，线1方向，线2起点，线2方向）
PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w){PDD u = p - q; db t = cross(w, u) / cross(v, w); return {p.x + t * v.x, p.y + t * v.y};}
//将直线顺时针转b：参数--（方向向量，角度）
PDD rotate(PDD a, db b) {return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};}
//求直线的中垂线：参数--（起点，终点）
pair<PDD, PDD> get_mid_c_line(PDD a, PDD b){return {(a + b) / 2, rotate(b - a, PI / 2)};}
//两点距离：参数--（点1，点2）
db get_dis(PDD a, PDD b) {db dx = a.x - b.x, dy = a.y - b.y; return sqrt(dx * dx + dy * dy);}
//求三点确定的圆：参数--（点1，点2，点3）
Circle get_circle(PDD a, PDD b, PDD c){auto u = get_mid_c_line(a, b), v = get_mid_c_line(a, c); PDD t = get_line_intersection(u.x, u.y, v.x, v.y); return {t, get_dis(t, a)};}

int main()
{
    scanf("%d", &n);
    for(int i = 1; i <= n; i ++) scanf("%lf %lf", &p[i].x, &p[i].y);

    //随机增量，将点随机打乱顺序 -- 考虑一个点作为覆盖圆上的概率大致是 3/n
    random_shuffle(p + 1, p + 1 + n);
    Circle o = {p[1], 0};//初始将第一个点本身作为最小圆
    for(int i = 1; i <= n; i ++)
        if(cmp(o.r, get_dis(p[i], o.c)) < 0)
        {
            o = {p[i], 0};
            for(int j = 1; j < i; j ++)
                if(cmp(o.r, get_dis(p[j], o.c)) < 0)
                {
                    o = {(p[i] + p[j]) / 2, get_dis(p[i], p[j]) / 2};
                    for(int k = 1; k < j; k ++)
                        if(cmp(o.r, get_dis(p[k], o.c)) < 0)
                            o = get_circle(p[i], p[j], p[k]);
                }

        }
    printf("%.10lf\n", o.r);
    printf("%.10lf %.10lf\n", o.c.x, o.c.y);

    return 0;
}