文章目录

前言
精度
点/向量相关
- 表示
- 向量基本运算
- 角度相关
- - 向量夹角
  - 旋转
直线/线段相关
- 表示
- 点与线
- - 求点到直线垂足
  - 求点关于直线的对称点
  - 点与直线的位置关系
  - 点与直线的距离
- 线与线
- - 直线与直线的位置关系
  - - 共线与垂直
    - 判断线段与线段是否相交
    - 求直线与直线的交点
  - 角平分线
  - 中垂线
多边形
- 表示
- 求多边形面积
- 判断多边形凹凸性
- 判定点在多边形内
- - 判定点在凸包内
圆
- 表示
- 位置关系
- - 点与圆的位置关系
  - 线与圆的位置关系
  - 圆与圆的位置关系
- 三角形和圆
- - 三角形内切圆
  - 三角形外接圆
- 交点相关
- - 圆与直线交点
  - 圆与圆的交点
- 切点/线相关
- - 点到圆的切点
  - 圆与圆公切线的切点
- 面积相关
- - 圆、扇形和弓形
  - 圆与多边形面积交
  - 圆与圆的面积交
二维凸包
- 定义
- 维护
- 代码
旋转卡壳
- 代码
半平面交
- 半平面的表示
- 流程
- 细节
- 代码
其他
- 皮克定理
- 距离
- - 曼哈顿距离和切比雪夫距离的转化
source

前言

强烈推荐AOJ的入门几何题库。
本文是一篇入门全解，需要的前置知识只有基础的高中的解析几何和向量相关。
计算几何的一些原则：

精度即生命。炸精度是一个糟糕的计算几何实现经常出现的问题（就像我一开始尝试自己造轮子常常遇到的），为了避免过大的精度误差，应该尽量减少正反三角函数、除法、根号的使用次数（但似乎根号函数的精度还是不错的）。（实在不行，long double 或许能救我们一命，但是会变慢）
简洁即美德。这一条是建立在不过分影响精度的前提之上的。很多时候对于某个问题，存在多种情况，但对于一种较简单情况给出的某种解法归纳发现在其他情况下也是正确的，就可以大大减少分类讨论（比如线段求交问题）。此外，代码合理的模块化也能使代码变得更加简洁。

所以，尽管我们自己往往也能乱搞出对于某个问题的一种做法，但是为了精度更高、更加简洁的实现，学习计算几何是必要的。
后面的函数里可能会直接调用前面讲完的一些函数，这也是代码模块化简化代码的一个体现。

注意：任何时候用除法都要考虑考虑会不会除0！

精度

计算几何常用浮点数，由于浮点数本身难以避免、只能减小的精度误差，常常不使用通常的，=、>、< 运算。
设 $e p s$ 为一个极小量。（通常在 $10^{-12},10^{-8}]$ 之间）
那么三种运算符再浮点运算下就可以表示为：

$\to\operatorname{abs}(a-b)<eps$
$\to a<b-eps$
$\to a>b+eps$

这不是死记硬背的。其意义就是假设当 $a = b$ 时，表达式不应该成立；那么类似的，<= 和 >= 应该在 $a = b$ 时成立，所以应该把对应的 $e p s$ 之前的符号反过来。

点/向量相关

表示

由于在计算几何中点和向量常常是等价的，所以归为一类。
只需要记录横纵坐标即可。

//definition
struct V{double x,y;V():x(0),y(0){}V(double a,double b):x(a),y(b){}
};
V ans[10];//declared for other functions
int tot;
inline void input(V &a){a.x=read();a.y=read();}
void print(const V &a,int op=1){printf("%.10lf %.10lf",a.x,a.y);putchar(op?10:32);}
//op:endl or space

向量基本运算

基本的加、减、乘、除、点积、叉积、模长、中点、垂直向量、向量单位化、三角形面积。

拓展：
高维点积：定义不变，但在投影方面几何意义只对二维和三维空间有效。
高维叉积：叉积在k维下的定义是三行k列的行列式，第一行为各维度的单位向量，二三行为两个点各自在对应维度的坐标。

注意：

严谨来说，向量叉积是一个向量而不是一个数，它的模长等于 $x_1y_2-x_2y_1$ ，方向垂直与纸面，遵循右手定则，这也是叉积可以表示负面积的原因。
这里的垂直向量没有关注方向。

//basic operation
inline V operator + (const V &a,const V &b){return (V){a.x+b.x,a.y+b.y};}
inline V operator - (const V &a,const V &b){return (V){a.x-b.x,a.y-b.y};}
inline V operator * (const double &x,const V &a){return (V){a.x*x,a.y*x};}
inline V operator * (const V &a,const double &x){return (V){a.x*x,a.y*x};}
inline V operator / (const V &a,const double x){return (V){a.x/x,a.y/x};}
inline bool operator == (const V &a,const V &b){return abs(a.x-b.x)<eps&&abs(a.y-b.y)<eps;}
inline bool operator != (const V &a,const V &b){return !(a==b);}
inline double operator * (const V &a,const V &b){return a.x*b.x+a.y*b.y;}
inline double operator ^ (const V &a,const V &b){return a.x*b.y-a.y*b.x;}
inline double len(const V &a){return sqrt(a.x*a.x+a.y*a.y);}
inline V mid(const V &a,const V &b){return (V){(a.x+b.x)/2,(a.y+b.y)/2};}
inline V cui(const V &a){return (V){a.y,-a.x};}//not take direction into account
inline V danwei(const V &a){return a/len(a);}
inline double tri_S(const V &a,const V &b,const V &c){return abs((b-a)^(c-a))/2;}//always be non-negative
inline bool operator < (const V &a,const V &b){return a.x<b.x-eps||(abs(a.x-b.x)<eps&&a.y<b.y-eps);
}

角度相关

向量夹角

利用点积易得。
注意这个角的值域是0到π。

inline double angle(const V &a,const V &b) { return acos(a * b / len(a) / len(b)); }
inline bool dun(const V &a,const V &b,const V &c){return ((b-a)*(c-a))<-eps;}//angle:BAC
inline bool rui(const V &a,const V &b,const V &c){return ((b-a)*(c-a))>eps;}
inline bool zhi(const V &a,const V &b,const V &c){return abs(((b-a)*(c-a)))<eps;}

旋转

若原来向量为 $(x, y)$ ，旋转角为 $θ\theta$ ，那么旋转后的向量就是 $(xcos⁡θ−ysin⁡θ,xsin⁡θ+ycos⁡θ)(x\cos\theta-y\sin\theta,x\sin\theta+y\cos\theta)$ 。
可以用虚数乘法， $(x+yi)×(cos⁡θ,sin⁡θi)(x+yi)\times (\cos\theta,\sin\theta \space i)$ ，结合几何意义简单证明。

inline V rotate(const V &o,double t){double s=sin(t),c=cos(t);return (V){o.x*c-o.y*s,o.x*s+o.y*c};
}

直线/线段相关

表示

直线常用方向向量加上线上任意一点表示。
线段常用两个端点表示。
把方向向量转化为单位向量会比较方便（下面直线的方向向量默认为单位向量）。
（下文如非特别说明，在大多数时候不在特别区分线段和直线，两者的算法是大同小异的。）

//definition
struct line{V d,a,b;
};
inline line trans(double a,double b,double c,double d){//given coordinatesV dd(c-a,d-b),x(a,b),y(c,d);dd=dd/len(dd);return (line){dd,x,y};
}
inline line trans(const V &a,const V &b){//given pointsV dd(b.x-a.x,b.y-a.y);dd=dd/len(dd);return (line){dd,a,b};
}
inline void input(line &l){double a=read(),b=read(),c=read(),d=read();l=trans(a,b,c,d);return;
}

