Simpson法相当好用啊！神奇的骗分算法！
1 /************************************************************** 2 Problem: 1502 3 User: zhuohan123 4 Language: C++ 5 Result: Accepted 6 Time:228 ms 7 Memory:1312 kb 8 ****************************************************************/ 9 10 #include <iostream> 11 #include <cstdio> 12 #include <cstring> 13 #include <cmath> 14 #include <algorithm> 15 using namespace std; 16 const double eps=1e-8,pi=3.141592653589793238; 17 int n; 18 struct point 19 { 20 double x,y; 21 point(){} 22 point(double X,double Y){x=X,y=Y;} 23 }; 24 struct line 25 { 26 point s,e; 27 line(){} 28 line(point S,point E){s=S,e=E;} 29 double y(double x){return (e.y*s.x-e.x*s.y-e.y*x+s.y*x)/(s.x-e.x);} 30 }l[510];int lnum; 31 struct circle{double x,r;}c[510]; 32 double h[510]; 33 double f(double x) 34 { 35 double s=0; 36 for(int i=1;i<n;i++) 37 { 38 if(abs(x-c[i].x)+eps<c[i].r)s=max(s,sqrt(c[i].r*c[i].r-(x-c[i].x)*(x-c[i].x))); 39 if(l[i].s.x<x+eps&&l[i].e.x>x-eps)s=max(s,l[i].y(x)); 40 } 41 return s*2; 42 } 43 const double dev=1e-6; 44 inline double simpson(double l,double r,double fl,double fm,double fr){return (fl+4*fm+fr)/6*(r-l);} 45 double integral(double l,double fl,double m,double fm,double r,double fr,double pre) 46 { 47 double lm=(l+m)/2,rm=(m+r)/2,flm=f(lm),frm=f(rm); 48 double intl=simpson(l,m,fl,flm,fm),intr=simpson(m,r,fm,frm,fr); 49 return abs(intl+intr-pre)<dev?intl+intr:integral(l,fl,lm,flm,m,fm,intl)+integral(m,fm,rm,frm,r,fr,intr); 50 } 51 int main(int argc, char *argv[]) 52 { 53 double alp;scanf("%d%lf",&n,&alp);n++; 54 for(int i=1;i<=n;i++)scanf("%lf",&h[i]); 55 for(int i=1;i<n;i++)scanf("%lf",&c[i].r);c[n].r=0; 56 double s=1e10,e=-1e10; 57 for(int i=1;i<=n;i++) 58 { 59 h[i]+=h[i-1]; 60 c[i].x=h[i]/tan(alp); 61 s=min(s,c[i].x-c[i].r); 62 e=max(e,c[i].x+c[i].r); 63 } 64 for(int i=1;i<n;i++) 65 { 66 double dx=c[i+1].x-c[i].x,dr=c[i].r-c[i+1].r; 67 if(abs(dx)<abs(dr)+eps)continue ; 68 l[++lnum]=line(point(c[i].x+c[i].r/dx*dr,sqrt(c[i].r*c[i].r-(c[i].r/dx*dr)*(c[i].r/dx*dr))) 69 ,point(c[i+1].x+c[i+1].r/dx*dr,sqrt(c[i+1].r*c[i+1].r-(c[i+1].r/dx*dr)*(c[i+1].r/dx*dr)))); 70 } 71 double m=(s+e)/2,fm=f(m); 72 printf("%.2lf ",integral(s,0,m,fm,e,0,simpson(s,e,0,fm,0))); 73 return 0; 74 }