Link
Luogu - https://www.luogu.org/problemnew/show/SP8073
BZOJ - https://www.lydsy.com/JudgeOnline/problem.php?id=2178
我就不讲Simpson了
Green Formula
对平面单联通区域
D
\rm D
D 及其正向轮廓线
L
\rm L
L (分段光滑),且
P
,
Q
P,Q
P,Q 在
D
\rm D
D 上有一阶连续偏导
∬
D
(
∂
Q
∂
x
−
∂
P
∂
y
)
d
x
d
y
=
∮
L
P
d
x
+
Q
d
y
\mathrm{\iint\limits_{_D}}\left(\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}\right)dxdy=\mathrm{\oint\limits_{_L}}Pdx+Qdy
D∬(∂x∂Q−∂y∂P)dxdy=L∮Pdx+Qdy
取
Q
=
x
,
P
=
−
y
Q=x,P=-y
Q=x,P=−y 得到
∬
D
2
d
x
d
y
=
∮
L
(
−
y
d
x
+
x
d
y
)
\iint\limits_{_{\rm D}}2dxdy=\oint\limits_{_{\rm L}}(\pmb-ydx+xdy)
D∬2dxdy=L∮(−−−ydx+xdy)
那么
∬
D
d
x
d
y
=
1
2
∮
L
(
−
y
d
x
+
x
d
y
)
\iint\limits_{\rm D}dxdy=\frac{1}{2}\oint\limits_{\rm L}(\pmb-ydx+xdy)
D∬dxdy=21L∮(−−−ydx+xdy)
圆
C
C
C 上弧
L
L
L 参数方程
{
x
=
x
0
+
r
cos
t
y
=
y
0
+
r
sin
t
,
t
∈
[
θ
1
,
θ
2
]
\begin{cases}\begin{array}{rcl} x&=&x_0+r\cos t&\\ y&=&y_0+r\sin t& \end{array},\quad t\in[\theta_1,\theta_2]\end{cases}
{xy==x0+rcosty0+rsint,t∈[θ1,θ2]
有
S
=
∮
L
(
−
y
d
x
+
x
d
y
)
=
∫
θ
1
θ
2
−
(
y
0
+
r
sin
t
)
d
(
x
0
+
r
cos
t
)
+
(
x
0
+
r
cos
t
)
d
(
y
0
+
r
sin
t
)
=
⋯
=
∫
θ
1
θ
2
(
r
2
+
x
0
r
cos
t
+
y
0
r
sin
t
)
  
d
t
\begin{array}{rcl} S=\displaystyle\oint\limits_L(\pmb-ydx+xdy)&=&\displaystyle\int\limits_{\theta_1}^{\theta_2}\pmb-(y_0+r\sin t)d(x_0+r\cos t)+(x_0+r\cos t)d(y_0+r\sin t)\\ &=&\cdots\\ &=&\displaystyle\int\limits_{\theta_1}^{\theta_2}(r^2+x_0r\cos t+y_0r\sin t)\;dt \end{array}
S=L∮(−−−ydx+xdy)===θ1∫θ2−−−(y0+rsint)d(x0+rcost)+(x0+rcost)d(y0+rsint)⋯θ1∫θ2(r2+x0rcost+y0rsint)dt
(中间公式被吞过,懒得再打一遍了……反正就那样嘛,不懂微分怎么算的话百度微分运算法则)
然后最后对每个圆算一下产生贡献的弧,这个差分一下就可以辣
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
const int MAXN = 1005;
const double eps = 1e-12;
const double tpi = acos(-1.0);
int n;
double Ans;
struct Circle {
double x, y, r;
inline double f(const double& t) {
return (t * r + sin(t) * x - cos(t) * y) * r;
}
inline double oint(const double& theta_l, const double theta_r) {
return r * (r * (theta_r - theta_l) + x * (sin(theta_r) - sin(theta_l)) - y * (cos(theta_r) - cos(theta_l)));
// return f(theta_r) - f(theta_l);
}
inline bool operator == (const Circle& o) const {
return (x==o.x)&&(y==o.y)&&(r==o.r);
}
inline bool operator < (const Circle& o) const {
return (x==o.x)?((y==o.y)?(r<o.r):(y<o.y)):(x<o.x);
}
}C[MAXN];
inline double defabs(const double&x) {
return (x<0)?(-x):x;
}
int tot;
pair<double, bool> Pt[(MAXN<<1)+5];
inline void Calc(const int& x) {
tot = 0;
register int Sm = 0;
register double DistX, DistY, Dist, Angle, Cosine, Theta, Thetb;
for (register int i = 1; i <= n; ++i) {
if (i == x) continue;
DistX = C[i].x - C[x].x;
DistY = C[i].y - C[x].y;
Dist = sqrt(DistX * DistX + DistY * DistY);
if (C[x].r + Dist <= C[i].r + eps) return;
if (C[i].r + Dist <= C[x].r + eps || Dist >= C[i].r + C[x].r - eps) continue;
Angle = atan2(DistY, DistX),
Cosine = acos((Dist*Dist+C[x].r*C[x].r-C[i].r*C[i].r)/(2.0*Dist*C[x].r));
Theta = Angle - Cosine;
Thetb = Angle + Cosine;
if (Theta < -tpi) Theta += tpi * 2;
if (Thetb >= tpi) Thetb -= tpi * 2;
if (Theta > Thetb) ++Sm;
Pt[++tot] = make_pair(Theta, 0);
Pt[++tot] = make_pair(Thetb, 1);
}
Pt[0] = make_pair(-tpi, 0);
Pt[++tot] = make_pair(tpi, 0);
sort(Pt + 1,Pt + tot);
for (register int i = 1; i <= tot; ++i) {
if (!Sm) Ans += C[x].oint(Pt[i-1].first, Pt[i].first);
(Pt[i].second)?(--Sm):(++Sm);
}
}
int main() {
scanf("%d", &n);
for (register int i = 0, j = 1; j <= n; ++j) {
++i;
scanf("%lf%lf%lf", &C[i].x , &C[i].y, &C[i].r);
}
sort(C + 1, C + 1 + n);
n = unique(C + 1, C + 1 + n) - C - 1;
for (register int i = 1; i <= n; ++i) {
Calc(i);
}
printf("%.3f", Ans / 2);
return 0;
}
没有名字的一般计算几何方法
好!
当然也可以求面积交
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
const int MAXN = 1005;
const double eps = 1e-12;
const double tpi = acos(-1.0);
int n;
double Ans;
struct Circle {
double x, y, r;
inline double f(const double& t) {
return (t * r + sin(t) * x - cos(t) * y) * r;
}
inline double oint(const double& theta_l, const double theta_r) {
return r * (r * (theta_r - theta_l) + x * (sin(theta_r) - sin(theta_l)) - y * (cos(theta_r) - cos(theta_l)));
// return f(theta_r) - f(theta_l);
}
inline bool operator == (const Circle& o) const {
return (x==o.x)&&(y==o.y)&&(r==o.r);
}
inline bool operator < (const Circle& o) const {
return (x==o.x)?((y==o.y)?(r<o.r):(y<o.y)):(x<o.x);
}
}C[MAXN];
inline double defabs(const double&x) {
return (x<0)?(-x):x;
}
int tot;
pair<double, bool> Pt[(MAXN<<1)+5];
inline void Calc(const int& x) {
tot = 0;
register int Sm = 0;
register double DistX, DistY, Dist, Angle, Cosine, Theta, Thetb;
for (register int i = 1; i <= n; ++i) {
if (i == x) continue;
DistX = C[i].x - C[x].x;
DistY = C[i].y - C[x].y;
Dist = sqrt(DistX * DistX + DistY * DistY);
if (C[x].r + Dist <= C[i].r + eps) return;
if (C[i].r + Dist <= C[x].r + eps || Dist >= C[i].r + C[x].r - eps) continue;
Angle = atan2(DistY, DistX),
Cosine = acos((Dist*Dist+C[x].r*C[x].r-C[i].r*C[i].r)/(2.0*Dist*C[x].r));
Theta = Angle - Cosine;
Thetb = Angle + Cosine;
if (Theta < -tpi) Theta += tpi * 2;
if (Thetb >= tpi) Thetb -= tpi * 2;
if (Theta > Thetb) ++Sm;
Pt[++tot] = make_pair(Theta, 0);
Pt[++tot] = make_pair(Thetb, 1);
}
Pt[0] = make_pair(-tpi, 0);
Pt[++tot] = make_pair(tpi, 0);
sort(Pt + 1,Pt + tot);
for (register int i = 1; i <= tot; ++i) {
if (!Sm) Ans += C[x].oint(Pt[i-1].first, Pt[i].first);
(Pt[i].second)?(--Sm):(++Sm);
}
}
int main() {
double Tem = 0;
scanf("%d", &n);
for (register int i = 1; i <= n; ++i) {
scanf("%lf%lf%lf", &C[i].x , &C[i].y, &C[i].r);
Tem += C[i].r * C[i].r * tpi;
}
sort(C + 1, C + 1 + n);
n = unique(C + 1, C + 1 + n) - C - 1;
for (register int i = 1; i <= n; ++i) {
Calc(i);
}
printf("%.3f", Tem - Ans / 2);
return 0;
}