计算几何初步

一.点与向量.

首先，在计算几何中的向量都是坐标表示，而点也是用坐标表示，所以大家都是直接把点和向量看成一个东西的.

然后是各种运算，对于两个向量$\vec{a}=(x_1,y_1),\vec{b}=(x_2,y_2)$，且它们对应的角度分别为$\alpha ,\beta$，那么有如下运算：
$$\begin{gathered} \vec{a}\pm \vec{b}=(x_1\pm x_2,y_1\pm y_2) \\ \vec{a}\cdot \vec{b}=|\vec{a}||\vec{b}|\cos (\beta -\alpha ) \\ \vec{a}\times \vec{b}=|\vec{a}||\vec{b}|\sin (\beta -\alpha ) \end{gathered}$$

由于结构体重载的时候并没有$\cdot$和$\times$所以我们用*和%代替.

同时我们用~a表示$|\vec{a}{|}^2$.

我们还定义了向量的单位化Unit函数和法向量Norm函数，以及一个判定一二还是三四象限的Quad函数.

共线判定：由于向量叉积的意义为$\vec{a}$和$\vec{b}$为两条边构成的平行四边形面积，所以我们还可以用$\vec{a}\times \vec{b}=0$来判定向量共线，即Check_coll(a,b)函数.

左右判定：由于向量叉积求的是有向面积，我们可以用$\vec{b}\times \vec{a}>0$来判定$\vec{a}$在$\vec{b}$的逆时针方向度数为$(0,\pi )$，我们称$\vec{a}$在$\vec{b}$左边，即Check_left(a,b)函数.

向量旋转：若将向量$\vec{a}=(x_1,y_1)$旋转$\theta$度得到向量$\vec{b}=(x_2,y_2)$，那么有一个公式：
$$\left[\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right]\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right]=\left[\begin{array}{c}{x}_2\\ y_2\end{array}\right]$$

即Get_rtt(a,th)函数.

极角：以$x$轴非负半轴为始边逆时针旋转到向量$\vec{a}$，这个角度即极角可以直接用atan2(a.y,a.x)这个内置的函数求得，不过值域为$(-\pi ,\pi ]$，即Get_rad(a)函数.

夹角：若$\vec{a}$和$\vec{b}$的角度分别为$\alpha ,\beta$，则Get_rad(a,b)求得是$\beta -\alpha$，即Get_rad(a,b)函数求的是$\vec{a}$和$\vec{b}$之间的夹角.

Check_cmpxy(a,b)是比较$a<b$，以$y$为第一关键字，$x$为第二关键字.

重载的<符号是以极角为关键字的，减掉$ltl$是为了写下面的凸包.

代码如下：

const double eps=1e-8,pi=acos(-1);

int sgn(double x){return fabs(x)<=eps?0:x>eps?1:-1;}
double Rad(double th){if (th<-pi-eps) th+=2*pi;if (th>pi+eps) th-=2*pi;return th;}

struct vec{
  double x,y;
  vec(double X=0,double Y=0){x=X;y=Y;}
  vec operator + (const vec &p)const{return vec(x+p.x,y+p.y);}
  vec operator - (const vec &p)const{return vec(x-p.x,y-p.y);}
  vec operator - ()const{return vec(-x,-y);}
  double operator * (const vec &p)const{return x*p.x+y*p.y;}
  double operator % (const vec &p)const{return x*p.y-y*p.x;}
  double operator ~ ()const{return x*x+y*y;}
  bool operator == (const vec &p)const{return fabs(x-p.x)<=eps&&fabs(y-p.y)<=eps;}
  bool operator != (const vec &p)const{return fabs(x-p.x)>eps||fabs(y-p.y)>eps;}
  vec Unit(){double t=sqrt(x*x+y*y);return vec(x/t,y/t);}
  vec Norm(){return vec(-y,x);}
  bool Quad(){return y>eps||fabs(y)<=eps&&x>=-eps;}
}ltl;
vec operator * (const vec &a,const double &p){return vec(a.x*p,a.y*p);}
vec operator * (const double &p,const vec &a){return vec(a.x*p,a.y*p);}
vec operator / (const vec &a,const double &p){return vec(a.x/p,a.y/p);}
typedef vec point;
typedef vector<point> polygon;
//* -> 向量点积
//% -> 向量叉积
//~ -> 向量长度平方
//Unit -> 单位化
//Orth -> 法向量，原向量逆时针旋转pi/2
//Quad -> 一二/三四象限 1/0
//polygon -> 用若干个点存储多边形

