多个圆面积的并（自适应辛普森）
问题引入
题目大意：给出n个圆，求这n个圆面积的并
输入格式：
第一行输入一个正整数 n 表示圆的个数 （1<=n<=1000）
接下来n行，每行三个正整数 x , y, r，表示圆心的坐标和圆的半径
输出格式：
输出一个实数表示答案，保留三位小数
这道题可以用自适应辛普森乱搞，这里先简单入门一下什么是自适应辛普森（Simpson）积分
自适应辛普森
定积分：
定积分就是求f(x)在区间[a, b]中的图像面积，如下图：
即阴影部分面积为：
但是大部分函数的原函数并不存在，或者求解过程极其复杂，因此可以使用自适应辛普森法递归求解近似值。
辛普森公式：
学习Simpson公式最好是要对积分有一定了解，但是不了解也没关系（像我这样高数忘光的人只能生搬硬套）
如果感兴趣可以自己去推导
自适应辛普森
有了Simpson公式之后，我们便可以把区间划分为一段段的小区间，然后在计算求和
通过递归函数integral(double L,double R)，每一次递归将区间 [L, R] 分成 [L,mid]和[mid,R], 其中mid=(L+R)/2，再分别计算 f(x) 在区间[L,R]、[L,mid] 、[mid, R] 的积分Sum, Sl, Sr， 当Sum和Sl+Sr的误差满足题目的精度要求时，便可退出递归。
double F(double x){ //相应的积分函数 return (c*x+d)/(a*x+b); } double Simpson(double a,double b){ //辛普森公式 double c=(a+b)/2; return (b-a)*(F(a)+F(b)+4*F(c))/6; } double integral(double L,double R){ //递归求解 double mid=(L+R)/2; double S=Simpson(L,R),Sl=Simpson(L,mid),Sr=Simpson(mid,R); if(fabs(Sl+Sr-S)<eps) return Sl+Sr; else return integral(L,mid)+integral(mid,R); }
例题:
poj4525
https://www.luogu.com.cn/problem/P4525
ac代码：
#include<iostream> #include<cstdio> #include<algorithm> #include<cstring> #include<cmath> #include<vector> #include<queue> #include<stack> #include<map> #define INF 0x3f3f3f3f using namespace std; typedef long long LL; typedef unsigned long long ULL; const int maxn=2e5+6; const int mod=1e9+7; const double pi=acos(-1.0); const double eps=1e-11; double a,b,c,d,l,r; double F(double x){ return (c*x+d)/(a*x+b); } double Simpson(double a,double b){ double c=(a+b)/2; return (b-a)*(F(a)+F(b)+4*F(c))/6; } double integral(double L,double R){ double mid=(L+R)/2; double S=Simpson(L,R),Sl=Simpson(L,mid),Sr=Simpson(mid,R); if(fabs(Sl+Sr-S)<eps) return Sl+Sr; else return integral(L,mid)+integral(mid,R); } int main(){ scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&l,&r); printf("%.6lf\n",integral(l,r)); return 0; }
HDU1724 Ellipse
http://acm.hdu.edu.cn/showproblem.php?pid=1724
ac代码：
#include<iostream> #include<cstdio> #include<algorithm> #include<cstring> #include<cmath> #include<vector> #include<queue> #include<stack> #include<map> #define INF 0x3f3f3f3f using namespace std; typedef long long LL; typedef unsigned long long ULL; const int maxn=2e5+6; const int mod=1e9+7; const double pi=acos(-1.0); const double eps=1e-11; double a,b,l,r; double F(double x){ return b*sqrt(1-x*x/(a*a)); } double Simpson(double a,double b){ double c=(a+b)/2; return (b-a)*(F(a)+F(b)+4*F(c))/6; } double integral(double L,double R){ double mid=(L+R)/2; double S=Simpson(L,R),Sl=Simpson(L,mid),Sr=Simpson(mid,R); if(fabs(Sl+Sr-S)<eps) return Sl+Sr; else return integral(L,mid)+integral(mid,R); } int main(){ int t; scanf("%d",&t); while(t--){ scanf("%lf%lf%lf%lf",&a,&b,&l,&r); printf("%.3lf\n",2*integral(l,r)); } return 0; }
大概了解了自适应辛普森后，我们回到求解多个圆面积的交的问题上
例如要求上图5个圆相交形成的红色部分的面积，我们可以将这些圆看成一个整体的函数，然后就可以用辛普森乱搞了。
首先对于被其他圆包含的圆，对答案并不会产生贡献，因此我们可以先预处理，遍历一遍所有圆，将被其他圆包含的圆去除
（圆已先按半径大小从小到大排序）
void init(){ //预处理，删除被包含的圆 int cnt=0; int i,j; for(i=1;i<=n;i++){ for(j=i+1;j<=n;j++){ if(c[i].r-c[j].r+dis(c[i].o,c[j].o)<eps) break; } if(j>n) c[++cnt]=c[i]; } n=cnt; }
之后将圆按最左端点的x坐标（圆心x坐标减去半径即可，最右端则加上半径）从小到大排序
bool cmp(cir a,cir b){ return a.o.x-a.r<b.o.x-b.r; }
接下来就可以辛普森乱搞了
遍历所有圆，每找出一段连通的圆，便做一次自适应辛普森，累加答案即可
（如图可以分成3段）
void solve(){ //寻找每段联通的圆，进行计算 double l,r; int i=1,j; while(i<=n){ l=c[i].o.x-c[i].r,r=c[i].o.x+c[i].r; for(j=i+1;j<=n;j++){ if(c[j].o.x-c[j].r>r) break; //圆j最左端的x坐标比前面圆的最右端还大，说明不满足上述要求 r=max(r,c[j].o.x+c[j].r); } st=i,ed=j-1,i=j; ans+=integral(l,r); } }
接下来就是直接套上自适应辛普森的板子，这里最关键的就是已知 x 的值，怎么求它对应的 f(x)
若 x 值为 a，f(x)其实就求直线x=a, 被圆覆盖的长度，如下图红色实线部分
所以我们可以求出直线x=a与所有圆的交点，做贪心覆盖，求其长度
double F(double x){ //求F(x) v.clear(); for(int i=st;i<=ed;i++){ //处理出交点的纵坐标，保存为线段 if(x>c[i].o.x+c[i].r||x<c[i].o.x-c[i].r) continue; double l=x-c[i].o.x; l=sqrt(c[i].r*c[i].r-l*l); v.push_back((seg){c[i].o.y-l,c[i].o.y+l}); } if(v.empty()) return 0.00; sort(v.begin(),v.end()); //开始贪心覆盖 double s=0,last=v[0].l; for(int i=0;i<v.size();i++){ if(v[i].r>last){ s+=v[i].r-max(last,v[i].l); last=v[i].r; } } return s; }
题目链接：https://www.luogu.com.cn/problem/SP8073
ac代码：
#include<iostream> #include<cstdio> #include<algorithm> #include<cstring> #include<cmath> #include<vector> #include<queue> #include<stack> #include<map> #define INF 1e16 using namespace std; typedef long long LL; typedef unsigned long long ULL; const int maxn=2e5+6; const int mod=1e9+7; const double pi=acos(-1.0); const double eps=1e-7; struct Point{ double x,y; Point(double x=0,double y=0):x(x),y(y) {} }; struct seg{ double l,r; friend bool operator <(seg a,seg b){ if(a.l==b.l) return a.r<b.r; return a.l<b.l; } }; struct cir{ Point o; double r; friend bool operator <(cir a,cir b){ return a.r<b.r; } }c[maxn]; double dis(Point a,Point b){ return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); } int n; vector<seg>v; double ans=0; int st,ed; double F(double x){ //求F(x) v.clear(); for(int i=st;i<=ed;i++){ //处理出交点的纵坐标，保存为线段 if(x>c[i].o.x+c[i].r||x<c[i].o.x-c[i].r) continue; double l=x-c[i].o.x; l=sqrt(c[i].r*c[i].r-l*l); v.push_back((seg){c[i].o.y-l,c[i].o.y+l}); } if(v.empty()) return 0.00; sort(v.begin(),v.end()); //开始贪心覆盖 double s=0,last=v[0].l; for(int i=0;i<v.size();i++){ if(v[i].r>last){ s+=v[i].r-max(last,v[i].l); last=v[i].r; } } return s; } double Simpson(double a,double b){ //自适应辛普森 double c=(a+b)/2; return (b-a)*(F(a)+F(b)+4*F(c))/6; } double integral(double L,double R){ double mid=(L+R)/2; double S=Simpson(L,R),Sl=Simpson(L,mid),Sr=Simpson(mid,R); if(fabs(Sl+Sr-S)<eps) return Sl+Sr; else return integral(L,mid)+integral(mid,R); } void init(){ //预处理，删除被包含的圆 int cnt=0; int i,j; for(i=1;i<=n;i++){ for(j=i+1;j<=n;j++){ if(c[i].r-c[j].r+dis(c[i].o,c[j].o)<eps) break; } if(j>n) c[++cnt]=c[i]; } n=cnt; } bool cmp(cir a,cir b){ return a.o.x-a.r<b.o.x-b.r; } void solve(){ //寻找每段联通的圆，进行计算 double l,r; int i=1,j; while(i<=n){ l=c[i].o.x-c[i].r,r=c[i].o.x+c[i].r; for(j=i+1;j<=n;j++){ if(c[j].o.x-c[j].r>r) break; r=max(r,c[j].o.x+c[j].r); } st=i,ed=j-1,i=j; ans+=integral(l,r); } } int main(){ scanf("%d",&n); double l=INF,r=-INF; for(int i=1;i<=n;i++){ scanf("%lf%lf%lf",&c[i].o.x,&c[i].o.y,&c[i].r); l=min(l,c[i].o.x-c[i].r); r=max(r,c[i].o.x+c[i].r); } sort(c+1,c+1+n); init(); sort(c+1,c+1+n,cmp); solve(); printf("%.3lf\n",ans); return 0; }
