计算几何算法模板

1. 二维几何

1.1常用函数模板

const double eps = 1e-8;
const double pi = acos(-1);
typedef pair<double, double> PDD;
#define x first
#define y second
struct Line
{
    PDD st, ed; // 存起点和终点
}lines[maxn];
int sign(double x)  // 判断正负
{
    if(fabs(x) < eps) return 0;
    if(x < 0) return -1;
    return 1;
}
int dcmp(double x, double 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, double t)  // 向量乘系数
{
    return {a.x * t, a.y * t};
}
double operator*(PDD a, PDD b)  // 点积
{
    return {a.x * b.x + a.y * b.y};
}
PDD operator/ (PDD a, double b)  // 实现除法
{
    return {a.x / b, a.y / b};
}
double dot(double x1, double y1, double x2, double y2)  // 求点积
{
    return x1 * x2 + y1 * y2;
}
double cross(double x1, double y1, double x2, double y2)  // 求叉积
{
    return x1 * y2 - x2 * y1;
}
double cross(PDD a, PDD b)
{
    return a.x * b.y - a.y * b.x;
}
double area(PDD a, PDD b, PDD c)  // 求面积
{
    return cross(b - a, c - a);
}
double get_length(PDD a)  // 取模
{
    return sqrt(dot(a.x, a.y, a.x, a.y));
}
double get_angle(PDD a, PDD b)  // 计算向量夹角
{
    return acos(a * b / get_length(a) / get_length(b));
}
double project(PDD a, PDD b, PDD c) // 求ac在ab上的投影长度
{
    return (b - a) * (c - a) / get_length(b - a);
}
PDD norm(PDD a)  // 求单位向量
{
    return a / get_length(a);
}
PDD rotate(PDD a, double b)  // 绕原点逆时针旋转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)  // 求两直线交点
{
    auto u = p - q;
    double t = cross(w, u) / cross(v, w);
    return {p.x + v.x * t, p.y + v.y * t};
}
PDD get_line_intersection(Line a, Line b)  // 求两直线交点
{
    return get_line_intersection(a.st, a.ed - a.st, b.st, b.ed - b.st);
}
double distance_to_line(PDD p, PDD a, PDD b)  // 点p到直线ab的距离
{
    auto v1 = b - a, v2 = p - a;
    return fabs(cross(v1, v2) / get_length(v1));
}
//判断b 和 c 的交点是否在直线a右侧
bool on_right(Line a, Line b, Line c)
{
    auto o = get_line_intersection(b, c);
    return sign(area(a.st, a.ed, o)) <= 0;
}
bool on_segment(PDD p, PDD a, PDD b)  // 判断点p是否在线段ab上
{
    // 首先，三点共线，然后点积小于等于0
    return !sign((p - a) * (p - b)) && sign((p - a) & (p - b)) <= 0;
}
double distance_to_segment(PDD p, PDD a, PDD b) // 点p到线段ab的距离
{
    if(a == b) return get_length(p - a);
    auto v1 = b - a, v2 = p - a, v3 = p - b;
    if(sign(v1 * v2) < 0) return get_length(v2);
    if(sign(v1 * v3) < 0) return get_length(v3);
    return distance_to_line(p, a, b);
}
// 求直线和圆的交点，存在pa和pb中，利用向量求法
double get_line_circle_intersection(PDD a, PDD b, PDD &pa, PDD &pb)
{
    auto e = get_line_intersection(a, b - a, r, rotate(b - a, pi / 2));
    auto min_d = get_dist(r, e);
    if(!on_segment(e, a, b)) min_d = min(get_dist(r, a), get_dist(r, b));
    if(dcmp(R, min_d) <= 0) return min_d;
    auto len = sqrt(R * R - get_dist(r, e) * get_dist(r, e));
    pa = e + norm(a - b) * len;
    pb = e + norm(b - a) * len;
    return min_d;
}
// 求向量a和向量b围成的扇形的面积
double get_sector(PDD a, PDD b)
{
    double angle = acos(a & b / get_len(a) / get_len(b));
    if(sign(a * b) < 0) angle = -angle;
    return R * R * angle / 2;
}

1.2 距离转换

1.3 Pick定理

Pick定理：给定顶点均为整点的简单多边形，皮克定理说明了其面积 $A$ 和内部格点数目$i$,边上格点数目b的关系：
$A=i+b/2-1$。

1.4 多边形

