BZOJ1502 simpson积分初步

#include<cstdio>
#include<algorithm>
#include<cmath>
#define inf 1e12
using namespace std;
const int N=1011;
typedef double db;
const db eps=1e-5;
db alpha;
int n,num;
struct point{
    db x,y;
    point(db X=0.00,db Y=0.00){
        x=X,y=Y;
    }
};
struct circle{
    point c;
    db r;
    circle(point C=point(0.00,0.00),db R=0.00){
        c=C,r=R;
    }
}c[N];
struct line{
    point s,t;
    db k,b;
    line(point S=point(0.00,0.00),point T=point(0.00,0.00)){
        s=S,t=T;
        if(s.x>s.x)swap(s,t);
        k=(s.y-t.y)/(s.x-t.x);
        b=s.y-k*s.x;
    }
    inline db f(db x){return k*x+b;}
}l[N];
inline db sqr(db x){
    return x*x;
}
inline int dcmp(db x){
    if(fabs(x)<eps)
        return 0;
    return x<0?-1:1;
}
inline db f(db x){
    db res=0;
    for(register int i=1;i<=n;++i){
        db d=fabs(x-c[i].c.x);
        if(dcmp(d-c[i].r)>0)continue;
        db len=2*sqrt(sqr(c[i].r)-sqr(d));
        res=max(res,len);
    }
    for(register int i=1;i<=num;++i)
        if(l[i].s.x<=x&&x<=l[i].t.x)
            res=max(res,2*l[i].f(x));
    return res;
}
inline db calc(db l,db r){
    return (f(l)+f(r)+4.00*f((l+r)*0.5))*(r-l)/6.00;
}
inline db simpson(db l,db r,db now){
    db mid=(l+r)*0.5;
    db x=calc(l,mid),y=calc(mid,r);
    if(dcmp(now-x-y)==0)return now;
    return simpson(l,mid,x)+simpson(mid,r,y);
}
inline void solve(){
    db L=inf,R=-inf;
    for(register int i=1;i<=n+1;++i)
        L=min(L,c[i].c.x-c[i].r),R=max(R,c[i].c.x+c[i].r);
    for(register int i=1;i<=n;++i){
        db d=c[i+1].c.x-c[i].c.x;
        if(dcmp(d-fabs(c[i].r-c[i+1].r))<0)continue;
        db sina=(c[i].r-c[i+1].r)/d,cosa=sqrt(1-sqr(sina));
        l[++num]=line(point(c[i].c.x+c[i].r*sina,c[i].r*cosa),point(c[i+1].c.x+c[i+1].r*sina,c[i+1].r*cosa));
    }
    printf("%.2lf\n",simpson(L,R,calc(L,R)));
}
int main(){
    scanf("%d%lf",&n,&alpha);
    db h;
    for(register int i=1;i<=n+1;++i){
        scanf("%lf",&h);
        c[i]=circle(point(h/tan(alpha)+c[i-1].c.x,0),0);
    }
    for(register int i=1;i<=n;++i)
        scanf("%lf",&c[i].r);
    solve();
    return 0;
}