bool Check_coll(vec a,vec b){return fabs(a%b)<=eps;}
bool Check_left(vec a,vec b){return b%a>eps;}
vec Get_rtt(vec a,double th){double s=sin(th),c=cos(th);return vec(c*a.x-s*a.y,s*a.x+c*a.y);}
double Get_rad(vec a){return atan2(a.y,a.x);}
double Get_rad(vec a,vec b){return Rad(Get_rad(b)-Get_rad(a));}
bool Check_cmpxy(vec a,vec b){return fabs(a.y-b.y)<=eps?a.x<b.x:a.y<b.y;}
//Check_coll -> 共线判定
//Check_left -> 左右判定，a在b左边返回1
//Get_rev -> a逆时针旋转th
//Get_rad -> 得到向量a对应角度，值域为(-pi,\pi]
//Check_cmpxy -> 按照y为第一关键字，x为第二关键字比较的小于

bool operator < (vec a,vec b){
  a=a-ltl;b=b-ltl;
  if (Check_coll(a,b)&&a*b>=-eps) return ~a<~b;
  return a.Quad()^b.Quad()?a.Quad()<b.Quad():Check_left(b,a);
}
//< -> 减去ltl后，从x轴负半轴开始的极角顺序

还有一份坐标均为整数的板子：

struct vec{
  LL x,y;
  vec(LL X=0,LL Y=0){x=X;y=Y;}
  vec operator + (const vec &p)const{return vec(x+p.x,y+p.y);}
  vec operator - (const vec &p)const{return vec(x-p.x,y-p.y);}
  vec operator - ()const{return vec(-x,-y);}
  LL operator * (const vec &p)const{return x*p.x+y*p.y;}
  LL operator % (const vec &p)const{return x*p.y-y*p.x;}
  LL operator ~ ()const{return x*x+y*y;}
  bool operator == (const vec &p)const{return x==p.x&&y==p.y;}
  bool operator != (const vec &p)const{return x^p.x||y^p.y;}
  vec Norm(){return vec(-y,x);}
  bool Quad(){return y>0||!y&&x>=0;}
}ltl;
vec operator * (const vec &a,const LL &p){return vec(a.x*p,a.y*p);}
vec operator * (const LL &p,const vec &a){return vec(a.x*p,a.y*p);}
vec operator / (const vec &a,const LL &p){return vec(a.x/p,a.y/p);}
typedef vec point;
typedef vector<point> polygon;

bool Check_coll(vec a,vec b){return !abs(a%b);}
bool Check_left(vec a,vec b){return b%a>0;}
bool Check_cmpxy(vec a,vec b){return a.y^b.y?a.y<b.y:a.x<b.x;}

bool operator < (vec a,vec b){
  a=a-ltl;b=b-ltl;
  if (Check_coll(a,b)&&a*b>=0) return ~a<~b;
  return a.Quad()^b.Quad()?a.Quad()<b.Quad():Check_left(b,a);
}

二.直线.

struct line{
  point a[2];
  line(point A=point(),point B=point()){a[0]=A;a[1]=B;}
  point &operator [] (const int &p){return a[p];}
  vec v(){return a[1]-a[0];}
};
//v -> 得到向量

point Get_cross(line a,line b){return a[0]+b.v()%(a[0]-b[0])/(a.v()%b.v())*a.v();}
int Get_rela(line a,line b){return Check_coll(a.v(),b.v())?!Check_coll(a[0]-b[0],a.v()):2;}
point Get_pro(point a,line b){return b[0]+b.v()*(a-b[0])*b.v()/~b.v();}
point Get_ref(point a,line b){return 2*Get_pro(a,b)-a;}
//Get_cross -> 计算交点
//Get_rela -> 判断关系0/1/2 重合/平行/有交
//Get_pro -> 求点a在直线b上的投影
//Get_ref -> 求点a关于直线b的对称点

bool Check_cross(line a,line b){
  if (Check_coll(a.v(),b.v())){
	if (!Check_coll(a[0]-b[0],a.v())) return 0;
	if (max(a[0].x,a[1].x)<min(b[0].x,b[1].x)-eps) return 0;
	if (min(a[0].x,a[1].x)>max(b[0].x,b[1].x)+eps) return 0;
	if (max(a[0].y,a[1].y)<min(b[0].y,b[1].y)-eps) return 0;
	if (min(a[0].y,a[1].y)>max(b[0].y,b[1].y)+eps) return 0;
	return 1;
  }
  if (sgn((a[0]-b[0])%(a[0]-a[1]))*sgn((a[0]-b[1])%(a[0]-a[1]))>0) return 0;
  if (sgn((b[0]-a[0])%(b[0]-b[1]))*sgn((b[0]-a[1])%(b[0]-b[1]))>0) return 0;
  return 1;
}

//Check_cross 判断线段是否相交

三.多边形.

求多边形面积：

double Get_area(polygon a){
  int n=a.size();
  if (n<3) return 0;
  double res=0;
  a.push_back(a[0]);
  for (int i=0;i<n;++i) res+=a[i]%a[i+1];
  return fabs(res)/2;
}

