BZOJ1502 月下柠檬树

Link

Solution

Simpson积分,拟合二次函数曲线。 公式为\(\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;
}

Powered by Hexo and Hexo-theme-hiker

Copyright © 2013 - 2017 Hey,Lucida here. All Rights Reserved.

Lucida Lu 保留所有权利