洛谷P4207 [NOI2005] 月下柠檬树

题目描述

李哲非常非常喜欢柠檬树，特别是在静静的夜晚，当天空中有一弯明月温柔地照亮地面上的景物时，他必会悠闲地坐在他亲手植下的那棵柠檬树旁，独自思索着人生的哲理。

李哲是一个喜爱思考的孩子，当他看到在月光的照射下柠檬树投在地面上的影子是如此的清晰，马上想到了一个问题：树影的面积是多大呢？

李哲知道，直接测量面积是很难的，他想用几何的方法算，因为他对这棵柠檬树的形状了解得非常清楚，而且想好了简化的方法。

李哲将整棵柠檬树分成了 𝑛n 层，由下向上依次将层编号为 1,2,...,𝑛1,2,...,n。从第 11 到 𝑛−1n−1 层，每层都是一个圆台型，第 𝑛n 层(最上面一层)是圆锥型。对于圆台型， 其上下底面都是水平的圆。对于相邻的两个圆台，上层的下底面和下层的上底面重合。第 𝑛n 层(最上面一层)圆锥的底面就是第 𝑛−1n−1 层圆台的上底面。所有的底面 的圆心(包括树顶)处在同一条与地面垂直的直线上。李哲知道每一层的高度为 ℎ1,ℎ2,...,ℎ𝑛h1​,h2​,...,hn​，第 11 层圆台的下底面距地面的高度为 ℎ0h0​，以及每层的下底面的圆的半径 𝑟1,𝑟2,...,𝑟𝑛r1​,r2​,...,rn​。李哲用熟知的方法测出了月亮的光线与地面的夹角为 alphaalpha。

为了便于计算，假设月亮的光线是平行光，且地面是水平的，在计算时忽略树干所产生的影子。李哲当然会算了，但是他希望你也来练练手。

输入格式

第一行包含一个整数 𝑛n 和一个实数 alphaalpha，表示柠檬树的层数和月亮 的光线与地面夹角（单位为弧度）。

第二行包含 𝑛+1n+1 个实数 ℎ0,ℎ1,ℎ2,...,ℎ𝑛h0​,h1​,h2​,...,hn​ 表示树离地的高度和每层的高度。

第三行包含 𝑛n 个实数 𝑟1,𝑟2,...,𝑟𝑛r1​,r2​,...,rn​ 表示柠檬树每层下底面的圆的半径。

同一行相邻的两个数之间用一个空格分隔。

输入的所有实数的小数点后可能包含一至十位有效数字。

输出格式

输出一个实数，表示树影的面积，四舍五入保留两位小数。

输入输出样例

输入 #1

2  0.7853981633 
10.0  10.00  10.00 
4.00  5.00

输出 #1

171.97

说明/提示

对于 10%10% 的数据，𝑛≤1n≤1；

对于 30%30% 的数据，𝑛≤2n≤2；

对于 60%60% 的数据，𝑛≤20n≤20；

对于 100%100% 的数据，1≤𝑛≤5001≤n≤500，0.3<alpha<𝜋20.3<alpha<2π​，0<ℎ𝑖≤1000<hi​≤100，0<𝑟𝑖≤1000<ri​≤100。

Code：

