计算几何全家桶

已于 2022-06-15 14:59:49 修改
阅读量3.4k 收藏 51
点赞数 23
分类专栏： # 计算几何 模板 精华汇总 文章标签： 计算几何
于 2021-12-30 15:29:32 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/BUG_Creater_jie/article/details/122226874
版权

文章目录

前言

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

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

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

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

精度

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

a = b → abs ⁡ ( a − b ) < e p s a=b \to\operatorname{abs}(a-b)<eps a=b→abs(a−b)<eps
a < b → a < b − e p s a<b \to a<b-eps a<b→a<b−eps
a > b → a > b + e p s a>b \to a>b+eps a>b→a>b+eps

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

点/向量相关

表示

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

//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列的行列式，第一行为各维度的单位向量，二三行为两个点各自在对应维度的坐标。

注意：

  1. 严谨来说，向量叉积是一个向量而不是一个数，它的模长等于 x 1 y 2 − x 2 y 1 x_1y_2-x_2y_1 x1​y2​−x2​y1​，方向垂直与纸面，遵循右手定则，这也是叉积可以表示负面积的原因。
  2. 这里的垂直向量没有关注方向。
//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 ) (x,y) (x,y)，旋转角为 θ \theta θ，那么旋转后的向量就是 ( x cos ⁡ θ − y sin ⁡ θ , x sin ⁡ θ + y cos ⁡ θ ) (x\cos\theta-y\sin\theta,x\sin\theta+y\cos\theta) (xcosθ−ysinθ,xsinθ+ycosθ)。
可以用虚数乘法， ( x + y i ) × ( cos ⁡ θ , sin ⁡ θ   i ) (x+yi)\times (\cos\theta,\sin\theta \space i) (x+yi)×(cosθ,sinθ 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 coordinates
  V 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 points
  V 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){//directed
  return 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 clockwise
    else return 2;//clockwise
  }
  else{
    if(((o-l.a)*(o-l.b))<eps) return 5;//on segment
    else if(((o-l.a)*l.d)<-eps) return 3;//online back
    else 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 segment
  if(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;}
判断线段与线段是否相交

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

  1. 快速排斥实验：定义一条线段所在的矩形为边与坐标轴平行、能完全包含该线段的最小的矩形，如果两条线段所在的矩形没有公共部分，显然不可能有交点；否则，进入下一步跨立实验。
  2. 跨立实验：如果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 test
  return ((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 l1,l2 l1,l2，线上各有一点 x 1 , x 2 x1,x2 x1,x2，方向向量为 d 1 , d 2 d1,d2 d1,d2。
假设交点 p = x 1 + k d 1 p=x1+kd1 p=x1+kd1，就有 ( p − x 2 ) × d 2 = 0 (p-x2)\times d2=0 (p−x2)×d2=0。
叉积是满足分配率的，拆开化简得到：
k = ( x 2 − x 1 ) × d 2 d 1 × d 2 k=\frac{(x2-x1)\times d2}{d1\times d2} k=d1×d2(x2−x1)×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:BAC
  V 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 O O，则多边形面积为：
S = 1 2 ∑ i = 1 n O A i → × O A → i % n + 1 S=\frac{1}{2}\sum_{i=1}^n\overrightarrow{OA_{i}}\times \overrightarrow{OA}_{i\%n+1} S=21​i=1∑n​OAi​ ​×OA i%n+1​
为了方便， O O O 直接取原点即可。
时间复杂度 O ( n ) O(n) 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 ) (a,b,c) (a,b,c)，作为一个三元组。
判断多边形凹凸的充要条件：对于多边形所有的 n n n 个这样三元组， c c c 和 a b ab ab 的位置关系（指顺时针或逆时针）必然全部都是相同的。
叉积判断符号即可。
注意如果没有保证相邻三个点不共线的话需要特判。
时间复杂度 O ( n ) O(n) 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 line
    int nop=o>0?1:-1;
    if(!op) op=nop;
    else if(op!=nop) return false;
  }
  return true;
}