1.4.1三角形

(1) 求面积

1. 叉积
2. 海伦公式

p = (a + b + c) / 2;
S = sqrt(p * (p - a) * (p - b) * (p - c));

(2) 三角形四心

1. 外心，外接圆圆心
   三条边中垂线交点，到三角形三个顶点的距离相等。
2. 内心，内切圆圆心
   角平分线交点，到三边距离相等。
3. 垂心
   三条垂线交点。
4. 重心
   三条中线交点(到三角形三顶点距离的平方和最小的点，三角形内到三边距离之积最大的点)
   (3) 求多边形面积
   可以从第一个顶点出发把凸多边形分成n - 2个三角形，然后把面积加起来。

double polygon_area(PDD p[], int n)
{
    double s = 0;
    for(int i = 1; i < n - 1; i ++){
        s += cross(p[i] - p[0], p[i + 1] - p[i]);
    }
    return s / 2;
}

1.5 极角序

在平面选取一个顶点O，称为极点；做射线X，那么OX就是极轴。平面上任意一点和O点作向量。通常选定逆时针方向为正，一般我们以X轴为极轴，那么极角就是平面向量与X轴的夹角。

极角排序：
给定极轴，在平面上分布着若干点，将这些点相对于极轴的大小排序，就叫做极角排序。注意，极角相同时，按照X的升序处理。

// 方法1：利用atan2函数，速度快，精度低
#include<cmath>
bool cmp(PDD a, PDD b)
{
    if(dcmp(atan2(a.y, a.x), atan2(b.y, b.x)) == 0){
        return a.x < b.x;
    }
    return atan2(a.y, a.x) < atan2(b.y, b.x);
}

// 方法2：利用向量叉积排序
bool cmp(PDD a, PDD b)
{
    if(sign(cross(a, b)) == 0) return a.x < b.x;
    return sign(cross(a, b)) > 0;
}

// 自定义极点
PII o; 
bool cmp(PDD a, PDD b)
{
    auto p = a - o, q = b - o;
    if(sign(cross(p, q)) == 0) return p.x < q.x;
    return sign(cross(p, q)) > 0;
}

1.6 二维凸包

PDD q[N];
int stk[N];
bool used[N];
// 求凸包周长
double 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)
        {
            if (area(q[stk[top - 2]], q[stk[top - 1]], q[i]) < 0)
                used[stk[ -- top]] = false;
            else top -- ;
        }
        stk[top ++ ] = i;
        used[i] = true;
    }

    used[0] = false;
    for (int i = n - 1; i >= 0; i -- )
    {
        if (used[i]) continue;
        while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0)
            top -- ;
        stk[top ++ ] = i;
    }
    top --;
    double res = 0;
    for (int i = 1; i <= top; i ++ )
        res += get_dist(q[stk[i - 1]], q[stk[i]]);
    return res;
}

1.7 半平面交

Line lines[maxn];
int n, m, cnt;
PDD pg[maxn], ans[maxn];
int q[maxn];

bool cmp(const Line& a, const Line& b)
{
    double thetaA = get_angle(a), thetaB = get_angle(b);
    if(!dcmp(thetaA, thetaB)) return area(a.st, a.ed, b.ed) < 0;
    return thetaA < thetaB;
}
// 求交集面积
double half_plane_intersection()
{
    sort(lines, lines + cnt, cmp);
    int hh = 0, tt = -1;
    for(int i = 0; i < cnt; i ++){
        if(i && !dcmp(get_angle(lines[i]), get_angle(lines[i - 1]))) continue;
        while(hh + 1 <= tt && on_right(lines[i], lines[q[tt - 1]], lines[q[tt]])) tt --;
        while(hh + 1 <= tt && on_right(lines[i], lines[q[hh]], lines[q[hh + 1]])) hh ++;
        q[++ tt] = i;
    }
    while(hh + 1 <= tt && on_right(lines[q[hh]], lines[q[tt - 1]], lines[q[tt]])) tt --;
    while(hh + 1 <= tt && on_right(lines[q[tt]], lines[q[hh]], lines[q[hh + 1]])) hh ++;
    q[++ tt] = q[hh];
    int k = 0;
    for(int i = hh; i < tt; i ++){
        ans[k ++] = get_line_intersection(lines[q[i]], lines[q[i + 1]]);
    }
    double res = 0;
    for(int i = 1; i + 1 < k; i ++){
        res += area(ans[0], ans[i], ans[i + 1]);
    }
    return res / 2;
}

