为方便整理,暂时只更新无注释版本(最后更新于2022.5.1)
0x80 前置知识
- 弧长计算公式:
$L=n^{\circ} \times r \times \frac{\pi}{180^{\circ}}, \quad L=\alpha \times r_{\circ}$ 其中$n$是圆心角度数(角度制), $r$是半径, $L$是圆心角弧长, $α$是圆心角度数(弧度制)。
pi = acos(-1)
double R_to_D(double rad)
{
return 180 / pi * rad;
}
double D_to_R(double D)
{
return pi / 180 * D;
}
2. 多边形
2.1 三角形
2.1.1 面积
(1) 叉积
(2) 海伦公式
p = (a + b + c) / 2;
S = sqrt(p(p - a) * (p - b) * (p - c));
2.1.2 三角形四心
(1) 外心,外接圆圆心
三边中垂线交点。到三角形三个顶点的距离相等
(2) 内心,内切圆圆心
角平分线交点,到三边距离相等
(3) 垂心
三条垂线交点
(4) 重心
三条中线交点(到三角形三顶点距离的平方和最小的点,三角形内到三边距离之积最大的点)
2.2 普通多边形
通常按逆时针存储所有点
2.2.1 定义
(1) 多边形
由在同一平面且不再同一直线上的多条线段首尾顺次连接且不相交所组成的图形叫多边形
(2) 简单多边形
简单多边形是除相邻边外其它边不相交的多边形
(3) 凸多边形
过多边形的任意一边做一条直线,如果其他各个顶点都在这条直线的同侧,则把这个多边形叫做凸多边形
任意凸多边形外角和均为360°
任意凸多边形内角和为(n−2)180°
2.2.2 常用函数
(1) 求多边形面积(不一定是凸多边形)
我们可以从第一个顶点除法把凸多边形分成n − 2个三角形,然后把面积加起来。
double polygon_area(Point p[], int n)
{
double s = 0;
for (int i = 1; i + 1 < n; i ++ )
s += cross(p[i] - p[0], p[i + 1] - p[i]);
return s / 2;
}
(2) 判断点是否在多边形内(不一定是凸多边形)
a. 射线法,从该点任意做一条和所有边都不平行的射线。交点个数为偶数,则在多边形外,为奇数,则在多边形内。
b. 转角法
(3) 判断点是否在凸多边形内
只需判断点是否在所有边的左边(逆时针存储多边形)。
2.3 皮克定理
皮克定理是指一个计算点阵中顶点在格点上的多边形面积公式该公式可以表示为:
S = a + b/2 - 1
其中a表示多边形内部的点数,b表示多边形边界上的点数,S表示多边形的面积。
0x81 基础操作
const double dinf = 1e18, eps = 1e-8;
int sgn(double x)
{
if(fabs(x) < eps) return 0;
if(x < 0) return -1;
return 1;
}
struct Point
{
double x, y;
Point() {}
Point(double x, double y) : x(x), y(y) {}
Point operator + (Point B) { return Point(x + B.x, y + B.y); }
Point operator - (Point B) { return Point(x - B.x, y - B.y); }
Point operator * (double k) { return Point(x * k, y * k); }
Point operator / (double k) { return Point(x / k, y / k); }
bool operator == (Point B) { return sgn(x - B.x) == 0 && sgn(y - B.x) == 0; }
bool operator < (Point B) { return sgn(x - B.x) < 0 || sgn(x - B.x) == 0 && sgn(y - B.y) < 0; }
};
typedef Point Vector;
double Dot(Vector A, Vector B) {return A.x * B.x + A.y * B.y;}
double Cross(Vector A, Vector B) {return A.x * B.y - A.y * B.x;}
double area(Point A, Point B, Point C) {return Cross(B - A, C - A);}
int dcmp(double x, double y)
{
if(fabs(x - y) < eps) return 0;
if(x < y) return -1;
return 1;
}
double get_dist(Point A, Point B)
{
double dx = A.x - B.x;
double dy = A.y - B.y;
return sqrt(dx * dx + dy * dy);
}
// 直线AB与点P的关系
int relation(Point A, Point B, Point C)
{
// 1 left -1 right 0 in
int c = sgn(Cross((B - A), (C - A)));
if(c < 0) return 1;
else if(c > 0) return -1;
return 0;
}
// 线段AB与点P的关系
int relation(Point A, Point B, Point C)
{
// 1 left -1 right/out 0 in
int c = sgn(Cross((B - A), (C - A)));
if(c < 0) return 1;
else if(c > 0 || !on_segment(C, A, B)) return -1;
return 0;
}
double get_length(Point A)
{
return sqrt(Dot(A, A));
}
// 计算向量夹角
double get_angle(Point A, Point B)
{
return acos(Dot(A, B) / get_length(A) / get_length(B));
}
// 计算两个向量构成的平行四边形有向面积
double area(Point A, Point B, Point C)
{
return Cross(B - A, C - A);
}
// 向量A顺时针旋转C的角度:
Point rotate(Point A, double angle)
{
return Point(A.x * cos(angle) + A.y * sin(angle), - A.x * sin(angle) + A.y * cos(angle));
}
(1) 判断点在直线上 A x B = 0
(2) 两直线相交
// Cross(v, w) == 0则两直线平行或者重合
Point get_intersection(Line a, Line b)
{
Vector u = a.P - b.P;
double t = Cross(b.v, u) / Cross(a.v, b.v);
return a.P + a.v * t;
}
(3) 点到直线的距离
double distance_to_line(Point P, Point A, Point B)
{
Vector v1 = B - A, v2 = P - A;
return fabs(Cross(v1, v2) / get_length(v1));
}
(4) 点到线段的距离
double distance_to_segment(Point P, Point A, Point B)
{
if (A == B) return get_length(P - A);
Vector v1 = B - A, v2 = P - A, v3 = P - B;
if (sgn(Dot(v1, v2)) < 0) return get_length(v2);
if (sgn(Dot(v1, v3)) > 0) return get_length(v3);
return distance_to_line(P, A, B);
}
(5) 点在直线上的投影长度
double get_project(Point P, Point A, Point B)
{
Vector v1 = B - A, v2 = P - A;
return fabs(Dot(v1, v2) / get_length(v1));
}
(6) 点是否在线段上
bool on_segment(Point P, Point A, Point B)
{
return sgn(Cross(P - A, P - B)) == 0 && sgn(Dot(P - A, P - B)) <= 0;
}
(7) 判断两线段是否相交
bool segment_intersection(Point a1, Point a2, Point b1, Point b2)
{
double c1 = Cross(a2 - a1, b1 - a1), c2 = Cross(a2 - a1, b2 - a1);
double c3 = Cross(b2 - b1, a2 - b1), c4 = Cross(b2 - b1, a1 - b1);
return sgn(c1) * sgn(c2) <= 0 && sgn(c3) * sgn(c4) <= 0;
}
(8) 求向量的单位向量
Point norm(Point a)
{
return a / get_length(a);
}
0x82 凸包
Andrew法求凸包(栈中元素即为构成凸包的顶点)
//求凸包周长
int n, top;
Point q[N];
int stk[N];
bool st[N];
double andrew()
{
sort(q, q + n);
for (int i = 0; i < n; i ++ )
{
while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) <= 0)
{
if (area(q[stk[top - 1]], q[stk[top]], q[i]) < 0)
st[stk[top -- ]] = false;
else top -- ;
}
stk[ ++ top] = i;
st[i] = true;
}
st[0] = false;
for (int i = n - 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;
}
0x83 求半平面交
struct Line
{
Point P;
Vector v;
double deg;
Line() {}
Line(Point P, Vector v) : P(P), v(v) { deg = atan2(v.y, v.x); }
bool operator < (const Line &L) const { return deg < L.deg; }
double value(double x)
{
return P.y + (x - P.x) * v.y / v.x;
}
};
bool on_left(Line L, Point P) { return sgn(Cross(L.v, P - L.P)) > 0; }
Point get_intersection(Line a, Line b)
{
Vector u = a.P - b.P;
double t = Cross(b.v, u) / Cross(a.v, b.v);
return a.P + a.v * t;
}
vector<Point> half_plane_intersection(vector<Line> &L)
{
int n = L.size();
sort(L.begin(), L.end());
int first, last;
vector<Point> p(n);
vector<Line> q(n);
vector<Point> ans;
q[first = last = 0] = L[0];
for(int i = 1; i < n; i ++ )
{
while(first < last && !on_left(L[i], p[last - 1])) last --;
while(first < last && !on_left(L[i], p[first])) first ++;
q[ ++ last] = L[i];
if(fabs(Cross(q[last].v, q[last - 1].v)) < eps)
{
last--;
if(on_left(q[last], L[i].P)) q[last] = L[i];
}
if(first < last) p[last - 1] = get_intersection(q[last - 1], q[last]);
}
while(first < last && !on_left(q[first], p[last - 1])) last --;
if(last - first <= 1) return ans;
p[last] = get_intersection(q[last], q[first]);
for(int i = first; i <= last; i ++ ) ans.push_back(p[i]);
return ans;
}
double Polygon_area(vector<Point> &p)
{
int n = p.size();
double area = 0;
for(int i = 0; i < n; i ++ ) area += Cross(p[i], p[(i + 1) % n]);
return area / 2;
}
0x84 最小圆覆盖
Point get_intersection(Line a, Line b)
{
Vector u = a.P - b.P;
double t = Cross(b.v, u) / Cross(a.v, b.v);
return a.P + a.v * t;
}
Line get_line(Point a, Point b) {
return {(a + b) / 2, rotate(b - a, pi / 2)};
}
struct Circle
{
Point P;
double r;
Circle(Point P = Point(), double r = 0) : P(P), r(r) {}
Point point(double a) //通过圆心角求坐标
{
return Point(P.x + cos(a) * r, P.y + sin(a) * r);
}
};
Circle get_circle(Point a, Point b, Point c)
{
auto u = get_line(a, b), v = get_line(a, c);
auto p = get_intersection(u, v);
return {p, get_dist(p, a)};
}
int n;
Point q[N];
signed main()
{
//ios::sync_with_stdio(0); cin.tie(0); cin.tie(0);
cin >> n;
for(int i = 0; i < n; i ++ ) cin >> q[i].x >> q[i].y;
random_shuffle(q, q + n);
/*
一个斜着的的椭圆形,它的长轴距离x轴的正方向有a度
把所有点先向右旋转a度,然后再缩小P倍数, 然后求最小圆的覆盖就行
double a, p;
cin >> a >> p;
for (int i = 0; i < n; i ++ )
{
q[i] = rotate(q[i], a / 180 * pi);
q[i].x /= p;
}
*/
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_circle(q[i], q[j], q[k]);
}
}
printf("%.10lf\n", c.r);
printf("%.10lf %.10lf\n", c.P.x, c.P.y);
return 0;
}
0x85 三维凸包
// double eps = 1e-10
double rand_eps()
{
return ((double)rand() / RAND_MAX - 0.5) * eps;
}
struct Point_3D
{
double x, y, z;
void shake()
{
x += rand_eps(), y += rand_eps(), z += rand_eps();
}
Point_3D(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}
bool operator == (const Point_3D &A) const
{
return x == A.x && y == A.y && z == A.z;
}
double len()
{
return sqrt(x * x + y * y + z * z);
}
};
typedef Point_3D Vector_3D;
Vector_3D operator+(const Vector_3D &A, const Vector_3D &B) { return Vector_3D(A.x + B.x, A.y + B.y, A.z + B.z); }
Vector_3D operator-(const Point_3D &A, const Point_3D &B) { return Vector_3D(A.x - B.x, A.y - B.y, A.z - B.z); }
Vector_3D operator*(const Vector_3D &A, double p) { return Vector_3D(A.x * p, A.y * p, A.z * p); }
Vector_3D operator/(const Vector_3D &A, double p) { return Vector_3D(A.x / p, A.y / p, A.z / p); }
double Dot_3D(const Vector_3D &A, const Vector_3D &B) { return A.x * B.x + A.y * B.y + A.z * B.z; }
Vector_3D Cross_3D(const Vector_3D &A, const Vector_3D &B) { return Vector_3D(A.y * B.z - A.z * B.y, A.z * B.x - A.x * B.z, A.x * B.y - A.y * B.x); }
//转单位向量
Vector_3D to_unitVector(Vector_3D u)
{
double k = sqrt(Dot_3D(u, u));
return Vector_3D(u / k);
}
struct Plane
{
int v[3];
Plane(int a, int b, int c)
{
v[0] = a;
v[1] = b;
v[2] = c;
}
Vector_3D Normal(const vector<Point_3D> &P) const
{
return Cross(P[v[1]] - P[v[0]], P[v[2]] - P[v[0]]);
}
int CanSee(const vector<Point_3D> &P, int i) const
{
return Dot_3D(P[i] - P[v[0]], Normal(P)) > 0;
}
double area(const vector<Point_3D> &P) // 求面积
{
return Normal(P).len() / 2;
}
};
vector<Plane> get_convex_3d(const vector<Point_3D> & P)
{
int n = P.size();
vector vis(n, vector<int>(n));
vector<Plane> cur;
cur.push_back(Plane(0, 1, 2));
cur.push_back(Plane(2, 1, 0));
for(int i = 3; i < n; i++)
{
vector<Plane> next;
// 计算每条边的“左面”的可见性
for(int j = 0; j < cur.size(); j++)
{
Plane & f = cur[j];
int res = f.CanSee(P, i);
if(!res) next.push_back(f);
for(int k = 0; k < 3; k ++ )
{
vis[f.v[k]][f.v[(k + 1) % 3]] = res;
}
}
for(int j = 0; j < cur.size(); j ++ )
for(int k = 0; k < 3; k ++ )
{
int a = cur[j].v[k], b = cur[j].v[(k + 1) % 3];
if(vis[a][b] != vis[b][a] && vis[a][b]) next.push_back(Plane(a, b, i));
}
cur = next;
}
return cur;
}
int n;
signed main()
{
ios::sync_with_stdio(0); cin.tie(0); cin.tie(0);
cin >> n;
vector<Point_3D> q(n);
for (int i = 0; i < n; i ++ )
{
cin >> q[i].x >> q[i].y >> q[i].z;
q[i].shake();
}
vector<Plane> res = get_convex_3d(q);
double ans = 0;
for (auto c : res) ans += c.area(q);
cout << fixed << setprecision(6) << ans << '\n';
return 0;
}
0x86 旋转卡壳
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)
st[stk[ -- top]] = false;
else top -- ;
}
stk[top ++ ] = i;
st[i] = true;
}
st[0] = false;
for (int i = n - 1; i >= 0; i -- )
{
if (st[i]) continue;
while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0)
top -- ;
stk[top ++ ] = i;
}
top -- ;
}
// 求最远点距
double rotating_calipers()
{
if (top <= 2) return get_dist(q[0], q[n - 1]);
double res = 0;
for (int i = 0, j = 2; i < top; i ++ )
{
Point 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;
}
//求最小矩形覆盖
int dcmp(double x, double y)
{
if(fabs(x - y) < eps) return 0;
if(x < y) return -1;
return 1;
}
Point rotate(Point A, double angle)
{
return Point(A.x * cos(angle) + A.y * sin(angle), - A.x * sin(angle) + A.y * cos(angle));
}
double get_length(Point A)
{
return sqrt(Dot(A, A));
}
double get_project(Point a, Point b, Point c)
{
return Dot((b - a), (c - a)) / get_length(b - a);
}
Point norm(Point a)
{
return a / get_length(a);
}
int n, top;
double min_area = dinf;
Point ans[5];
Point q[N];
int stk[N];
bool st[N];
// andrew() 同上
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(get_project(d, e, q[stk[b]]), get_project(d, e, q[stk[b + 1]])) < 0) b = (b + 1) % top;
if (!i) c = a;
while (dcmp(get_project(d, e, q[stk[c]]), get_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_length(e - d);
auto w = Dot((y - z), (e - d)) / get_length(e - d);
if (h * w < min_area)
{
min_area = h * w;
ans[0] = d + norm(e - d) * get_project(d, e, y);
ans[3] = d + norm(e - d) * get_project(d, e, z);
auto u = norm(rotate(e - d, -pi / 2));
ans[1] = ans[0] + u * h;
ans[2] = ans[3] + u * h;
}
}
}
0x87 三角刨分
//求圆心在坐标点的圆与多边形的面积交
double get_circle_line_intersection(Point a, Point b, Point& pa, Point& pb)
{
auto e = get_line_intersection(Line(a, b - a), Line(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(Point a, Point b)
{
auto angle = acos((Dot(a, b)) / get_length(a) / get_length(b));
if (sgn(Cross(a, b)) < 0) angle = -angle;
return R * R * angle / 2;
}
double get_circle_triangle_area(Point a, Point b)
{
auto da = get_dist(r, a), db = get_dist(r, b);
if (dcmp(R, da) >= 0 && dcmp(R, db) >= 0) return Cross(a, b) / 2;
if (!sgn(Cross(a, b))) return 0;
Point 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 Cross(a, pb) / 2 + get_sector(pb, b);
if (dcmp(R, db) >= 0) return get_sector(a, pa) + Cross(pa, b) / 2;
return get_sector(a, pa) + Cross(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);
}
0x88 扫描线
//求矩形面积并(输入左上和右下端点)
int n;
PII l[N], r[N];
PII q[N];
int range_area(int a, int b)
{
int cnt = 0;
for(int i = 0; i < n; i ++ )
if(l[i].first <= a && r[i].first >= b)
q[cnt ++ ] = {l[i].second, r[i].second};
if(!cnt) return 0;
sort(q, q + cnt);
int res = 0;
int st = q[0].first, ed = q[0].second;
for (int i = 1; i < cnt; i ++ )
if (q[i].first <= ed) ed = max(ed, q[i].second);
else
{
res += ed - st;
st = q[i].first, ed = q[i].second;
}
res += ed - st;
return res * (b - a);
}
signed main()
{
ios::sync_with_stdio(0); cin.tie(0); cin.tie(0);
cin >> n;
vector<int> xs;
for(int i = 0; i < n; i ++ )
{
cin >> l[i].first >> l[i].second >> r[i].first >> r[i].second;
xs.emplace_back(l[i].first), xs.emplace_back(r[i].first);
}
sort(xs.begin(), xs.end());
int res = 0;
for(int i = 0; i + 1 < xs.size(); i ++ )
if(xs[i] != xs[i + 1])
res += range_area(xs[i], xs[i + 1]);
cout << res << '\n';
return 0;
}
//求三角形面积并
bool on_segment(Point P, Point A, Point B) //在线段内
{
return sgn(Dot((P - A), (P - B))) <= 0;
}
struct Line
{
Point P;
Vector v;
double deg;
Line() {}
Line(Point P, Vector v) : P(P), v(v) { deg = atan2(v.y, v.x); }
bool operator < (const Line &L) const { return deg < L.deg; }
double value(double x)
{
return P.y + (x - P.x) * v.y / v.x;
}
};
Point get_intersection(Point p, Point v, Point q, Point w)
{
if (!sgn(Cross(v, w))) return {dinf, dinf};
auto u = p - q;
auto t = Cross(w, u) / Cross(v, w);
auto o = p + v * t;
if (!on_segment(o, p, p + v) || !on_segment(o, q, q + w))
return {dinf, dinf};
return o;
}
int n;
Point tr[N][3];
Point q[N];
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_intersection(t[j], t[(j + 1) % 3] - t[j], {a, -dinf}, {0, dinf * 2});
if (dcmp(o.x, dinf))
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;
}
signed main()
{
//ios::sync_with_stdio(0); cin.tie(0); cin.tie(0);
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_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, dinf)) xs.push_back(o.x);
}
sort(xs.begin(), xs.end());
double res = 0;
for (int i = 0; i + 1 < xs.size(); i ++ )
if (dcmp(xs[i], xs[i + 1]))
res += range_area(xs[i], xs[i + 1]);
cout << fixed << setprecision(2) << res << '\n';
return 0;
}
0x89 自适应辛普森积分
// 求任意函数在一段区间内的面积
// const double eps = 1e-12
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;
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);
}
// printf("%.6lf\n", asr(a, b, simpson(a, b))); //a, b为积分上下限
// 求圆的面积并
typedef pair<double, double> PDD;
int n;
struct Circle
{
PDD r;
double R;
}c[N];
PDD q[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.first), R = c[i].R;
if (dcmp(X, R) < 0)
{
auto Y = sqrt(R * R - X * X);
q[cnt ++ ] = {c[i].r.second - Y, c[i].r.second + Y};
}
}
if (!cnt) return 0;
sort(q, q + cnt);
double res = 0, st = q[0].first, ed = q[0].second;
for (int i = 1; i < cnt; i ++ )
if (q[i].first <= ed) ed = max(ed, q[i].second);
else
{
res += ed - st;
st = q[i].first, ed = q[i].second;
}
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);
}
signed main()
{
ios::sync_with_stdio(0); cin.tie(0); cin.tie(0);
cin >> n;
//double l = 2000, r = -2000;
for (int i = 0; i < n; i ++ )
{
cin >> c[i].r.first >> c[i].r.second >> c[i].R;
l = min(l, c[i].r.first - c[i].R), r = max(r, c[i].r.first + c[i].R);
}
cout << fixed << setprecision(3) << asr(l - 100, r + 100, simpson(l, r)) << '\n';
return 0;
}
0x90 罗德里格斯旋转方程$(Rodrigues)$
罗德里格斯旋转公式,用于表示空间中任一向量$\vec{v}$,沿任一旋转轴$\vec{k}$,旋转任一角度$\theta$后,得到的结果:
$\vec{v}_{r o t}=\vec{v} \cos \theta+(1-\cos \theta)(\vec{k} \cdot \vec{v}) \cdot \vec{k}+\sin \theta * \vec{k} \times \vec{v}$
//罗格里格斯公式 (r为弧度制, 向量v沿着单位向量u旋转r角度)
Vector_3D cal(Vector_3D u, Vector_3D v, double r)
{
return v * cos(r) + u * Dot_3D(u, v) * (1 - cos(r)) + Cross_3D(u, v) * sin(r);
}