判定点在多边形内

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

  1. 如果当前边经过了该点，直接返回相应结果
  2. 由于可能出现射线经过多边形顶点导致一个点被计算两次的情况，强制令这个点在下面的边被统计。实现上，当 max ⁡ ( y 1 , y 2 ) = y 0 \max(y1,y2)=y0 max(y1,y2)=y0 时，统计该次答案，而当 min ⁡ ( y 1 , y 2 ) = y 0 \min(y1,y2)=y0 min(y1,y2)=y0 时，不统计该次答案。
  3. 出现水平边时，直接忽略即可，结合第二条的处理原则，就依然是正确的。
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 edge
    if(abs(u.y-v.y)<eps) continue;//ignore horizontal
    if(max(u.y,v.y)<o.y-eps) continue;//equal will not continue
    if(min(u.y,v.y)>o.y-eps) continue;//equal will continue
    double 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 ) a_1=(0,0) a1​=(0,0)。（如果不满足整体平移即可）
lowerbound 找到极角夹在询问点两侧的相邻点对，判断是否在二者连成的边内部即可。
注意需要对于在 a 1 , a n a_1,a_n a1​,an​ 外侧的点需要特判。

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;//intersect
	else if(abs(d-c.r)<eps) return 2;//tangent
	else 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 cross
  else if(abs(d-c1.r-c2.r)<eps) return 3;//circumscribed
  else if(max(c1.r,c2.r)<min(c1.r,c2.r)+d -eps) return 2;//intersect
  else if(abs(max(c1.r,c2.r)-min(c1.r,c2.r)-d)<eps) return 1;//inscribed
  else return 0;//include
}

三角形和圆

三角形内切圆

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

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 a,b,c a,b,c，则 R = a b c ( 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)}} R=(a+b+c)(a+b−c)(a+c−b)(b+c−a) ​abc​。

交点相关

圆与直线交点

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

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 point
  if(abs(c.r-dd)<eps){//one cross point
    ans[++tot]=x;return;
  }
  double dis=calc(c.r,dd);
  V a=x+dis*l.d,b=x-dis*l.d;//two cross points
  ans[++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 point
  V 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 angle
  V 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 lines
    double 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 lines
  	double 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 sector
  double 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 bow
	return 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 OAB
  if(on(trans(a,b),c.o)) return 0.0;
  int f=(((a-c.o)^(b-c.o))>0)?1:-1;//direction
  line 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 triangle
  else 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 sector
    if(!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=0
  else if(op<=1) return min(cir_S(c1),cir_S(c2));//completly include
  else{//partly include
  	double 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 ( n log ⁡ n ) O(n\log n) O(nlogn)。
这里介绍较为简单的 Graham 扫描法。

首先，找出所有点中按照 X − Y X-Y X−Y 排序最小的点 O O O。
然后以 O 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_CLOCKWISE
  double 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);
#endif
  n=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 n 个点，求出所有点中的最远点对。
n ≤ 1 0 5 n\le 10^5 n≤105

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

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

给出一些定义：

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

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

考虑换一个角度，每次枚举边，然后用到边距离最远的点和边的两端点的距离来更新答案。（每次更新答案的点其实都是对踵点）
显然最优答案一定会被枚举到。
不难发现，如果我们逆时针枚举边，最远点的位置也是在逆时针旋转。
那么我们利用类似 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_CLOCKWISE
  	double 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_JUDGE
  freopen("a.in","r",stdin);
  freopen("a.out","w",stdout);
#endif
	n=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 → × X P → > 0 \overrightarrow{d}\times \overrightarrow{XP}>0 d ×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);
}

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

  1. 每加入一个半平面时， 若队尾与队尾前一个的交点在当前半平面外，则不断弹出队尾。
  2. 类似的，若队首和后一个的交点在当前半平面外，也不断弹出队首。
  3. 加入当前半平面，并计算其与上一个半平面的交点。
  4. 插入所有半平面后，若队尾与上一个的交点在队首半平面外，那么队尾也是无用的，弹出队尾所有没有用的半平面。
  5. 最后，插入队尾与队首的交点，得到的封闭的凸包就是答案。
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;
}

细节

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

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

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

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

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

