「NOI2005」月下柠檬树 - Simpson自适应

题目大意

有一棵由圆台组成的树，平行光线以∠$\alpha$射向这棵树，求地上影子的面积。

题目大意2

求本人心理阴影面积。
记得这道题是从初一开始就膜的神题，今天总算正式A了。

---

初步想法

重要条件：平行光线。
已知：平行光线对二维图形的投影形状与大小完全一样。
平行光线对于图形的唯一影响就是位置发生了改变。

那么对于三维物体，拆分成无数个二维图形。
对于本题，将圆台拆分成无数个圆，投影到地上它们有相同的公切线。
也就是说，仅仅是每一个圆台的上下顶起着确定图形范围的作用，整个影子是由圆与它们的公切线组成。
样例对应的投影就是下图：

那么本题就转化成了求上图的面积。

---

算法确定

求图形面积并，求矩形面积并我们可以用矩形切割或线段树。
但是求非规则图形面积并怎么办？
我们可以用Simpson自适应。
Simpson积分公式：

$$\int_a^b {f(x)} \, d x\approx \frac{b-a}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]$$

Simpson自适应的学习笔记见这里。
对于此题，我们求出相邻两个圆的切线，然后调用Simpson自适应，但是Simpson自适应需要调用函数值，这图像根本不是函数啊。
但是我们可以用图形的最高点代替其函数值。
对于此题，我们可以判断$x$是否会与线段、圆相交，然后算出其在线段、圆上对应的$y$值，取一个最大值即为答案。
对于如何计算公切线，如何计算点值等可以使用中学几何的方法画图解决，此处不细讲，可参考代码。

---

代码

#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;
inline const int Get_Int() {
	int num=0,bj=1;
	char x=getchar();
	while(x<'0'||x>'9') {
		if(x=='-')bj=-1;
		x=getchar();
	}
	while(x>='0'&&x<='9') {
		num=num*10+x-'0';
		x=getchar();
	}
	return num*bj;
}
const int maxn=1005;
const double eps=1e-8;
int dcmp(double x) {
	if(fabs(x)<=eps)return 0;
	if(x>eps)return 1;
	return -1;
}
struct Point {
	double x,y;
	Point() {}
	Point(double _x,double _y):x(_x),y(_y) {}
	Point operator - (const Point& a) const {
		return Point(a.x-x,a.y-y);
	}
};
typedef Point Vector;
struct Line {
	Point p1,p2;
	double k,b;
	Line() {}
	Line(Point _p1,Point _p2):p1(_p1),p2(_p2) {
		Vector _v=_p1-_p2;
		k=_v.y/_v.x;
		b=_p1.y-k*_p1.x;
	}
	double f(double x) {
		return k*x+b;
	}
} l[maxn*2];
int cnt=0,n;
double alpha,h[maxn];
struct Circle {
	double x,r;
	Circle() {}
	Circle(double _x,double _r):x(_x),r(_r) {}
} a[maxn];
void Tanget(Circle a,Circle b) { //求圆公切线 
	double _sin=(a.r-b.r)/(b.x-a.x),_cos=sqrt(1-_sin*_sin),_tan=_sin/_cos;
	l[++cnt]=Line(Point(a.x+a.r*_sin,a.r*_cos),Point(b.x+b.r*_sin,b.r*_cos));
}
double f(double x) {
	double Max=0;
	for(int i=1; i<=cnt; i++)
		if(x>=l[i].p1.x&&x<=l[i].p2.x)Max=max(Max,l[i].f(x)); //直线上 
	for(int i=1; i<=n; i++)
		if(x>=a[i].x-a[i].r&&x<=a[i].x+a[i].r)Max=max(Max,sqrt(a[i].r*a[i].r-(x-a[i].x)*(x-a[i].x))); //圆上 
	return Max;
}
double Cal(double Left,double Right) {
	double mid=(Left+Right)/2;
	return (Right-Left)/6*(f(Left)+4*f(mid)+f(Right));
}
double Simpson(double Left,double Right,double ans) {
	double mid=(Left+Right)/2,ans1=Cal(Left,mid),ans2=Cal(mid,Right);
	if(dcmp(ans-ans1-ans2)==0)return ans;
	return Simpson(Left,mid,ans1)+Simpson(mid,Right,ans2);
}
int main() {
	ios::sync_with_stdio(false);
	cin>>n>>alpha;
	for(int i=1; i<=n+1; i++) {
		cin>>h[i];
		h[i]+=h[i-1];
	}
	for(int i=1; i<=n; i++)cin>>a[i].r;
	double Tan=tan(alpha),Left=1e15,Right=-1e15;
	a[n+1]=Circle(h[n+1]/Tan,0); //顶点半径为0 
	Left=min(Left,a[n+1].x),Right=max(Right,a[n+1].x);
	for(int i=n; i>=1; i--) {
		a[i].x=h[i]/Tan;
		Left=min(Left,a[i].x-a[i].r),Right=max(Right,a[i].x+a[i].r);
		if(dcmp(a[i+1].x-a[i].x-fabs(a[i+1].r-a[i].r))>=0)Tanget(a[i],a[i+1]); //不包含则求公切线 
	}
	printf("%0.2lf\n",Simpson(Left,Right,Cal(Left,Right))*2); //辛普森仅计算了y>0部分 
	return 0;
}