[NOI2005] bzoj 1502 月下柠檬树 - 辛普森积分

贴板子，注意边界细节

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<cmath>
#define gc getchar()
#define N 510
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
#define db double
using namespace std;
inline int inn()
{
    int x,ch;while((ch=gc)<'0'||ch>'9');
    x=ch^'0';while((ch=gc)>='0'&&ch<='9')
        x=(x<<1)+(x<<3)+(ch^'0');return x;
}
const db eps=1e-6;int n,lc;db h[N];
inline db gabs(db x) { return x<0?-x:x; }
struct circle{
    db x,r;circle(db _x=0.0,db _r=0.0):x(_x),r(_r){}
    inline db f(db x0) { return 2*sqrt(max(0.0,r*r-(x-x0)*(x-x0))); }
}c[N];
struct line{
    db x1,y1,x2,y2;
    line(db _x1=0,db _y1=0,db _x2=0,db _y2=0)
    {   x1=_x1,y1=_y1,x2=_x2,y2=_y2;    }
    inline db f(db x0)
    {
        if(x0<x1||x0>x2) return 0.0;
        return 2*(y1+(x0-x1)/(x2-x1)*(y2-y1));
    }
}l[N];
inline db F(db x0,db ans=0.0)
{
    for(int i=1;i<=n;i++) ans=max(ans,c[i].f(x0));
    for(int i=1;i<=lc;i++) ans=max(ans,l[i].f(x0));
    return ans;
}
inline db simpson(db a,db b) { return (F(a)+F(b)+4*F((a+b)/2))/6*(b-a); }
inline db S(db a,db b,db A,db Eps=eps)
{
    db c=(a+b)/2,L=simpson(a,c),R=simpson(c,b);
    if(gabs(L+R-A)<15*Eps) return L+R+(L+R-A)/15;
    return S(a,c,L,Eps/2)+S(c,b,R,Eps/2);
}
inline db S(db a,db b) { return S(a,b,simpson(a,b)); }
int main()
{
    db alp,dr,dx,len,R=0,L=100000000.0;
    scanf("%d%lf",&n,&alp);
    for(int i=1;i<=n+1;i++) scanf("%lf",&h[i]);
    for(int i=1;i<=n;i++) scanf("%lf",&c[i].r);
    for(int i=1;i<=n+1;i++) h[i]+=h[i-1];
    for(int i=1;i<=n+1;i++) c[i].x=h[i]/tan(alp);
    for(int i=1;i<=n;i++)
    {
        dr=c[i+1].r-c[i].r,dx=c[i+1].x-c[i].x,len=sqrt(dx*dx-dr*dr);
        if(dx<=gabs(c[i].r-c[i+1].r)) continue;else lc++;
        l[lc].x1=c[i].x-dr/dx*c[i].r,l[lc].x2=c[i+1].x-dr/dx*c[i+1].r,
        l[lc].y1=len/dx*c[i].r,l[lc].y2=len/dx*c[i+1].r;
//      debug(l[lc].x1)sp,debug(l[lc].x2)sp,debug(l[lc].y1)sp,debug(l[lc].y2)ln;
    }
    for(int i=1;i<=n+1;i++) L=min(L,c[i].x-c[i].r);
    for(int i=1;i<=n+1;i++) R=max(R,c[i].x+c[i].r);
    return !printf("%.2f\n",S(L,R));
}