1.8 最小圆覆盖

给定二维平面上的$n$个点，找一个最小的能包含所有点的圆。

struct Line
{
    PDD st, ed;
};
struct Circle
{
    PDD p;
    double r;
};
Line get_line(Line a)  // 求直线a的中垂线
{
    return {(a.st + a.ed) / 2, rotate(a.ed - a.st, pi / 2)};
}
Circle get_clrcle(PDD a, PDD b, PDD c)  // 根据三点确定一个圆
{
    auto u = get_line({a, b}), v = get_line({a, c});
    auto p = get_line_intersection(u.st, u.ed, v.st, v.ed);
    return {p, get_dist(p, a)};
}
Circle get_the_min_circle()  // 求最小圆
{
    random_shuffle(q, q + n);
    Circle c({q[0], 0});
    for(int i = 1; i < n; i ++){
        if(dcmp(c.r, get_dist(c.p, q[i])) < 0){
            c = {q[i], 0};
            for(int j = 0; j < i; j ++){
                if(dcmp(c.r, get_dist(c.p, q[j])) < 0){
                    c = {(q[i] + q[j]) / 2, get_dist(q[i], q[j]) / 2};
                    for(int k = 0; k < j; k ++){
                        if(dcmp(c.r, get_dist(c.p, q[k])) < 0){
                            c = get_clrcle(q[i], q[j], q[k]);
                        }
                    }
                }
            }
        }
    }
    return c;
}

1.9 最小矩形覆盖

给定二维平面上的一组点，求覆盖这组点的最小的矩形的面积。

int n;
PDD q[maxn];
PDD ans[maxn];
double min_area = inf;
int stk[maxn], top;
bool used[maxn];

void get_convex()  // 求凸包
{
    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)
        {
            if (area(q[stk[top - 2]], q[stk[top - 1]], q[i]) < 0)
                used[stk[ -- top]] = false;
            else top -- ;
        }
        stk[top ++ ] = i;
        used[i] = true;
    }

    used[0] = false;
    for (int i = n - 1; i >= 0; i -- )
    {
        if (used[i]) continue;
        while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0)
            top -- ;
        stk[top ++ ] = i;
    }
    top --;
}
void rotating_calipers()  // 旋转卡壳
{
    for(int i = 0, a = 2, b = 1, c = 2; i < top; i ++){
        auto 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(d, e, q[stk[b]]), project(d, e, q[stk[b + 1]])) < 0) b = (b + 1) % top;
        if (!i) c = a;
        while (dcmp(project(d, e, q[stk[c]]), project(d, e, q[stk[c + 1]])) > 0) c = (c + 1) % top;
        auto x = q[stk[a]], y = q[stk[b]], z = q[stk[c]];
        auto h = area(d, e, x) / get_len(e - d);
        auto w = ((y - z) & (e - d)) / get_len(e - d);
        if (h * w < min_area)
        {
            min_area = h * w;
            ans[0] = d + norm(e - d) * project(d, e, y);
            ans[3] = d + norm(e - d) * project(d, e, z);
            auto u = norm(rotate(e - d, -pi / 2));
            ans[1] = ans[0] + u * h;
            ans[2] = ans[3] + u * h;
        }
    }
}
int main()
{
    cin >> n;
    for(int i = 0; i < n; i ++) cin >> q[i].x >> q[i].y;
    get_convex();
    rotating_calipers();
    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) < 0)
            k = i;
    printf("%.5lf\n", min_area);
    for (int i = 0; i < 4; i ++, k ++ )
    {
        auto x = ans[k % 4].x, y = ans[k % 4].y;
        if (!sign(x)) x = 0;
        if (!sign(y)) y = 0;
        printf("%.5lf %.5lf\n", x, y);
    }
    return 0;
}

1.10 旋转卡壳

利用旋转卡壳求解二维平面最远点对

