[NOI2005] 月下柠檬树

阅读本题解前，您需要掌握三角函数的有关知识

本文中有大篇幅的数学公式，如有头晕，目眩，癫痫等不适症状，请立即停止阅读本文并及时就医

前置前置芝士 微积分

微积分的基本定理：

假设 $F'(x) = f(x)$ ,那么有：

$$\int_a^b f(x) dx = F(b)- F(a)$$

(假设 $a < b$)

因为指数函数的微分公式 $f(x) = x ^ n \Rightarrow f'(x) = nx^{n-1}$

所以，显然

$$\int_a^b (Ax^2+Bx+C)dx$$

$$= \int_a^b (\frac{Ax^3}{3} + \frac{Bx^2}{2} + \frac{Cx}{1})$$

$$= \frac{Ab^3}{3} - \frac{Aa^3}{3} + \frac{Bb ^2}{2} - \frac{Ba ^ 2}{2} + Cb - Ca$$

前置芝士 Simpson 公式

我看到题解中没有一篇认真讲过辛普森公式，于是就来水一篇题解。

基本思想

将函数近似看做二次函数

由于这个原因，这个公式及其的简单且不精确。

公式推导

$$\int_a^b f(x) dx$$

上方是最最基础的微积分求面积的公式；

刚刚的思想中讲到了思路是将 $f(x)$ 看做 $ax^2 + bx + c$ ,于是有:

$$= \int_a^b (Ax^2+Bx+C)dx$$

$$\approx \int_a^b (\frac{Ax^2}{3} + \frac{Bx}{2} + \frac{C}{1})$$

$$= \frac{Ab^3}{3} - \frac{Aa^3}{3} + \frac{Bb ^2}{2} - \frac{Ba ^ 2}{2} + Cb - Ca$$

$$= \frac{2Ab^3 - 2Aa^3 + 3Bb^2 - 3Ba^2 + 6b - 6CaC}{6}$$

$$= \frac{2A(a^3 - b^3) +3B(b^2 - a^2) + 6C(b-a)}{6}$$

$$= \frac{2A(b-a)(a^2 + b^2 + ab) + 3B(a + b)(b-a) + 6C(b-a)}{6}$$

$$= \frac{(b-a)(2Aa^2 + 2Ab^2 + 2Aab + 3Ba + 3Bb + 6C)}{6}$$

$$= (b-a)\frac{(Aa^2 + Ba + C) + (Ab^2 + Bb + C) + (Aa^2 + Ab^2 + 2Ba + 2Bb + 2Aab + 4C)}{6}$$

$$= (b-a)\frac{f(a) + f(b) + 4(\frac{Aa^2}{4} + \frac{Ab^2}{4} + \frac{Ba}{2} + \frac{Bb}{2} + \frac{2Aab}{4} + C)}{6}$$

$$= (b-a)\frac{f(a) + f(b) + 4(\frac{Aa^2 + Aab}{4} + \frac{Ab^2 + Aab}{4} + \frac{Ba + Bb}{2} + C)}{6}$$

$$= (b-a)\frac{f(a) + f(b) + 4(A\frac{a(a + b)}{4} + A\frac{b(a + b)}{4} + B\frac{a +b}{2} + C)}{6}$$

$$= (b-a)\frac{f(a) + f(b) + 4(A\frac{(a + b) ^ 2}{4} + B\frac{a +b}{2} + C)}{6}$$

$$= (b-a)\frac{f(a) + f(b) + 4f(\frac{a + b}{2})}{6}$$

所以，我们得到：

$$\int_a^b f(x) dx \approx (b-a)\frac{f(a) + f(b) + 4f(\frac{a + b}{2})}{6}$$

处理精度问题

为了控制精度，以及不必要的问题，保留两位小数时 eps 一般取 $10^{-6} \sim 10^{-8}$

有了 Simpson 公式，一个自然的想法是把积分区间拆成多个小区间后求和。

但是分成区间的个数和长度因积分区间和精度要求甚至被积函数而异。

换句话说，分的区间数太少不满足精度要求，太多了会 TLE 。

「自适应」就是求积分时能够自动控制切割的区间大小和长度，使精度能满足要求。

具体地： solve(l , r) 表示求 $\int_r^l f(x)dx$

