BZOJ1502: [NOI2005]月下柠檬树

因为光线平行，射过去的面积其实是一个个圆和他们之间的公切线
算面积可以套自适应Simpson积分，好像大概是这个东西$S(l,r)=\dfrac{r-l}{6}(h(l)+h(r)+4h(mid))$

code：

#include<set>
#include<map>
#include<deque>
#include<queue>
#include<stack>
#include<cmath>
#include<ctime>
#include<bitset>
#include<string>
#include<vector>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<climits>
#include<complex>
#include<iostream>
#include<algorithm>
#define ll long long
using namespace std;

const int maxn = 510;
const double eps = 1e-7;
inline double sqr(const double x){return x*x;}

int n; double alp;
struct cir{double xi,ri;}r[maxn]; int rn;
inline bool cover(cir x,cir y)
{
    if(x.ri<y.ri) swap(x,y);
    double dis=fabs(x.xi-y.xi);
    return dis<x.ri-y.ri+eps;
}
struct Seg{int i,j;double alpha,l,r,k,b;}seg[maxn]; int sn;

double f(double x)
{
    double ans=0;
    for(int i=1;i<=rn;i++)
        if(fabs(x-r[i].xi)<r[i].ri+eps) ans=max(ans,sqrt(sqr(r[i].ri)-sqr(x-r[i].xi)));
    for(int i=1;i<=sn;i++)
    {
        int L=seg[i].i,R=seg[i].j;
        double q=seg[i].alpha;
        if(r[L].ri>=r[R].ri)
        {
            double t1=r[L].xi+r[L].ri*sin(q),t2=r[R].xi+r[R].ri*sin(q);
            if(x<=t1||x>=t2) continue;
            double k=(x-t1)/(t2-t1);
            ans=max(ans,r[L].ri*cos(q)-(r[L].ri-r[R].ri)*cos(q)*k);
        }
        else
        {
            double t1=r[L].xi-r[L].ri*sin(q),t2=r[R].xi-r[R].ri*sin(q);
            if(x<=t1||x>=t2) continue;
            double k=(t2-x)/(t2-t1);
            ans=max(ans,r[R].ri*cos(q)-(r[R].ri-r[L].ri)*cos(q)*k);
        }
        //if(seg[i].l<=x&&x<=seg[i].r) ans=max(ans,seg[i].k*x+seg[i].b);
    }
    //printf("%lf %lf\n",x,ans);
    return ans*2.0;
}
double cal(double L,double R)
{
    double fl=f(L),fr=f(R),fm=f((L+R)/2.0);
    return (R-L)/6.0*(fl+fr+fm*4.0);
}
double solve(double L,double R,double now)
{
    double mid=(L+R)/2.0;
    double t1=cal(L,mid),t2=cal(mid,R);
    if(fabs(t1+t2-now)<eps) return now;
    return solve(L,mid,t1)+solve(mid,R,t2);
}

int main()
{
    scanf("%d%lf",&n,&alp); alp=1.0/tan(alp);
    double las=0;
    for(int i=0;i<=n;i++)
    {
        double x; scanf("%lf",&x); x=x*alp;
        las+=x; r[++rn].xi=las;
    }
    for(int i=1;i<rn;i++) scanf("%lf",&r[i].ri);
    r[rn].ri=eps;

    //for(int i=1;i<=rn;i++) printf("%lf %lf\n",r[i].xi,r[i].ri);
    for(int i=1;i<rn;i++)
    {
        if(cover(r[i],r[i+1])) continue;
        ++sn; double q=asin(fabs(r[i+1].ri-r[i].ri)/(r[i+1].xi-r[i].xi));
        seg[sn].i=i,seg[sn].j=i+1,seg[sn].alpha=q;
        /*seg[sn].k=tan(q);
        seg[sn].l=r[i].xi-r[i].ri*sin(q),seg[sn].r=r[i+1].xi-r[i+1].ri*sin(q);
        double h=r[i].ri*cos(q);
        seg[sn].b=h-seg[sn].k*seg[sn].l;
        printf("%lf %lf %lf %lf\n",seg[sn].l,seg[sn].r,seg[sn].k,seg[sn].b);*/
    }
    double ml=0,mr=0;
    for(int i=1;i<=rn;i++) mr=max(mr,r[i].xi+r[i].ri),ml=min(ml,r[i].xi-r[i].ri);
    printf("%.2lf\n",solve(ml,mr,0));

    return 0;
}