单组询问点在任意多边形内的判定：

bool Check_in(polygon a,point x){
  int n=a.size();
  double res=0;
  a.push_back(a[0]);
  for (int i=0;i<n;++i){
	point u=a[i]-x,v=a[i+1]-x;
	if (Check_coll(u,v)&&u*v<=eps) return 1;
	res+=Get_rad(u,v);
  }
  return fabs(res)>eps;
}

多组询问点在凸多边形内的判定：

bool Check_in(const polygon &a,point x){
  ltl=a[0];
  if (x==ltl) return 1;
  if (Check_left(x-ltl,a.back()-ltl)||Check_left(a[1]-ltl,x-ltl)) return 0;
  int p=lower_bound(a.begin(),a.end()-1,x)-a.begin();
  return !Check_left(a[p]-x,x-a[p-1])&&Check_cmpxy(ltl,x);
}

四.凸包.

凸包：

polygon Graham(polygon a){
  int n=a.size();
  polygon res;
  ltl=a[0];
  for (int i=1;i<n;++i)
	if (Check_cmpxy(a[i],ltl)) ltl=a[i];
  sort(a.begin(),a.end());
  for (int i=0;i<n;++i){
	for (;res.size()>=2&&!Check_left(a[i]-res.back(),res.back()-res[res.size()-2]);res.pop_back());
	res.push_back(a[i]);
  }
  return res;
}

旋转卡壳（平面最远点对）：

double Get_dia(polygon a){
  int n=a.size();
  double res=0;
  a.push_back(a[0]);
  for (int i=0,j=0;i<n;++i){
	for (;(a[i+1]-a[i])%(a[j]-a[i])<(a[i+1]-a[i])%(a[j+1]-a[i])-eps;j=(j+1)%n);
	res=max(res,max(~(a[j]-a[i]),~(a[j]-a[i+1])));
	res=max(res,max(~(a[j+1]-a[i]),~(a[j+1]-a[i+1])));
  }
  return sqrt(res);
}

Minkowski和：

polygon Minkowski(polygon a,polygon b){
  int n=a.size(),m=b.size();
  polygon res;
  a.push_back(a[0]);b.push_back(b[0]);
  res.push_back(ltl=a[0]+b[0]);
  for (int i=0,j=0;2333;){
	Check_left(a[i+1]-a[i],b[j+1]-b[j])?j=(j+1)%m:i=(i+1)%n;
	point t=a[i]+b[j];
	if (res.size()>=2&&Check_coll(t-res.back(),res.back()-res[res.size()-2])) res.pop_back();
	res.push_back(t);
	if (!i&&!j) break;
  }
  res.pop_back();
  return res;
}

半平面交：

bool cmp(line a,line b){return Check_coll(a.v(),b.v())&&a.v()*b.v()>eps?Check_left(a[0]-b[0],b.v()):Check_cmp(a.v(),b.v());}

polygon HPI(vector<line>a){
  int n=a.size(),hd,tl;
  polygon res,qk(n+1);
  vector<line>q(n+1);
  sort(a.begin(),a.end(),cmp);
  q[hd=tl=0]=a[0];
  for (int i=1;i<n;++i){
	if (!Check_cmp(q[tl].v(),a[i].v())) continue;
	for (;hd<tl&&Check_left(a[i].v(),qk[tl-1]-a[i][0]);--tl);
	for (;hd<tl&&Check_left(a[i].v(),qk[hd]-a[i][0]);++hd);
	qk[tl]=Get_cross(a[i],q[tl]);q[++tl]=a[i];
  }
  for (;hd<tl&&Check_left(q[hd].v(),qk[tl-1]-q[hd][0]);--tl);
  for (;hd<tl&&Check_left(q[tl].v(),qk[hd]-q[tl][0]);++hd);
  qk[tl]=Get_cross(q[hd],q[tl]);
  return polygon(qk.begin()+hd,qk.begin()+tl+1);
}

五.圆.

圆：

struct circle{
  point o;
  double r;
  circle(point O=point(),double R=0){o=O;r=R;}
  point Point(double th){return point(o.x+cos(th)*r,o.y+sin(th)*r);}
};
//Point -> 从x非负半轴开始逆时针th对应的点

三点定圆（如果三点共线返回两个最远点为直径端点的圆）：

circle Get_circle(point a,point b,point c){
  if (Check_coll(a-b,b-c)){
    if (Check_cmpxy(a,b)) swap(a,b);
    if (Check_cmpxy(a,c)) swap(a,c);
    if (Check_cmpxy(b,c)) swap(b,c);
    return circle((a+c)/2,sqrt(~(a-c))/2);
  }
  point t=(a+b)/2;
  line t0=line(t,t+(a-b).Norm());
  t=(b+c)/2;
  line t1=line(t,t+(b-c).Norm());
  t=Get_cross(t0,t1);
  return circle(t,sqrt(~(a-t)));
}

