BZOJ1502: [NOI2005]月下柠檬树

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

额，本蒟蒻的第一道自适应辛普森积分的题。

我觉得这位大佬写的还蛮好的：博客

然后也第一次知道辛普森复杂度为 O(玄学) 。。。

还有，几何真有趣用！

附代码：

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cmath>
#define MAXN 510
#define eps (1e-7)
using namespace std;
int n;
double alpha,h[MAXN],r[MAXN];
struct line{
    double k,b,left,right;
}a[MAXN];
inline int read(){
    int date=0,w=1;char c=0;
    while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}
    while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}
    return date*w;
}
inline double get_length(double z,double x){
    return sqrt(z*z-x*x);
}
void get_line(int x,int y){
    if(fabs(r[x]-r[y])<eps){
        a[x].left=h[x];a[x].right=h[y];
        a[x].k=0;a[x].b=r[x];
        return;
    }
    double deltax=h[y]-h[x],deltay=fabs(r[x]-r[y]),left,right;
    if(r[x]>r[y]){
        a[x].left=h[x]+r[x]*deltay/deltax;
        a[x].right=h[y]+(a[x].left-h[x])*r[y]/r[x];
        left=get_length(r[x],a[x].left-h[x]);
        right=get_length(r[y],a[x].right-h[y]);
        a[x].k=(left-right)/(a[x].left-a[x].right);
        a[x].b=left-a[x].left*a[x].k;
    }
    else{
        a[x].right=h[y]-r[y]*deltay/deltax;
        a[x].left=h[x]-(h[y]-a[x].right)*r[x]/r[y];
        left=get_length(r[x],a[x].left-h[x]);
        right=get_length(r[y],a[x].right-h[y]);
        a[x].k=(left-right)/(a[x].left-a[x].right);
        a[x].b=left-a[x].left*a[x].k;
    }
}
inline double f(double x){
    double ans=0;
    for(int i=1;i<=n;i++){
        if(h[i]-r[i]<x&&x<h[i]+r[i])ans=max(ans,get_length(r[i],x-h[i]));
        if(a[i].left<=x&&x<=a[i].right)ans=max(ans,a[i].k*x+a[i].b);
    }
    return ans;
}
inline double simpson(double l,double r){
    double mid=(l+r)/2.0;
    return (f(l)+f(mid)*4.0+f(r))*(r-l)/6.0;
}
double query(double l,double r,double ans){
    double mid=(l+r)/2.0,left=simpson(l,mid),right=simpson(mid,r);
    if(fabs(left+right-ans)<=eps*15.0)return left+right+(left+right-ans)/15.0;
    return query(l,mid,left)+query(mid,r,right);
}
void work(){
    double left=h[1]-r[1],right=h[n];
    for(int i=1;i<=n;i++){
        left=min(left,h[i]-r[i]);
        right=max(right,h[i]+r[i]);
    }
    printf("%.2lf\n",query(left,right,simpson(left,right))*2.0);
}
void init(){
    n=read();scanf("%lf",&alpha);
    alpha=1.00/tan(alpha);
    scanf("%lf",&h[1]);
    h[1]*=alpha;
    for(int i=2;i<=n+1;i++){
        scanf("%lf",&h[i]);
        h[i]*=alpha;
        h[i]+=h[i-1];
    }
    for(int i=1;i<=n;i++)scanf("%lf",&r[i]);
    n++;
    r[n]=0.0;
    for(int i=1;i<n;i++)get_line(i,i+1);
}
int main(){
    init();
    work();
    return 0;
}