基本模板
/* WXCU_ACM Yuxiao */
const double PI = acos(-1.0), eps = 1e-8;
int sign(double x) {
return fabs(x) < eps ? 0 : (x < 0) ? -1 : 1;
}
int dcmp(double x, double y) {
return sign(x - y);
}
double f(double x) {
return sin(x) / x;
}
double simpson(double l, double r) { // 求 f(x) 在l ~ 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);
}
/* 平面几何：点 线 */
struct Point {
double x, y;
Point() {}
Point(double x, double y) : x(x), y(y) {}
Point operator+ (Point t) {
return Point(x + t.x, y + t.y);
}
Point operator- (Point t) {
return Point(x - t.x, y - t.y);
}
Point operator* (double k) {
return Point(x * k, y * k);
}
Point operator/ (double k) {
return Point(x / k, y / k);
}
bool operator== (Point t) {
return sign(x - t.x) == 0 && sign(y - t.y) == 0;
}
bool operator< (Point t) { // 用于凸包逆时针排序
return sign(x - t.x) < 0 || (sign(x - t.x) == 0 && sign(y - t.y) < 0);
}
};
typedef Point Vector;
double cross(Point a, Point b) {
return a.x * b.y - a.y * b.x;
}
double dot(Vector a, Vector b) { // 点积
return a.x * b.x + a.y * b.y;
}
double length(Vector a) { // 向量a的长度
return sqrt(dot(a, a));
}
double area(Point a, Point b, Point c) { // 叉积
return cross(b - a, c - a);
}
double angle(Vector a, Vector b) { // 两向量夹角
return acos(dot(a, b) / length(a) / length(b));
}
double distance(Point a, Point b) { // 两点距离
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
Point P_rotate(Point a, double k) { // 点a绕原点逆时针旋转k度
return {a.x * cos(k) + a.y * sin(k), -a.x * sin(k) + a.y * cos(k)};
}
Vector V_rotate(Vector a, double k) { // 向量a绕原点逆时针旋转k度
return {a.x * cos(k) - a.y * sin(k), a.x * sin(k) + a.y * cos(k)};
}
Vector normal(Vector a) { // 向量A的单位法向量
return {-a.y / length(a), a.x / length(a)};
}
bool parallel(Vector a, Vector b) { // 向量平行或重合
return sign(cross(a, b)) == 0;
}
struct Line {
Point x, y;
Line() {}
Line(Point x, Point y) : x(x), y(y) {}
// 根据一个点和倾角angle确定直线 0 <= angle <= PI
Line(Point p, double angle) {
x = p;
if (sign(angle - PI / 2) == 0) {
y = (x + Point(0, 1));
} else {
y = (x + Point(1, tan(angle)));
}
}
// ax + by + c = 0
Line(double a, double b, double c) {
if (sign(a) == 0) {
x = Point(0, -c / b);
y = Point(1, -c / b);
} else if (sign(b) == 0) {
x = Point(-c / a, 0);
y = Point(-c / a, 1);
} else {
x = Point(0, -c / b);
y = Point(1, (-c - a) / b);
}
}
};
typedef Line Segment;
double angle(Line a) { // 直线倾角
return atan2(a.y.y - a.x.y, a.y.x - a.x.x);
}
bool lcmp(Line a, Line b) { // 用于半平面交逆时针排序
double A = angle(a), B = angle(b);
if (!dcmp(A, B)) {
return area(a.x, a.y, b.y) < 0;
} else {
return A < B;
}
}
double PL_distance(Point p, Line a) { // 点到直线的距离
return fabs(cross(p - a.x, a.y - a.x)) / distance(a.x, a.y);
}
double PS_distance(Point p, Segment a) { // 点到线段的距离
if (sign(dot(p - a.x, a.y - a.x)) < 0 || sign(dot(p - a.y, a.x - a.y)) < 0) {
return min(distance(p, a.x), distance(p, a.y));
}
return PL_distance(p, a);
}
int PL_relation(Point p, Line a) { // 点和直线的关系：-1在左侧，0在直线上，1在右侧
int c = sign(cross(p - a.x, a.y - a.x));
return !c ? 0 : c < 0 ? -1 : 1;
}
int L_relation(Line a, Line b) { // 两直线关系：0平行，1重合，2相交
if (sign(cross(a.y - a.x, b.y - b.x)) == 0) {
return PL_relation(a.x, b) == 0 ? 1 : 0;
} else {
return 2;
}
}
bool PS_relation(Point p, Segment a) { // 点和线段的关系：0点p不在线段v上；1点p在线段v上
return sign(cross(p - a.x, a.y - a.x)) == 0 && sign(dot(p - a.x, p - a.y)) <= 0;
}
Point PL_project(Point p, Line a) { // 点在直线上的投影
double k = dot(a.y - a.x, p - a.x) / length(a.y - a.x) / length(a.y - a.x);
return a.x + (a.y - a.x) * k;
}
Point line_intersection(Point p, Point v, Point q, Point w) { // 两直线的交点
Point u = p - q;
double t = cross(w, u) / cross(v, w);
return p + v * t;
}
Point line_intersection(Line a, Line b) { // 两直线的交点
return line_intersection(a.x, a.y - a.x, b.x, b.y - b.x);
}
Line mid_line(Point a, Point b) { // 线段的中垂线
return {(a + b) / 2, V_rotate((b - a), PI / 2)};
}
/* 平面几何：圆*/
struct Circle {
Point c; double r;
Circle() {}
Circle(Point c, double r) : c(c), r(r) {}
Circle(double x, double y, double R) {
c = Point(x, y);
r = R;
}
};
Circle get_circle(Point a, Point b, Point c) { // 三点确定一条直线
Line u = mid_line(a, b), v = mid_line(a, c);
Point p = line_intersection(u.x, u.y, v.x, v.y);
return {p, distance(p, a)};
}
int PC_relation(Point p, Circle c) { // 点和圆的关系: 0点在圆内, 1点在圆上, 2点在圆外.
int d = sign(distance(p, c.c) - c.r);
return !d ? 0 : d < 0 ? -1 : 1;
}
int LC_relation(Line a, Circle c) { // 直线和圆的关系：0直线在圆内, 1直线和圆相切, 2直线在圆外
int d = sign(PL_distance(c.c, a) - c.r);
return !d ? 0 : d < 0 ? -1 : 1;
}
int SC_relation(Segment a, Circle c) { // 线段和圆的关系：0线段在圆内, 1线段和圆相切, 2线段在圆外
int d = sign(PS_distance(c.c, a) - c.r);
return !d ? 0 : d < 0 ? -1 : 1;
}
int LC_intersection(Line a, Circle c, Point &x, Point &y) { // 直线和圆的交点 x, y是交点 返回值是交点个数
if (LC_relation(a, c) == 2) {
return 0; // 无交点
}
Point q = PL_project(c.c, a);
double d = PL_distance(c.c, a);
double k = sqrt(c.r * c.r - d * d);
if (sign(k) == 0) { // 1个交点，直线和圆相切
x = q, y = q;
return 1;
}
Point n = (a.y - a.x) / length(a.y - a.x);
x = q + n * k, y = q - n * k;
return 2;
}
/* 平面几何：多边形*/
struct Polygon {
int n; // 多边形的定点数
Point p[N]; // 多边形的点
Line a[N]; // 多边形的边
};
int PP_relation(Point np, Polygon pg) { // 判断点和任意多边形的关系: 3点上; 2边上; 1内部; 0外部
auto p = pg.p;
auto n = pg.n;
for (int i = 0; i < n; i ++) { // 点在多边形的顶点上
if (p[i] == np) {
return 3;
}
}
for (int i = 0; i < n; i ++) { // 点在多边形的边上
Line a = Line(p[i], p[(i + 1) % n]);
if (PS_relation(np, a)) {
return 2;
}
}
int num = 0;
for (int i = 0; i < n; i ++) {
int j = (i + 1) % n;
int c = sign(cross(np - p[j], p[i] - p[j]));
int u = sign(p[i].y - np.y);
int v = sign(p[j].y - np.y);
if (c > 0 && u < 0 && v >= 0) {
num ++;
} else if (c < 0 && u >= 0 && v < 0) {
num --;
}
}
return num != 0; // 1内部，0外部
}
double Polygon_area(Polygon pg) { // 多边形面积
auto p = pg.p;
auto n = pg.n;
double res = 0;
for (int i = 0; i < n; i ++) {
res += cross(p[i], p[(i + 1) % n]);
}
return res / 2;
}
Point Polygon_center(Polygon pg) {
auto p = pg.p;
auto n = pg.n;
Point res = {0, 0};
if (Polygon_area(pg) == 0) {
return res;
}
for (int i = 0; i < n; i ++) {
res = res + (p[i] + p[(i + 1) % n]) * cross(p[i], p[(i + 1) % n]);
}
return res / Polygon_area(pg) / 6;
}
/* 立体（三维）几何：点 线*/
struct Point3 {
double x, y, z;
Point3() {}
Point3(double x, double y, double z) : x(x), y(x), z(z) {}
Point3 operator+ (Point3 t) {
return Point3(x + t.x, y + t.y, z + t.z);
}
Point3 operator- (Point3 t) {
return Point3(x - t.x, y - t.y, z - t.z);
}
Point3 operator* (double k) {
return Point3(x * k, y * k, z * k);
}
Point3 operator/ (double k) {
return Point3(x / k, y / k, z / k);
}
bool operator== (Point3 t) {
return sign(x - t.x) == 0 && sign(y - t.y) == 0 && sign(z - t.z) == 0;
}
};
typedef Point3 Vector3;
Vector3 cross(Point3 a, Point3 b) {
return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
}
double dot(Vector3 a, Vector3 b) { // 点积
return a.x * b.x + a.y * b.y + a.z * b.z;
}
double length(Vector3 a) { // 向量a的长度
return sqrt(dot(a, a));
}
Vector3 area(Point3 a, Point3 b, Point3 c) { // 叉积
return cross(b - a, c - a);
}
double angle(Vector3 a, Vector3 b) { // 两向量夹角
return acos(dot(a, b) / length(a) / length(b));
}
double distance(Point3 a, Point3 b) { // 两点距离
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y) + (a.z - b.z) * (a.z - b.z));
}
struct Line3 {
Point3 x, y;
Line3() {}
Line3(Point3 x, Point3 y) : x(x), y(y) {}
};
typedef Line3 Segment3;
double PL_distance(Point3 p, Line3 a) { // 点到直线的距离
return length(cross(a.y - a.x, p - a.x)) / distance(a.x, a.y);
}
double PS_distance(Point3 p, Segment3 a) { // 点到线段的距离
if (sign(dot(p - a.x, a.y - a.x)) < 0 || sign(dot(p - a.y, a.x - a.y)) < 0) {
return min(distance(p, a.x), distance(p, a.y));
}
return PL_distance(p, a);
}
bool PL_relation(Point3 p, Line3 a) { // 点在直线上 0在，1不在
return sign(length(cross(a.x - p, a.y - p))) == 0 && sign(dot(a.x - p, a.y - p)) == 0;
}
Point3 PL_project(Point3 p, Line3 a) { // 点在直线上的投影
double k = dot(a.y - a.x, p - a.x) / length(a.y - a.x) / length(a.y - a.x);
return a.x + (a.y - a.x) * k;
}
andrew二维凸包
double andrew(int cnt) { // 二维凸包
sort(q, q + cnt);
cnt = unique(q, q + cnt) - q;
int top = 0;
for (int i = 0; i < cnt; i ++) {
while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) >= 0) {
st[stk[top --]] = false;
}
stk[++ top] = i;
st[i] = true;
}
st[0] = false;
for (int i = cnt - 1; i >= 0; i --) {
if (st[i]) {
continue;
}
while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) >= 0) {
top --;
}
stk[++ top] = i;
}
double res = 0;
for (int i = 2; i <= top; i ++) {
res += get_dist(q[stk[i - 1]], q[stk[i]]);
}
return res;
}
最小圆覆盖
random_shuffle(q, q + n); // 最小圆覆盖
Circle c = {q[0], 0};
for (int i = 1; i < n; i ++) {
if (dcmp(c.r, distance(c.c, q[i])) < 0) {
c = {q[i], 0};
for (int j = 0; j < i; j ++) {
if (dcmp(c.r, distance(c.c, q[j])) < 0) {
c = {(q[i] + q[j]) / 2, distance(q[i], q[j]) / 2};
for (int k = 0; k < j; k ++) {
if (dcmp(c.r, distance(c.c, q[k])) < 0) {
c = get_circle(q[i], q[j], q[k]);
}
}
}
}
}
}
圆的面积并
/*
Created by YuXiao
Nothing to say.
*/
#define YuXiao() cin.tie(0), cout.tie(0), ios::sync_with_stdio(false)
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <vector>
#include <iomanip>
#define x first
#define y second
using namespace std;
typedef pair<double, double> PDD;
const int N = 1010;
const double eps = 1e-8, INF = 1e9;
int n;
PDD q[N];
vector<PDD> segs;
struct Circle
{
PDD r;
double R;
} c[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)
{
int cnt = 0;
for (int i = 0; i < n; i ++)
{
auto X = fabs(x - c[i].r.x), R = c[i].R;
if (dcmp(X, R) < 0)
{
auto Y = sqrt(R * R - X * X);
q[cnt ++] = {c[i].r.y - Y, c[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;
}
return res + ed - st;
}
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)
{
auto mid = (l + r) / 2;
auto left = simpson(l, mid), right = simpson(mid, r);
if (fabs(s - left - right) < eps) return left + right;
return asr(l, mid, left) + asr(mid, r, right);
}
void merge_range()
{
sort(segs.begin(), segs.end());
vector<PDD> new_segs;
double l = -INF, r = -INF;
for (auto &seg : segs)
{
if (dcmp(seg.x, r) <= 0) r = max(r, seg.y);
else
{
if (dcmp(l, -INF)) new_segs.push_back({l, r});
l = seg.x, r = seg.y;
}
}
new_segs.push_back({l, r});
segs = new_segs;
}
int main()
{
YuXiao();
cin >> n;
for (int i = 0; i < n; i ++)
{
cin >> c[i].r.x >> c[i].r.y >> c[i].R;
segs.push_back({c[i].r.x - c[i].R, c[i].r.x + c[i].R});
}
merge_range();
double res = 0;
for (auto &seg : segs) res += asr(seg.x, seg.y, simpson(seg.x, seg.y));
cout << fixed << setprecision(3) << res << endl;
return 0;
}
扫描线
/*
Created by YuXiao
Nothing to say.
*/
#define YuXiao() cin.tie(0), cout.tie(0), ios::sync_with_stdio(false)
#include <iostream>
#include <algorithm>
#include <cstring>
#include <vector>
#define x first
#define y second
using namespace std;
typedef long long LL;
typedef pair<int, int> PII;
const int N = 1010;
int n;
PII l[N], r[N];
PII q[N];
vector<int> xs;
LL 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);
LL 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);
}
int main()
{
YuXiao();
cin >> n;
for (int i = 0; i < n; i ++)
{
cin >> l[i].x >> l[i].y >> r[i].x >> r[i].y;
xs.push_back(l[i].x), xs.push_back(r[i].x);
}
sort(xs.begin(), xs.end());
LL res = 0;
for (int i = 0; i < xs.size() - 1; i ++)
if (xs[i] != xs[i + 1])
res += range_area(xs[i], xs[i + 1]);
cout << res << endl;
return 0;
}
三维凸包
/* 立体（三维）几何：立方体*/
struct CH3D {
struct Plane {
int a, b, c; // 凸包一个面上三个点的编号
bool ok; // 该面是否属于凸包上的面
};
int n, num; // 初始顶点数、凸包表面的三角形数
Point3 P[N]; // 凸包表面的三角形数
Plane F[N * 8];
int g[N][N];
void read(int i) {
cin >> P[i].x >> P[i].y >> P[i].z;
}
double volume4(Point3 a, Point3 b, Point3 c, Point3 d) { // 四面体有向体积 * 6
return dot(cross(b - a, c - a), d - a);
}
double cmp3(Point3 &a, Plane &f) { // 为正：点在面同向
Point3 p1 = P[f.b] - P[f.a];
Point3 p2 = P[f.c] - P[f.a];
Point3 p3 = a - P[f.a];
return dot(cross(p1, p2), p3);
}
void deal(int p, int a, int b) {
int f = g[a][b];
Plane add;
if (F[f].ok) {
if (cmp3(P[p], F[f]) > eps) {
dfs(p, f);
} else {
add.a = b, add.b = a, add.c = p, add.ok = true;
g[p][b] = g[a][p] = g[b][a] = num;
F[num ++] = add;
}
}
}
void dfs(int p, int now) { // 递归搜索所有应该从凸包内删除的面
F[now].ok = false;
deal(p, F[now].b, F[now].a);
deal(p, F[now].c, F[now].b);
deal(p, F[now].a, F[now].c);
}
bool same(int s, int t) {
Point3 &a = P[F[s].a];
Point3 &b = P[F[s].b];
Point3 &c = P[F[s].c];
return fabs(volume4(a, b, c, P[F[t].a])) < eps &&
fabs(volume4(a, b, c, P[F[t].b])) < eps &&
fabs(volume4(a, b, c, P[F[t].c])) < eps;
}
void get_convex3() { // 三维凸包
num = 0;
Plane add;
/* -------------------------保证四个点不共面---------------------------- */
bool flag = true;
for (int i = 1; i < n; i ++) {
if (!(P[0] == P[i])) {
swap(P[1], P[i]);
flag = false;
break;
}
}
if (flag) { return ; } flag = true;
for (int i = 2; i < n; i ++) {
if (length((cross(P[1] - P[0], P[i] - P[0]))) > eps) {
swap(P[2], P[i]);
flag = false;
break;
}
}
if (flag) { return ; } flag = true;
for (int i = 3; i < n; i ++) {
if (fabs(dot(cross(P[1] - P[0], P[2] - P[0]), P[i] - P[0])) > eps) {
swap(P[3], P[i]);
flag = false;
break;
}
}
if (flag) { return ; }
/* ---------------------------------------------------------------------- */
for (int i = 0; i < 4; i ++) {
add.a = (i + 1) % 4, add.b = (i + 2) % 4, add.c = (i + 3) % 4, add.ok = true;
if (cmp3(P[i], add) > 0) {
swap(add.b, add.c);
}
g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = num;
F[num ++] = add;
}
for (int i = 4; i < n; i ++) {
for (int j = 0; j < num; j ++) {
if (F[j].ok && cmp3(P[i], F[j]) > eps) {
dfs(i, j);
break;
}
}
}
int tmp = num; num = 0;
for (int i = 0; i < tmp; i ++) {
if (F[i].ok) {
F[num ++] = F[i];
}
}
}
double convex3_area() { // 表面积
double res = 0;
if (n == 3) {
Vector3 p = cross(P[2] - P[0], P[1] - P[0]);
return length(p) / 2;
}
for (int i = 0; i < num; i ++) {
res += length(area(P[F[i].a], P[F[i].b], P[F[i].c]));
}
return res / 2;
}
double volume() { // 体积
double res = 0;
Point3 tmp = Point3(0, 0, 0);
for (int i = 0; i < num; i ++) {
res += volume4(tmp, P[F[i].a], P[F[i].b], P[F[i].c]);
}
return fabs(res / 6);
}
int triangle() { // 表面三角形个数
return num;
}
int polygon() { // 表面多边形个数
int res = 0;
for (int i = 0; i < num; i ++) {
bool flag = true;
for (int j = 0; j < i; j ++) {
if (same(i, j)) {
flag = false;
break;
}
}
res += flag;
}
return res;
}
Point3 center() { // 重心
Point3 res = Point3(0, 0, 0);
Point3 o = Point3(0, 0, 0);
double all = 0;
for (int i = 0; i < num; i ++) {
double v = volume4(o, P[F[i].a], P[F[i].b], P[F[i].c]);
res = res + (((o + P[F[i].a] + P[F[i].b] + P[F[i].c]) / 4) * v);
all += v;
}
res = res / all;
return res;
}
double PP_distance(Point3 p, int i) { // 点到面的距离
double a = fabs(volume4(P[F[i].a], P[F[i].b], P[F[i].c], p));
double b = length(cross(P[F[i].b] - P[F[i].a], P[F[i].c] - P[F[i].a]));
return a / b;
}
};
三角形面积并
/*
Created by YuXiao
Nothing to say.
*/
// O(n^3logn)
#define YuXiao() cin.tie(0), cout.tie(0), ios::sync_with_stdio(false)
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <vector>
#include <iomanip>
#define x first
#define y second
using namespace std;
typedef pair<double, double> PDD;
const int N = 110;
const double eps = 1e-8, INF = 1e6;
int n;
PDD tr[N][3];
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, double t) { return {a.x * t, a.y * t}; }
double operator* (PDD a, PDD b) { return a.x * b.y - a.y * b.x; }
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; }
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;
}
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;
auto 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_area(double 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;
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_area(a, 1) + line_area(b, 0)) * (b - a) / 2;
}
int main()
{
YuXiao();
cin >> n;
vector<double> xs;
for (int i = 0; i < n; i ++)
{
for (int j = 0; j < 3; j ++)
{
cin >> tr[i][j].x >> tr[i][j].y;
xs.push_back(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 ++)
{
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)) xs.push_back(o.x);
}
sort(xs.begin(), xs.end());
double res = 0;
for (int i = 0; i < xs.size() - 1; i ++)
if (dcmp(xs[i], xs[i + 1]))
res += range_area(xs[i], xs[i + 1]);
cout << fixed << setprecision(2) << res << endl;
return 0;
}
三角剖分
/*
Created by YuXiao
Nothing to say.
*/
#define YuXiao() cin.tie(0), cout.tie(0), ios::sync_with_stdio(false)
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <iomanip>
#define x first
#define y second
using namespace std;
typedef pair<double, double> PDD;
const int N = 60;
const double eps = 1e-8, PI = acos(-1);
double R;
int n;
PDD q[N], 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, double t) { return {a.x * t, a.y * t}; }
double operator* (PDD a, PDD b) { return a.x * b.y - a.y * b.x; }
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; }
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;
}
double area(PDD a, PDD b, PDD c)
{
return (b - a) * (c - a);
}
double get_len(PDD a)
{
return sqrt(a & a);
}
double get_dist(PDD a, PDD b)
{
return get_len(b - a);
}
double project(PDD a, PDD b, PDD c)
{
return ((c - a) & (b - a)) / get_len(b - a);
}
PDD rotate(PDD a, double 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);
}
bool on_segment(PDD p, PDD a, PDD b)
{
return !sign((p - a) * (p - b)) && sign((p - a) & (p - b)) <= 0;
}
PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w)
{
auto u = p - q;
auto t = w * u / (v * w);
return p + v * t;
}
double get_circle_line_intersection(PDD a, PDD b, PDD &pa, PDD &pb)
{
auto e = get_line_intersection(a, b - a, r, rotate(b - a, PI / 2));
auto mind = get_dist(r, e);
if (!on_segment(e, a, b)) mind = min(get_dist(r, a), get_dist(r, b));
if (dcmp(R, mind) <= 0) return mind;
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 mind;
}
double get_sector(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;
}
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 mind = get_circle_line_intersection(a, b, pa, pb);
if (dcmp(R, mind) <= 0) return get_sector(a, b);
if (dcmp(R, da) >= 0) return a * pb / 2 + get_sector(pb, b);
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);
}
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);
}
int main()
{
YuXiao();
while (cin >> R >> n)
{
for (int i = 0; i < n; i ++) cin >> q[i].x >> q[i].y;
cout << fixed << setprecision(2) << work() << endl;
}
return 0;
}
半平面交
double half_plane_intersection() { // 半平面交
sort(l, l + cnt, lcmp);
int hh = 0, tt = -1;
for (int i = 0; i < cnt; i ++) {
if (i && !dcmp(angle(l[i]), angle(l[i - 1]))) continue;
while (hh + 1 <= tt && PL_relation(line_intersection(l[q[tt - 1]], l[q[tt]]), l[i]) == 1) tt --;
while (hh + 1 <= tt && PL_relation(line_intersection(l[q[hh]], l[q[hh + 1]]), l[i]) == 1) hh ++;
q[++ tt] = i;
}
while (hh + 1 <= tt && PL_relation(line_intersection(l[q[tt - 1]], l[q[tt]]), l[hh]) == 1) tt --;
while (hh + 1 <= tt && PL_relation(line_intersection(l[q[hh]], l[q[hh + 1]]), l[tt]) == 1) hh ++;
q[++ tt] = q[hh];
int k = 0;
for (int i = hh; i < tt; i ++) ans[k ++] = line_intersection(l[q[i]], l[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;
}
评论