[BZOJ1502]月下柠檬树(自适应辛普森积分)

1502: [NOI2005]月下柠檬树

Time Limit: 5 Sec  Memory Limit: 64 MB
Submit: 1387  Solved: 739
[Submit][Status][Discuss]

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

HINT

Source

[Submit][Status][Discuss]

这个题的正解是各种分类讨论求面积，但是可以用辛普森积分骗分。
关于辛普森积分的推导可以看 https://www.zhihu.com/question/47728235
其实这个东西很好卡掉，但是对于很多数据还是可以较为精确地出解的。

首先发现要求的东西是轴对称的，这意味着我们可以只求x坐标上面的函数区域没。这个图像的解析式还是很难求，但是对于每个横坐标上对应的函数值可以比较轻松地求出，这个时候就可以用Simpson积分了。
普通辛普森积分根本不可能出解，因为仅仅一个二次函数不可能拟合这么复杂的图像。但是要分成很多份的话时间复杂度又过高，这个时候就需要引入自适应：如果要对[a,b]积分，可以先对[a,b]拟合，再对[a,(a+b)/2]和[(a+b)/2,b]拟合，如果两个相差很小则直接退出，否则继续递归。
个人认为这个方法在这里可以成功很大程度是因为图像中有很多直线围城的区域，而自适应Simpson积分在直线积分方面的表现应该是很不错的。

具体做法看 https://www.cnblogs.com/DaD3zZ-Beyonder/p/5676841.html
上面讲的很清楚了，直接作为模板即可。

#include<cmath>
#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=l; i<=r; i++)
typedef double db;
using namespace std;

const int N=1010;
double eps=1e-5,inf=1e12;
int n,num;
db alpha;
struct P{ db x,y; P(db X=0,db Y=0){ x=X; y=Y; } };
struct Ci{ db r; P c; Ci (P C=(P(0,0)),db R=0):c(C),r(R){}; }C[N];

struct Li{
    P s,t; db k,b;
    Li (P S=P(0,0),P T=P(0,0)){
        s=S,t=T; if (s.x>t.x) swap(s,t);
        k=(s.y-t.y)/(s.x-t.x); b=s.y-k*s.x;
    }
    db f(db x){ return k*x+b; }
}l[N];

int dcmp(db x){ if (fabs(x)<eps) return 0; return (x<0) ? -1 : 1; } 
db F(db x){
    db res=0;
    rep(i,1,n){
        db d=fabs(x-C[i].c.x);
        if (dcmp(d-C[i].r)>0) continue;
        db len=2*sqrt(C[i].r*C[i].r-d*d);
        res=max(res,len);
    }
    rep(i,1,num) if (x>l[i].s.x && x<=l[i].t.x) res=max(res,2*l[i].f(x));
    return res;
}

db calc(db l,db r){ db mid=(l+r)/2; return (F(l)+F(r)+F(mid)*4)*(r-l)/6; }
db Simpson(db l,db r,db now){
    db mid=(l+r)/2,x=calc(l,mid),y=calc(mid,r);
    if (!dcmp(now-x-y)) return now; else return Simpson(l,mid,x)+Simpson(mid,r,y);
}

void solve(){
    db L=inf,R=-inf;
    rep(i,1,n+1) L=min(L,C[i].c.x-C[i].r),R=max(R,C[i].c.x+C[i].r);
    rep(i,1,n){
        db d=C[i+1].c.x-C[i].c.x;
        if (dcmp(d-fabs(C[i].r-C[i+1].r))<0) continue;
        db sina=(C[i].r-C[i+1].r)/d,cosa=sqrt(1-sina*sina);
        l[++num]=Li(P(C[i].c.x+C[i].r*sina,C[i].r*cosa),P(C[i+1].c.x+C[i+1].r*sina,C[i+1].r*cosa));
    }
    printf("%.2lf
",Simpson(L,R,calc(L,R)));
}

int main(){
    scanf("%d%lf",&n,&alpha); db h,r;
    rep(i,1,n+1)
        scanf("%lf",&h),C[i]=Ci((P((h/tan(alpha))+C[i-1].c.x,0)),0);
    rep(i,1,n) scanf("%lf",&r),C[i].r=r;
    solve();
    return 0;
}