bzoj 1502 计算几何

计算几何就是麻烦，
我只想学学辛普森函数和自适应辛普森算法，
计算几何其他内容。。。就理性放弃吧。。。。。。

$S(l,r) = (f(l)+ 4* f(mid) + f(r))/6$

---

---

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>

#define pow2(x) ((x)*(x))

const int MAXN = 600 ,INF = 1<<30;const double  eps = 0.0000001;

int n ; double alpha;

double r0[MAXN] = {0};
double h[MAXN] = {0};

struct circle{double x,r;}c[MAXN] = {0};
struct line{double x1,x2,y1,y2;}l[MAXN] = {0};

int dcmp(const double &x)
{
    if(fabs(x) < eps)return 0;
    else if(x > 0)   return 1;
    else            return -1;
}


double f(double x)
{
    double h = 0;

    for(int i = 1 ;i <= n ;i++)
    {
        if(c[i].x - c[i].r <= x && x <= c[i].x + c[i].r)
        {
            h = std::max(h , sqrt(pow2(c[i].r) - pow2(c[i].x - x)));
        }
        if(i < n && l[i].x1 <= x && x <= l[i].x2)
        {
            h = std::max(h , l[i].y1+(x-l[i].x1)*(l[i].y2-l[i].y1)/(l[i].x2-l[i].x1));
        }
    }
    return h;
}

double simpson(const double &fl,const double &fm,const double &fr,const double &l,const double &r){return (r - l)*(fl + 4*fm + fr)/6;}

double fun_simpson(const double fll,const double fm ,const double frr ,const double ll ,const double rr)
{
    double m = (ll + rr)/2 , lm = (ll + m)/2, mr = (m + rr)/2, flm = f(lm),fmr = f(mr); 

    double s = simpson(fll,fm,frr,ll,rr) ;
    double sl = simpson(fll,flm,fm,ll,m) ; 
    double sr = simpson(fm,fmr,frr,m,rr);

    if(dcmp(s - (sl + sr)) != 0)
    {
        sl = fun_simpson(fll,flm,fm,ll,m);
        sr = fun_simpson(fm,fmr,frr,m,rr);

        s = sl + sr;
    }

    return s;
}

int main()
{
    double ans = 0, ll = INF, rr = -INF , mid;

#ifndef ONLINE_JUDGE    
    freopen("bzoj1502.in","r",stdin);
    freopen("bzoj1502.out","w",stdout);
#endif  


    scanf("%d",&n); n = n + 1;
    scanf("%lf",&alpha);


    for(int i = 1 ; i <= n ;i++)
    {
      scanf("%lf", &h[i]);  
      h[i] += h[i-1];
    }
    for(int i = 1 ; i < n ;i++)
      scanf("%lf", &r0[i]);

    for(int i = 1 ; i <= n ;i++)
    {
        c[i].r = r0[i];
        c[i].x = h[i] / tan(alpha);

        ll = std::min(ll ,c[i].x - c[i].r);
        rr = std::max(rr ,c[i].x + c[i].r);
    }

    for(int i = 1 ; i < n ;i++)
    {
        l[i].x1 = c[i].x + c[i].r*(c[i].r-c[i+1].r)/(c[i+1].x-c[i].x);
        l[i].y1 = sqrt(pow2(c[i].r) - pow2(l[i].x1-c[i].x));
        l[i].x2 = c[i+1].x + c[i+1].r*(c[i].r-c[i+1].r)/(c[i+1].x-c[i].x);
        l[i].y2 = sqrt(pow2(c[i+1].r) - pow2(l[i].x2-c[i+1].x));
    }

    mid = (ll + rr)/2;

    ans = fun_simpson(f(ll),f(mid),f(rr),ll,rr)*2;

    printf("%.2lf",ans);

#ifndef ONLINE_JUDGE    
    fclose(stdin);
    fclose(stdout);
#endif  
    return 0;
}