[NOI2005]月下柠檬树[计算几何（simpson）]

#include<cmath>
#include<cstdio>
#include<algorithm>
#define pf(x) ((x)*(x))
using namespace std;
const int N=2000+5;
const double eps=1e-6;
typedef pair<double,double> point;
typedef pair<double,double> circle;
struct line{
    point s,t;
    double k,b;
    line(){}
    line(point _s,point _t){
        s=_s;t=_t;
        k=(s.second-t.second)/(s.first-t.first);
        b=s.second-s.first*k;
    }
    const double f(const double x){
        return k*x+b;
    }
};
int n,n1;double alpha,H[N];
point p;line L[N];circle C[N];
double lb=2e9,rb;
double sina,cosa,tana;
inline void add(const circle &a,const circle &b){
    n1++;
    sina=(a.second-b.second)/(b.first-a.first);
    cosa=sqrt(1-pf(sina));
    tana=sina/cosa;
    L[n1].s=make_pair(a.first+a.second*sina,a.second*cosa);
    L[n1].t=make_pair(b.first+b.second*sina,b.second*cosa);
    L[n1].k=-tana;
    L[n1].b=L[n1].s.second-L[n1].s.first*L[n1].k;
}
inline const double F(const double x){
    double re=0;
    for(int i=1;i<=n1;i++) if(x>=L[i].s.first&&x<=L[i].t.first) re=max(re,L[i].f(x));
    for(int i=1;i<=n;i++) if(x>=C[i].first-C[i].second&&x<=C[i].first+C[i].second) re=max(re,sqrt(pf(C[i].second)-pf(x-C[i].first)));
    return re;
}
inline const double simpson(const double l,const double r){
    double mid=(l+r)/2;
    return (F(l)+F(r)+4*F(mid))*(r-l)/6;
}
inline double asr(double l,double r,double eps,double last){
    double mid=(l+r)/2;
    double L=simpson(l,mid),R=simpson(mid,r);
    if(fabs(L+R-last)<=15*eps) return L+R+(L+R-last)/15;
    return asr(l,mid,eps/2,L)+asr(mid,r,eps/2,R);
}
inline int cmp(const double x){
    if(fabs(x)<eps) return 0;
    return x>0?1:-1;
}
int main(){
    scanf("%d%lf",&n,&alpha);
    for(int i=1;i<=n+1;i++) scanf("%lf",&H[i]),H[i]+=H[i-1];
    for(int i=1;i<=n;i++) scanf("%lf",&C[i].second);
    double ta=tan(alpha);
    p=make_pair(H[n+1]/ta,0);rb=max(rb,p.first);
    double x,r,l,h;
    C[n].first=H[n]/ta;
    x=C[n].first;r=C[n].second;
    lb=min(lb,x-r);
    rb=max(rb,x+r);
    if(x+r<p.first){
        l=pf(r)/(p.first-x);// 射影定理 
        h=sqrt(pf(r)-pf(l));
        L[++n1]=line(make_pair(x+l,h),p);
    }
    for(int i=n-1;i;i--){
        C[i].first=H[i]/ta;
        x=C[i].first;r=C[i].second;
        lb=min(lb,x-r);
        rb=max(rb,x+r);
        if(cmp(C[i+1].first-x-fabs(C[i+1].second-r))>0)//内含 
            add(C[i],C[i+1]);
    }
    printf("%.2lf\n",2*asr(lb,rb,eps,simpson(lb,rb)));
    return 0;
}