【洛谷】P4207 [NOI2005]月下柠檬树

题解

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

代码

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <ctime>
//#define ivorysi
#define MAXN 505
#define eps 1e-7
#define pb push_back
using namespace std;
typedef long long int64;
typedef unsigned int u32;
typedef double db;
inline db o(db x) {
    return x * x;
}
bool dcmp(db a,db b) {
    return fabs(a - b) <= eps;
}
struct Point {
    db x,y;
    Point(){}
    Point(db _x,db _y) {
	x = _x;y = _y;
    }
    friend Point operator + (const Point &a,const Point &b) {
	return Point(a.x + b.x,a.y + b.y);
    }
    friend Point operator - (const Point &a,const Point &b) {
	return Point(a.x - b.x,a.y - b.y);
    }
    friend Point operator * (const Point &a,const db &d) {
	return Point(a.x * d,a.y * d);
    }
    friend db operator * (const Point &a,const Point &b) {
	return a.x * b.y - a.y * b.x;
    }
    db norm() {
	return sqrt(x * x + y * y);
    }
    friend db dot(const Point &a,const Point &b) {
	return a.x * b.x + a.y * b.y; 
    }
}P[MAXN];
struct Seg {
    Point a,b;
    db d;
    Seg(){}
    Seg(Point _a,Point _b) {
	a = _a;b = _b;
	d = (b.y - a.y) / (b.x - a.x);
    }
    friend Point Cross_Point(const Seg &s,const Seg &t) {
	db S1 = (s.a - t.a) * (t.b - t.a);
	db S2 = (s.b - t.b) * (t.a - t.b);
	return s.a + (s.b - s.a) * (S1 / (S1 + S2));
    }
}S[MAXN];
db alpha;
db H[MAXN],R[MAXN];
int N;
db F(db x) {
    db res = 0;
    for(int i = 1 ; i <= N ; ++i) {
	if(i != N) {
	    if(P[i + 1].x - P[i].x > fabs(R[i] - R[i + 1])) {
		if(S[i].a.x + eps <= x && x <= S[i].b.x - eps) {
		    res = max(res,S[i].a.y + (x - S[i].a.x) * S[i].d);
		}
	    }
	}
	if(x >= P[i].x - R[i] && x <= P[i].x + R[i]) {
	    res = max(res,sqrt(o(R[i]) - o(P[i].x - x)));
	}
    }
    return res;
}
db calc(db L,db R) {
    db mid = (L + R) / 2;
    return (R - L) * (F(L) + F(R) + 4 * F((L + R) / 2)) / 6;
}
db Simpson(db L,db R) {
    db mid = (L + R) / 2;
    db Slr = calc(L,R),Sl = calc(L,mid),Sr = calc(mid,R);
    if(dcmp(Slr,Sl + Sr)) return Sl + Sr;
    else return Simpson(L,mid) + Simpson(mid,R);
}
void Init() {
    scanf("%d%lf",&N,&alpha);
    ++N;
    for(int i = N ; i >= 1; --i) scanf("%lf",&H[i]);
    R[1] = 0;
    for(int i = N ; i >= 2; --i) scanf("%lf",&R[i]);
    for(int i = 1 ; i <= N ; ++i) H[i] /= tan(alpha);
    P[1] = Point(0,0);
    for(int i = 2 ; i <= N ; ++i) {
	P[i] = Point(P[i - 1].x + H[i - 1],0);
    }
}
void Solve() {
    db r = 0,l = 0;
    for(int i = 1 ; i <= N ; ++i) {
	l = min(P[i].x - R[i],l);
	r = max(r,P[i].x + R[i]);
    }
    for(int i = 1 ; i < N ; ++i) {
	if(P[i + 1].x - P[i].x > fabs(R[i] - R[i + 1])) {
	    db C = (R[i] - R[i + 1]) / (P[i + 1].x - P[i].x);
	    db s1 = P[i].x + R[i] * C,s2 = P[i + 1].x + R[i + 1] * C;
	    db h1 = sqrt(o(R[i]) - o(P[i].x - s1)),h2 = sqrt(o(R[i + 1]) - o(P[i + 1].x - s2));
	    S[i] = Seg(Point(s1,h1),Point(s2,h2));
	}
    }
    printf("%.2lf
",Simpson(l,r) * 2);
}
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    Init();
    Solve();
}