点与线

求点到直线垂足

把线上的代表点移动投影的距离。
这里的投影是带符号的，符号体现移动的方向。

inline V cui(const line &l,const V &o){//directedreturn l.a+((o-l.a)*l.d)*l.d;
}

求点关于直线的对称点

求出垂足后用中点公式即可。

inline V duichen(const line &l,const V &o){V u=cui(l,o);return (V){2*u.x-o.x,2*u.y-o.y};
}

点与直线的位置关系

共线可以通过叉积等于0来判断。
不共线时，可以通过叉积的正负判断是在直线的哪一侧。
共线时，如果是线段或射线，还可以利用点积的正负判断与端点的位置关系。
（对于线段而言，可以直接通过模长判断是否在线段上）

inline bool on_line(const line &l,const V &o){return abs((l.d^(o-l.a)))<eps;}
inline bool on_seg(const line &l,const V &o){return abs(len(o-l.a)+len(o-l.b)-len(l.a-l.b))<eps;
}
inline int pos(const line &l,const V &o){if(!on_line(l,o)){if((l.d^(o-l.a))>eps) return 1;//counter clockwiseelse return 2;//clockwise}else{if(((o-l.a)*(o-l.b))<eps) return 5;//on segmentelse if(((o-l.a)*l.d)<-eps) return 3;//online backelse return 4;//online front}
}

点与直线的距离

直线时，可以利用用叉积等面积法求出距离。
线段时，若垂足在线段上，还是点到直线的距离；否则是到两端点距离的最小值。
垂足是否在线段上可以利用是否有钝角来判断。

inline double dis(const line &l,const V &o,int op=0){//op=0:dis on line,op=1:dis on segmentif(op&&(dun(l.a,o,l.b)||dun(l.b,o,l.a))) return min(len(o-l.a),len(o-l.b));else return abs((o-l.a)^(o-l.b))/len(l.a-l.b);
}

线与线

直线与直线的位置关系

共线与垂直

如果方向向量叉积是0则共线（注意还要分平行和重合），如果点积是0则垂直。

inline bool gongxian(const line &a,const line &b){return abs(a.d^b.d)<eps;}
inline bool cuizhi(const line &a,const line &b){return abs(a.d*b.d)<eps;}

判断线段与线段是否相交

需要两个步骤：快速排斥实验和跨立实验。

快速排斥实验：定义一条线段所在的矩形为边与坐标轴平行、能完全包含该线段的最小的矩形，如果两条线段所在的矩形没有公共部分，显然不可能有交点；否则，进入下一步跨立实验。
跨立实验：如果l1的两端点在l2两侧，同时l2的两端点也在l1两侧（也可以在线段上），则得出两线段有公共交点，否则则没有公共交点。

实现上，快速排斥实验可以用坐标 min、max 的比较轻松解决，跨立实验则可以用前面讲的叉积求点与线的关系来实现，注意叉积为0（端点在线段上）是符合条件的。

直观上感觉似乎跨立实验已经很充要了，但是由于判断的时候叉积等于0认为合法，会把没有公共部分的共线线段判为有交点，所以快速排斥实验是必要的。（也可以用特判共线来取代）

inline bool jiao(const line &u,const line &v){if(min(u.a.x,u.b.x)>max(v.a.x,v.b.x)+eps||max(u.a.x,u.b.x)<min(v.a.x,v.b.x)-eps||min(u.a.y,u.b.y)>max(v.a.y,v.b.y)+eps||max(u.a.y,u.b.y)<min(v.a.y,v.b.y)-eps) return false;//rapid rejection testreturn ((u.a-v.a)^v.d)*((u.b-v.a)^v.d)<eps&&((v.a-u.a)^u.d)*((v.b-u.a)^u.d)<eps;//straddle test
}

求直线与直线的交点

首先应保证两直线不共线。
设两直线为 $l 1, l 2$ ，线上各有一点 $x 1, x 2$ ，方向向量为 $d 1, d 2$ 。
假设交点 $p = x 1 + k d 1$ ，就有 $(p−x2)×d2=0(p-x2)\times d2=0$ 。
叉积是满足分配率的，拆开化简得到：
$k=(x2−x1)×d2d1×d2k=\frac{(x2-x1)\times d2}{d1\times d2}$
带回去即可。

inline V jiaodian(line u,line v){double k=((v.a-u.a)^v.d)/(u.d^v.d);return u.a+(k*u.d);
}

角平分线

求出两边的单位向量d1,d2，d1+d2就是角平分线的方向向量。

inline line pingfen(const V &a,const V &b,const V &c){//angle:BACV d1=(b-a)/len(b-a),d2=(c-a)/len(c-a),d=(d1+d2)/len(d1+d2);return (line){d,a,a+d};
}

中垂线

找到中点和与原线段垂直的方向向量即可。

inline line zhongcui(const V &a,const V &b){V x=mid(a,b),d=danwei(cui(b-a));return (line){d,x,x+d};
}

多边形

表示

用一个V数组顺时针或逆时针存储多边形的所有顶点即可。
由于封装的话老打点会比较烦，就直接传数组和顶点数n了。

求多边形面积

取任意一个点 $O$ ，则多边形面积为：
$S=12∑i=1nOAi→×OA→i%n+1S=\frac{1}{2}\sum_{i=1}^n\overrightarrow{OA_{i}}\times \overrightarrow{OA}_{i\%n+1}$
为了方便， $O$ 直接取原点即可。
时间复杂度 $O (n)$ 。

inline double S(const V  *a,const int n){double res(0);for(int i=1;i<=n;i++) res+=(a[i]^a[i%n+1]);return res/2;
}

判断多边形凹凸性

取相邻的三个点 $(a, b, c)$ ，作为一个三元组。
判断多边形凹凸的充要条件：对于多边形所有的 $n$ 个这样三元组， $c$ 和 $a b$ 的位置关系（指顺时针或逆时针）必然全部都是相同的。
叉积判断符号即可。
注意如果没有保证相邻三个点不共线的话需要特判。
时间复杂度 $O (n)$ 。

bool is_convex(const V *a,const int n){a[0]=a[n];a[n+1]=a[1];int op=0;for(int i=1;i<=n;i++){double o=(a[i+1]-a[i])^(a[i]-a[i-1]);if(abs(o)<eps) continue;//three neighbouring points on the same lineint nop=o>0?1:-1;if(!op) op=nop;else if(op!=nop) return false;}return true;
}

判定点在多边形内

一个看起来很好写但写起来很恶心的东西。
推荐使用点数判别法。
向左水平引出一条射线，计算其与多边形产生多少个交点，奇内偶外，然后还有亿点点细节：

如果当前边经过了该点，直接返回相应结果
由于可能出现射线经过多边形顶点导致一个点被计算两次的情况，强制令这个点在下面的边被统计。实现上，当 $max⁡(y1,y2)=y0\max(y1,y2)=y0$ 时，统计该次答案，而当 $min⁡(y1,y2)=y0\min(y1,y2)=y0$ 时，不统计该次答案。
出现水平边时，直接忽略即可，结合第二条的处理原则，就依然是正确的。

int in_poly(const V *a,const int n,const V o){int res(0);for(int i=1;i<=n;i++){V u=a[i],v=a[i+1];if(on_seg(trans(u,v),o)) return 1;//on the edgeif(abs(u.y-v.y)<eps) continue;//ignore horizontalif(max(u.y,v.y)<o.y-eps) continue;//equal will not continueif(min(u.y,v.y)>o.y-eps) continue;//equal will continuedouble x=u.x+(o.y-u.y)/(v.y-u.y)*(v.x-u.x);if(x<o.x) res^=1;}return res?2:0;//odd:in even:out
}

判定点在凸包内

首先，需要保证 $a_1=(0,0)$ 。（如果不满足整体平移即可）
lowerbound 找到极角夹在询问点两侧的相邻点对，判断是否在二者连成的边内部即可。
注意需要对于在 $a_1,a_n$ 外侧的点需要特判。

bool cmp2(V a,V b){return (a^b)>0;
}
bool in(const V *c,const V &o){if(((c[n]^o)>eps)||((o^c[2])>eps)) return 0;int pl=lower_bound(c+1,c+1+n,o,cmp2)-c-1;return ((c[pl+1]-c[pl])^(o-c[pl]))>-eps;
}

圆

表示

记录圆心和半径即可表示。

struct cir{V o;double r;cir(double a=0,double b=0,double c=0):o(a,b),r(c){}
};
inline void input(cir &cc){cc.o.x=read();cc.o.y=read();cc.r=read();}

位置关系

点与圆的位置关系

判断距离即可。

inline bool incir(const cir &c,const V &p){return len(p-c.o)<c.r-eps;}
inline bool oncir(const cir &c,const V &p){return (len(p-c.o)-c.r)<eps;}
inline bool outcir(const cir &c,const V &p){return len(p-c.o)>c.r+eps;}

线与圆的位置关系

按照线圆距和半径的关系判断即可。

inline double dis(const cir &c,const line &l){return dis(l,c.o);}
inline int pos(const cir &c,const line &l){double d=dis(c,l);if(d<c.r-eps) return 1;//intersectelse if(abs(d-c.r)<eps) return 2;//tangentelse if(d>c.r+eps) return 3;//disjoint
}

圆与圆的位置关系

分为外离、外切、相交、内切、内含五种。
公切线的数量对应为4、3、2、1、0。
通过圆心距和半径之间的关系判断即可。

inline double dis(const cir &c1,const cir &c2){return len(c1.o-c2.o);}
inline int pos(const cir &c1,const cir &c2){double d=dis(c1,c2);if(d>c1.r+c2.r+eps) return 4;//do not crosselse if(abs(d-c1.r-c2.r)<eps) return 3;//circumscribedelse if(max(c1.r,c2.r)<min(c1.r,c2.r)+d -eps) return 2;//intersectelse if(abs(max(c1.r,c2.r)-min(c1.r,c2.r)-d)<eps) return 1;//inscribedelse return 0;//include
}

三角形和圆

三角形内切圆

圆心就是两条平分线的交点即可。
半径可以用公式 $2Sa+b+c\dfrac{2S}{a+b+c}$ ，也可以调用点到直线距离。

inline cir incir(const V &a,const V &b,const V &c){V x(jiaodian(pingfen(a,b,c),pingfen(b,a,c)));return cir(x,dis(trans(a,b),x));
}

三角形外接圆

我使用的是直接求中垂线的交点。
但是这么做有点炸精度…开longdouble才能过。

inline cir outcir(const V &a,const V &b,const V &c){V x(jiaodian(zhongcui(a,b),zhongcui(a,c)));return cir(x,len(x-a));
}
inline void input(cir &cc){in(cc.o);cc.r=read();}

法二：
设三边边长为 $a, b, c$ ，则 $R=abc(a+b+c)(a+b−c)(a+c−b)(b+c−a)R=\dfrac{abc}{\sqrt{(a+b+c)(a+b-c)(a+c-b)(b+c-a)}}$ 。

交点相关

圆与直线交点

求出圆心到直线距离后勾股定理即可。

inline double calc(double xie,double zhi){return sqrt(xie*xie-zhi*zhi);}//Pythagorean Theorem
inline void cross_line_cir(const cir &c,const line &l){tot=0;V x=chui(l,c.o);double dd=len(x-c.o);if(c.r<dd-eps) return;//none cross pointif(abs(c.r-dd)<eps){//one cross pointans[++tot]=x;return;}double dis=calc(c.r,dd);V a=x+dis*l.d,b=x-dis*l.d;//two cross pointsans[++tot]=a;ans[++tot]=b;return;
}

圆与圆的交点

首先得保证存在交点。
考虑两圆心和任意一个交点组成的三角形，可以发现这个三角形的三边都是已知的。
我们就可以用余弦定理求出夹角，然后转上去即可。

inline void cross_cir_cir(const cir &c1,const cir c2){int op=pos(c1,c2);if(op==4||op==0) return;//none cross pointV L=c2.o-c1.o;double a=c2.r,b=c1.r,c=len(L);double t=acos((b*b+c*c-a*a)/(2*b*c));//find the angleV x=c1.o+c1.r*rotate((L)/len(L),t),y=duichen(trans(c1.o,c2.o),x);if(x<y) print(x,0),print(y,1);else print(y,0),print(x,1);
}

切点/线相关

点到圆的切点

切点、圆心、给定点自然而然形成了一个直角三角形。
求出角后转一下即可。

inline void qiedian(const cir &c,const V &p){tot=0;line L=trans(c.o,p);double t=acos(c.r/len(c.o-p));V a=c.o+(c.r*rotate(L.d,t)),b=duichen(L,a);ans[++tot]=a;if(a!=b) ans[++tot]=b;return;
}

圆与圆公切线的切点

分为外侧切线和内测切线。
两者大致的思路类似，都是试图求出切点圆心连线与两圆心连线的夹角，然后转上去。
符号的正负很和谐，不必特殊考虑。

inline void common_qie(const cir &c1,const cir &c2){tot=0;int op=pos(c1,c2);if(op){//outside tangent linesdouble d=c1.r-c2.r,t=acos(d/len(c1.o-c2.o));V u=c1.o+(c1.r*danwei(rotate(c2.o-c1.o,t))),v=c1.o+(c1.r*danwei(rotate(c2.o-c1.o,-t)));ans[++tot]=u;if(u!=v) ans[++tot]=v;   }if(op>2){//inside tangent linesdouble d=len(c2.o-c1.o)/(c1.r+c2.r)*c1.r,t=acos(c1.r/d);V u=c1.o+(c1.r*danwei(rotate(c2.o-c1.o,t))),v=c1.o+(c1.r*danwei(rotate(c2.o-c1.o,-t)));ans[++tot]=u;if(u!=v) ans[++tot]=v;}
}

面积相关

圆、扇形和弓形

比较基础，套初中公式即可。

inline double cir_S(const cir &c){return pi*c.r*c.r;}
inline double shan(const cir &c,const V &a,const V &b){//S of sectordouble t=ang(a-c.o,b-c.o);return t/2*c.r*c.r;
}
inline double gong(const cir &c,const V &a,const V &b){//S of bowreturn shan(c,a,b)-tri_S(c.o,a,b);
}

圆与多边形面积交

和求多边形面积类似的，问题可以转化为求多边形 n 对相邻顶点与圆心组成的三角形与圆的有向面积交。
然后需要亿点点大力分讨…（具体可见代码和注释，不太难理解）
一个实现的 trick 是先把有向面积的符号求出来，后面求面积的绝对值就行了。

inline double O_tri(const cir &c,V a,V b){
//directed S of Intersection between circle O and triangle OABif(on(trans(a,b),c.o)) return 0.0;int f=(((a-c.o)^(b-c.o))>0)?1:-1;//directionline l=trans(a,b);int f1=incir(c,a),f2=incir(c,b);if(f1&&f2)  return f*tri_S(a,b,c.o);//both incircle:the S of triangleelse if(!f1&&!f2){//both outcircle:a sector(possibly minus a bow)V u=(c.o+(c.r*danwei(a-c.o))),v=(c.o+(c.r*danwei(b-c.o)));double S=shan(c,u,v);    if(dis(l,c.o,1)<c.r-eps){cross_line_cir(c,l);assert(tot==2);S-=gong(c,ans[1],ans[2]);}return f*S;}else{//one in and one out:a traingle and a sectorif(!f1) swap(a,b);double sa=Sin(b-a,c.o-a),su=sa/c.r*len(c.o-a),A=asin(sa),U=asin(su);if((b-a)*(c.o-a)<-eps) A=pi-A;double t=pi-A-U,dis=sin(t)/sa*c.r;    V u=a+(dis*danwei(b-a)),v=c.r*danwei(b-c.o);return f*(tri_S(c.o,a,u)+shan(c,u,v));}
}
double common_cir_poly(const V *a,const cir &c,int n){double res(0);for(int i=1;i<=n;i++) res+=O_tri(c,a[i],a[i+1]);return res;
}

圆与圆的面积交

首先把较为简单的外离和内含的情况判掉，重点考虑相交。
先考虑比较简单的情形。交点与两圆心连线夹角都是锐角时，相交部分是两个弓形，利用余弦定理求出相交部分的圆心角，然后用扇形面积减三角形面积即可。
然后发现如果变成了钝角，三角形面积会变成一个负面积，减完等于加上三角形面积的绝对值，结果还是对的，十分和谐，不需要特判。

inline double common_cir_cir(const cir &c1,const cir &c2){int op=pos(c1,c2);if(op>=3) return 0;//common area=0else if(op<=1) return min(cir_S(c1),cir_S(c2));//completly includeelse{//partly includedouble L=len(c1.o-c2.o);double t1=2*acos((L*L+c1.r*c1.r-c2.r*c2.r)/(2*L*c1.r));double t2=2*acos((L*L+c2.r*c2.r-c1.r*c1.r)/(2*L*c2.r));return c1.r*c1.r*t1/2-c1.r*c1.r*sin(t1)/2+c2.r*c2.r*t2/2-c2.r*c2.r*sin(t2)/2;}
}

二维凸包

解析集合的第一课。
关键特征：周长最小。此时一定是凸包。

定义

凸包：在平面上能包含所有给定点的最小凸多边形叫做凸包。

性质：凸包的周长是所有能包含给定点的多边形中最小的。

维护

凸包的主流求法有 Andrew 和 Graham，两者的复杂度瓶颈都在于排序，为 $O(nlog⁡n)O(n\log n)$ 。
这里介绍较为简单的 Graham 扫描法。

首先，找出所有点中按照 $X - Y$ 排序最小的点 $O$ 。
然后以 $O$ 作为原点，把所有其它点按照极角排序，极角相同时按距离升序排序（实测这里也可以降序，但是绝对不能无序）。
可以用快排重载 cmp 函数的方法实现，利用叉积判断。

bool cmp(V a,V b){double d=(a-p[1])^(b-p[1]);if(abs(d)>eps) return d>0;else return len(a-p[1])<len(b-p[1]);
}

排好序后，开一个栈维护当前凸包中的点。
如果新点破坏了凸性，则不断退栈。
最终栈内的元素就是凸包中的点。

是否破坏凸性可以用叉积判断，如果新点和栈顶形成的向量向右拐了（叉积小于0），则说明破坏了凸性。

注意：这里判断退栈的条件最好带等！，一方面，共线时在边上的顶点没有什么意义，另一方面，当有重合点时，不带等会导致程序错误。

void ConvexHull(V *p,int &n){int top=0;sort(p+1,p+1+n);sort(p+2,p+1+n,cmp);top=0;for(int i=1;i<=n;i++){while((top>1&&((zhan[top]-zhan[top-1])^(p[i]-zhan[top]))<=0)) --top;zhan[++top]=p[i];}memcpy(p,zhan,sizeof(zhan));n=top;return;
}

代码

P2742 [USACO5.1]圈奶牛Fencing the Cows /【模板】二维凸包

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define debug(...) fprintf(stderr,__VA_ARGS__)
inline ll read(){ll x(0),f(1);char c=getchar();while(!isdigit(c)){if(c=='-')f=-1;c=getchar();}while(isdigit(c)){x=(x<<1)+(x<<3)+c-'0';c=getchar();}return x*f;
}//basic declare
//#define double long double
const double eps=1e-10;
const double pi=acos(-1.0);//---about vectors (or points)//definition
struct V{double x,y;V():x(0),y(0){}V(double a,double b):x(a),y(b){}
};
V ans[10];//declared for other functions
int tot;
inline void input(V &a){scanf("%lf%lf",&a.x,&a.y);}
void print(const V &a,int op=1){printf("%.10lf %.10lf",a.x,a.y);putchar(op?10:32);}
//op:endl or space//basic operation
inline V operator + (const V &a,const V &b){return (V){a.x+b.x,a.y+b.y};}
inline V operator - (const V &a,const V &b){return (V){a.x-b.x,a.y-b.y};}
inline V operator * (const double &x,const V &a){return (V){a.x*x,a.y*x};}
inline V operator * (const V &a,const double &x){return (V){a.x*x,a.y*x};}
inline V operator / (const V &a,const double x){return (V){a.x/x,a.y/x};}
inline bool operator == (const V &a,const V &b){return abs(a.x-b.x)<eps&&abs(a.y-b.y)<eps;}
inline bool operator != (const V &a,const V &b){return !(a==b);}
inline double operator * (const V &a,const V &b){return a.x*b.x+a.y*b.y;}
inline double operator ^ (const V &a,const V &b){return a.x*b.y-a.y*b.x;}
inline double len(const V &a){return sqrt(a.x*a.x+a.y*a.y);}
inline V mid(const V &a,const V &b){return (V){(a.x+b.x)/2,(a.y+b.y)/2};}
inline V chui(const V &a){return (V){a.y,-a.x};}//not take direction into account
inline V danwei(const V &a){return a/len(a);}
inline double tri_S(const V &a,const V &b,const V &c){return abs((b-a)^(c-a))/2;}//always be non-negative
inline bool operator < (const V &a,const V &b){return a.x<b.x-eps||(abs(a.x-b.x)<eps&&a.y<b.y-eps);
}
inline double ang(const V &a,const V &b){return acos((a*b)/len(a)/len(b));}
inline V rotate(const V &o,double t){//COUNTER_CLOCKWISEdouble s=sin(t),c=cos(t);return (V){o.x*c-o.y*s,o.x*s+o.y*c};
}const int N=1e5+100;
const int M=505;
int n,m;
V p[N],zhan[N];
bool cmp(V a,V b){double d=(a-p[1])^(b-p[1]);if(abs(d)>eps) return d>0;else return len(a-p[1])<len(b-p[1]);
}
void ConvexHull(V *p,int &n){int top=0;sort(p+1,p+1+n);sort(p+2,p+1+n,cmp);top=0;for(int i=1;i<=n;i++){while((top>1&&((zhan[top]-zhan[top-1])^(p[i]-zhan[top]))<=0)) --top;zhan[++top]=p[i];}memcpy(p,zhan,sizeof(zhan));n=top;return;
}
signed main(){
#ifndef ONLINE_JUDGE//freopen("a.in","r",stdin);//freopen("a.out","w",stdout);
#endifn=read();for(int i=1;i<=n;i++) input(p[i]);ConvexHull(p,n);zhan[n+1]=p[1];double res(0);for(int i=1;i<=n;i++) res+=len(zhan[i+1]-zhan[i]);//print(zhan[i]);printf("%.2lf\n",res);return 0;
}
/*
3 5
0 -2
-5 3
0 -7
*/

旋转卡壳

前置知识：凸包
个人感觉很像 two-pointers 算法。
能够在优秀的线性时间复杂度内完成总多求最值（周长、面积…）的神奇操作。

给出情境：

给出平面内的 $n$ 个点，求出所有点中的最远点对。
$n≤105n\le 10^5$

首先有一个结论：最远点对一定都是点集的凸包的顶点。
较为显然，证明可以考虑把凸包内的点延长到凸包一条边上，边两边的顶点一定有一个更优。

那么我们就转化成了求凸包上的最远点对，这个问题也叫做凸包的直径问题。

给出一些定义：

凸包的切线：若一条直线过凸包上的一点或一边，且整个凸包都在直线的同侧或在线上，那么我们就称这条直线为凸包的切线。
对踵点：如果经过凸包的两个顶点，可以作两条平行的凸包的切线，那么就称这两个点是一对对踵点。

不难发现，最远点对一定是一对对踵点。
~~然而个人感觉旋转卡壳这个知识点完全不需要这个概念。~~

考虑换一个角度，每次枚举边，然后用到边距离最远的点和边的两端点的距离来更新答案。（每次更新答案的点其实都是对踵点）
显然最优答案一定会被枚举到。
不难发现，如果我们逆时针枚举边，最远点的位置也是在逆时针旋转。
那么我们利用类似 two-pointers 的思想就可以线性的求出答案。
问题得以解决。

实现的细节上，我比较喜欢的方法是一开始先扫一遍暴力找到指针的起始位置，而不是倍长（野蛮）或者每次移动指针都更新答案（玄学）。

代码

P1452 [USACO03FALL]Beauty Contest G /【模板】旋转卡壳

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define debug(...) fprintf(stderr,__VA_ARGS__)
inline ll read(){ll x(0),f(1);char c=getchar();while(!isdigit(c)){if(c=='-')f=-1;c=getchar();}while(isdigit(c)){x=(x<<1)+(x<<3)+c-'0';c=getchar();}return x*f;
}//basic declare
//#define double long double
const double eps=1e-10;
const double pi=acos(-1.0);//---about vectors (or points)//definition
struct V{double x,y;V():x(0),y(0){}V(double a,double b):x(a),y(b){}
};
V ans[10];//declared for other functions
int tot;
inline void input(V &a){scanf("%lf%lf",&a.x,&a.y);}
void print(const V &a,int op=1){printf("%.10lf %.10lf",a.x,a.y);putchar(op?10:32);}
//op:endl or space//basic operation
inline V operator + (const V &a,const V &b){return (V){a.x+b.x,a.y+b.y};}
inline V operator - (const V &a,const V &b){return (V){a.x-b.x,a.y-b.y};}
inline V operator * (const double &x,const V &a){return (V){a.x*x,a.y*x};}
inline V operator * (const V &a,const double &x){return (V){a.x*x,a.y*x};}
inline V operator / (const V &a,const double x){return (V){a.x/x,a.y/x};}
inline bool operator == (const V &a,const V &b){return abs(a.x-b.x)<eps&&abs(a.y-b.y)<eps;}
inline bool operator != (const V &a,const V &b){return !(a==b);}
inline double operator * (const V &a,const V &b){return a.x*b.x+a.y*b.y;}
inline double operator ^ (const V &a,const V &b){return a.x*b.y-a.y*b.x;}
inline double len(const V &a){return sqrt(a.x*a.x+a.y*a.y);}
inline V mid(const V &a,const V &b){return (V){(a.x+b.x)/2,(a.y+b.y)/2};}
inline V chui(const V &a){return (V){a.y,-a.x};}//not take direction into account
inline V danwei(const V &a){return a/len(a);}
inline double tri_S(const V &a,const V &b,const V &c){return abs((b-a)^(c-a))/2;}//always be non-negative
inline bool operator < (const V &a,const V &b){return a.x<b.x-eps||(abs(a.x-b.x)<eps&&a.y<b.y-eps);
}
inline double ang(const V &a,const V &b){return acos((a*b)/len(a)/len(b));}
inline V rotate(const V &o,double t){//COUNTER_CLOCKWISEdouble s=sin(t),c=cos(t);return (V){o.x*c-o.y*s,o.x*s+o.y*c};
}const int N=1e5+100;
const int M=505;
int n,m;
V p[N],zhan[N];
bool cmp(V a,V b){double d=(a-p[1])^(b-p[1]);if(abs(d)>eps) return d>0;else return len(a-p[1])<len(b-p[1]);
}
void graham(V *p,int &n){int top=0;sort(p+1,p+1+n);sort(p+2,p+1+n,cmp);top=0;for(int i=1;i<=n;i++){while((top>1&&((zhan[top]-zhan[top-1])^(p[i]-zhan[top]))<=0)) --top;zhan[++top]=p[i];}memcpy(p,zhan,sizeof(zhan));n=top;return;
}
inline ll calc(const V &a,const V &b){return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+eps;
}
ll rotate_calipers(V *p,int n){ll res(0);int pl=1;for(int i=2;i<=n;i++){if(((p[2]-p[1])^(p[pl]-p[2]))<((p[2]-p[1])^(p[i]-p[2]))) pl=i;}for(int i=1;i<=n;i++){while(((p[i+1]-p[i])^(p[pl]-p[i+1]))<((p[i+1]-p[i])^(p[pl+1]-p[i+1]))){pl=(pl+1)%n;res=max(res,max(calc(p[i],p[pl]),calc(p[i+1],p[pl])));}res=max(res,max(calc(p[i],p[pl]),calc(p[i+1],p[pl])));}return res;
}
signed main(){
#ifndef ONLINE_JUDGEfreopen("a.in","r",stdin);freopen("a.out","w",stdout);
#endifn=read();for(int i=1;i<=n;i++) input(p[i]);graham(p,n);p[0]=p[n];p[n+1]=p[1];//for(int i=1;i<=n;i++) print(p[i]);//putchar('\n');printf("%lld\n",rotate_calipers(p,n));return 0;
}
/*
3 5
0 -2
-5 3
0 -7
*/

半平面交

感觉应用很灵活的一个算法，一切有两个变量的线性规划问题都可以转化为半平面交。
有时可能要注意取等问题（指射箭），~~但多数时候似乎并不用~~。

半平面的表示

首先，我们需要一个简洁的方式表达半平面。
不禁想到，我们可以利用叉积的符号来判断点是在一条线的顺时针还是逆时针方向。那么类似的，我们可以用一个点x和一个方向向量d表示一个半平面，一个点p属于这个半平面，当且仅当 $d→×XP→>0\overrightarrow{d}\times \overrightarrow{XP}>0$ （也就是点在方向向量的左侧），是否取等取决于半平面本身的开闭。

inline bool in(const line &l,const V o){return (l.d^(o-l.a))>eps;}
inline bool out(const line &l,const V o){return (l.d^(o-l.a))<-eps;}

流程

有了一个合适的表示方法后，我们就可以开始尝试解决半平面交问题了，顾名思义，这个问题就是：

给出若干个半平面，求它们的交。

（显然，得到的一定是一个凸包。）
怎么做呢？
首先，我们把所有点半平面按照极角排序，极角可以利用 atan2(y,x) 函数很方便的求出。
当两个向量共线时，显然只有靠左的那个是有用的，我们让有用的放在前面。（这样在后面的算法中遇到没用的时候可以直接跳过，避免求平行线的交点而出现错误）

bool cmp(const line l1,const line l2){double t1=atan2(l1.d.y,l1.d.x),t2=atan2(l2.d.y,l2.d.x);if(abs(t1-t2)>eps) return t1<t2;return in(l2,l1.a);
}

考虑维护一个双端队列来维护当前对答案有用的半平面，并记录每个半平面和前一个半平面的交点。
维护分为如下步骤：

每加入一个半平面时，若队尾与队尾前一个的交点在当前半平面外，则不断弹出队尾。
类似的，若队首和后一个的交点在当前半平面外，也不断弹出队首。
加入当前半平面，并计算其与上一个半平面的交点。
插入所有半平面后，若队尾与上一个的交点在队首半平面外，那么队尾也是无用的，弹出队尾所有没有用的半平面。
最后，插入队尾与队首的交点，得到的封闭的凸包就是答案。

void half_plane(){sort(l+1,l+1+tot,cmp);st=1,ed=0;for(int i=1;i<=tot;i++){if(st<=ed&&abs(q[ed].d^l[i].d)<eps) continue;//和上一个共线，必然是无用向量 while(ed-st+1>=2&&out(l[i],pt[ed])) --ed;//先弹队尾 while(ed-st+1>=2&&out(l[i],pt[st+1])) ++st;q[++ed]=l[i];if(ed-st+1>=2) pt[ed]=jiaodian(q[ed-1],q[ed]);}while(ed-st+1>=2&&out(q[st],pt[ed])) --ed;//弹出队尾无用向量 if(ed-st+1>=2) pt[ed+1]=jiaodian(q[st],q[ed]),++ed;//加入队尾与队首的交点 m=ed-st;for(int i=1;i<=m;i++) a[i]=pt[st+i];a[m+1]=a[1];return;
}

细节

可不可以先弹队首再弹队尾？

不可以！
考虑这两种写法不一样的情况，就是当我判断在半平面外的那个交点恰好就是队首和队尾的交点时。
也就是：
在这里插入图片描述
（图片来自：OI-wiki 半平面交）

显然这个时候弹队首会出错。

如何判断交集为空或无边界？

考虑在正无穷远处加入四条与坐标轴平行的边，把整个平面框成一个大矩形，在进行半平面交算法。
这样就不会出现无边界情况，要判断也可以通过凸包内是否出现了那四条边来进行。
而且，这样之后，凸包元素不足三个就变成了答案无解的充要条件。
非常优美简洁。

代码

P4196 [CQOI2006]凸多边形 /【模板】半平面交

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define debug(...) fprintf(stderr,__VA_ARGS__)
inline ll read(){ll x(0),f(1);char c=getchar();while(!isdigit(c)){if(c=='-')f=-1;c=getchar();}while(isdigit(c)){x=(x<<1)+(x<<3)+c-'0';c=getchar();}return x*f;
}//basic declare
//#define double long double
const double eps=1e-10;
const double pi=acos(-1.0);//---about vectors (or points)//definition
struct V{double x,y;V():x(0),y(0){}V(double a,double b):x(a),y(b){}
};
inline void input(V &a){scanf("%lf%lf",&a.x,&a.y);}
void print(const V &a,int op=1){printf("%g %g",a.x,a.y);putchar(op?10:32);}
//op:endl or space//basic operation
inline V operator + (const V &a,const V &b){return (V){a.x+b.x,a.y+b.y};}
inline V operator - (const V &a,const V &b){return (V){a.x-b.x,a.y-b.y};}
inline V operator * (const double &x,const V &a){return (V){a.x*x,a.y*x};}
inline V operator * (const V &a,const double &x){return (V){a.x*x,a.y*x};}
inline V operator / (const V &a,const double x){return (V){a.x/x,a.y/x};}
inline bool operator == (const V &a,const V &b){return abs(a.x-b.x)<eps&&abs(a.y-b.y)<eps;}
inline bool operator != (const V &a,const V &b){return !(a==b);}
inline double operator * (const V &a,const V &b){return a.x*b.x+a.y*b.y;}
inline double operator ^ (const V &a,const V &b){return a.x*b.y-a.y*b.x;}
inline double len(const V &a){return sqrt(a.x*a.x+a.y*a.y);}
inline V mid(const V &a,const V &b){return (V){(a.x+b.x)/2,(a.y+b.y)/2};}
inline V chui(const V &a){return (V){a.y,-a.x};}//not take direction into account
inline V danwei(const V &a){return a/len(a);}
inline double ang(const V &a,const V &b){return acos((a*b)/len(a)/len(b));}
inline V rotate(const V &o,double t){//COUNTER_CLOCKWISEdouble s=sin(t),c=cos(t);return (V){o.x*c-o.y*s,o.x*s+o.y*c};
}
struct line{V d,a;
};
void print(const line &l,int op=1){printf("point: ");print(l.a,0);printf("dir: ");print(l.d,op);
}
inline line trans(const V &a,const V &b){return (line){b-a,a};
}
inline line trans(const double &a,const double &b,const double &c,const double &d){
//given one point and directionreturn (line){(V){c,d},(V){a,b}};
}
inline V jiaodian(const line &u,const line &v){if(u.d==v.d){print(u);print(v);exit(0);}double k=((v.a-u.a)^v.d)/(u.d^v.d);return u.a+k*u.d;
}inline bool in(const line &l,const V o){return (l.d^(o-l.a))>0;}
inline bool out(const line &l,const V o){return (l.d^(o-l.a))<-eps;}const int N=1005;
int n,m;int tot,st,ed;
V a[N],pt[N];
line l[N],q[N];bool cmp(const line l1,const line l2){double t1=atan2(l1.d.y,l1.d.x),t2=atan2(l2.d.y,l2.d.x);if(abs(t1-t2)>eps) return t1<t2;return in(l2,l1.a);
}const double inf=1e9;
void half_plane(){l[++tot]=trans(0,0,1,0);//插入边界l[++tot]=trans(X,0,0,1);l[++tot]=trans(0,Y,-1,0);l[++tot]=trans(0,0,0,-1);sort(l+1,l+1+tot,cmp);st=1,ed=0;for(int i=1;i<=tot;i++){if(st<=ed&&abs(q[ed].d^l[i].d)<eps) continue;//和上一个共线，必然是无用向量 while(ed-st+1>=2&&out(l[i],pt[ed])) --ed;//先弹队尾 while(ed-st+1>=2&&out(l[i],pt[st+1])) ++st;q[++ed]=l[i];if(ed-st+1>=2) pt[ed]=jiaodian(q[ed-1],q[ed]);}while(ed-st+1>=2&&out(q[st],pt[ed])) --ed;//弹出队尾无用向量 if(ed-st+1>=2) pt[ed+1]=jiaodian(q[st],q[ed]),++ed;//加入队尾与队首的交点 m=ed-st;for(int i=1;i<=m;i++) a[i]=pt[st+i];a[m+1]=a[1];return;
}
signed main(){
#ifndef ONLINE_JUDGEfreopen("a.in","r",stdin);freopen("a.out","w",stdout);
#endifn=read();for(int i=1;i<=n;i++){m=read();for(int j=1;j<=m;j++) input(a[j]);a[m+1]=a[1];for(int j=1;j<=m;j++) l[++tot]=trans(a[j],a[j+1]);}half_plane();if(m<3){printf("0.000");return 0;}double res(0);for(int i=1;i<=m;i++) res+=(a[i]^a[i+1])/2;if(res<0) res=-res;printf("%.3lf\n",res);return 0;
}
/*
4
2 8
4 2
6 5
-8 2
*/

其他

皮克定理

条件：整点任意多边形。
内容：面积=整点数+1/2*边点数-1。
（为什么减一？可以考虑极限情况，多边形退化成单格点时，面积是0）

拓展：对于平行四边形的格点依然成立；对于三角形的格点，面积要乘二，即2面积=整点数+1/2边点数-1，可以理解成变成三角形后面积减少了一半所以要乘二才能相等。

距离

规定：以下的 $d$ 表示维度， $x_{i,d}$ 表示点 $i$ 在第 $d$ 维的坐标。

欧几里得距离 ： $∑d(xi,d−xj,d)2\sqrt{\sum_d (x_{i,d}-x_{j,d})^2 }$
曼哈顿距离： $∑d∣xi,d−xj,d∣\sum_d|x_{i,d-}x_{j,d}|$
切比雪夫距离： $max_d|x_{i,d-}x_{j,d}|$

曼哈顿距离和切比雪夫距离的转化

曼哈顿距离下的点 $(x, y)$ 等价于切比雪夫距离下的点 $(x + y, x - y)$ 。
反过来，切比雪夫下的点 $(x, y)$ 等价于切比雪夫距离下的点 $(x+y2,x−y2)(\dfrac{x+y}{2},\dfrac{x-y}{2})$ 。
证明：暴力拆分曼哈顿的定义或者结合单位正方形。

source

//basic declare
//#define double long double
const double eps=1e-10;
const double pi=acos(-1.0);//---about vectors (or points)//definition
struct V{double x,y;V():x(0),y(0){}V(double a,double b):x(a),y(b){}
};
V ans[10];//declared for other functions
int tot;
inline void input(V &a){a.x=read();a.y=read();}
void print(const V &a,int op=1){printf("%.10lf %.10lf",a.x,a.y);putchar(op?10:32);}
//op:endl or space//basic operation
inline V operator + (const V &a,const V &b){return (V){a.x+b.x,a.y+b.y};}
inline V operator - (const V &a,const V &b){return (V){a.x-b.x,a.y-b.y};}
inline V operator * (const double &x,const V &a){return (V){a.x*x,a.y*x};}
inline V operator * (const V &a,const double &x){return (V){a.x*x,a.y*x};}
inline V operator / (const V &a,const double x){return (V){a.x/x,a.y/x};}
inline bool operator == (const V &a,const V &b){return abs(a.x-b.x)<eps&&abs(a.y-b.y)<eps;}
inline bool operator != (const V &a,const V &b){return !(a==b);}
inline double operator * (const V &a,const V &b){return a.x*b.x+a.y*b.y;}
inline double operator ^ (const V &a,const V &b){return a.x*b.y-a.y*b.x;}
inline double len(const V &a){return sqrt(a.x*a.x+a.y*a.y);}
inline V mid(const V &a,const V &b){return (V){(a.x+b.x)/2,(a.y+b.y)/2};}
inline V chui(const V &a){return (V){a.y,-a.x};}//not take direction into account
inline V danwei(const V &a){return a/len(a);}
inline double tri_S(const V &a,const V &b,const V &c){return abs((b-a)^(c-a))/2;}//always be non-negative
inline bool operator < (const V &a,const V &b){return a.x<b.x-eps||(abs(a.x-b.x)<eps&&a.y<b.y-eps);
}//operations about angle and vectors
inline double ang(const V &a,const V &b){return acos((a*b)/len(a)/len(b));}
inline bool dun(const V &a,const V &b,const V &c){return ((b-a)*(c-a))<-eps;}//angle:BAC
inline bool rui(const V &a,const V &b,const V &c){return ((b-a)*(c-a))>eps;}
inline bool zhi(const V &a,const V &b,const V &c){return abs(((b-a)*(c-a)))<eps;}
inline V rotate(const V &o,double t){//COUNTER_CLOCKWISEdouble s=sin(t),c=cos(t);return (V){o.x*c-o.y*s,o.x*s+o.y*c};
}//---about lines//definition
struct line{V d,a,b;
};
inline line trans(double a,double b,double c,double d){//given coordinatesV dd(c-a,d-b),x(a,b),y(c,d);dd=dd/len(dd);return (line){dd,x,y};
}
inline line trans(const V &a,const V &b){//given pointsV dd(b.x-a.x,b.y-a.y);dd=dd/len(dd);return (line){dd,a,b};
}
inline void input(line &l){double a=read(),b=read(),c=read(),d=read();l=trans(a,b,c,d);return;
}//functions about points and lines
inline V chui(const line &l,const V &o){//directedreturn l.a+((o-l.a)*l.d)*l.d;
}
inline V duichen(const line &l,const V &o){V u=chui(l,o);return (V){2*u.x-o.x,2*u.y-o.y};
}
inline bool on_line(const line &l,const V &o){return abs((l.d^(o-l.a)))<eps;}
inline bool on_seg(const line &l,const V &o){return abs(len(o-l.a)+len(o-l.b)-len(l.a-l.b))<eps;
}
inline int pos(const line &l,const V &o){if(!on_line(l,o)){if((l.d^(o-l.a))>0) return 1;//counter clockwiseelse return 2;//clockwise}else{if(((o-l.a)*(o-l.b))<=0) return 5;//on segmentelse if(((o-l.a)*l.d)<0) return 3;//online backelse return 4;//online front}
}
inline double dis(const line &l,const V &o,int op=0){//op=0:dis on line,op=1:dis on segmentif(op&&(dun(l.a,o,l.b)||dun(l.b,o,l.a))) return min(len(o-l.a),len(o-l.b));else return abs((o-l.a)^(o-l.b))/len(l.a-l.b);
}//functions about lines and lines
inline bool gongxian(const line &a,const line &b){return abs(a.d^b.d)<eps;}
inline bool chuizhi(const line &a,const line &b){return abs(a.d*b.d)<eps;}
inline bool jiao(const line &u,const line &v){if(min(u.a.x,u.b.x)>max(v.a.x,v.b.x)+eps||max(u.a.x,u.b.x)<min(v.a.x,v.b.x)-eps||min(u.a.y,u.b.y)>max(v.a.y,v.b.y)+eps||max(u.a.y,u.b.y)<min(v.a.y,v.b.y)-eps) return false;//rapid rejection testreturn ((u.a-v.a)^v.d)*((u.b-v.a)^v.d)<eps&&((v.a-u.a)^u.d)*((v.b-u.a)^u.d)<eps;//straddle test
}
inline double Sin(const V &a,const V &b){return abs(a^b)/len(a)/len(b);}
inline V jiaodian(const line &u,const line &v){if(on_line(v,u.a)) return u.a;if(on_line(v,u.b)) return u.b;double s1=tri_S(u.a,v.a,v.b),s2=tri_S(u.b,v.a,v.b);return u.a+((s1/(s1+s2))*(u.b-u.a));
}
inline line pingfen(const V &a,const V &b,const V &c){//angle:BACV d1=danwei(b-a),d2=danwei(c-a),d=danwei(d1+d2);return (line){d,a,a+d};
}
inline line zhongchui(const V &a,const V &b){V x=mid(a,b),d=danwei(chui(b-a));return (line){d,x,x+d};
}//---about polygons//definition
//V a[N];
//int n;//funtions
double S(const V  *a,const int n){double res(0);for(int i=1;i<=n;i++) res+=(a[i]^a[i%n+1]);return res/2;
}
bool is_convex(const V *a,const int n){int op=0;for(int i=1;i<=n;i++){double o=(a[i+1]-a[i])^(a[i]-a[i-1]);if(abs(o)<eps) continue;//three neighbouring points on the same lineint nop=o>0?1:-1;if(!op) op=nop;else if(op!=nop) return false;}return true;
}
int in_poly(const V *a,const int n,const V o){int res(0);for(int i=1;i<=n;i++){V u=a[i],v=a[i+1];if(on_seg(trans(u,v),o)) return 1;//on the edgeif(abs(u.y-v.y)<eps) continue;//ignore horizontalif(max(u.y,v.y)<o.y-eps) continue;//equal will not continueif(min(u.y,v.y)>o.y-eps) continue;//equal will continuedouble x=u.x+(o.y-u.y)/(v.y-u.y)*(v.x-u.x);if(x<o.x) res^=1;}return res?2:0;//odd:in even:out
}//---about circles//definition
struct cir{V o;double r;cir(double a=0,double b=0,double c=0):o(a,b),r(c){}cir(V O,double c=0):o(O),r(c){}
};
inline void input(cir &cc){input(cc.o);cc.r=read();}//functions about position relation
inline bool incir(const cir &c,const V &p){return len(p-c.o)<c.r-eps;}
inline bool oncir(const cir &c,const V &p){return (len(p-c.o)-c.r)<eps;}
inline bool outcir(const cir &c,const V &p){return len(p-c.o)>c.r+eps;}
inline double dis(const cir &c,const line &l){return dis(l,c.o);}
inline int pos(const cir &c,const line &l){double d=dis(c,l);if(d<c.r-eps) return 1;//intersectelse if(abs(d-c.r)<eps) return 2;//tangentelse if(d>c.r+eps) return 3;//disjoint
}
inline double dis(const cir &c1,const cir &c2){return len(c1.o-c2.o);}
inline int pos(const cir &c1,const cir &c2){double d=dis(c1,c2);if(d>c1.r+c2.r+eps) return 4;//do not crosselse if(abs(d-c1.r-c2.r)<eps) return 3;//circumscribedelse if(max(c1.r,c2.r)<min(c1.r,c2.r)+d -eps) return 2;//intersectelse if(abs(max(c1.r,c2.r)-min(c1.r,c2.r)-d)<eps) return 1;//inscribedelse return 0;//include
}//functions about circles and triangles
inline cir incir(const V &a,const V &b,const V &c){V x(jiaodian(pingfen(a,b,c),pingfen(b,a,c)));return cir(x,dis(trans(a,b),x));
}
inline cir outcir(const V &a,const V &b,const V &c){V x(jiaodian(zhongchui(a,b),zhongchui(a,c)));return cir(x,len(x-a));
}//functions about cross points
inline double calc(double xie,double zhi){return sqrt(xie*xie-zhi*zhi);}//Pythagorean Theorem
inline void cross_line_cir(const cir &c,const line &l){tot=0;V x=chui(l,c.o);double dd=len(x-c.o);if(c.r<dd-eps) return;//none cross pointif(abs(c.r-dd)<eps){//one cross pointans[++tot]=x;return;}double dis=calc(c.r,dd);V a=x+dis*l.d,b=x-dis*l.d;//two cross pointsans[++tot]=a;ans[++tot]=b;return;
}
inline void cross_cir_cir(const cir &c1,const cir c2){int op=pos(c1,c2);if(op==4||op==0) return;//none cross pointV L=c2.o-c1.o;double a=c2.r,b=c1.r,c=len(L);double t=acos((b*b+c*c-a*a)/(2*b*c));//find the angleV x=c1.o+c1.r*rotate((L)/len(L),t),y=duichen(trans(c1.o,c2.o),x);if(x<y) print(x,0),print(y,1);else print(y,0),print(x,1);
}//functions about tangent points and lines
inline void qiedian(const cir &c,const V &p){tot=0;line L=trans(c.o,p);double t=acos(c.r/len(c.o-p));V a=c.o+(c.r*rotate(L.d,t)),b=duichen(L,a);ans[++tot]=a;if(a!=b) ans[++tot]=b;return;
}
inline void common_qie(const cir &c1,const cir &c2){tot=0;int op=pos(c1,c2);if(op){//outside tangent linesdouble d=c1.r-c2.r,t=acos(d/len(c1.o-c2.o));V u=c1.o+(c1.r*danwei(rotate(c2.o-c1.o,t))),v=c1.o+(c1.r*danwei(rotate(c2.o-c1.o,-t)));ans[++tot]=u;if(u!=v) ans[++tot]=v;   }if(op>2){//inside tangent linesdouble d=len(c2.o-c1.o)/(c1.r+c2.r)*c1.r,t=acos(c1.r/d);V u=c1.o+(c1.r*danwei(rotate(c2.o-c1.o,t))),v=c1.o+(c1.r*danwei(rotate(c2.o-c1.o,-t)));ans[++tot]=u;if(u!=v) ans[++tot]=v;}
}//functions about area
inline double cir_S(const cir &c){return pi*c.r*c.r;}
inline double shan(const cir &c,const V &a,const V &b){//S of sectordouble t=ang(a-c.o,b-c.o);return t/2*c.r*c.r;
}
inline double gong(const cir &c,const V &a,const V &b){//S of bowreturn shan(c,a,b)-tri_S(c.o,a,b);
}
inline double O_tri(const cir &c,V a,V b){
//directed S of Intersection between circle O and triangle OABif(on_line(trans(a,b),c.o)) return 0.0;int f=(((a-c.o)^(b-c.o))>0)?1:-1;//directionline l=trans(a,b);int f1=incir(c,a),f2=incir(c,b);if(f1&&f2)  return f*tri_S(a,b,c.o);//both incircle:the S of triangleelse if(!f1&&!f2){//both outcircle:a sector(possibly minus a bow)V u=(c.o+(c.r*danwei(a-c.o))),v=(c.o+(c.r*danwei(b-c.o)));double S=shan(c,u,v);    if(dis(l,c.o,1)<c.r-eps){cross_line_cir(c,l);assert(tot==2);S-=gong(c,ans[1],ans[2]);}return f*S;}else{//one in and one out:a traingle and a sectorif(!f1) swap(a,b);double sa=Sin(b-a,c.o-a),su=sa/c.r*len(c.o-a),A=asin(sa),U=asin(su);if((b-a)*(c.o-a)<-eps) A=pi-A;double t=pi-A-U,dis=sin(t)/sa*c.r;    V u=a+(dis*danwei(b-a)),v=c.r*danwei(b-c.o);return f*(tri_S(c.o,a,u)+shan(c,u,v));}
}
double common_cir_poly(const V *a,const cir &c,int n){double res(0);for(int i=1;i<=n;i++) res+=O_tri(c,a[i],a[i+1]);return res;
}
inline double common_cir_cir(const cir &c1,const cir &c2){int op=pos(c1,c2);if(op>=3) return 0;//common area=0else if(op<=1) return min(cir_S(c1),cir_S(c2));//completly includeelse{//partly includedouble L=len(c1.o-c2.o);double t1=2*acos((L*L+c1.r*c1.r-c2.r*c2.r)/(2*L*c1.r));double t2=2*acos((L*L+c2.r*c2.r-c1.r*c1.r)/(2*L*c2.r));return c1.r*c1.r*t1/2-c1.r*c1.r*sin(t1)/2+c2.r*c2.r*t2/2-c2.r*c2.r*sin(t2)/2;}
}