NOI2005 月下柠檬树

题目:

  李哲非常非常喜欢柠檬树,特别是在静静的夜晚,当天空中有一弯明月温柔地照亮地面上的景物时,他必会悠闲地坐在他亲手植下的那棵柠檬树旁,独自思 索着人生的哲理。
  李哲是一个喜爱思考的孩子,当他看到在月光的照射下柠檬树投在地面上的 影子是如此的清晰,马上想到了一个问题:树影的面积是多大呢?
  李哲知道,直接测量面积是很难的,他想用几何的方法算,因为他对这棵柠 檬树的形状了解得非常清楚,而且想好了简化的方法。
  李哲将整棵柠檬树分成了 n 层,由下向上依次将层编号为 1,2,…,n。从第 1 到 n-1 层,每层都是一个圆台型,第 n 层(最上面一层)是圆锥型。对于圆台型, 其上下底面都是水平的圆。对于相邻的两个圆台,上层的下底面和下层的上底面 重合。第 n 层(最上面一层)圆锥的底面就是第 n-1 层圆台的上底面。所有的底面 的圆心(包括树顶)处在同一条与地面垂直的直线上。李哲知道每一层的高度为 h1,h2,…,hn,第 1 层圆台的下底面距地面的高度为 h0,以及每层的下底面的圆的 半径 r1,r2,…,rn。李哲用熟知的方法测出了月亮的光线与地面的夹角为 alpha。
  为了便于计算,假设月亮的光线是平行光,且地面是水平的,在计算时忽略 树干所产生的影子。李哲当然会算了,但是他希望你也来练练手。

思路:

  这题是在模拟赛里做到的,之前完全没听过自适应辛普森积分
  关于自适应辛普森积分:
  先给出辛普森公式:
  推导可见这篇文章
  当然这个式子的精度是不太靠谱的,因为 越小精度越高,所以如果一个区间 的精度不令人满意的话,就可以把 拆成 继续计算。伪代码如下:

def f(x) = ...
def simpson(l, r) = (r - l) * (f(l) + f(r) + 4 * f((l + r) / 2)) / 6
def calc(l, r) :
    m = (l + r) / 2
    if(abs(simpson(l, m) + simpson(m, r) - simpson(l, r)) < eps)
        return simpson(l, r)
    else
        return calc(l, m) + calc(m, r)

  所以现在的问题就是如果求出 ,也就是某一条直线在影子上的截距。
  画图发现,影子是由每个圆盘投影出的圆和两圆之间的公切线组成的封闭图形。所以可以之间用数学方法求出每个圆的坐标和半径、每条公切线的解析式和区间。查询时暴力枚举每一个即可,复杂度 。本地试了一下 是能跑出来的。

代码

#include <bits/stdc++.h>
#define P2(x) ((x) * (x))
#define rep(i, a, b) for(int i(a), i##_END_(b); i <= i##_END_; i++)
#define drep(i, a, b) for(int i(a), i##_END_(b); i >= i##_END_; i--)
#define File(_) freopen(#_ ".in", "r", stdin), freopen(#_ ".out", "w", stdout)
using namespace std;
template<class T> inline bool tomax(T &a, T b) {return a < b ? a = b, 1 : 0;}
template<class T> inline bool tomin(T &a, T b) {return b < a ? a = b, 1 : 0;}
typedef long long ll;
typedef double db;
const db PI = acos(-1), EPS = 1e-8;
const int N = 505;

db h[N], r[N];

struct Point {
    db x, y;
};
struct Line {
    db k, b, l, r;
    db calc(db x) {
        return k * x + b;
    }
    bool in(db x) {
        return l <= x && x <= r;
    }
};
struct Circle {
    db p, r;
    db calc(db x) {
        return sqrt(P2(r) - P2(x - p));
    }
    bool in(db x) {
        return fabs(x - p) <= r;
    }
    bool in(Circle a) {
        return a.r + fabs(a.p - p) <= r;
    }
};

Line getLine(Point a, Point b) {
    db k = (a.y - b.y) / (a.x - b.x);
    db t = a.y - a.x * k;
    return (Line) {k, t, a.x, b.x};
}
Line getLine(Circle a, Circle b) {
    db d = b.p - a.p;
    db o = asin((b.r - a.r) / d);
    Point p1 = (Point) {a.p - a.r * sin(o), a.r * cos(o)};
    Point p2 = (Point) {b.p - b.r * sin(o), b.r * cos(o)};
    return getLine(p1, p2);
}

Circle cr[N];
Line ln[N];
int cntLn, n;

db f(db x) {
    db ans = 0;
    rep(i, 1, n) 
        if(cr[i].in(x)) tomax(ans, cr[i].calc(x));
    rep(i, 1, cntLn)
        if(ln[i].in(x)) tomax(ans, ln[i].calc(x));
    return ans * 2;
}
db simpson(db a, db b) {
    db c = (a + b) / 2;
    return (b - a) * (f(a) + f(b) + f(c) * 4) / 6;
}
db asr(db a, db b, db ans) {
    db c = (a + b) / 2;
    db l = simpson(a, c), r = simpson(c, b);
    if(fabs(l + r - ans) < EPS) 
        return l + r;
    return asr(a, c, l) + asr(c, b, r);
}

int main() {
    File(lemon);
    db cot;
    scanf("%d%lf", &n, &cot);
    cot = 1 / tan(cot);
    rep(i, 1, n + 1) {
        scanf("%lf", h + i);
        h[i] += h[i - 1];
    }
    rep(i, 1, n) scanf("%lf", r + i);
    rep(i, 1, n + 1) cr[i] = (Circle) {h[i] * cot, r[i]};
    rep(i, 1, n) 
        if(!cr[i].in(cr[i + 1]) || !cr[i + 1].in(cr[i]))
            ln[++cntLn] = getLine(cr[i], cr[i + 1]);
    db mx = -1e9, mn = 1e9;
    rep(i, 1, n + 1) {
        tomax(mx, h[i] * cot + r[i]);
        tomin(mn, h[i] * cot - r[i]);
    }
    printf("%.2lf\n", asr(mn, mx, simpson(mn, mx)));
    return 0;
}