bzoj1502 simpson求面积
题目：http://www.lydsy.com/JudgeOnline/problem.php?id=1502
题解：
simpson积分求面积，s = (f(a)+f(b)+4*f(c))/6*Δx ，c=(a+b)/2。
题中的树投影下来是一些圆和相邻圆的公切线组成的一个封闭图形，并且上下对称，所以可以只求上半部分。
simpson求面积时，若f(x)代价很大，要尽量减少其的重复调用。其次尽量优化f(x)函数的计算。
写完后还要自己出极限随机数据，将eps调到能接受的最大值。
1 #include <cstdio> 2 #include <cmath> 3 #include <iostream> 4 #define maxn 510 5 #define eps 1e-6 6 using namespace std; 7 8 int sg( double x ) { 9 return (x>-eps)-(x<eps); 10 } 11 struct Vector { 12 double x, y; 13 Vector(){} 14 Vector( double x, double y ) : x(x), y(y) {} 15 Vector operator+( Vector b ) { return Vector(x+b.x,y+b.y); } 16 Vector operator-( Vector b ) { return Vector(x-b.x,y-b.y); } 17 Vector operator*( double b ) { return Vector(x*b,y*b); } 18 double len() { 19 return sqrt( x*x+y*y ); 20 } 21 }; 22 typedef Vector Point; 23 24 struct Line { 25 Point A, B; 26 double k, b; 27 double lf, rg; 28 Line(){} 29 Line( Point A, Point B ):A(A),B(B) { 30 if( A.x>B.x ) swap(A,B); 31 lf = A.x; 32 rg = B.x; 33 k = (B.y-A.y)/(B.x-A.x); 34 b = A.y-A.x*k; 35 } 36 inline bool in( double x ) { 37 return sg(x-lf)>=0 && sg(rg-x)>=0; 38 } 39 inline double f( double x ) { return k*x+b; } 40 void out() { 41 printf( "(%lf,%lf) (%lf,%lf) \n", A.x, A.y, B.x, B.y ); 42 } 43 }; 44 struct Circle { 45 Point C; 46 double r; 47 double lf, rg; 48 Circle(){} 49 void make() { 50 lf = C.x-r; 51 rg = C.x+r; 52 } 53 inline bool in( double x ) { 54 return sg(x-lf)>=0 && sg(rg-x)>=0; 55 } 56 inline double f( double x ) { 57 return sqrt( r*r-(C.x-x)*(C.x-x) ); 58 } 59 void out() { 60 printf( "(%lf,%lf) r=%lf\n", C.x, C.y, r ); 61 } 62 }; 63 64 int n; 65 double alpha; 66 Line lines[maxn]; 67 Circle cirs[maxn]; 68 69 Line cir_tang( Circle & a, Circle & b ) { 70 double ang = acos( (a.r-b.r)/(b.C.x-a.C.x) ); 71 Point A = a.C + Vector(cos(ang),sin(ang))*a.r; 72 Point B = b.C + Vector(cos(ang),sin(ang))*b.r; 73 return Line(A,B); 74 } 75 76 double F( double x ) { 77 double y = -1.0; 78 for( int i=0; i<=n; i++ ) { 79 if( cirs[i].lf-eps<=x && x<=cirs[i].rg+eps ) 80 y = max( y, cirs[i].f(x) ); 81 if( i&& lines[i].lf-eps<=x && x<=lines[i].rg+eps ) 82 y = max( y, lines[i].f(x) ); 83 } 84 return y; 85 } 86 inline double simpson( double a, double b, double fa, double fb, double fc ) { 87 return (fa+4*fc+fb)*(b-a)/6; 88 } 89 double adapt( double a, double b, double c, double fa, double fb, double fc ) { 90 double d = (a+c)/2, e = (c+b)/2; 91 double fd = F(d), fe = F(e); 92 double sa = simpson(a,c,fa,fc,fd); 93 double sb = simpson(c,b,fc,fb,fe); 94 double ss = simpson(a,b,fa,fb,fc); 95 if( fabs(sa+sb-ss)<15*eps ) return sa+sb; 96 return adapt(a,c,d,fa,fc,fd)+adapt(c,b,e,fc,fb,fe); 97 } 98 99 int main() { 100 scanf( "%d%lf", &n, &alpha ); 101 double k = 1/tan(alpha); 102 double hsum = 0.0, h; 103 double lf=1e200, rg=-1e200; 104 for( int i=0; i<=n; i++ ) { 105 scanf( "%lf", &h ); 106 hsum += h; 107 cirs[i].C.x = hsum*k; 108 cirs[i].C.y = 0.0; 109 } 110 for( int i=0; i<n; i++ ) { 111 scanf( "%lf", &cirs[i].r ); 112 cirs[i].make(); 113 lf = min( lf, cirs[i].lf ); 114 rg = max( rg, cirs[i].rg ); 115 } 116 cirs[n].r = 0.0; 117 rg = max( rg, cirs[n].C.x ); 118 for( int i=1; i<=n; i++ ) 119 lines[i] = cir_tang( cirs[i-1], cirs[i] ); 120 printf( "%.2lf\n", adapt(lf,rg,(lf+rg)/2,F(lf),F(rg),F((lf+rg)/2) )*2 ); 121 }