1. 取 mid = l+r >> 1;
2. 用 Simpson 公式近似计算 $f(x)$ 在区间 $[l,mid]$ 和区间 $[mid,r]$ 内的积分，分别为 $l_s$ 和 $r_s$ ，及 $[l,r]$ 的近似积分 $sum$ ;
3. 如果 $l_s + r_s$ 与 $sum$ 的误差允许(即在 eps 之内)，则返回 $sum$ ;
4. 否则递归 solve(l , mid) 和 solve(mid , r) ，返回这两个递归计算结果的和;

摘自 这篇博文。

代码

double Function (double x) {
	//Your Function!
}

double Solve (double l , double r) {
	double mid = (l + r) / 2;
	return (Function(l) + Function(r) + Function(mid) * 4) * (r - l) / 6;
}

double Simpson (double l , double r , double eps , double res) {
        double mid = (l + r) / 2;
        double x1 = Solve(l , mid) , x2 = Solve(mid , r);
        if(abs(x1 + x2 - res) <= 15 * eps) return x1 + x2 + (x1 + x2 - res) / 15;
        else return Simpson(l , mid , eps / 2 , x1) + Simpson(mid , r , eps / 2 , x2);
}

思路

对于一个圆，投影下去还是一个圆；

定义树的轴在投影平面上经过的点为原点，定一个正方向，建立平面直角坐标系，
能发现，对于一个半径为 $r$，高度为 $h$ 的圆，投影到平面上是圆心坐标为 $(\cot(α)h,0)$，半径为 $r$ 的圆；

想象有一个水平的平面，竖直向上移，可以把树切出一堆圆，对于这些圆，把它们投影求个并就是答案；

对于每个圆台，它一堆圆的并就是先求上下两个面的圆的投影，再对投影求外公切线，围成的图形；

这个图形挺特殊，所以可以对不规则的函数下方面积考虑使用自适应 Simpson；

然后，这道题目就做完了。

代码实现

#include <bits/stdc++.h>
using namespace std;

const double eps = 1e-6;

const int N = 510;

struct Point {
	double x , y;
}S[N] , E[N] , A[N] , B[N];

double Alpha;
int n;
double h;

double sqr (double x) {
	return x * x;
}

void Cal (Point &s , Point &e , Point a , Point b) {
	if(abs(a.y - b.y) < eps) return s = a , e = b , void();
	double k = (b.y - a.y) / (b.x - a.x);
	s.x = a.x - k * a.y;
	s.y = sqrt(sqr(a.y) - sqr(a.x - s.x));
	e.x = b.x - k * b.y;
	e.y = sqrt(sqr(b.y) - sqr(b.x - e.x));
}

double Function (double x) {
	double y = 0;
	for (int i = 1; i <= n + 1; i++) {
		if(abs(x - A[i].x) <= A[i].y) {
			y = max(y , sqrt(sqr(A[i].y) - sqr(x - A[i].x)));
		}
	}
	for (int i = 1; i <= n; i++) {
		if(A[i + 1].x - A[i].x - abs(A[i].y - A[i + 1].y) > eps && S[i].x <= x && x <= E[i].x) {
			y = max(y , S[i].y + (x - S[i].x) * (E[i].y - S[i].y) / (E[i].x - S[i].x));
		}
	}
	return y;
}

double Solve (double l , double r) {
	double mid = (l + r) / 2;
	return (Function(l) + Function(r) + Function(mid) * 4) * (r - l) / 6;
}

double Simpson (double l , double r) {
	double mid = (l + r) / 2;
	double x1 = Solve(l , mid) , x2 = Solve(mid , r) , x3 = Solve(l , r);
	if(abs(x1 + x2 - x3) < eps) return x1 + x2;
	else return Simpson(l , mid) + Simpson(mid , r);
}

int main() {
	cin >> n >> Alpha;
	Alpha = 1 / tan(Alpha);
	
	for (int i = 1; i <= n + 1; i++) {
		cin >> A[i].x;
		h += A[i].x;
		A[i].x = h * Alpha;
	}
	
	for (int i = 1; i <= n; i++) {
		cin >> A[i].y;
	}
	
	A[n + 1].y = A[n + 2].y = 0;
	
	double l , r;
	
	l = A[1].x;
	r = A[n + 1].x;
	
	for (int i = 1; i <= n; i++) {
		l = min(l , A[i].x - A[i].y);
		r = max(r , A[i].x + A[i].y);
		if (A[i + 1].x - A[i].x - abs(A[i + 1].y - A[i].y) > eps) {
			Cal(S[i] , E[i] , A[i] , A[i + 1]);
		}
	}
	
	printf("%.2lf\n" ,2 * Simpson(l , r));
	
	return 0;
}