int stk[maxn], top;
bool used[maxn];

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)
        {
            if (area(q[stk[top - 2]], q[stk[top - 1]], q[i]) < 0)
                used[stk[ -- top]] = false;
            else top -- ;
        }
        stk[top ++ ] = i;
        used[i] = true;
    }

    used[0] = false;
    for (int i = n - 1; i >= 0; i -- )
    {
        if (used[i]) continue;
        while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0)
            top -- ;
        stk[top ++ ] = i;
    }
    top -- ;
}
int rotating_calipers()
{
    if(top <= 2) return get_dist(q[0], q[n - 1]);
    int res = 0;
    for(int i = 0, j = 2; i < top; i ++){
        auto d = q[stk[i]], e = q[stk[i + 1]];
        while(area(d, e, q[stk[j]]) < area(d, e, q[stk[j + 1]])) j = (j + 1) % top;
        res = max(res, max(get_dist(d, q[stk[j]]), get_dist(e, q[stk[j]])));
    }
    return res;
}

1.11 三角剖分

// 求任意多边形和圆形面积交集的面积
// 利用三角剖分，将圆心和每条边的两个点连线，每次求三角形和原型交集的面积，累加求和即可。
// 要进行五种情况的分类讨论

// 求三角形oab和圆面积交集的面积，圆心处于原点
PDD r;
double R;
double get_circle_triangle_area(PDD a, PDD b)
{
    auto da = get_dist(r, a), db = get_dist(r, b);
    if(dcmp(R, da) >= 0 && dcmp(R, db) >= 0) return a * b / 2;
    if(!sign(a * b)) return 0;
    PDD pa, pb;
    auto min_d = get_line_circle_intersection(a, b, pa, pb);
    if(dcmp(R, min_d) <= 0) return get_sector(a, b);
    if(dcmp(R, da) >= 0) return get_sector(pb, b) + a * pb / 2;
    if(dcmp(R, db) >= 0) return get_sector(a, pa) + pa * b / 2;
    return get_sector(a, pa) + pa * pb / 2 + get_sector(pb, b);
}
// 求面积，时间复杂度O(n)
double work()
{
    double res = 0;
    for(int i = 0; i < n; i ++){
        res += get_circle_triangle_area(q[i], q[(i + 1) % n]);
    }
    return fabs(res);
}

1.12 扫描线

在二维平面中，求多个规则矩形面积的并。

// 矩形个数不多时，利用区间合并和扫描线。
// 按照x坐标大小进行排序，然后在每两个相邻的x之间，枚举全部矩形，矩形一定是跨过或者刚好在区间内，对这些矩形的y值进行区间合并，求出y方向的总长度，然后求面积再求和即可。
#include<iostream>
#include<cstring>
#include<vector>
#include<algorithm>
using namespace std;
const int maxn = 1010;
#define int long long
#define x first
#define y second
typedef pair<int,int> PII;
PII l[maxn], r[maxn];
int n;
PII q[maxn];
int range_area(int a, int b)
{
    int cnt = 0;
    for(int i = 0; i < n; i ++){
        if(l[i].x <= a && r[i].x >= b){
            q[cnt ++] = {l[i].y, r[i].y};
        }
    }
    if(!cnt) return 0;
    sort(q, q + cnt);
    int res = 0;
    int 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;
        }
    }
    res += ed - st;
    return res * (b - a);
}
signed main()
{
    cin >> n;
    vector<int> v;
    for(int i = 0; i < n; i ++){
        cin >> l[i].x >> l[i].y >> r[i].x >> r[i].y;
        v.push_back(l[i].x);
        v.push_back(r[i].x);
    }
    sort(v.begin(), v.end());
    v.erase(unique(v.begin(), v.end()), v.end());
    int res = 0;
    for(int i = 0; i < v.size() - 1; i ++){
        res += range_area(v[i], v[i + 1]);
    }
    cout << res << endl;
    return 0;
}

