bzoj 2178 圆的面积并【辛普森积分】

最新推荐文章于 2019-04-22 14:58:00 发布

Neo__Z

最新推荐文章于 2019-04-22 14:58:00 发布

阅读量460

点赞数

CC 4.0 BY-SA版权

分类专栏：辛普森积分 --------计算几何-------- bzoj 文章标签：圆的面积并辛普森积分 bzoj

本文链接：https://blog.youkuaiyun.com/cdsszjj/article/details/78738891

bzoj 同时被 3 个专栏收录

145 篇文章

订阅专栏

--------计算几何--------

43 篇文章

订阅专栏

辛普森积分

1 篇文章

订阅专栏

本文介绍了一种利用辛普森积分法求解由多个圆构成的区域面积的方法。通过将直线与圆的交点长度总和作为积分函数，并采用递归方式逼近精确面积，实现了高效准确的计算。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

解题思路：

辛普森积分可以求过三点 $(l,f(l)),(mid,f(mid)),(r,f(r))$ ，其中 $mid=(l+r)/2$ 的抛物线在 $(l,r)$ 一段的积分。其公式为：

$\int _l^rf(x)=\frac{(r-l)}{6}(f(l)+4f(mid)+f(r))$

在这道题中，我们可以将 $f(i)$ 取成直线 $x=i$ 与所有圆的截线总长(重叠部分只算一次)，就可以算 $(l,r)$ 一段所截的面积了。
由于不是标准的抛物线形式，所以我们可以递归求解，直到 $ansl+ansr$ 和 $ans$ 的误差在eps以内。这道题eps要取到1e-13。
还有个优化就是先去掉被其他圆包含的圆。

#include<cstdio>
#include<iostream>
#include<cstring>
#include<string>
#include<algorithm>
#include<cmath>
#include<ctime>
#include<vector>
#include<set>
#include<complex>
#define ll long long
using namespace std;

int getint()
{
    int i=0,f=1;char c;
    for(c=getchar();(c<'0'||c>'9')&&c!='-';c=getchar());
    if(c=='-')f=-1,c=getchar();
    for(;c>='0'&&c<='9';c=getchar())i=(i<<3)+(i<<1)+c-'0';
    return i*f;
}

const int N=1005;
const double eps=1e-13,INF=1e9;
struct point
{
    double x,y;
    point(){}
    point(double _x,double _y):
        x(_x),y(_y){}
    inline friend point operator - (const point &a,const point &b)
    {return point(a.x-b.x,a.y-b.y); }
    inline double dis() {return sqrt(x*x+y*y);} 
    inline friend bool operator <(const point &a,const point &b)
    {
        if(fabs(a.x-b.x)<eps)return a.y>b.y;
        return a.x<b.x;
    }
}l[N];
struct circle
{
    point O;double r;
    inline friend bool operator < (const circle &a,const circle &b)
    {return a.r<b.r;}
}C[N];
int n;

void pre()
{
    sort(C+1,C+n+1);
    int tmp=0;
    for(int i=1,j;i<=n;i++)
    {
        for(j=i+1;j<=n;j++)
            if((C[i].O-C[j].O).dis()+C[i].r-C[j].r<eps)break;
        if(j==n+1)C[++tmp]=C[i];
    }
    n=tmp;
}

double Cut(double x)
{
    int cnt=0;double len;
    for(int i=1;i<=n;i++)
        if(abs(C[i].O.x-x)<C[i].r)
        {
            len=C[i].O.x-x;
            len=sqrt(C[i].r*C[i].r-len*len);
            l[++cnt]=point(C[i].O.y-len,C[i].O.y+len);
        }
    sort(l+1,l+cnt+1);
    int top=0;
    for(int i=1;i<=cnt;i++)
        if(!top||l[i].y>l[top].y)l[++top]=l[i];
    len=0;
    for(int i=1;i<=top;i++)
        if(i==1||l[i].x>l[i-1].y)len+=l[i].y-l[i].x;
        else len+=l[i].y-l[i-1].y;
    return len;
}

double Simpson(double l,double mid,double r,double fl,double fmid,double fr)
{
    double lenl=Cut((l+mid)/2),lenr=Cut((mid+r)/2);
    double ans=(r-l)*(fl+4*fmid+fr)/6,ansl=(mid-l)*(fl+4*lenl+fmid)/6,ansr=(r-mid)*(fmid+4*lenr+fr)/6;
    if(abs(ans-ansl-ansr)<eps)return ans;
    return Simpson(l,(l+mid)/2,mid,fl,lenl,fmid)+Simpson(mid,(mid+r)/2,r,fmid,lenr,fr);
}

int main()
{
    //freopen("lx.in","r",stdin);
    //freopen("lx.out","w",stdout);
    n=getint();
    double l=INF,r=-INF;
    for(int i=1;i<=n;i++)
    {
        C[i].O.x=getint(),C[i].O.y=getint(),C[i].r=getint();
        l=min(l,C[i].O.x-C[i].r),r=max(r,C[i].O.x+C[i].r);
    }   
    pre();
    printf("%0.3f",Simpson(l,(l+r)/2,r,0,Cut((l+r)/2),0));
    return 0;
}