1502: [NOI2005]月下柠檬树

#include<cmath>
#include<cstdio>
#include<string>
#include<cstring>
#include<algorithm>
using namespace std;
const double eps=1e-5;
struct circle
{
double x,y;
double r;
}a[10001];
struct circle2
{
double x,y;
double r;
double ll,rr;
}sx[10001];
struct line
{
double x1,y1,x2,y2;
}lin[10001];
double h[10001];
bool flag[10001];
int p,n;
inline double fabs(double x)
{
if(x<0)
x=-x;
return x;
}
inline double max(double x,double y)
{
if(x>y)
return x;
return y;
}
inline double geth(double x)
{
int i;
double len=0,dis;
for(i=1;i<=n;i++)
{
if(sx[i].ll>=x||sx[i].rr<=x)
continue;
dis=sqrt(sx[i].r*sx[i].r-(sx[i].x-x)*(sx[i].x-x));
len=max(len,dis*(double)2);
}
for(i=1;i<=p;i++)
{
if(lin[i].x1>x||lin[i].x2<x)
continue;
dis=(lin[i].y1+(lin[i].y2-lin[i].y1)/(lin[i].x2-lin[i].x1)*(x-lin[i].x1));
len=max(len,dis*(double)2);
}
return len;
}
inline double cal(double len,double lh,double mh,double rh)
{
return (lh+mh*(double)4+rh)*len/(double)6;
}
inline double simpson(double l,double mid,double r,double lh,double mh,double rh,double s)
{
double mid1=(l+mid)/(double)2,mid2=(mid+r)/(double)2;
double mh1=geth(mid1),mh2=geth(mid2);
double s1=cal(mid-l,lh,mh1,mh),s2=cal(r-mid,mh,mh2,rh);
if(fabs(s1+s2-s)<eps)
return s1+s2;
return simpson(l,mid1,mid,lh,mh1,mh,s1)+simpson(mid,mid2,r,mh,mh2,rh,s2);
}
int main()
{
memset(flag,false,sizeof(flag));
double alpha;
scanf("%d%lf",&n,&alpha);
int i,j;
for(i=0;i<=n;i++)
{
scanf("%lf",&h[i]);
if(i!=0)
h[i]=h[i]+h[i-1];
a[i+1].x=h[i]/tan(alpha);
a[i+1].y=0;
}
double x;
for(i=1;i<=n;i++)
{
scanf("%lf",&x);
a[i].r=x;
}
for(i=1;i<=n;i++)
{
if(a[i+1].x-a[i].x<fabs(a[i+1].r-a[i].r))
continue;
p++;
lin[p].x1=a[i].x-a[i].r*(a[i+1].r-a[i].r)/(a[i+1].x-a[i].x);
lin[p].y1=a[i].y+sqrt(a[i].r*a[i].r-(a[i].x-lin[p].x1)*(a[i].x-lin[p].x1));
lin[p].x2=a[i+1].x-a[i+1].r*(a[i+1].r-a[i].r)/(a[i+1].x-a[i].x);
lin[p].y2=a[i+1].y+sqrt(a[i+1].r*a[i+1].r-(a[i+1].x-lin[p].x2)*(a[i+1].x-lin[p].x2));
}
for(i=1;i<=n;i++)
{
sx[i].r=a[i].r;
sx[i].x=a[i].x;
sx[i].y=a[i].y;
sx[i].ll=a[i].x-a[i].r;
sx[i].rr=a[i].x+a[i].r;
}
double l=sx[1].x-sx[1].r,r=lin[p].x2,mid,lh,rh,mh;
for(i=1;i<=n;i++)
{
l=min(l,sx[i].ll);
r=max(r,sx[i].rr);
}
mid=(l+r)/(double)2;
double ans=0;
lh=geth(l);
rh=geth(r);
mh=geth(mid);
ans=simpson(l,mid,r,lh,mh,rh,cal(r-l,lh,mh,rh));
printf("%.2lf\n",ans);
return 0;
}