// 点数非常大时，需要用线段树进行维护
#include<iostream>
#include<cstring>
#include<algorithm>
#include<vector>
using namespace std;
const int maxn = 10010;
typedef struct node
{
	int l, r;
	int cnt;
	double len;
}Node;
typedef struct node1
{
	double x, y1, y2;
	int k;
	double operator < (const struct node1 &w) const{
		return x < w.x;
	}
}Segment;
Segment a[maxn * 2];
Node tr[maxn * 8];
int n;
vector<double> ys;
int find(double y)
{
	return lower_bound(ys.begin(), ys.end(), y) - ys.begin();
}
void pushup(Node &u, Node &l, Node &r)
{
	if(u.cnt) u.len = ys[u.r + 1] - ys[u.l];
	else if(u.l != u.r)
	{
		u.len = l.len + r.len;
	}else u.len = 0;
}
void pushup(int u)
{
	pushup(tr[u], tr[u << 1], tr[u << 1 | 1]);
}
void build(int u, int l, int r)
{
	tr[u] = {l, r, 0, 0};
	if(l != r)
	{
	    int mid = l + r >> 1;
	    build(u << 1, l, mid), build(u << 1 | 1, mid + 1, r);
	    pushup(u);
	}
}
void modify(int u, int l, int r, int d)
{
	if(tr[u].l >= l && tr[u].r <= r)
	{
		tr[u].cnt += d;
		pushup(u);
	}
	else
	{
		int mid = tr[u].l + tr[u].r >> 1;
		if(l <= mid) modify(u << 1, l, r, d);
		if(r > mid) modify(u << 1 | 1, l, r, d);
		pushup(u);
	}
}
int main()
{
	int t = 1;
	while(cin >> n, n)
	{
		ys.clear();
		for(int i = 0, j = 0; i < n; i ++){
			double x1, y1, x2, y2; cin >> x1 >> y1 >> x2 >> y2;
			a[j ++] = {x1, y1, y2, 1};
			a[j ++] = {x2, y1, y2, -1};
			ys.push_back(y1); ys.push_back(y2);
		}
		sort(ys.begin(), ys.end());
		ys.erase(unique(ys.begin(), ys.end()), ys.end());
		build(1, 0, ys.size() - 2);
		sort(a, a + n * 2);
		double res = 0;
		for(int i = 0; i < n * 2; i ++)
		{
			if(i > 0) res += tr[1].len * (a[i].x - a[i - 1].x);
			modify(1, find(a[i].y1), find(a[i].y2) - 1, a[i].k);
		}
		printf("Test case #%d\n", t ++);
		printf("Total explored area: %.2lf\n\n", res);
	}
	return 0;
}

求三角形并的面积

给出n个三角形，利用扫描线求三角形并集的面积。
按照三角形的所有顶点和交点进行划分区间，每个区间内部是一些梯形，最后转化为求这些梯形面积的和，可以考虑在左右区间上进行区间合并，加快计算。
时间复杂度：$O(n^{3}logn)$

#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
using namespace std;
const int maxn = 110;
#define x first
#define y second
typedef pair<double, double> PDD;
const double eps = 1e-8;
const double pi = acos(-1);
const double inf = 1e12;
PDD tr[maxn][3];
PDD q[maxn];
int n;
int sign(double x)
{
    if(fabs(x) < eps) return 0;
    if(x < 0) return -1;
    return 1;
}
int dcmp(double x, double 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, double t)
{
    return {a.x * t, a.y * t};
}
PDD operator / (PDD a, double t)
{
    return {a.x / t, a.y / t};
}
double operator * (PDD a, PDD b)
{
    return a.x * b.y - a.y * b.x;
}
double operator & (PDD a, PDD b)
{
    return a.x * b.x + a.y * b.y;
}
bool on_segment(PDD p, PDD a, PDD b)
{
    return sign((p - a) & (p - b)) <= 0;
}
PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w)
{
    if(!sign(v * w)) return {inf, inf};
    auto u = p - q;
    double t = w * u / (v * w);
    auto o = p + v * t;
    if(!on_segment(o, p, p + v) || !on_segment(o, q, q + w)) return {inf, inf};
    return o;
}
double line_range(double a, int side)  // 求区间上线段并集长度，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;
        if(!dcmp(t[0].x, a) && !dcmp(t[1].x, a)){  // 特判一下左边和直线重合的情况
            if(side){
                q[cnt ++] = {t[0].y ,t[1].y};
            }
        }
        else if(!dcmp(t[2].x, a) && !dcmp(t[1].x, a)){  // 特判一下右边和直线重合的情况
            if(!side){
                q[cnt ++] = {t[2].y, t[1].y};
            }
        }else{
            double d[3];
            int u = 0;
            for(int j = 0; j < 3; j ++){
                auto o = get_line_intersection(t[j], t[(j + 1) % 3] - t[j], {a, -inf}, {0, inf * 2});
                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);  // 求一遍区间合并
    double 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;
            st = q[i].x, ed = q[i].y;
        }
    }
    res += ed - st;
    return res;
}
double range_area(double a, double b)
{
    return (line_range(a, 1) + line_range(b, 0)) * (b - a) / 2;
}
int main()
{
    cin >> n;
    vector<double> v;
    for(int i = 0; i < n; i ++){
        for(int j = 0; j < 3; j ++){
            cin >> tr[i][j].x >> tr[i][j].y;
            v.push_back(tr[i][j].x);
        }
        sort(tr[i], tr[i] + 3);  // 排序后，方便求交点
    }
    for(int i = 0; i < n; i ++){  // 求每两条边的交点
        for(int j = 0; j < n; j ++){
            for(int x = 0; x < 3; x ++){
                for(int y = 0; y < 3; y ++){
                    auto 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)) v.push_back(o.x);  // 如果存在交点
                }
            }
        }
    }
    sort(v.begin(), v.end());
    v.erase(unique(v.begin(), v.end()), v.end());
    double res = 0;
    for(int i = 0; i < v.size() - 1; i ++){
        res += range_area(v[i], v[i + 1]);
    }
    printf("%.2lf\n", res);
    return 0;
}

