BZOJ 1502 月下柠檬树 自适应辛普森积分

题意：给你一颗由n个圆台组成的柠檬树，再给你一个光线，问你这棵树的投影的面积。

思路：由于光线是平行的，所以圆的投影还是原来的圆，样例的投影如下图：

其中DE，CF为原来圆台的母线，现在为两个圆的切线。最后无脑辛普森就好。。。具体的就是对于一个下标x，暴力找到x这个位置的所有圆和切线的最大值。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cstdlib>
using namespace std;

const int maxn = 500 + 10;
struct Line{
    double x1, y1;
    double x2, y2;
    double k;
    double b;
    Line() {}
    Line(double x1, double y1, double x2, double y2, double k, double b) :
        x1(x1), y1(y1), x2(x2), y2(y2), k(k), b(b){}
};
int n;
double rad;
int line_cnt;
Line line[maxn];
double h[maxn], r[maxn];

double F(double x){
    double ans = 0;
    for(int i = 0; i <= n; i++){
        if(x >= h[i] - r[i] && x <= h[i] + r[i]){
            ans = max(ans, sqrt(r[i] * r[i] - fabs(h[i] - x) * fabs(h[i] - x)));
        }
    }
    for(int i = 0; i < line_cnt; i++){
        if(x >= line[i].x1 && x <= line[i].x2){
            ans = max(ans, line[i].k * x + line[i].b);
        }
    }
    return ans;
}

double simpson(double a, double b){
    double c = a + (b - a) / 2;
    return (F(a) + 4 * F(c) + F(b)) * (b - a) / 6;
}

double asr(double a, double b, double eps, double A){
    double c = a + (b - a) / 2;
    double L = simpson(a, c), R = simpson(c, b);
    if(fabs(L + R - A) <= eps) return L + R + (L + R - A) / 15;
    return asr(a, c, eps / 2, L) + asr(c, b, eps / 2, R);
}

double asr(double a, double b, double eps){
    return asr(a, b, eps, simpson(a, b));
}

void solve(){
    scanf("%d%lf", &n, &rad);
    double cur_h, sum_h = 0;
    for(int i = 0; i <= n; i++){
        scanf("%lf", &cur_h);
        h[i] = (cur_h + sum_h) / tan(rad);
        sum_h += cur_h;
    }
    double l = 1e9, r1 = 0;
    for(int i = 0; i <= n; i++){
        if(i != n) scanf("%lf", &r[i]);
        else r[i] = 0;
        l = min(l, h[i] - r[i]);
        r1 = max(r1, h[i] + r[i]);
    }
    line_cnt = 0;
    for(int i = 0; i < n; i++){
        double sina, cosa;
        if(fabs(r[i] - r[i+1]) >= fabs(h[i] - h[i+1])) continue;
        sina = (r[i] - r[i+1]) / (h[i+1] - h[i]);
        cosa = sqrt(1 - sina * sina);
        double x1 = h[i] + r[i] * sina;
        double y1 = r[i] * cosa;
        double x2 = h[i+1] + r[i+1] * sina;
        double y2 = r[i+1] * cosa;
        double k = (y1 - y2) / (x1 - x2);
        double b = y1 - k * x1;
        line[line_cnt++] = Line(x1, y1, x2, y2, k, b);
    }
    printf("%.2lf\n", asr(l, r1, 1e-5) * 2);
}
int main()
{
    solve();
    return 0;
}