NOI2005 月下柠檬树

题目

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

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

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

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

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

输入描述

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

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

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

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

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

输出描述

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

样例输入

2  0.7853981633 
10.0  10.00  10.00 
4.00  5.00

样例输出

171.97

题解

原来自适应simpson积分是个很简单的东西！
我们尝试分析一下影子，圆的投影还是圆，圆锥的尖投影成一个点，而圆台的棱是圆的公切线，我们把圆心投影出来，发现平面上圆心的距离是两两高度差/tan(alpha)
这是一个轴对称图形，我们只需要求一侧就好了
用simpson积分的公式
\(S = \frac{b - a}{6}(f(a) + 4 * f(\frac{a + b}{2}) + f(b))\)
计算区间就好了，啥，你说肯定不对……
确实肯定不对，然而你可以递归，如何判断这个区间计算的和正确答案相差无几呢，就是左右分别积分出来的值和整个区间的积分相差小于精度的时候就可以认为积分正确了，返回即可，否则左右递归计算积分

代码实现

#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int inf = 0x3f3f3f3f;
const double eps = 1e-6;

const int N = 510;

int n;

double angle;

struct Point {
    double x, y;
}circle[N], line[N << 1];

void get_point() {//求得切线的左右两个坐标。
    for(int i = 1; i < n; i++) {
        double x1 = circle[i].x, x2 = circle[i + 1].x, r1 = circle[i].y, r2 = circle[i + 1].y;
        double ce = r1 * (r1 - r2) / (x2 - x1), xa = x1 + ce, ya = sqrt(r1 * r1 - ce * ce);
        double gd = r2 * (r1 - r2) / (x2 - x1), xb = x2 + gd, yb = sqrt(r2 * r2 - gd * gd);
        line[2 * i - 1] = {xa, ya}, line[2 * i] = {xb, yb};
    }
}

double f(double x) {
    double ans = 0.0;
    for(int i = 1; i <= n; i++) {//在圆里面，
        if(x < circle[i].x + circle[i].y && x > circle[i].x - circle[i].y) {
            ans = max(ans, sqrt(circle[i].y * circle[i].y - (x - circle[i].x) * (x - circle[i].x)));
        }
    }
    for(int i = 1; i <= 2 * (n - 1); i += 2) {//被切线包围，
        if (x >= line[i].x && x <= line[i + 1].x) {
            ans = max(ans, (line[i + 1].y - line[i].y) / (line[i + 1].x - line[i].x) * (x - line[i].x) + line[i].y);
        }
    }
    //这些值一定是取最大嘛，因为它可能同时符合圆，切线的要求，所以我们得取覆盖面积最大的。
    return ans;
}

double sim(double l, double r) {
    return (r - l) * (f(l) + f(r) + f((l + r) / 2.0) * 4.0) / 6.0;
}

double asr(double l, double r, double eps, double ans) {
    double mid = (l + r) / 2.0;
    double ansl = sim(l, mid), ansr = sim(mid, r);
    if(fabs(ansl + ansr - ans) < 15.0 * eps) return ansl + ansr + (ansl + ansr - ans) / 15.0;
    return asr(l, mid, eps / 2.0, ansl) + asr(mid, r, eps / 2.0, ansr);
}

int main() {
    // freopen("in.txt", "r", stdin);
    // freopen("out.txt", "w", stdout);
    // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    scanf("%d %lf", &n, &angle);
    angle = tan(angle);
    n++;
    for(int i = 1; i <= n; i++) {   
        scanf("%lf", &circle[i].x);
        circle[i].x /= angle;
        circle[i].x += circle[i - 1].x;
    }
    for(int i = 1; i < n; i++) {
        scanf("%lf", &circle[i].y);
    }
    circle[n].y = 0;
    get_point();
    double l = circle[1].x - circle[1].y, r = circle[n].x;//确定积分区间。
    for(int i = 1; i <= n; i++) {
        l = min(l, circle[i].x - circle[i].y);
        r = max(r, circle[i].x + circle[i].y);
    }
    printf("%.2f", 2.0 * asr(l, r, eps, sim(l, r)));
	return 0;
}