1.13 自适应辛普森积分

给定一个积分函数和上下限，求积分值。

const double eps = 1e-12;
double l, r;
double f(double x)
{
    // 根据题意来写
    return sin(x) / x;
}
double simpson(double l, double r)
{
    double mid = (l + r) / 2;
    return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}
double asr(double l, double r, double s)
{
    double mid = (l + r) / 2;
    double 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()
{
    cin >> l >> r;
    printf("%.6lf\n", asr(l, r, simpson(l, r)));
}

求圆的的并的面积

#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int maxn = 1010;
const double eps = 1e-8;
const double pi = acos(-1);
#define x first
#define y second
typedef pair<double, double> PDD;
typedef struct node
{
    PDD r;
    double R;
}Circle;
Circle a[maxn];
PDD q[maxn];
int n;
int dcmp(double x, double y)
{
    if(fabs(x - y) < eps) return 0;
    if(x < y) return -1;
    return 1;
}
double f(double x)  // 所有圆和X = x这条直线交集的长度
{
    int cnt = 0;
    for(int i = 0; i < n; i ++){
        auto X = fabs(a[i].r.x - x), R = a[i].R;
        if(dcmp(X, R) < 0){
            auto Y = sqrt(R * R - X * X);
            q[cnt ++] = {a[i].r.y - Y, a[i].r.y + Y};
        }
    }
    if(!cnt) return 0;
    sort(q, q + cnt); // 区间合并
    double 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;
            st = q[i].x, ed = q[i].y;
        }
    }
    res += ed - st;
    return res;
}
double simpson(double l, double r)
{
    auto mid = (l + r) / 2;
    return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}
double asr(double l, double r, double s)
{
    double mid = (l + r) / 2;
    auto 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);
}
double get_area()
{
    cin >> n;
    double l = 2000, r = -2000; // 积分范围，根据题意来定
    for(int i = 0; i < n; i ++){
        cin >> a[i].r.x >> a[i].r.y >> a[i].R;
        l = min(l, a[i].r.x - a[i].R), r = max(r, a[i].r.x + a[i].R);
    }
    printf("%.3lf\n", asr(l - 100, r + 100, simpson(l, r)));
}

2. 三维计算几何

常用函数模板:

double rand_eps()
{
    return ((double)rand() / RAND_MAX - 0.5) * eps;
}

struct Point  // 三维坐标系中的点
{
    double x, y, z;
    void shake()  // 给坐标值加一个微小抖动
    {
        x += rand_eps(), y += rand_eps(), z += rand_eps();
    }
    Point operator- (Point t)
    {
        return {x - t.x, y - t.y, z - t.z};
    }
    double operator& (Point t) 
    {
        return x * t.x + y * t.y + z * t.z;
    }
    Point operator* (Point t)
    {
        return {y * t.z - t.y * z, z * t.x - x * t.z, x * t.y - y * t.x};
    }
    double 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]]);
    }
    double area()  // 求面积
    {
        return norm().len() / 2;
    }
    bool above(Point a) // 判断点a是否在平面上方
    {
        return ((a - q[v[0]]) & norm()) >= 0;
    }
}plane[N], np[N];

2.1 三维凸包

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 ++){
            bool t = plane[j].above(q[i]);
            if(!t) np[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 ++){
            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]){
                    np[cnt ++] = {a, b, i};
                }
            }
        }
        m = cnt;
        for(int j = 0; j < m; j ++) plane[j] = np[j];
    }
}

// 求面积
double res = 0;
for(int i = 0; i < m; i ++){
    res += plane[i].area();
}