BZOJ1502 月下柠檬树

发表于 | 更新于

Link

Solution

Simpson积分，拟合二次函数曲线。 公式为(r−l)(f(l)+f(r)+4f(mid))6\frac {(r-l)(f(l)+f(r)+4f(mid))}{6} 涉及圆的公切线的求法，要用到一点点平面几何。

Code

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
100
101
102
103
104
105
106
107
//Code by Lucida
#include<bits/stdc++.h>
#define red(x) scanf("%d",&x)
#define fred(x) scanf("%lf",&x)
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
typedef double ld;
using std::swap;
const int MAXN=500+1;
const ld eps=1e-6,pi=acos(-1.0),INF=1e100;
int fcmp(ld x)
{
    if(-eps<x && x<eps) return 0;
    return x<0?-1:1;
}
template <class T> T sqr(T x){return x*x;}
template <class T> T abs(T x){return x>0?x:-x;}
struct vec
{
    ld x,y;
    vec(ld _x=0,ld _y=0):x(_x),y(_y){}
    vec operator -(){return vec(-x,-y);}
    ld norm(){return sqrt(x*x+y*y);}
    void operator +=(vec a){x+=a.x,y+=a.y;}
    void operator -=(vec a){x-=a.x,y-=a.y;}
    void operator /=(ld a){x/=a,y/=a;}
    void operator *=(ld a){x*=a,y*=a;}
    bool operator ==(vec a){return x==a.x && y==a.y;}
	bool operator <(vec a){if(x!=a.x) return x<a.x;return y<a.y;}
};
typedef vec point;
vec operator +(vec a,vec b){return vec(a.x+b.x,a.y+b.y);}
vec operator -(vec a,vec b){return vec(a.x-b.x,a.y-b.y);}
vec operator *(vec a,ld t){return vec(a.x*t,a.y*t);}
vec operator /(vec a,ld t){return vec(a.x/t,a.y/t);}
ld inner(vec a,vec b){return a.x*b.x+a.y*b.y;}
ld outer(vec a,vec b){return a.x*b.y-a.y*b.x;}
struct circle
{
	ld x,r;
	circle(){}
	circle(ld _x,ld _r):x(_x),r(_r){}
	bool cover(circle a){return fcmp(abs(x-a.x)+a.r-r)<=0;}
}c[MAXN];
struct line
{
	point a;vec v;
	line(){}
	line(point _a,vec _v):a(_a),v(_v){}
}le[MAXN];
int n;
ld f(ld x)
{
	ld res=0;
	for(int i=1;i<=n;i++)
	{
		if(fcmp(c[i].x-c[i].r-x)<0 && fcmp(c[i].x+c[i].r-x)>0)
			chkmx(res,sqrt(sqr(c[i].r)-sqr(x-c[i].x)));
		point p=le[i].a,q=le[i].a+le[i].v;
		//if(p>q) swap(p,q);
		if(fcmp(p.x-x)<0 && fcmp(x-q.x)<0)
			chkmx(res,p.y+(x-p.x)/(q.x-p.x)*(q.y-p.y));
	}
	return res;
}
ld Simpson(ld l,ld r,ld lv,ld rv,ld mv){return (r-l)*(lv+rv+4*mv)/6;}
ld calc(ld l,ld r,ld lv,ld rv)
{
	ld mid=(l+r)/2,mv=f(mid),res;
	if(fcmp((res=Simpson(l,r,lv,rv,mv))-Simpson(l,mid,lv,mv,f((l+mid)/2))-Simpson(mid,r,mv,rv,f((mid+r)/2))))
		res=calc(l,mid,lv,mv)+calc(mid,r,mv,rv);
	return res;
}

line comtan(circle a,circle b)
{
	if(a.cover(b) || b.cover(a))
		return line(point(0,0),vec(1,0));
	if(a.x>b.x) swap(a,b);//necs
	point p,q;
	p.x=(a.r-b.r)/abs(a.x-b.x)*a.r,p.y=sqrt(sqr(a.r)-sqr(p.x));
	q.x=b.r/a.r*p.x,q.y=b.r/a.r*p.y;
	p.x+=a.x,q.x+=b.x;
	return line(p,q-p);
}
int main()
{
	static ld h[MAXN],ra[MAXN];
	//freopen("input","r",stdin);
	ld alpha;
	red(n),fred(alpha);
	alpha=1.0/tan(alpha);
	n++;
	for(int i=1;i<=n;i++) fred(h[i]),h[i]*=alpha,h[i]+=h[i-1];
	for(int i=1;i<n;i++) fred(ra[i]);ra[n]=0;
	ld l=INF,r=-INF;
	for(int i=1;i<=n;i++)
	{
		c[i]=circle(h[i],ra[i]);
		chkmn(l,h[i]-ra[i]);
		chkmx(r,h[i]+ra[i]);
	}
	for(int i=1;i<n;i++)
		le[i]=comtan(c[i],c[i+1]);
	printf("%.2lf\n",calc(l,r,f(l),f(r))*2);
	return 0;
}