bzoj2178 圆的面积并（Simpson积分）

最新推荐文章于 2020-06-03 07:32:15 发布

原创最新推荐文章于 2020-06-03 07:32:15 发布 · 416 阅读

0 ·

本内容遵循CC 4.0 BY-SA版权协议

标签

#c++

Simpson 专栏收录该内容

2 篇文章

订阅专栏

博客围绕用辛普森积分求解圆形面积并展开。指出求图形面积本质可用Simpson解决，定义函数f(k)，通过勾股定理和扫描线方法求该函数。还提到因原写法常数大，改为新写法，以及自适应辛普森积分的细节，如删被包含圆、选积分范围。

不会有向面积表示很难受
但是辛普森积分大法好。
不过
为什么我的 $S i m p s o n$ 的常数是别人的 $10$ 倍呢？

首先，所有求图形面积的题，其实本质都可以用 $S i m p s o n$ 来解决。
我们定义 $f (k)$ 表示 $x = k$ 这条直线处于圆形中的长度的总和。
不难发现，如果我们将这个函数进行积分的话，直接就能求出来圆形的面积并了

那么问题就转化成了如果求这个函数，其实不难发现，我们可以通过勾股定理求出来和每个相交的圆的两个交点，然后类似扫描线，按照 $y$ 来排序，每次从下到上，将下面的交点设成 $1$ ，然后上面的弄成 $- 1$ ，那我们只需要计算前缀和大于1的区域即可。
由于我之前的写法常数太大，所以我后来改成了 $s u n c o n g b o$ 大爹的写法。

inline double f(double now)
{
	cnt=0;
	for (register int i=1;i<=n;++i)
	{
		tmp = fabs(a[i].x.x-now);
		if (tmp-a[i].r<=eps)
		{
			d = sqrtl(a[i].r*a[i].r-tmp*tmp);
			++cnt;
			v[cnt]=mk(a[i].x.y-d,a[i].x.y+d);
	    }
    }
    sort(v+1,v+1+cnt);
    double last=0;
	int sum=0;
    double ans=0;	
    int i=1;
    /*for (register int i=1;i<=cnt;++i)
    {
    	if (sum) ans=ans+v[i].first-last;
    	last=v[i].first;
    	sum=sum+v[i].second;
    }*/
    while(i<=cnt)
	{
		double k = v[i].first,mx = v[i].second; ++i;
		while(i<=cnt && dcmp(v[i].first-mx)<=0)
		{
			mx = max(mx,v[i].second);
			++i;
		}
		ans += mx-k;
	}
    return ans;
}

剩下的就是普通的自适应辛普森积分了

不过这里有一些细节
1.要删掉被其他圆完全包含的圆
2.为了避免左右端点对答案的影响，我们选择 $[m n + e p s, m x - e p s]$ 这个范围进行积分。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<vector>
#include<queue>
#include<cmath>
#include<map>
#include<set>
#define pb push_back
#define mk make_pair
#define ll long long
#define lson ch[x][0]
#define rson ch[x][1]
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")
#pragma GCC optimize("inline")
#pragma GCC optimize("-fgcse")
#pragma GCC optimize("-fgcse-lm")
#pragma GCC optimize("-fipa-sra")
#pragma GCC optimize("-ftree-pre")
#pragma GCC optimize("-ftree-vrp")
#pragma GCC optimize("-fpeephole2")
#pragma GCC optimize("-ffast-math")
#pragma GCC optimize("-fsched-spec")
#pragma GCC optimize("unroll-loops")
#pragma GCC optimize("-falign-jumps")
#pragma GCC optimize("-falign-loops")
#pragma GCC optimize("-falign-labels")
#pragma GCC optimize("-fdevirtualize")
#pragma GCC optimize("-fcaller-saves")
#pragma GCC optimize("-fcrossjumping")
#pragma GCC optimize("-fthread-jumps")
#pragma GCC optimize("-funroll-loops")
#pragma GCC optimize("-fwhole-program")
#pragma GCC optimize("-freorder-blocks")
#pragma GCC optimize("-fschedule-insns")
#pragma GCC optimize("inline-functions")
#pragma GCC optimize("-ftree-tail-merge")
#pragma GCC optimize("-fschedule-insns2")
#pragma GCC optimize("-fstrict-aliasing")
#pragma GCC optimize("-fstrict-overflow")
#pragma GCC optimize("-falign-functions")
#pragma GCC optimize("-fcse-skip-blocks")
#pragma GCC optimize("-fcse-follow-jumps")
#pragma GCC optimize("-fsched-interblock")
#pragma GCC optimize("-fpartial-inlining")
#pragma GCC optimize("no-stack-protector")
#pragma GCC optimize("-freorder-functions")
#pragma GCC optimize("-findirect-inlining")
#pragma GCC optimize("-frerun-cse-after-loop")
#pragma GCC optimize("inline-small-functions")
#pragma GCC optimize("-finline-small-functions")
#pragma GCC optimize("-ftree-switch-conversion")
#pragma GCC optimize("-foptimize-sibling-calls")
#pragma GCC optimize("-fexpensive-optimizations")
#pragma GCC optimize("-funsafe-loop-optimizations")
#pragma GCC optimize("inline-functions-called-once")
#pragma GCC optimize("-fdelete-null-pointer-checks")
#include<ctime>