#include<cstdio>
#include<cmath>
#include<algorithm>
const double PI=acos(-1);
#define EPS 1e-10
struct Vector{
	double x,y;
	inline Vector(){x=y=0;}
	inline Vector(double _x,double _y){x=_x;y=_y;}
	inline double len()const{return __builtin_sqrt(x*x+y*y);}
	inline double len2()const{return x*x+y*y;}
	inline double atan()const{return __builtin_atan2(y,x);}
	inline void operator += (const Vector &a){x+=a.x;y+=a.y;}
	inline void operator -= (const Vector &a){x-=a.x;y-=a.y;}
	inline void operator *= (const double &a){x*=a;y*=a;}
	inline void operator /= (const double &a){x/=a;y/=a;}
	inline Vector operator + (const Vector &a)const{return Vector(x+a.x,y+a.y);}
	inline Vector operator - (const Vector &a)const{return Vector(x-a.x,y-a.y);}
	inline Vector operator * (const double &a)const{return Vector(x*a,y*a);}
	inline Vector operator / (const double &a)const{return Vector(x/a,y/a);}
	inline double operator * (const Vector &a)const{return x*a.x+y*a.y;}
	inline double operator ^ (const Vector &a)const{return x*a.y-y*a.x;}
};
struct Circle{
	Vector o;
	double r;
};
struct Line{
	Vector way,o;
	inline Line(){}
	inline Line(double a,double b,double c){//ax+by+c=0
		if(std::abs(a)>=EPS) o=Vector(-c/a,0);
		else o=Vector(0,-c/b);
		way=Vector(b,-a);
	}
	inline Line(const Vector &a,const Vector &b){o=a;way=b-a;}
};
inline Vector rotate(const Vector &v,double a){//counterclockwise, will not change the len
	double x=v.x,y=v.y,cos=__builtin_cos(a),sin=__builtin_sin(a);
	return Vector(x*cos-y*sin,x*sin+y*cos);
}
inline double getAngle(const Vector &a,const Vector &b){
	return b.atan()-a.atan();
}
inline Vector getIntersection(const Vector &v,const Circle &c,int type){
	double a=v.atan();
	return Vector(c.o.x+__builtin_cos(a+PI/2)*c.r*type,c.o.y+__builtin_sin(a+PI/2)*c.r*type);
}
inline Line getCommonTangent(const Circle &a,const Circle &b,Vector *pa=NULL,Vector *pb=NULL,int type=1){
//type: 1= on left of the vector (b-a) -1= on right
	Vector d=b.o-a.o;
	double ang=__builtin_acos((a.r-b.r)/d.len());
	d=rotate(d,type*(ang-PI/2));
	Vector _pa=getIntersection(d,a,type),_pb=getIntersection(d,b,type);
	if(pa) *pa=_pa;
	if(pb) *pb=_pb;
	return Line(_pa,_pb);
}
inline double getYByX(const Circle &c,double x){
	if(x>c.o.x+c.r) return 0;
	return __builtin_sqrt(c.r*c.r-(x-c.o.x)*(x-c.o.x));
}
inline double getYByX(const Line &l,double x){
	return l.o.y+l.way.y/l.way.x*(x-l.o.x);
}
#define N 100006
int n;
double alpha;
Circle c[N];
Line tangent[N];
Vector va[N],vb[N];
double h[N];
inline void init(){
	n++;
	double cot=1/__builtin_tan(alpha);
	c[1].o=Vector(0,0);
	for(int i=2;i<=n;i++) c[i].o=c[i-1].o+Vector(h[i-1]*cot,0);
	for(int i=1;i<n;i++) tangent[i]=getCommonTangent(c[i],c[i+1],vb+i,va+i+1);
}
inline double calc(double x){
	double ans=0;
	for(int i=1;i<=n;i++) ans=std::max(ans,getYByX(c[i],x));
	for(int i=1;i<n;i++)if(vb[i].x<=x&&x<=va[i+1].x) ans=std::max(ans,getYByX(tangent[i],x));
//		printf("x= %lf   ans= %lf\n",x,ans);
	return ans;
}
inline double simpson(double l,double r){
	return (r-l)*(calc(l)+4*calc((l+r)/2)+calc(r))/6;
}
double asr(double l,double r,double eps,double ans){
	double mid=(l+r)/2;
	double fl=simpson(l,mid),fr=simpson(mid,r);
	if(std::abs(fl+fr-ans)<=15*eps) return fl+fr+(fl+fr-ans)/15;
	return asr(l,mid,eps/2,fl)+asr(mid,r,eps/2,fr);
}
int main(){
	scanf("%d%lf",&n,&alpha);
	for(int i=0;i<=n;i++) scanf("%lf",h+i);
	for(int i=1;i<=n;i++) scanf("%lf",&c[i].r);
	init();
		#ifdef LOCAL
		puts("\n\n");
		for(int i=1;i<=n;i++) printf("%lf %lf r=%lf    %lf %lf  %lf %lf\n",c[i].o.x,c[i].o.y,c[i].r,va[i].x,va[i].y,vb[i].x,vb[i].y);
		puts("\n\n");
		#endif
	double l=1e10,r=-1e10;
	for(int i=1;i<=n;i++) l=std::min(l,c[i].o.x-c[i].r),r=std::max(r,c[i].o.x+c[i].r);
		#ifdef LOCAL
		printf("l= %lf   r= %lf\n\n\n",l,r);
		#endif
	printf("%.2lf\n",asr(l,r,1e-3,simpson(l,r))*2);
	return 0;
}