[NOI2005]月下柠檬树 luogu P4207 计算几何 自适应辛普森积分法

题目链接
这个题做的有点点小难受… 被初中数学知识给打败了
还有一个小错误调了蛮久的

首先自己观察一下投影 可以拿周围的东西尝试模拟一下

发现圆投影上去的东西是一个和原来全等的圆
然后斜线就是某一层的高度除以倾斜角的正切

借用16bit戦争的一张图

搞出来就是一个这样的东西
圆是和原来全等的 那么我们要考虑的都是切线 还有切点的位置

接下来到了初中数学内容 怎么求两个圆的切线

因为这个题我们建系可以让圆心都在$x$轴上，做起来会比较方便

求直线解析式的话利用相似的只是求出倾斜角正切就很好办了

搞出切线方程并且确定了位置之后

对于这种类似图形的面积并我们当然是采用$Simpson$积分直接做

然后就没了。

Codes

#include<bits/stdc++.h>

using namespace std;

typedef double db;

const int N = 500 + 10;
const db eps = 1e-6, Pi = acos(-1);
db k[N], b[N], L[N], R[N];
db x[N], h[N], r[N], alpha;
int n; 

db pf(db X) {
    return X * X;
}

db f(db X) {
    db y = 0;
    for(int i = 1; i <= n; ++ i) {
        if(L[i] <= X && X <= R[i])
            y = max(y, k[i] * X + b[i]);
        if(x[i] - r[i] <= X && X <= x[i] + r[i])
            y = max(y, sqrt(pf(r[i]) - pf(x[i] - X)));
    }
    return y;
}

db integral(db l, db r) {
    db mid = (l + r) / 2;
    return (r - l) / 6 * (f(l) + f(r) + f(mid) * 4);
}

db Simpson(db l, db r, db Eps) {
    db mid = (l + r) / 2, res = integral(l, r);
    db res1 = integral(l, mid), res2 = integral(mid, r);
    if(fabs(res1 + res2 - res) <= Eps * 15) return res1 + res2 + (res1 + res2 - res) / 15;
    return Simpson(l, mid, Eps / 2) + Simpson(mid, r, Eps / 2);
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("4207.in", "r", stdin);
    freopen("4207.out", "w", stdout);
#endif
    db st = 1e9, ed = -1e9;
    scanf("%d%lf", &n, &alpha);
    alpha = 1 / tan(alpha);
    for(int i = 1; i <= n + 1; ++ i)
        scanf("%lf", &h[i]), h[i] += h[i - 1], x[i] = alpha * h[i];
    for(int i = 1; i <= n; ++ i)
        scanf("%lf", &r[i]);    
    for(int i = 1; i <= n + 1; ++ i)
        st = min(st, x[i] - r[i]), ed = max(ed, x[i] + r[i]);
    k[n] = -r[n] / sqrt(pf(x[n + 1] - x[n]) - pf(r[n])); 
    db K = -1 / k[n], B = x[n] / k[n]; b[n] = -k[n] * x[n + 1];
    L[n] = (B - b[n]) / (k[n] - K); R[n] = x[n + 1];
    for(int i = 1; i < n; ++ i) {
        db K, B1, B2, x_1 = x[i], x_2 = x[i + 1];
        if(r[i + 1] > r[i]) {
            db len = r[i] * (x_2 - x_1) / (r[i + 1] - r[i]);
            k[i] = r[i] / sqrt(pf(len) - pf(r[i]));
            b[i] = -k[i] * (x_1 - len);
        }
        else {
            db len = r[i + 1] * (x_2 - x_1) / (r[i] - r[i + 1]);
            k[i] = -r[i + 1] / sqrt(pf(len) - pf(r[i + 1]));
            b[i] = -k[i] * (x_2 + len);
        }
        K = -1 / k[i], B1 = x_1 / k[i], B2 = x_2 / k[i];
        L[i] = (B1 - b[i]) / (k[i] - K), R[i] = (B2 - b[i]) / (k[i] - K);
    }
//  for(int i = 1; i <= n; ++ i)
//      cout << k[i] << ' ' << b[i] << ' ' << L[i] << ' ' << R[i] << endl;
    printf("%.2lf", 2 * Simpson(st, ed, eps));
    return 0;
}