直线与圆的交点：

vector<point>Get_cross(line a,circle b){
  point t=Get_pro(b.o,a);
  double d2=~(t-b.o),d=sqrt(d2);
  if (b.r<d-eps) return vector<point>();
  if (fabs(b.r-d)<=eps) vector<point>(1,t);
  vector<point>res;
  d=sqrt(b.r*b.r-d2);
  res.push_back(t+a.v().Unit()*d);
  res.push_back(t-a.v().Unit()*d);
  return res;
}

两圆的交点（必须保证两个圆不重合）：

vector<point>Get_cross(circle a,circle b){
  double d2=~(a.o-b.o),d=sqrt(d2);
  if (a.r+b.r<d-eps||fabs(a.r-b.r)>d+eps) return vector<point>();
  if (fabs(a.r+b.r-d)<=eps||fabs(fabs(a.r-b.r)-d)<=eps)
	return vector<point>(1,a.o+(b.o-a.o).Unit()*a.r);
  vector<point>res;
  double x=acos((d2+a.r*a.r-b.r*b.r)/(2*d*a.r)),y=Get_rad(b.o-a.o);
  res.push_back(a.Point(y+x));
  res.push_back(a.Point(y-x));
  return res;
}

过定点求圆的切线（终点在圆上）：

vector<line>Get_tan(point a,circle b){
  double d=~(a-b.o);
  if (b.r*b.r>d+eps) return vector<line>();
  if (fabs(b.r*b.r-d)<=eps) return vector<line>(1,line(a+(b.o-a).Norm(),a));
  vector<line>res;
  d=b.r/sqrt(d);
  double th=acos(d);
  res.push_back(line(a,b.o+Get_rtt(a-b.o,th)*d));
  res.push_back(line(a,b.o+Get_rtt(a-b.o,-th)*d));
  return res;
}

求两圆内外公切线（必须保证两个圆不重合，返回的切线起点在$a$ 上）：

vector<line>Get_extan(circle a,circle b){
  if (a.r<=eps&&b.r<=eps) return vector<line>(1,line(a.o,b.o));
  vector<line>res;
  if (fabs(a.r-b.r)<=eps){
	vec t=(b.o-a.o).Norm().Unit()*a.r;
	res.push_back(line(a.o+t,b.o+t));
	res.push_back(line(a.o-t,b.o-t));
	return res;
  }
  int flag=0;
  if (b.r<a.r-eps) swap(a,b),flag=1;
  res=Get_tan(a.o,circle(b.o,b.r-a.r));
  for (int vs=res.size(),i=0;i<vs;++i){
	vec t=(res[i][1]-b.o).Unit()*a.r;
	res[i]=line(res[i][0]+t,res[i][1]+t);
	if (flag) swap(res[i][0],res[i][1]);
  }
  if (!flag&&res.size()==1) swap(res[0][0],res[0][1]);
  return res;
}

vector<line>Get_intan(circle a,circle b){
  if (a.r<=eps&&b.r<=eps) return vector<line>(1,line(a.o,b.o));
  vector<line>res=Get_tan(a.o,circle(b.o,a.r+b.r));
  for (int vs=res.size(),i=0;i<vs;++i){
	vec t=(res[i][1]-b.o).Unit()*a.r;
	res[i]=line(res[i][0]-t,res[i][1]-t);
  }
  if (res.size()==1) swap(res[0][0],res[0][1]);
  return res;
}

最小圆覆盖：

circle Get_mincircle(vector<point>a){
  int n=a.size();
  for (int i=0;i<5;++i) random_shuffle(a.begin(),a.end());
  circle res=circle(a[0],0);
  for (int i=1;i<n;++i){
    if (~(a[i]-res.o)<=res.r*res.r+eps) continue;
    res=circle(a[i],0);
    for (int j=0;j<i;++j)
      if (~(a[j]-res.o)>res.r*res.r+eps){
        res=circle((a[i]+a[j])/2,sqrt(~(a[i]-a[j]))/2);
        for (int k=0;k<j;++k)
          if (~(a[k]-res.o)>res.r*res.r+eps) res=Get_circle(a[i],a[j],a[k]);
      }
  }
  return res;
}

六.simpson积分.

simpson积分法：

double Get_integ(double l,double r){return (r-l)*(Get_f(l)+4*Get_f((l+r)/2)+Get_f(r))/6;}

double Simpson(double l,double r,double s){
  double mid=(l+r)/2,a=Get_integ(l,mid),b=Get_integ(mid,r);
  return fabs(a+b-s)<=eps?s+(a+b-s)/15:Simpson(l,mid,a)+Simpson(mid,r,b);
}