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 }
View Code

 

posted @ 2015-02-05 17:40  idy002  阅读(272)  评论(0)  编辑  收藏  举报