Geometry

这里把大部分板子都放过来了。

除了三维计算几何和圆反演的内容(这两块恐怕要等我上大学之后)。。。

讲的会比较简略,可能会把 Whk 的内容搬过来。

然后告诫一句:这个东西 OI 基本考不到,花太多时间在上面是不值得的。。。

为什么我会写这个东西?无聊而已。。。而且我想多学一点东西,不仅仅是为了做题。

这些东西都是《算法竞赛入门到进阶》上的,一部分是杜老师的,反正都是干货(毕竟在 OI 里似乎没有什么讲头,而且非常码农)。

那张闵可夫斯基和的图片是 shadowice1984 的。

要注明的话东西比较多,又不是写论文,所以我就懒了 QAQ。

Luogu 只有三道模板,二维凸包、旋转卡壳和半平面交。三维的还有一部分模板但我不会。。。想要测模板可以自己去交。。。

一个事情就是这里的大部分模板没有测过(除了上述模板题以及其相关的函数),所以很有可能有锅 QAQ,但我也懒得对拍了,所以非常抱歉,对拍就咕了。

例题会补上的(咕咕咕)。

一、常数定义

eps 用来保证精度,Pi 后面会用。

代码:

1
2
3
// Some cosntant.
const ld Pi=acos(-1),eps=1e-10;
const ll N=1e5,inf=1e15;

二、Sgn 和 Cmp 函数

为了方便保证精度写的函数。

代码:

1
2
3
// The Sgn function & Cmp function.
inline ll Sgn(ld x) {return x<-eps?-1:x>eps;}
inline ll Cmp(ld x,ld y) {return Sgn(x-y);}

三、点与向量

点可以表示向量,所以就不 typedef 了。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// Point & Vector.
struct Point{
ld x,y;
inline Point() {}
inline Point(ld x_,ld y_):x(x_),y(y_){}
inline Point operator+(Point p) {return Point(x+p.x,y+p.y);}
inline Point operator-(Point p) {return Point(x-p.x,y-p.y);}
inline Point operator*(ld d) {return Point(x*d,y*d);}
inline Point operator/(ld d) {return Point(x/d,y/d);}
inline bool operator==(Point p) {return !Sgn(x-p.x)&&!Sgn(y-p.y);}
inline bool operator<(Point p) {
return Sgn(x-p.x)<0||(!Sgn(x-p.x)&&Sgn(y-p.y)<0);
}
};
Point p[N+5],tmp_p[N+5];

四、计算两点间距离

勾股定理。。。

代码:

1
2
3
4
// To calculate the distance of two Points.
inline ld Dist(Point a,Point b) {
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

五、向量点积类操作

坐标下的点积运算不用说了吧?

$$\begin{aligned}\mathbf{a}\cdot\mathbf{b} & =|\mathbf{a}||\mathbf{b}|\cos \theta\\ &=(x_a,y_a)\cdot(x_b,y_b)\\ &= x_ax_b-y_ay_b\end {aligned}$$

可以用这个判断垂直(以及夹角的锐钝),计算夹角,通过平方还可以计算向量长度。

实际上就是:

$$\mathbf{a}\cdot\mathbf{b}=0\Leftrightarrow \mathbf{a}\perp \mathbf{b}$$

$$\cos \theta =\dfrac{\mathbf{a\cdot b}}{|\mathbf{a}||\mathbf{b}|}$$

$$|\mathbf{a}|=\sqrt{\mathbf{a^2}}$$

然后就是表示运算的这个点($\cdot$)是不能省略的。

1
2
3
4
5
// Dot product.
inline ld Dot(Point a,Point b) {return a.x*b.x+a.y+b.y;}
inline ld Len(Point a) {return sqrt(Dot(a,a));}
inline ld Len2(Point a) {return Dot(a,a);}
inline ld Angle(Point a,Point b) {return acos(Dot(a,b)/Len(a)/Len(b));}

六、向量叉积类操作

叉积众所周知:

$$\begin{aligned}\mathbf{a}\times \mathbf{b}&= |\mathbf{a}||\mathbf{b}|\sin \theta\\ &= (x_a,y_a)\times (x_b,y_b)\\ &= x_ay_b-x_by_a\end{aligned}$$

这个东西的功用多一点,可以算两个向量所夹的平行四边形的面积(就是叉积值,所以下面那个函数表示求一个三角形的面积的两倍,但是一般情况下这个值是有正负的,表示顺时针还是逆时针),还可以判断平行之类的(叉积值等于 0)。

一些 Trivial 的操作一并附在下面了,就顺带讲掉。

一个是将向量旋转 $\theta$。

众所周知啊,复数相乘在复平面上表现为模长相乘,辅角相加。

设向量 $\mathbf{a}=x+yi$,那么旋转 $\theta$ 相当于乘上一个单位向量 $\mathbf{e}$,其辅角为 $\theta$,显然 $\mathbf{e}=\cos \theta+i\sin\theta$。

那么就有:

$$\begin{aligned}\mathbf{a}\times \mathbf{e}&=(x+yi)(\cos \theta+i\sin\theta)\\&=x\cos\theta-y\sin\theta+(x\sin\theta+y\cos\theta)i\end{aligned}$$

坐标不用我讲了吧?

以及一个求法向量的操作,其实就是将向量逆时针旋转 $\dfrac{\pi}{2}$,然后再变成单位向量即可。

逆时针旋转 $\dfrac{\pi}{2}$ 即:

$$(x,y)\rightarrow(-y,x)$$

代码:

1
2
3
4
5
6
7
8
// Cross product.
inline ld Cross(Point a,Point b) {return a.x*b.y-a.y*b.x;}
inline ld Area2(Point a,Point b,Point c) {return Cross(b-a,c-a);}
inline bool Parallel(Point a,Point b) {return !Sgn(Cross(a,b));}
inline Point Rotate(Point a,ld rad) {
return Point(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
inline Point Normal(Point a) {return Point(-a.y/Len(a),a.x/Len(a));}

七、直线与线段

这一块的东西巨大多繁琐。

  1. 定义结构体中的直线由两个点表示(在后面半平面的表示中含义会有所不同),不需要极角序的话可以不用记录角度。结构体中后面的几个析构函数分别是将点角表示和解析式表示转化为两点表示,最后的重载运算符是为了极角排序用的。

  2. 点与直线的关系,用叉积正负判断是在直线逆时针方向还是在顺时针方向还是在直线上。

  3. 点与线段的关系。首先根据叉积判断点在直线上,再根据点积正负判断点在线段上。

  4. 点与直线的距离,就是高的长度。其实就是构造一个平行四边形然后除以高。因为有正负之分所以取绝对值。

  5. 然后是点到线段的投影(因为是点所以最后还要加上起点)。实际上就是:

$$\begin{aligned}|\mathbf{a}|\cos\theta\mathbf{e}&=|\mathbf{a}|\dfrac{\mathbf{a}\cdot\mathbf{b}}{|\mathbf{a}||\mathbf{b}|}\cdot\dfrac{\mathbf{b}}{|\mathbf{b}|}\\&= \dfrac{\mathbf{a}\cdot\mathbf{b}}{\mathbf{b}^2}\cdot\mathbf{b}\end{aligned}$$

  1. 点关于直线的对称点。先求投影向量得到一个平行四边形,利用其对角线与边的加减的到对称点的向量表示。

  2. 点与线段的距离。实际上就是分类讨论一下,如果就是高的话直接求即可,反之取到两端点的距离的最小值即可。是否是高根据点积来得到夹角的钝锐性判断。

  3. 直线之间的关系,先判断叉积是否为 0(这样的话不是平行就是重合,根据点是否在直线上判断即可),反之肯定是相交。

  4. 求直线交点。可以联立求解,这里将一个几何方法:

    图中有直线 $AB$ 与 $CD$ 交于点 $P$,根据这个面积、底与高的简单关系和相似可知:

    $$\dfrac{DP}{CP}=\dfrac{S_{\triangle ABD}}{S_{\triangle ABC}}=\dfrac{\overrightarrow{AD}\times\overrightarrow{AB}}{\overrightarrow{AB}\times\overrightarrow{AC}}=\dfrac{x_D-x_P}{x_P-x_C}=\dfrac{y_D-y_P}{y_P-y_C}$$

    简单变换一下就成了:

    $$x_P=\dfrac{S_{\triangle ABD}x_C+S_{\triangle ABC}x_D}{S_{\triangle ABD}+S_{\triangle{ABC}}}$$

    $$y_P=\dfrac{S_{\triangle ABD}y_C+S_{\triangle ABC}y_D}{S_{\triangle ABD}+S_{\triangle ABC}}$$

    面积显然可以叉积求,这个东西的比值显然可以被叉积之间互相消掉,细节就不讲了。

    当然因为这里有做除法的危险操作,如果题目没有说明三点共线的情况不存在一定要特判!!!

    直线交点

  5. 线段相交。判断对于每一条线段,另外一条线段的两个点是否分别在这条线段的顺时针和逆时针方向。

  6. 线段交点。先判断是否相交,再判断直线交点。

代码:

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
// Line.
struct Line{
Point p1,p2;ld ang;
inline Line() {}
inline Line(Point p1_,Point p2_):p1(p1_),p2(p2_){ang=atan2(p2_.y,p2_.x);}
inline Line(Point p,double angle) {
p1=p;if(Sgn(angle-Pi/2)==0) {p2=(p1+Point(0,1));}
else {p2=(p1+Point(1,tan(angle)));}
}
inline Line(ld a,ld b,ld c) {
if(!Sgn(a)) {p1=Point(0,-c/b);p2=Point(1,-c/b);}
else if(!Sgn(b)) {p1=Point(-c/a,0);p2=Point(-c/a,1);}
else {p1=Point(0,-c/b);p2=Point(1,(-c-a)/b);}
}
inline bool operator<(Line &l) {return ang<l.ang;}
};
inline ll Point_Line_Relation(Point p,Line v) {
ll c=Sgn(Cross(p-v.p1,v.p2-v.p1));
if(c<0) return 1;if(c>0) return 2;return 0;
}
inline bool Point_on_Seg(Point p,Line v) {
return !Sgn(Cross(p-v.p1,v.p2-v.p1))&&Sgn(Dot(p-v.p1,p-v.p2))<=0;
}
inline ld Dis_Point_Line(Point p,Line v) {
return fabs(Cross(p-v.p1,v.p2-v.p1))/Dist(v.p1,v.p2);
}
inline Point Point_Line_Proj(Point p,Line v) {
ld k=Dot(v.p2-v.p1,p-v.p1)/Len2(v.p2-v.p1);
return v.p1+(v.p1-v.p2)*k;
}
inline Point Point_Line_Symmetry(Point p,Line v) {
Point q=Point_Line_Proj(p,v);
return Point(2*q.x-p.x,2*q.y-p.y);
}
inline ld Dis_Point_Seg(Point p,Line v) {
if(Sgn(Dot(p-v.p1,v.p2-v.p1))<0||Sgn(Dot(p-v.p2,v.p1-v.p2))<0) {
return min(Dist(p,v.p1),Dist(p,v.p2));
}
return Dis_Point_Line(p,v);
}
inline ll Line_Relation(Line v1,Line v2) {
if(Sgn(Cross(v1.p2-v1.p1,v2.p2-v2.p1))==0) {
if(Point_Line_Relation(v1.p1,v2)==0) return 1;
else return 0;
}
return 2;
}
inline Point Cross_Point(Line v1,Line v2) {
ld s1=Cross(v1.p2-v1.p1,v2.p1-v1.p1),s2=Cross(v1.p2-v1.p1,v2.p2-v1.p1);
return Point(v2.p1.x*s2-v2.p2.x*s1,v2.p1.y*s2-v2.p2.y*s1)/(s2-s1);
}
inline bool Cross_Seg(Line v1,Line v2) {
ld c1=Cross(v1.p2-v1.p1,v2.p1-v1.p1),c2=Cross(v1.p2-v1.p1,v2.p2-v1.p1);
ld d1=Cross(v2.p2-v2.p1,v1.p1-v2.p1),d2=Cross(v2.p2-v2.p1,v1.p2-v2.p1);
return Sgn(c1)*Sgn(c2)<=0&&Sgn(d1)*Sgn(d2)<=0;
}

八、多边形

  1. 判断一个点是否在多边形内部。

    有一种方法是从该点引出一条射线,然后判断射线与多边形的边相交的次数,奇数次说明在多边形内部,偶数次说明在外部。

    另一种方法是转角法,每个点连接,然后绕多边形一周计算转角,如果说在转过 $2\pi$ 说明在内部,$\pi$ 说明在边界上。但是使用反三角函数会有精度问题。

    我们利用转角法的思想,换一种实现方式。与每一条边检查如下的东西:

    $$c=Cross(P-j,i-j)$$

    $$u=i.y-P.y$$

    $$v=j.y-P.y$$

    根据 $u$,$v$ 检查过 $P$ 的水平线是否经过该边,根据 $c$ 判断是在顺时针方向还是在逆时针方向。

    然后我们用一个 $num$ 来计数,若 $num>0$ 说明在内部;否则在外部(在边上可以单独判断)。这里具体看代码。我不太能解释,但是感性理解就是对的(确信)。

    1. 求多边形面积。叉积啥的可以相互抵消,所以实际上这个面积就是以原点做三角剖分,再把叉积加起来除以二(因为叉积的值是三角形面积值的两倍)。

    2. 求多边形重心。给多边形作三角剖分,按照三角形的有向面积求加权平均,然后就能求出重心(然而我也不知道什么原理)。

代码:

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
// Polygon.
inline ll Point_in_Polygon(Point pt,Point *p,ll n) {
for(ll i=0;i<n;i++) {if(p[i]==pt) return 3;}
for(ll i=0;i<n;i++) {
Line v=Line(p[i],p[(i+1)%n]);
if(Point_on_Seg(pt,v)) return 2;
}
ll num=0;
for(ll i=0;i<n;i++) {
ll j=(i+1)%n;
ll c=Sgn(Cross(pt-p[j],p[i]-p[j]));
ll u=Sgn(p[i].y-pt.y),v=Sgn(p[j].y-pt.y);
if(c>0&&u<0&&v>=0) num++;
if(c<0&&u>=0&&v<0) num--;
}
return num!=0;
}
inline ld Polygon_Area(Point *p,ll n) {
ld area=0;for(ll i=0;i<n;i++) area+=Cross(p[i],p[(i+1)%n]);return area/2;
}
inline Point Polygon_Centre(Point *p,ll n) {
Point ans(0,0);if(Polygon_Area(p,n)==0) return ans;
for(ll i=0;i<n;i++) ans=ans+(p[i]+p[(i+1)%n])*Cross(p[i],p[(i+1)%n]);
return ans/Polygon_Area(p,n)/6;
}

九、二维凸包

这里介绍 Graham 扫描法的变种 Andrew 算法。

算法做两次扫描,先从最左边的点沿“下凸包”扫描到最右边,再从最右边的点沿“上凸包”扫描到最左边,两个凸包合起来就是完整的凸包。

然后,按照横坐标从小到大排序,然后按照这个拐弯来判断凸包,实际上还是挺朴素的。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// Convex hull.
inline ll Convex_Hull(Point *p,ll n,Point *ch) {
sort(p,p+n);n=unique(p,p+n)-p;ll v=0;
for(ll i=0;i<n;i++) {
while(v>1&&Sgn(Cross(ch[v-1]-ch[v-2],p[i]-ch[v-2]))<=0) v--;
ch[v++]=p[i];
}
ll j=v;
for(ll i=n-2;i>=0;i--) {
while(v>j&&Sgn(Cross(ch[v-1]-ch[v-2],p[i]-ch[v-2]))<=0) v--;
ch[v++]=p[i];
}
if(n>1) v--;
return v;
}

十、最近点对

懒了,丢个链接

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// Closest pair.
inline bool cmpy(Point a,Point b) {return Sgn(a.y-b.y)<0;}
inline ld Closest_Pair(ll left,ll right) {
ld dis=inf;if(left==right) return dis;
if(left+1==right) return Dist(p[left],p[right]);
ll mid=(left+right)/2;
ld d1=Closest_Pair(left,mid),d2=Closest_Pair(mid+1,right);
dis=min(d1,d2);
ll k=0;
for(ll i=left;i<=right;i++) {
if(fabs(p[mid].x-p[i].x)<=dis) tmp_p[k++]=p[i];
}
sort(tmp_p,tmp_p+k,cmpy);
for(ll i=0;i<k;i++) {
for(ll j=i+1;j<k;j++) {
if(tmp_p[j].y-tmp_p[i].y>=dis) break;
dis=min(dis,Dist(tmp_p[i],tmp_p[j]));
}
}
return dis;
}

十一、旋转卡壳

先求凸包。

这个东西可以求凸包直径。

用两条平行线卡住凸包,然后旋转旋转。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
// Rotating calipers.
inline ld Convex_Diameter(Point *p,ll n) {
if(n<=1) return 0;
ll i_=0,j_=0;
for(ll k=1;k<=n;k++) {i_=p[k]<p[i_]?k:i_;j_=p[j_]<p[k]?k:j_;}
ll i=i_,j=j_;
ld ret=Dist(p[i],p[j]);
for(;i!=i_||j!=j_;ret=max(ret,Dist(p[i],p[j]))) {
if(Cross(p[(i+1)%n]-p[i],p[(j+1)%n]-p[j])>=0) {j++;j%=n;}
else {i++;i%=n;}
}
return ret;
}

十二、半平面交

向量逆时针方向可以表示一个半平面。

所以可以直接用一个点表示向量经过这个点,再用一个向量表示方向。

然后极角排序,双端队列扫描半平面,最后能得到半平面交出的多边形。

时间复杂度 $O(n\log n)$。

然后注意一下这个交点不是原来的函数了,因为半平面的表示方法不一样。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
// Half plane intersection.
inline Point Cross_Point(Line v1,Line v2) {
Point u=v1.p1-v2.p1;ld t=Cross(v2.p2,u)/Cross(v1.p2,v2.p2);
return v1.p1+v1.p2*t;
}
inline bool Onleft(Line l,Point p) {return Sgn(Cross(l.p2,p-l.p2))>0;}
inline ll Half_Plan_Intersection(Line *l,Point *ans,ll n) {
sort(l,l+n);ll h=0,t=0;
static Point p[N+5];static Line q[N+5];q[0]=l[0];
for(ll i=1;i<n;i++) {
while(h<t&&!Onleft(l[i],p[t-1])) t--;
while(h<t&&!Onleft(l[i],p[h])) h++;
q[++t]=l[i];
if(fabs(Cross(q[t].p2,q[t-1].p2))<eps) {
t--;if(Onleft(q[t],l[i].p1)) q[t]=l[i];
}
if(h<t) p[t-1]=Cross_Point(q[t-1],q[t]);
}
while(h<t&&!Onleft(q[h],p[t-1])) t--;
if(t-h<=1) return 0;p[t]=Cross_Point(q[t],q[h]);
ll amt=0;
for(ll i=h;i<=t;i++) ans[amt++]=p[i];
return t-h+1;
}

十三、圆

圆的结构体定义一个圆心的位置和一个半径。

点或线与圆的关系很简单。。。求个距离与半径比较就完了。

然后是直线与圆的交点,根据圆心到直线距离得到弦长,然后还能确定方向和直线的单位向量,从而得到了两个点的坐标。

代码:

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
// Circle.
struct Circle{
Point c;ld r;
Circle(){}
Circle(Point _c,ld _r):c(_c),r(_r){}
};
inline ll Point_Circle_Relation(Line v,Circle c) {
ld dst=Dis_Point_Line(c.c,v);
if(Sgn(dst-c.r)<0) return 0;
if(Sgn(dst-c.r)==0) return 1;
return 2;
}
inline ll Line_Circle_Relation(Line v,Circle c) {
ld dst=Dis_Point_Seg(c.c,v);
if(Sgn(dst-c.r)<0) return 0;
if(Sgn(dst-c.r)==0) return 1;
return 2;
}
inline ll Line_Cross_Circle(Line v,Circle c,Point &pa,Point &pb) {
if(Line_Circle_Relation(v,c)==2) return 0;
Point q=Point_Line_Proj(c.c,v);
ld d=Dis_Point_Line(c.c,v);
ld k=sqrt(c.r*c.r-d*d);
if(Sgn(k)==0) {pa=q;pb=q;return 1;}
Point n=(v.p2-v.p1)/Len(v.p2-v.p1);
pa=q+n*k;pb=q-n*k;return 2;
}

十四、最小圆覆盖

不会几何法,所以,模拟退火!!!(但是没有用随机)

但我也不知道没有用随机是否能满足一般的高斯正态分布。。。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
// Minimum Cover Circle.
inline void Min_Cover_Circle(Point *p,ll n,Point &c,ld &r) {
ld T=1000.0,delta=0.998;
c=p[0];ll pos;
while(T>eps) {
pos=0;r=0;
for(ll i=0;i<n;i++) {
if(Dist(c,p[i])>r) {r=Dist(c,p[i]);pos=i;}
}
c.x+=(p[pos].x-c.x)/r*T;
c.y+=(p[pos].y-c.y)/r*T;
T*=delta;
}
}

十五、闵可夫斯基和

求出 $p$,$q$ 凸壳,然后用双指针合并。

简单来说走 $n+m-1$ 步,然后借助斜率选择凸壳上的点。

合并过程

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
// Minkowski Sum.
inline void Minkowski_Sum(Point *f,ll n,Point *g,ll m,Point *q,ll &sz) {
sz=0;ll i=1,j=1;
for(ll k=1;k<=n+m-1;k++) {
q[++sz]=f[i]+g[j];
if(i==n) j++;
else if(j==m) i++;
else {
ld diff=Cross(f[i]+g[j+1]-q[sz],f[i+1]+g[j]-q[sz]);
if(Sgn(diff)<=0) j++;
else i++;
}
}
}

十六、模板集合

终于写完了,这边把几何模板合在一起放在这里。

代码:

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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<vector>
#define ll long long
#define ld long double
using namespace std;

// Some cosntant.
const ld Pi=acos(-1),eps=1e-10;
const ll N=1e5,inf=1e15;

// The Sgn function & Cmp function.
inline ll Sgn(ld x) {return x<-eps?-1:x>eps;}
inline ll Cmp(ld x,ld y) {return Sgn(x-y);}

// Point & Vector.
struct Point{
ld x,y;
inline Point() {}
inline Point(ld x_,ld y_):x(x_),y(y_){}
inline Point operator+(Point p) {return Point(x+p.x,y+p.y);}
inline Point operator-(Point p) {return Point(x-p.x,y-p.y);}
inline Point operator*(ld d) {return Point(x*d,y*d);}
inline Point operator/(ld d) {return Point(x/d,y/d);}
inline bool operator==(Point p) {return !Sgn(x-p.x)&&!Sgn(y-p.y);}
inline bool operator<(Point p) {
return Sgn(x-p.x)<0||(!Sgn(x-p.x)&&Sgn(y-p.y)<0);
}
};
Point p[N+5],tmp_p[N+5];

// To calculate the distance of two Points.
inline ld Dist(Point a,Point b) {
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

// Dot product.
inline ld Dot(Point a,Point b) {return a.x*b.x+a.y+b.y;}
inline ld Len(Point a) {return sqrt(Dot(a,a));}
inline ld Len2(Point a) {return Dot(a,a);}
inline ld Angle(Point a,Point b) {return acos(Dot(a,b)/Len(a)/Len(b));}

// Cross product.
inline ld Cross(Point a,Point b) {return a.x*b.y-a.y*b.x;}
inline ld Area2(Point a,Point b,Point c) {return Cross(b-a,c-a);}
inline Point Rotate(Point a,ld rad) {
return Point(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
inline Point Normal(Point a) {return Point(-a.y/Len(a),a.x/Len(a));}
inline bool Parallel(Point a,Point b) {return !Sgn(Cross(a,b));}

// Line.
struct Line{
Point p1,p2;ld ang;
inline Line() {}
inline Line(Point p1_,Point p2_):p1(p1_),p2(p2_){ang=atan2(p2_.y,p2_.x);}
inline Line(Point p,double angle) {
p1=p;if(Sgn(angle-Pi/2)==0) {p2=(p1+Point(0,1));}
else {p2=(p1+Point(1,tan(angle)));}
}
inline Line(ld a,ld b,ld c) {
if(!Sgn(a)) {p1=Point(0,-c/b);p2=Point(1,-c/b);}
else if(!Sgn(b)) {p1=Point(-c/a,0);p2=Point(-c/a,1);}
else {p1=Point(0,-c/b);p2=Point(1,(-c-a)/b);}
}
inline bool operator<(Line &l) {return ang<l.ang;}
};
inline ll Point_Line_Relation(Point p,Line v) {
ll c=Sgn(Cross(p-v.p1,v.p2-v.p1));
if(c<0) return 1;if(c>0) return 2;return 0;
}
inline bool Point_on_Seg(Point p,Line v) {
return !Sgn(Cross(p-v.p1,v.p2-v.p1))&&Sgn(Dot(p-v.p1,p-v.p2))<=0;
}
inline ld Dis_Point_Line(Point p,Line v) {
return fabs(Cross(p-v.p1,v.p2-v.p1))/Dist(v.p1,v.p2);
}
inline Point Point_Line_Proj(Point p,Line v) {
ld k=Dot(v.p2-v.p1,p-v.p1)/Len2(v.p2-v.p1);
return v.p1+(v.p1-v.p2)*k;
}
inline Point Point_Line_Symmetry(Point p,Line v) {
Point q=Point_Line_Proj(p,v);
return Point(2*q.x-p.x,2*q.y-p.y);
}
inline ld Dis_Point_Seg(Point p,Line v) {
if(Sgn(Dot(p-v.p1,v.p2-v.p1))<0||Sgn(Dot(p-v.p2,v.p1-v.p2))<0) {
return min(Dist(p,v.p1),Dist(p,v.p2));
}
return Dis_Point_Line(p,v);
}
inline ll Line_Relation(Line v1,Line v2) {
if(Sgn(Cross(v1.p2-v1.p1,v2.p2-v2.p1))==0) {
if(Point_Line_Relation(v1.p1,v2)==0) return 1;
else return 0;
}
return 2;
}
inline Point Cross_Point(Line v1,Line v2) {
ld s1=Cross(v1.p2-v1.p1,v2.p1-v1.p1),s2=Cross(v1.p2-v1.p1,v2.p2-v1.p1);
return Point(v2.p1.x*s2-v2.p2.x*s1,v2.p1.y*s2-v2.p2.y*s1)/(s2-s1);
}
inline bool Cross_Seg(Line v1,Line v2) {
ld c1=Cross(v1.p2-v1.p1,v2.p1-v1.p1),c2=Cross(v1.p2-v1.p1,v2.p2-v1.p1);
ld d1=Cross(v2.p2-v2.p1,v1.p1-v2.p1),d2=Cross(v2.p2-v2.p1,v1.p2-v2.p1);
return Sgn(c1)*Sgn(c2)<=0&&Sgn(d1)*Sgn(d2)<=0;
}

// Polygon.
inline ll Point_in_Polygon(Point pt,Point *p,ll n) {
for(ll i=0;i<n;i++) {if(p[i]==pt) return 3;}
for(ll i=0;i<n;i++) {
Line v=Line(p[i],p[(i+1)%n]);
if(Point_on_Seg(pt,v)) return 2;
}
ll num=0;
for(ll i=0;i<n;i++) {
ll j=(i+1)%n;
ll c=Sgn(Cross(pt-p[j],p[i]-p[j]));
ll u=Sgn(p[i].y-pt.y),v=Sgn(p[j].y-pt.y);
if(c>0&&u<0&&v>=0) num++;
if(c<0&&u>=0&&v<0) num--;
}
return num!=0;
}
inline ld Polygon_Area(Point *p,ll n) {
ld area=0;for(ll i=0;i<n;i++) area+=Cross(p[i],p[(i+1)%n]);return area/2;
}
inline Point Polygon_Centre(Point *p,ll n) {
Point ans(0,0);if(Polygon_Area(p,n)==0) return ans;
for(ll i=0;i<n;i++) ans=ans+(p[i]+p[(i+1)%n])*Cross(p[i],p[(i+1)%n]);
return ans/Polygon_Area(p,n)/6;
}

// Convex hull.
inline ll Convex_Hull(Point *p,ll n,Point *ch) {
sort(p,p+n);n=unique(p,p+n)-p;ll v=0;
for(ll i=0;i<n;i++) {
while(v>1&&Sgn(Cross(ch[v-1]-ch[v-2],p[i]-ch[v-2]))<=0) v--;
ch[v++]=p[i];
}
ll j=v;
for(ll i=n-2;i>=0;i--) {
while(v>j&&Sgn(Cross(ch[v-1]-ch[v-2],p[i]-ch[v-2]))<=0) v--;
ch[v++]=p[i];
}
if(n>1) v--;
return v;
}

// Closest pair.
inline bool cmpy(Point a,Point b) {return Sgn(a.y-b.y)<0;}
inline ld Closest_Pair(ll left,ll right) {
ld dis=inf;if(left==right) return dis;
if(left+1==right) return Dist(p[left],p[right]);
ll mid=(left+right)/2;
ld d1=Closest_Pair(left,mid),d2=Closest_Pair(mid+1,right);
dis=min(d1,d2);
ll k=0;
for(ll i=left;i<=right;i++) {
if(fabs(p[mid].x-p[i].x)<=dis) tmp_p[k++]=p[i];
}
sort(tmp_p,tmp_p+k,cmpy);
for(ll i=0;i<k;i++) {
for(ll j=i+1;j<k;j++) {
if(tmp_p[j].y-tmp_p[i].y>=dis) break;
dis=min(dis,Dist(tmp_p[i],tmp_p[j]));
}
}
return dis;
}

// Rotating calipers.
inline ld Convex_Diameter(Point *p,ll n) {
if(n<=1) return 0;
ll i_=0,j_=0;
for(ll k=1;k<=n;k++) {i_=p[k]<p[i_]?k:i_;j_=p[j_]<p[k]?k:j_;}
ll i=i_,j=j_;
ld ret=Dist(p[i],p[j]);
for(;i!=i_||j!=j_;ret=max(ret,Dist(p[i],p[j]))) {
if(Cross(p[(i+1)%n]-p[i],p[(j+1)%n]-p[j])>=0) {j++;j%=n;}
else {i++;i%=n;}
}
return ret;
}

// Half plane intersection.
// inline Point Cross_Point(Line v1,Line v2) {
// Point u=v1.p1-v2.p1;ld t=Cross(v2.p2,u)/Cross(v1.p2,v2.p2);
// return v1.p1+v1.p2*t;
// }
inline bool Onleft(Line l,Point p) {return Sgn(Cross(l.p2,p-l.p2))>0;}
inline ll Half_Plan_Intersection(Line *l,Point *ans,ll n) {
sort(l,l+n);ll h=0,t=0;
static Point p[N+5];static Line q[N+5];q[0]=l[0];
for(ll i=1;i<n;i++) {
while(h<t&&!Onleft(l[i],p[t-1])) t--;
while(h<t&&!Onleft(l[i],p[h])) h++;
q[++t]=l[i];
if(fabs(Cross(q[t].p2,q[t-1].p2))<eps) {
t--;if(Onleft(q[t],l[i].p1)) q[t]=l[i];
}
if(h<t) p[t-1]=Cross_Point(q[t-1],q[t]);
}
while(h<t&&!Onleft(q[h],p[t-1])) t--;
if(t-h<=1) return 0;p[t]=Cross_Point(q[t],q[h]);
ll amt=0;
for(ll i=h;i<=t;i++) ans[amt++]=p[i];
return t-h+1;
}

// Circle.
struct Circle{
Point c;ld r;
Circle(){}
Circle(Point _c,ld _r):c(_c),r(_r){}
};
inline ll Point_Circle_Relation(Line v,Circle c) {
ld dst=Dis_Point_Line(c.c,v);
if(Sgn(dst-c.r)<0) return 0;
if(Sgn(dst-c.r)==0) return 1;
return 2;
}
inline ll Line_Circle_Relation(Line v,Circle c) {
ld dst=Dis_Point_Seg(c.c,v);
if(Sgn(dst-c.r)<0) return 0;
if(Sgn(dst-c.r)==0) return 1;
return 2;
}
inline ll Line_Cross_Circle(Line v,Circle c,Point &pa,Point &pb) {
if(Line_Circle_Relation(v,c)==2) return 0;
Point q=Point_Line_Proj(c.c,v);
ld d=Dis_Point_Line(c.c,v);
ld k=sqrt(c.r*c.r-d*d);
if(Sgn(k)==0) {pa=q;pb=q;return 1;}
Point n=(v.p2-v.p1)/Len(v.p2-v.p1);
pa=q+n*k;pb=q-n*k;return 2;
}

// Minimum Cover Circle.
inline void Min_Cover_Circle(Point *p,ll n,Point &c,ld &r) {
ld T=1000.0,delta=0.998;
c=p[0];ll pos;
while(T>eps) {
pos=0;r=0;
for(ll i=0;i<n;i++) {
if(Dist(c,p[i])>r) {r=Dist(c,p[i]);pos=i;}
}
c.x+=(p[pos].x-c.x)/r*T;
c.y+=(p[pos].y-c.y)/r*T;
T*=delta;
}
}

// Minkowski Sum.
inline void Minkowski_Sum(Point *f,ll n,Point *g,ll m,Point *q,ll &sz) {
sz=0;ll i=1,j=1;
for(ll k=1;k<=n+m-1;k++) {
q[++sz]=f[i]+g[j];
if(i==n) j++;
else if(j==m) i++;
else {
ld diff=Cross(f[i]+g[j+1]-q[sz],f[i+1]+g[j]-q[sz]);
if(Sgn(diff)<=0) j++;
else i++;
}
}
}

int main() {
return 0;
}