1. 程式人生 > >bzoj-2178 圓的面積並

bzoj-2178 圓的面積並

題意:

給出平面上的n個圓,求它們的面積並;

n<=1000;

題解:

這題似乎有很多種姿勢來解,我學了一種比較Simple的;

對於三次以下多項式函式的定積分,有一個Simpson公式:

∫[l,r]f(x)=(r-l)(f(l)+f(r)+4f(mid))/6

公式可以利用導數證明,但是對於三次以上或者其他函式是不成立的;

比如圓的引數方程,三角函式之類的奇怪東西;

雖說如此,不成立的時候,也可以利用這個公式來乾點什麼。。

利用這個公式來擬合曲線!

對於一段區間[l,r],我們求得f(l),f(r),f(mid);

然後直接帶入公式,就得到了過三個點的二次函式曲線的面積;

這並不一定是解,但是離解應該也比較接近;

所以驗證一下,再套[l,mid]與[mid,r]的公式,也可以得到一個這個區間的面積;

如果這兩個面積差不多那答案就差不多啦,如果差的很多呢?

遞迴計算!

分別對[l,mid],[mid,r]區間遞迴算值,一直遞迴到可以接受的地步;

這個過程將曲線劃分越來越細,然後對每一段擬合一個函式曲線;

像是一個適應曲線的過程(?),所以這個演算法被稱為Simpson自適應法;

這個演算法顯然不夠完美,比如面對非連續函式幾乎無法得到正確結果;

即使是連續函式也可以構造出很多的神奇資料卡掉演算法(BZOJ這道題並沒有卡就是了);

對策就是在驗證以及遞迴計算的時候隨機分段,或者分成多段;

這樣時間效率會慢一些但是被卡的機率降低了很多;

BZOJ2178有些卡常數,EPS到1e-13也真是喪心病狂;

程式碼:

#include<math.h>
#include<stdio.h>
#include<string.h>
#include<algorithm>
#define N 1100
#define pr pair<double,double>
#define Fabs(x) ((x)>0?(x):-(x))
using namespace std;
const double pi=acos(-1.0);
const double EPS=1e-13;
const double INF=1e100;
struct Point
{
	int x,y;
	friend double dis(Point a,Point b)
	{
		return sqrt((double)(a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
	}
};
struct Circle
{
	Point p;
	int r;
	void read(){scanf("%d%d%d",&p.x,&p.y,&r);}
	friend bool operator <(Circle a,Circle b)
	{
		if(a.p.x-a.r<b.p.x-a.r)
			return a.p.x+a.r<b.p.x+a.r;
		return a.p.x-a.r<b.p.x-a.r;
	}
	pr f(double x)
	{
		if(r<=fabs(p.x-x))	return pr(0,0);
		double t=r*r-(p.x-x)*(p.x-x);
		t=sqrt(t);
		return pr(p.y-t,p.y+t);
	}
}O[N];
bool ban[N];
pr p[N];
int n;
double Cut(double x)
{
	double ret=0,last=-INF;
	int cnt=0;
	for(int i=1;i<=n;i++)
	{
		p[++cnt]=O[i].f(x);
		if(p[cnt]==pr(0,0))
			cnt--;
	}
	sort(p+1,p+cnt+1);
	for(int i=1;i<=cnt;i++)
	{
		if(p[i].first>last)
			ret+=p[i].second-p[i].first,last=p[i].second;
		else if(p[i].second>last)
			ret+=p[i].second-last,last=p[i].second;
	}
	return ret;
}
double Simpson(double l,double r,double mid,double Cl,double Cr,double Cm)
{
	double tCl=Cut((l+mid)/2),tCr=Cut((mid+r)/2);
	double ans=(r-l)*(Cl+Cr+4*Cm)/6,lans=(mid-l)*(Cl+Cm+4*tCl)/6,rans=(r-mid)*(Cr+Cm+4*tCr)/6;
	if(Fabs(lans+rans-ans)<EPS)
		return ans;
	else
		return Simpson(l,mid,(l+mid)/2,Cl,Cm,tCl)+Simpson(mid,r,(mid+r)/2,Cm,Cr,tCr);
}
int main()
{
	int i,j,k;
	double l,r;
	scanf("%d",&n);
	l=INF,r=-INF;
	for(i=1;i<=n;i++)
	{
		O[i].read();
		l=min(l,(double)O[i].p.x-O[i].r);
		r=max(r,(double)O[i].p.x+O[i].r);
	}
	sort(O+1,O+n+1);
	for(i=1;i<=n;i++)
	{
		if(ban[i])	continue;
		for(j=i+1;j<=n;j++)
		{
			if(ban[j])	continue;
			if(dis(O[i].p,O[j].p)+O[j].r<=O[i].r)
				ban[j]=1;
		}
	}
	for(i=1;i<=n;i++)
	{
		if(ban[i])
		{
			swap(ban[i],ban[n]);
			swap(O[i--],O[n--]);
		}
	}
	printf("%.3lf\n",Simpson(l,r,(l+r)/2,0,0,Cut((l+r)/2)));
	return 0;
}