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调到能接受的最大值。

#include <cstdio>
#include <cmath>
#include <iostream>
#define maxn 510
#define eps 1e-6
using namespace std;

int sg( double x ) {
    return (x>-eps)-(x<eps);
}
struct Vector {
    double x, y;
    Vector(){}
    Vector( double x, double y ) : x(x), y(y) {}
    Vector operator+( Vector b ) { return Vector(x+b.x,y+b.y); }
    Vector operator-( Vector b ) { return Vector(x-b.x,y-b.y); }
    Vector operator*( double b ) { return Vector(x*b,y*b); }
    double len() {
        return sqrt( x*x+y*y );
    }
};
typedef Vector Point;

struct Line {
    Point A, B;
    double k, b;
    double lf, rg;
    Line(){}
    Line( Point A, Point B ):A(A),B(B) {
        if( A.x>B.x ) swap(A,B);
        lf = A.x;
        rg = B.x;
        k = (B.y-A.y)/(B.x-A.x);
        b = A.y-A.x*k;
    }
    inline bool in( double x ) {
        return sg(x-lf)>=0 && sg(rg-x)>=0;
    }
    inline double f( double x ) { return k*x+b; }
    void out() {
        printf( "(%lf,%lf) (%lf,%lf) 
", A.x, A.y, B.x, B.y );
    }
};
struct Circle {
    Point C;
    double r;
    double lf, rg;
    Circle(){}
    void make() {
        lf = C.x-r;
        rg = C.x+r;
    }
    inline bool in( double x ) {
        return sg(x-lf)>=0 && sg(rg-x)>=0;
    }
    inline double f( double x ) {
        return sqrt( r*r-(C.x-x)*(C.x-x) );
    }
    void out() {
        printf( "(%lf,%lf) r=%lf
", C.x, C.y, r );
    }
};

int n; 
double alpha;
Line lines[maxn];
Circle cirs[maxn];

Line cir_tang( Circle & a, Circle & b ) {
    double ang = acos( (a.r-b.r)/(b.C.x-a.C.x) );
    Point A = a.C + Vector(cos(ang),sin(ang))*a.r;
    Point B = b.C + Vector(cos(ang),sin(ang))*b.r;
    return Line(A,B);
}

double F( double x ) {
    double y = -1.0;
    for( int i=0; i<=n; i++ ) {
        if( cirs[i].lf-eps<=x && x<=cirs[i].rg+eps ) 
            y = max( y, cirs[i].f(x) );
        if( i&& lines[i].lf-eps<=x && x<=lines[i].rg+eps ) 
            y = max( y, lines[i].f(x) );
    }
    return y;
}
inline double simpson( double a, double b, double fa, double fb, double fc ) {
    return (fa+4*fc+fb)*(b-a)/6;
}
double adapt( double a, double b, double c, double fa, double fb, double fc ) {
    double d = (a+c)/2, e = (c+b)/2;
    double fd = F(d), fe = F(e);
    double sa = simpson(a,c,fa,fc,fd);
    double sb = simpson(c,b,fc,fb,fe);
    double ss = simpson(a,b,fa,fb,fc);
    if( fabs(sa+sb-ss)<15*eps ) return sa+sb;
    return adapt(a,c,d,fa,fc,fd)+adapt(c,b,e,fc,fb,fe);
}

int main() {
    scanf( "%d%lf", &n, &alpha );
    double k = 1/tan(alpha);
    double hsum = 0.0, h;
    double lf=1e200, rg=-1e200;
    for( int i=0; i<=n; i++ ) {
        scanf( "%lf", &h );
        hsum += h;
        cirs[i].C.x = hsum*k;
        cirs[i].C.y = 0.0;
    }
    for( int i=0; i<n; i++ )  {
        scanf( "%lf", &cirs[i].r );
        cirs[i].make();
        lf = min( lf, cirs[i].lf );
        rg = max( rg, cirs[i].rg );
    }
    cirs[n].r = 0.0;
    rg = max( rg, cirs[n].C.x );
    for( int i=1; i<=n; i++ ) 
        lines[i] = cir_tang( cirs[i-1], cirs[i] );
    printf( "%.2lf
", adapt(lf,rg,(lf+rg)/2,F(lf),F(rg),F((lf+rg)/2) )*2 );
}