bzoj 1502 月下柠檬树【Simpson积分】

投影到地面之后，会发现圆形在平行光下面积和形状是不会变的，也就是所要求的图形是若干个圆和把相邻两个圆连起来的公切线所组成的。
公切线和圆间距瞎求一下就行，注意要去掉被完全覆盖的圆
然后simpson即可
eps大概1e-6

#include<iostream> 
#include<cstdio> 
#include<cmath> 
using namespace std; 
const int N=1005;
const double eps=1e-6,inf=1e15;
double alp; 
int n,m,num; 
struct dian 
{ 
    double x,y; 
    dian (double X=0,double Y=0) 
	{
		x=X; y=Y;
	} 
}; 
struct yuan 
{ 
    double r; 
    dian c; 
    yuan(dian a=(dian){0,0},double R=0) 
	{
		c=a; r=R;
	} 
}a[N]; 
struct xian 
{ 
    dian s,t; 
    double k,b; 
    xian(dian S=(dian){0,0},dian T=(dian){0,0})  
    { 
        s=S,t=T; 
        if(s.x>t.x) swap(s,t); 
        k=(s.y-t.y)/(s.x-t.x); 
        b=s.y-k*s.x; 
    } 
    double f(double x) 
	{
		return k*x+b;
	} 
}l[N]; 
int cmp(double x) 
{
	if(fabs(x)<eps) 
		return 0; 
	return x<0? -1:1;
} 
double f(double x) 
{ 
    double re=0; 
    for(int i=1;i<=n;i++)
        { 
            double d=fabs(x-a[i].c.x);
            if(cmp(d-a[i].r)>0) 
				continue;
            double len=2*sqrt(a[i].r*a[i].r-d*d);
            re=max(re,len); 
        } 
    for(int i=1;i<=num;i++)
        if(x>=l[i].s.x && x<=l[i].t.x) 
			re=max(re,2*l[i].f(x));
    return re;
}
double sps(double l,double r,double now,double fl,double fr,double fm)
{//cout<<l<<" "<<r<<endl;
	double mid=(l+r)/2,ffl=f((l+mid)/2),ffr=f((mid+r)/2),p=(fl+fm+ffl*4)*(mid-l)/6,q=(fm+fr+ffr*4)*(r-mid)/6;
	if(cmp(now-p-q)==0)
		return now;
	else
		return sps(l,mid,p,fl,fm,ffl)+sps(mid,r,q,fm,fr,ffr);
}
int main() 
{ 
    scanf("%d%lf",&n,&alp);
    double h,r;
    for(int i=1;i<=n+1;i++) 
	{
        scanf("%lf",&h), 
        a[i]=(yuan){((dian){(h/tan(alp))+a[i-1].c.x,0}),0};
	}
    for(int i=1;i<=n;i++) 
        scanf("%lf",&r),a[i].r=r;
    double L=inf,R=-inf;
    for(int i=1;i<=n+1;i++) 
        L=min(L,a[i].c.x-a[i].r),R=max(R,a[i].c.x+a[i].r);
    for(int i=1;i<=n;i++) 
        { 
            double d=a[i+1].c.x-a[i].c.x; 
            if(cmp(d-fabs(a[i].r-a[i+1].r))<0) continue;
            double sina=(a[i].r-a[i+1].r)/d,cosa=sqrt(1-sina*sina);
            l[++num]=(xian){(dian){a[i].c.x+a[i].r*sina,a[i].r*cosa},(dian){a[i+1].c.x+a[i+1].r*sina,a[i+1].r*cosa}};
        } 
    // printf("%.2lf\n",Simpson(L,R,Calc(L,R)));
	double fl=f(L),fr=f(R),fm=f((L+R)/2);
	printf("%.2lf\n",sps(L,R,(fl+4*fm+fr)*(R-L)/6,fl,fr,fm));
    return 0;
}
/*
2 0.72953
9.61090 0.26021 4.47090
2.98979 2.00036
*/