using namespace std;

inline int read()
{
  int x=0,f=1;char ch=getchar();
  while (!isdigit(ch)) {if (ch=='-') f=-1;ch=getchar();}
  while (isdigit(ch)) {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
  return x*f;
}

const int maxn = 4010;
const double eps = 3e-13-5e-14;

int dcmp(double x)
{
	return x<-eps ? -1 : (x>eps ? 1 : 0);
}

struct Point
{
	double x,y;
};

Point operator + (Point a,Point b)
{
	return (Point){a.x+b.x,a.y+b.y};
}

Point operator - (Point a,Point b)
{
	return (Point){a.x-b.x,a.y-b.y};
}

struct R{
	double r;
	Point x;
};
double dis(Point a,Point b)
{
	Point now = a-b;
	return sqrt(now.x*now.x+now.y*now.y);
}

R a[maxn];
double mx=-1e9,mn=1e9;
pair<double,double> v[maxn];
int n;
double tmp;
double d;
int cnt;

inline double f(double now)
{
	cnt=0;
	for (register int i=1;i<=n;++i)
	{
		tmp = fabs(a[i].x.x-now);
		if (tmp-a[i].r<=eps)
		{
			d = sqrtl(a[i].r*a[i].r-tmp*tmp);
			++cnt;
			v[cnt]=mk(a[i].x.y-d,a[i].x.y+d);
	    }
    }
    sort(v+1,v+1+cnt);
    double last=0;
	int sum=0;
    double ans=0;	
    int i=1;
    /*for (register int i=1;i<=cnt;++i)
    {
    	if (sum) ans=ans+v[i].first-last;
    	last=v[i].first;
    	sum=sum+v[i].second;
    }*/
    while(i<=cnt)
	{
		double k = v[i].first,mx = v[i].second; ++i;
		while(i<=cnt && dcmp(v[i].first-mx)<=0)
		{
			mx = max(mx,v[i].second);
			++i;
		}
		ans += mx-k;
	}
    return ans;
}

double t1,t2,t3,tt1,tt2,tt3,mmid;

inline double count(double l,double r)
{
	mmid = (l+r)/2;
  	tt1 = f(l);
  	tt2 = f(r);
  	tt3 = f(mmid);
  	return (tt1+tt2+4*tt3)*(r-l)/6.0;
}

double solve(double l,double r)
{
	//if(dcmp(l-r)>0) return 0.0;
	double mid = (l+r)/2;
	t1 = count(l,r);
	t2 = count(l,mid);
	t3 = count(mid,r);
	if (fabs(t2+t3-t1)<=eps) return t1;
	else return solve(l,mid)+solve(mid+eps,r);
}

double l[maxn],r[maxn];
int m;
R b[maxn];
int tag[maxn];

int main()
{
//int st = clock();
  n=read();
  m=n;
  for (int i=1;i<=n;i++) 
  {	
  	a[i].x.x=read();
  	a[i].x.y=read();
  	a[i].r=read();
    l[i]=a[i].x.x-a[i].r;
    r[i]=a[i].x.x+a[i].r;
	mx=max(a[i].x.x+a[i].r,1.0*mx);
  	mn=min(a[i].x.x-a[i].r,1.0*mn);
  }
  for (int i=1;i<=m;i++)
  {
  	for (int j=1;j<=m;j++)
  	{
  		if (i==j) continue;
  		if (a[i].x.x==a[j].x.x && a[i].x.y==a[j].x.y && a[i].r==a[j].r) continue;
  		if (dcmp(dis(a[i].x,a[j].x) - (a[j].r-a[i].r))<=0) tag[i]=1;
    }
  }
  //cout<<clock()-st<<endl;;
  n=0;
  for (int i=1;i<=m;i++)
  {
  	if (tag[i]) continue;
  	b[++n]=a[i];
  }
  for (int i=1;i<=n;i++) a[i]=b[i];
  printf("%.3lf\n",solve(mn+eps,mx-eps));
  return 0;
}