BZOJ1502: [NOI2005]月下柠檬树

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const
int MAXN=550;
const
double eps=1e-6;
int n;
struct Point
{double x,y;}P[MAXN];

struct Circle
{Point o;double r;}C[MAXN];

struct Line
{
Point x,y;
double k,b;
inline void Bg()
{
k=(y.y-x.y)/(y.x-x.x);
b=y.y-k*y.x;
}
inline double D(double x)
{return k*x+b;}

}L[MAXN];
int Con=0;
#define abs(a) ((a)<0?(-(a)):(a))
inline double D(double x)
{
double d,ans=0;
for(int i=1;i<=n;i++)
{
d=abs(x-C[i].o.x);
if(d-C[i].r>-eps)
continue;
ans=max(ans,2*sqrt(C[i].r*C[i].r-d*d));
}
for(int i=1;i<=Con;i++)
if(x>=L[i].x.x&&x<=L[i].y.x)
ans=max(ans,2*L[i].D(x));
return ans;
}
#define Calc(a,b,c,d) (((b)+4*(c)+(d))*(a)/6)
double Simpson(double L,double M,double R,double fL,double fM,double fR,double sqr)
{
double Mid1=(L+M)/2,Mid2=(M+R)/2,D1=D(Mid1),D2=D(Mid2),Area1=Calc(M-L,fL,D1,fM),Area2=Calc(R-M,fM,D2,fR);
if(abs(Area1+Area2-sqr)<eps) return Area1+Area2;
return Simpson(L,Mid1,M,fL,D1,fM,Area1)+Simpson(M,Mid2,R,fM,D2,fR,Area2);
}

int main()
{
double Min=1<<29,Max=-Min,d,h,Ang,Sin,Cos;
scanf("%d%lf",&n,&Ang);
n++;
for(int i=1;i<=n;i++)
scanf("%lf",&h),C[i].o.y=0,C[i].o.x=h/tan(Ang)+C[i-1].o.x;
for(int i=1;i<n;i++)scanf("%lf",&C[i].r);
for(int i=1;i<=n;i++)
Min=min(Min,C[i].o.x-C[i].r),Max=max(Max,C[i].o.x+C[i].r);
for(int i=1;i<n;i++)
{
d=C[i+1].o.x-C[i].o.x;
if(d-abs(C[i+1].r-C[i].r)<-eps)continue;
Sin=(C[i].r-C[i+1].r)/d,Cos=sqrt(1-Sin*Sin);
L[++Con]=(Line){(Point){C[i].o.x+C[i].r*Sin,C[i].r*Cos},(Point){C[i+1].o.x+C[i+1].r*Sin,C[i+1].r*Cos}};
L[Con].Bg();
}
double L=Min,R=Max;
double M=(L+R)/2;
double fL=D(L),fM=D(M),fR=D(R);
printf("%.2lf\n",Simpson(L,M,R,fL,fM,fR,Calc(R-L,fL,fM,fR)));
return 0;
}