「NOI2005」月下柠檬树 - Simpson自适应 | Bill Yang's Blog

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

题目大意

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

题目大意2

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


初步想法

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

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

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


算法确定

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

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


代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
#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;
}
本篇文章有用吗?有用还不快点击!
0%