P4207 [NOI2005]月下柠檬树

题意

算是第一道辛普森积分的题吧，这里有篇洛谷题解相当详细，这里就不赘述了。

code:

#include<bits/stdc++.h>
using namespace std;
const int maxn=1010;
const double eps=1e-7;
int n;
double alpha;
struct Circle{double x,r;}cir[maxn];
struct Lines{double k,b,l,r;}line[maxn];
inline double calc(double a,double c){return sqrt(c*c-a*a);}
inline void work(int x)
{
	if(fabs(cir[x].r-cir[x+1].r)<=eps)
	{
		line[x].l=cir[x].x,line[x].r=cir[x+1].x;
		line[x].k=0,line[x].b=cir[x].r;
		return;
	}
	double dx=cir[x+1].x-cir[x].x,dr=fabs(cir[x].r-cir[x+1].r);
	if(cir[x].r<cir[x+1].r)
	{
		line[x].l=cir[x].x-cir[x].r*dr/dx;
		line[x].r=cir[x+1].x-(cir[x].x-line[x].l)*cir[x+1].r/cir[x].r;
		double ly=calc(cir[x].x-line[x].l,cir[x].r),ry=calc(cir[x+1].x-line[x].r,cir[x+1].r);
		line[x].k=(ry-ly)/(line[x].r-line[x].l);
		line[x].b=ly-line[x].l*line[x].k;
	}
	else 
	{
		line[x].l=cir[x].x+cir[x].r*dr/dx;
		line[x].r=cir[x+1].x+(line[x].l-cir[x].x)*cir[x+1].r/cir[x].r;
		double ly=calc(line[x].l-cir[x].x,cir[x].r),ry=calc(line[x].r-cir[x+1].x,cir[x+1].r);
		line[x].k=(ry-ly)/(line[x].r-line[x].l);
		line[x].b=ly-line[x].l*line[x].k;
	}
}
inline double F(double x)
{
	double res=0;
	for(int i=1;i<=n;i++)
		if(x>=cir[i].x-cir[i].r&&x<=cir[i].x+cir[i].r)
			res=max(res,calc(abs(cir[i].x-x),cir[i].r));
	for(int i=1;i<n;i++)
		if(x>=line[i].l&&x<=line[i].r)
			res=max(res,line[i].k*x+line[i].b);
	return res;
}
inline double simpson(double l,double r)
{
	double mid=(l+r)/2.0;
	return (F(l)+F(r)+4*F(mid))*(r-l)/6.0;
}
double solve(double l,double r,double ans,double eps)
{
	double mid=(l+r)/2.0,lans,rans;
	lans=simpson(l,mid),rans=simpson(mid,r);
	if(fabs(lans+rans-ans)<=15.0*eps)return lans+rans+(lans+rans-ans)/15.0;
	return solve(l,mid,lans,eps/2)+solve(mid,r,rans,eps/2);
}
int main()
{
	scanf("%d%lf",&n,&alpha);
	n++;
	for(int i=1;i<=n;i++)
	{
		scanf("%lf",&cir[i].x);
		cir[i].x*=1.0/tan(alpha);
		cir[i].x+=cir[i-1].x;
	}
	for(int i=1;i<n;i++)scanf("%lf",&cir[i].r);
	cir[n].r=0.0;
	for(int i=1;i<n;i++)work(i);
	double l=1e9,r=-1e9;
	for(int i=1;i<=n;i++)l=min(l,cir[i].x-cir[i].r),r=max(r,cir[i].x+cir[i].r);
	printf("%.2lf",2.0*solve(l,r,calc(l,r),1e-3));
	return 0;
}