luogu P4207 [NOI2005]月下柠檬树

luogu

为了方便,设树顶的那一点为一个半径为0的圆

首先我们把所有圆投影到平面上,然后得到的图形应该是一堆圆心在(x)轴上的圆+相邻两个圆公切线构成的梯形的并,我们可以只计算上面一半的面积然后乘2.因为上面一半的轮廓线不规则,考虑自适应Simpson求轮廓线和(x)轴围成的面积.现在的问题是求出某一个横坐标上与(x)轴垂直的直线和轮廓线的交点纵坐标,因为这个轮廓线是圆和公切线的并,那么我们只要枚举所有圆和公切线,然后对于这些图形在(x)处的纵坐标取最大值即可

至于圆的公切线,这可以把两个圆的半径同时减去较小的半径,然后就变成圆和点的切线,可以用一些(sincos)等算出来,最后平移会回原位置即可

#include<bits/stdc++.h>
#define LL long long
#define db double

using namespace std;
const int N=500+10;
const db pi=acos(-1),eps=1e-10;
int rd()
{
    int x=0,w=1;char ch=0;
    while(ch<'0'||ch>'9'){if(ch=='-') w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+(ch^48);ch=getchar();}
    return x*w;
}
int n;
db gg,h[N],r[N],li[N][2],ps[N][2];
db sq(db x){return x*x;}
db f(db x)
{
    db an=0;
    for(int i=1;i<n;++i)
	if(x>ps[i][0]-eps&&x<ps[i][1]+eps) an=max(an,li[i][0]*x+li[i][1]);
    for(int i=1;i<=n;++i)
    {
	    if(fabs(x-h[i])>r[i]-eps) continue;
    	an=max(an,sqrt(sq(r[i])-sq(fabs(x-h[i]))));
    }
    return an;
}
db g(db l,db r){return (f(l)+4*f((l+r)/2)+f(r))/6*(r-l);}
db cal(db l,db r,db s)
{
    db mid=(l+r)/2,s1=g(l,mid),s2=g(mid,r);
    if(fabs(s-s1-s2)<eps) return s;
    return cal(l,mid,s1)+cal(mid,r,s2);
}

int main()
{
//be careful of details
    n=rd()+1;
    scanf("%lf",&gg),gg=1/tan(gg);
    for(int i=1;i<=n;++i) scanf("%lf",&h[i]),h[i]=h[i]*gg+h[i-1];
    for(int i=1;i<n;++i) scanf("%lf",&r[i]);
    for(int i=1;i<n;++i)
    {
    	if(fabs(r[i]-r[i+1])<eps)
    	{
    	    li[i][0]=0,li[i][1]=r[i],ps[i][0]=h[i],ps[i][1]=h[i+1];
    	    continue;
    	}
    	db nx=0,ny=abs(r[i]-r[i+1]),xx=nx,yy=ny;
    	db a1=acos(ny/(h[i+1]-h[i])),ag=pi/2-a1,nk,nb;
    	if(r[i]>r[i+1])
    	{
    	    nx=xx*cos(ag)+yy*sin(ag),ny=yy*cos(ag)-xx*sin(ag);
    	    nk=-ny/(h[i+1]-h[i]-nx),nb=-nk*h[i+1]+r[i+1]*(sin(a1)-nk*sin(ag));
    	}
    	else
    	{
    	    nx=xx*cos(ag)-yy*sin(ag),ny=yy*cos(ag)+xx*sin(ag);
    	    nk=ny/(h[i+1]+nx-h[i]),nb=-nk*h[i]+r[i]*(sin(a1)+nk*sin(ag));
    	}
    	db dx=r[i]*sin(ag)*(r[i]<r[i+1]?-1:1),d2=r[i+1]*sin(ag)*(r[i]<r[i+1]?-1:1);
    	li[i][0]=nk,li[i][1]=nb;
    	if(fabs(min(r[i],r[i+1]))>eps) ps[i][0]=h[i]+dx,ps[i][1]=h[i+1]+d2;
    	else if(fabs(r[i])<eps) ps[i][0]=h[i],ps[i][1]=h[i+1]+nx;
    	else ps[i][0]=h[i]+nx,ps[i][1]=h[i+1];
    }
    db ql=1e9,qr=-1e9;
    for(int i=1;i<=n;++i) ql=min(ql,h[i]-r[i]),qr=max(qr,h[i]+r[i]);
    printf("%.2lf
",cal(ql,qr,g(ql,qr))*2);
    return 0;
}