代码

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_CLOCKWISE
  	double 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 direction
	return (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_JUDGE
  	freopen("a.in","r",stdin);
  	freopen("a.out","w",stdout);
#endif
	n=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 d d 表示维度， x i , d x_{i,d} xi,d​ 表示点 i i i 在第 d d d 维的坐标。

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

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

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

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_CLOCKWISE
  double 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 coordinates
  V 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 points
  V 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){//directed
  return 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 clockwise
    else return 2;//clockwise
  }
  else{
    if(((o-l.a)*(o-l.b))<=0) return 5;//on segment
    else if(((o-l.a)*l.d)<0) return 3;//online back
    else 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 segment
  if(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 test
  return ((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:BAC
  V 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 line
    int 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 edge
    if(abs(u.y-v.y)<eps) continue;//ignore horizontal
    if(max(u.y,v.y)<o.y-eps) continue;//equal will not continue
    if(min(u.y,v.y)>o.y-eps) continue;//equal will continue
    double 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;//intersect
	else if(abs(d-c.r)<eps) return 2;//tangent
	else 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 cross
  else if(abs(d-c1.r-c2.r)<eps) return 3;//circumscribed
  else if(max(c1.r,c2.r)<min(c1.r,c2.r)+d -eps) return 2;//intersect
  else if(abs(max(c1.r,c2.r)-min(c1.r,c2.r)-d)<eps) return 1;//inscribed
  else 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 point
  if(abs(c.r-dd)<eps){//one cross point
    ans[++tot]=x;return;
  }
  double dis=calc(c.r,dd);
  V a=x+dis*l.d,b=x-dis*l.d;//two cross points
  ans[++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 point
  V 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 angle
  V 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 lines
    double 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 lines
  	double 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 sector
  double 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 bow
	return 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 OAB
  if(on_line(trans(a,b),c.o)) return 0.0;
  int f=(((a-c.o)^(b-c.o))>0)?1:-1;//direction
  line 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 triangle
  else 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 sector
    if(!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=0
  else if(op<=1) return min(cir_S(c1),cir_S(c2));//completly include
  else{//partly include
  	double 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;
  }
}
wind__whisper
关注
  • 23
    点赞
  • 踩
  • 51
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
专栏目录
【计算几何全家桶】计算几何二维板子套路详解
SC_Linno的博客
11-16 758
title : 计算几何（二维） tags : ACM,数学,计算几何 date : 2021-10-23 author : Linno 向量 点乘与叉乘 叉乘一般用于求面积或者判断两点的相对位置：(x1,y1)×(x2,y2)=x1*y2-x2*y1; 点乘一般用于判断直线是否正交：(x1,y1)·(x2,y2)=x1*x2+y1*y2; inline double Len(Vt a){return sqrt(Dot(a,a));}//模长 inline Vt operator +(Vt a,Vt .
一文搞懂PointNet全家桶——强势的点云处理神经网络
热门推荐
limingmin2020的博客
03-21 2万+
PointNet、PointNet++、 F-PointNet一网打尽！
计算几何资料
09-25
ACM计算几何资料，从基本到全面的介绍计算几何题，有模板，全！
二维几何常用算法合集 繁凡的《计算几何全家桶（二）》1
08-03
1.平面扫描2 2.凸包3 3.旋转卡壳5 4.半平面交7 5.闵可夫斯基和 10 6.平面区域10 7.平面最近点对 12 1.平面扫描 2.凸包
[转]计算几何专题
weixin_30401605的博客
02-13 192
这两天在学习计算几何，随便说说自己的学习过程吧。 　　基本的叉积、点积和凸包等东西就不多说什么了，网上一搜一大堆，切一些题目基本熟悉了就差不多了。 　　一些基本的题目可以自己搜索，比如这个blog：http://blog.sina.com.cn/s/blog_49c5866c0100f3om.html 　　接下来，研究了半平面交，思想方法看07年朱泽园的国家队论文，模板代码参考自我校大牛韬哥...
二维几何基础大合集！《计算几何全家桶（一）》（基础运算、点、线、多边形、圆、网格）
繁凡さん的博客
09-27 6594
整理的算法模板合集： ACM模板 目录基本运算点与直线多边形 基本运算 #include<bits/stdc++.h> using namespace std; const int N = 5007, M = 50007, INF = 0x3f3f3f3f; const double DINF = 12345678910, eps = 1e-10; struct Point{ double x, y; Point(double x = 0, double y = 0):x(x
二维计算几何全家桶
最新发布
zcxxn的博客
02-09 554
计算几何——从入门到入土
冈萨雷斯《数字图像处理》全家桶
11-22
《冈萨雷斯《数字图像处理》全家桶》是一份涵盖了该经典教材多个版本及配套资源的综合包。这个包包含了以下文件： 1. **Digital Image Processing Using Matlab.zip**：这是一份基于Matlab的数字图像处理教程，可能...
我的所有优质博客全部开源啦（我自己原创的《ACM模板》《算法全家桶》《算法竞赛中的初等数论》 PDF免费下载）
繁凡さん的博客
11-13 1万+
你好呀ヾ(≧▽≦*)o 我是繁凡さん
计算几何
08-06
目录 ㈠ 点的基本运算 1. 平面上两点之间距离 1 2. 判断两点是否重合 1 3. 矢量叉乘 1 4. 矢量点乘 2 5. 判断点是否在线段上 2 6. 求一点饶某点旋转后的坐标 2 7. 求矢量夹角 2 ㈡ 线段及直线的基本运算 1. 点与线段的关系 3 2. 求点到线段所在直线垂线的垂足 4 3. 点到线段的最近点 4 4. 点到线段所在直线的距离 4 5. 点到折线集的最近距离 4 6. 判断圆是否在多边形内 5 7. 求矢量夹角余弦 5 8. 求线段之间的夹角 5 9. 判断线段是否相交 6 10.判断线段是否相交但不交在端点处 6 11.求线段所在直线的方程 6 12.求直线的斜率 7 13.求直线的倾斜角 7 14.求点关于某直线的对称点 7 15.判断两条直线是否相交及求直线交点 7 16.判断线段是否相交，如果相交返回交点 7 ㈢ 多边形常用算法模块 1. 判断多边形是否简单多边形 8 2. 检查多边形顶点的凸凹性 9 3. 判断多边形是否凸多边形 9 4. 求多边形面积 9 5. 判断多边形顶点的排列方向，方法一 10 6. 判断多边形顶点的排列方向，方法二 10 7. 射线法判断点是否在多边形内 10 8. 判断点是否在凸多边形内 11 9. 寻找点集的graham算法 12 10.寻找点集凸包的卷包裹法 13 11.判断线段是否在多边形内 14 12.求简单多边形的重心 15 13.求凸多边形的重心 17 14.求肯定在给定多边形内的一个点 17 15.求从多边形外一点出发到该多边形的切线 18 16.判断多边形的核是否存在 19 ㈣ 圆的基本运算 1 .点是否在圆内 20 2 .求不共线的三点所确定的圆 21 ㈤ 矩形的基本运算 1.已知矩形三点坐标，求第4点坐标 22 ㈥ 常用算法的描述 22 ㈦ 补充 1．两圆关系： 24 2．判断圆是否在矩形内： 24 3．点到平面的距离： 25 4．点是否在直线同侧： 25 5．镜面反射线： 25 6．矩形包含： 26 7．两圆交点： 27 8．两圆公共面积： 28 9. 圆和直线关系： 29 10. 内切圆： 30 11. 求切点： 31 12. 线段的左右旋： 31 13．公式： 32
多项式全家桶
olahiuj的博客
01-15 611
施工现场 多项式 在数学中，由若干个单项式相加组成的代数式叫做多项式（若有减法：减一个数等于加上它的相反数）。多项式中的每个单项式叫做多项式的项，这些单项式中的最高项次数，就是这个多项式的次数。其中多项式中不含字母的项叫做常数项。 为了方便接下来瞎逼逼，我们定义n次多项式A(x)=∑i=1nai⋅xiA(x)=\sum\limits_{i=1}^n a_i\cdot x^iA(x)=i=1∑n​...
三维几何基础大合集（三维点积叉积、点线面、凸包）《计算几何全家桶（三）》
繁凡さん的博客
10-03 2608
整理的算法模板合集： ACM模板 目录
计算几何模板
w750636248的专栏
07-26 685
#include #define MAXN 1000 #define offset 10000 #define eps 1e-8 #define PI acos(-1.0)//3.14159265358979323846 //判断一个数是否为0,是则返回true,否则返回false #define zero(x)(((x)>0?(x):-(x))<eps) //返回一个数的符号，正数返回1，负数返
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值