bzoj-2178 圓的面積並
阿新 • • 發佈:2019-01-23
題意:
給出平面上的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; }