[bzoj1502][NOI2005]月下柠檬树——自适应simpson积分

题目大意：

给定一个由圆台构成的树以及平行光线的夹角，求树的阴影面积。

思路：

一个圆台在平行光下投影必定是两个圆加上两条外公切线，题目所求即n个这种图形的面积并。
每个圆的面积和相邻两个圆的公切线可以利用三角函数来求解，考虑如何求解这个不规则图形的面积。
可以按照圆心连线为x轴建系，可以对于在x轴上的部分求解积分，不难发现这个不规则图形的大部分都是规则图形，于是直接自适应simpson积分即可。
注意几个细节，左右端点的选取要注意，以及尽量避免计算重复的值。

#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a;i<=b;++i)
typedef long long ll;

using namespace std;

void File(){
    freopen("bzoj1502.in","r",stdin);
    freopen("bzoj1502.out","w",stdout);
}

const double pi=acos(-1.0);
const int maxn=500+10;
int n;
double a,x[maxn],r[maxn],a1[maxn],a2[maxn],k[maxn],b[maxn];
bool is[maxn];

double geth(double p,int i){
    if(p<a1[i] || p>a2[i])return 0;
    return k[i]*p+b[i];
}

double fun(double p){
    double ret=0;
    REP(i,1,n+1)if(fabs(x[i]-p)<=r[i])
        ret=max(ret,sqrt(r[i]*r[i]-(x[i]-p)*(x[i]-p)));
    REP(i,1,n)ret=max(ret,geth(p,i));
    return ret;
}

double simpson(double L,double R,double lv,double mv,double rv){
    return (R-L)/6*(lv+4*mv+rv);
}

double solve(double L,double R,double eps,double lv,double mv,double rv){
    double mid=(L+R)/2,ss=fun((L+mid)/2),tt=fun((mid+R)/2);
    double s0=simpson(L,R,lv,mv,rv);
    double s1=simpson(L,mid,lv,ss,mv);
    double s2=simpson(mid,R,mv,tt,rv);
    if(fabs(s1+s2-s0)<=15*eps)return s1+s2+(s1+s2-s0)/15;
    return solve(L,mid,eps,lv,ss,mv)+solve(mid,R,eps,mv,tt,rv);
}

void init(){
    scanf("%d%lf",&n,&a);
    double h=0,t,ta=1.0/tan(a);
    REP(i,1,n+1){
        scanf("%lf",&t);
        h+=t;
        x[i]=ta*h;
    }
    REP(i,1,n)scanf("%lf",&r[i]);
    REP(i,1,n){
        double c=abs(r[i+1]-r[i])/(x[i+1]-x[i]);
        c=pi/2-asin(c);
        double si=sin(c),co=cos(c);
        double b1=si*r[i],b2=si*r[i+1];
        a1[i]=x[i]+(r[i]<r[i+1] ? -1 : 1)*co*r[i];
        a2[i]=x[i+1]+(r[i]<r[i+1] ? -1 : 1)*co*r[i+1];
        k[i]=(b1-b2)/(a1[i]-a2[i]);
        b[i]=b1-k[i]*a1[i];
    }
}

int main(){
    File();
    init();
    double L=999999999,R=-999999999;
    REP(i,1,n+1)L=min(L,x[i]-r[i]),R=max(R,x[i]+r[i]);
    printf("%.2lf\n",solve(L,R,1e-8,fun(L),fun((L+R)/2),fun(R))*2);
    return 0;
}