[BZOJ1502] [NOI2005]月下柠檬树

Description

李哲非常非常喜欢柠檬树，特别是在静静的夜晚，当天空中有一弯明月温柔地照亮地面上的景物时，他必会悠闲地

坐在他亲手植下的那棵柠檬树旁，独自思索着人生的哲理。李哲是一个喜爱思考的孩子，当他看到在月光的照射下

柠檬树投在地面上的影子是如此的清晰，马上想到了一个问题：树影的面积是多大呢？李哲知道，直接测量面积是

很难的，他想用几何的方法算，因为他对这棵柠檬树的形状了解得非常清楚，而且想好了简化的方法。李哲将整棵

柠檬树分成了n 层，由下向上依次将层编号为1,2,…,n。从第1到n-1 层，每层都是一个圆台型，第n 层(最上面一

层)是圆锥型。对于圆台型，其上下底面都是水平的圆。对于相邻的两个圆台，上层的下底面和下层的上底面重合

。第n 层(最上面一层)圆锥的底面就是第n-1 层圆台的上底面。所有的底面的圆心(包括树顶)处在同一条与地面垂

直的直线上。李哲知道每一层的高度为h1,h2,…,hn，第1 层圆台的下底面距地面的高度为h0，以及每层的下底面

的圆的半径r1,r2,…,rn。李哲用熟知的方法测出了月亮的光线与地面的夹角为alpha。

为了便于计算，假设月亮的光线是平行光，且地面是水平的，在计算时忽略树干所产生的影子。

李哲当然会算了，但是他希望你也来练练手

Input

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

第2行包含n+1个实数h0,h1,h2,…,hn，表示树离地的高度和每层的高度。

第3行包含n个实数r1,r2,…,rn，表示柠檬树每层下底面的圆的半径。

上述输入文件中的数据，同一行相邻的两个数之间用一个空格分隔。

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

1≤n≤500，0.3<alpha<π/2，0<hi≤100，0<ri≤100

Output

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

Sample Input

2 0.7853981633
10.0 10.00 10.00
4.00 5.00

Sample Output

171.97

Solution

注意到一个圆锥投影下来是一个圆加一个三角形，所以可以知道这个图形的样子。

然后辛普森积分就好了。

为啥我bzoj死活卡不过啊...

打表大法好

#include<bits/stdc++.h>
using namespace std;
 
void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}
 
void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');}

#define lf double

const int maxn = 5e2+10;
const lf eps = 5e-6;

int n;
lf h[maxn],rr[maxn],alpha,cot;
struct data {
	lf l,r,k,b,L,R;
}a[maxn];

lf sqr(lf x) {return x*x;}

lf f(lf x) {
	lf ans=0;
	for(int i=1;i<=n;i++) {
		lf l=a[i].L,r=a[i].R,t=0;
		if(l<=x&&x<=r) t=a[i].k*x+a[i].b;
		else if(x>=a[i].l-rr[i]) t=sqrt(sqr(rr[i])-sqr(a[i].l-x));
		else if(x<=a[i].r+rr[i+1]) t=sqrt(sqr(rr[i+1])-sqr(x-a[i].r));
		ans=t>ans?t:ans;
	}
	return ans;
}

lf _int(lf l,lf r) {return (r-l)*(f(l)+f(r)+4.0*f((l+r)*0.5))/6.0;}

lf simpson(lf l,lf r,lf res,lf fl,lf fr) {
	lf mid=(l+r)*0.5,fmid=f(mid);
	lf L=(mid-l)*(fl+fmid+4.0*f((l+mid)*0.5))/6.0;
	lf R=(r-mid)*(fmid+fr+4.0*f((mid+r)*0.5))/6.0;
	if(fabs(L+R-res)<=eps) return L+R;
	else return simpson(l,mid,L,fl,fmid)+simpson(mid,r,R,fmid,fr);
}

int main() {
	read(n);scanf("%lf",&alpha);
	cot=1.0/tan(alpha);
	for(int i=0;i<=n;i++) scanf("%lf",&h[i]),h[i]+=i?h[i-1]:0;
	for(int i=1;i<=n;i++) scanf("%lf",&rr[i]);
	lf l=1e9,r=-1e9;
	for(int i=1;i<=n;i++) {
		a[i].l=h[i-1]*cot,a[i].r=h[i]*cot;
		lf t;
		if(fabs(rr[i]-rr[i+1])<eps) {
			a[i].L=a[i].l,a[i].R=a[i].r,a[i].k=0,a[i].b=rr[i];
		} else {
			if(rr[i]<rr[i+1]) {
				lf x=rr[i]*(a[i].r-a[i].l)/(rr[i+1]-rr[i])+a[i].r-a[i].l;
				t=rr[i+1]/x;
				a[i].k=rr[i+1]/(x*sqrt(1.00-sqr(rr[i+1]/x)));
				a[i].b=-a[i].k*(a[i].r-x);
			} else {
				lf x=rr[i+1]*(a[i].r-a[i].l)/(rr[i]-rr[i+1])+a[i].r-a[i].l;
				t=rr[i]/x;
				a[i].k=-rr[i]/(x*sqrt(1.00-sqr(rr[i]/x)));
				a[i].b=-a[i].k*(a[i].l+x);
			}
			t=sqrt(1-sqr(t));
			a[i].L=(rr[i]*t-a[i].b)/a[i].k;
			a[i].R=(rr[i+1]*t-a[i].b)/a[i].k;
		}
		l=min(l,a[i].l-rr[i]);
		r=max(r,a[i].r+rr[i+1]);
	}
	printf("%.2lf\n",simpson(l,r,_int(l,r),f(l),f(r))*2);
	return 0;
}