一、常量&约定
为了方便地比较浮点数，我们需要设计一个三态比较函数，如下。
#define IL inline
const double pi=acos(-1.0);
const double eps=1e-8;
IL int sign(double x){
if(fabs(x)<eps)return 0;
return (x<0)?-1:1;
}
IL int dcmp(double x,double y){return sign(x-y);}
另外，除非特殊说明，所有涉及旋转/角度的操作，都默认逆时针为正。
二、点&向量
2.1 加、减、数乘
在计算几何中，由于三角函数运算的值不够精确，我们一般采用坐标运算，我们首先需要重载最基本的坐标运算。
struct Point{
double x,y;
Point(double X=0,double Y=0):x(X),y(Y){}
inline void in(){scanf("%lf%lf",&x,&y);}
inline void out(){printf("%.2lf %.2lf
",x,y);}
Point operator + (const Point &r) const{return Point(x+r.x,y+r.y);}
Point operator - (const Point &r) const{return Point(x-r.x,y-r.y);}
Point operator * (double p) const{return Point(x*p,y*p);}
Point operator / (double p) const{return Point(x/p,y/p);}
bool operator < (const Point &r) const{return x<r.x||(x==r.x&&y<r.y);}
bool operator == (const Point &r) const{return sign(x-r.x)==0&&sign(y-r.y)==0;}
};
typedef Point Vector;
2.2 点积
(oldsymbol{a}cdot oldsymbol{b}=|oldsymbol{a}||oldsymbol{b}|cos heta( heta=leftlangle oldsymbol{a},oldsymbol{b} ight angle)=x_1 x_2 +y_1 y_2)
IL double dot(Vector A,Vector B){return A.x*B.x+A.y*B.y;}
利用点积，我们可以比较方便地求出向量的模长，两点间的距离，向量的夹角，以及一个向量在另一个向量上的投影。
IL double length(Vector A){return sqrt(dot(A,A));}
IL double dis(Point A,Point B){return length(A-B);}
IL double angle(Vector A,Vector B){return acos(dot(A,B)/length(A)/length(B));}
IL double project(Point a,Point b,Point c){return dot(b-a,c-a)/length(b-a);}
2.3 叉积
(oldsymbol{a} imes oldsymbol{b}=|oldsymbol{a}||oldsymbol{b}|sin heta( heta=leftlangle oldsymbol{a},oldsymbol{b} ight angle)=x_1 y_2 -x_2 y_1)
IL double cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}
叉积的结果可以表示两个向量(oldsymbol a)，(oldsymbol b)共起点时，所构成平行四边形的有向面积，这个性质会在之后的许多操作中得以应用。此外，有向面积的正负可以帮助我们判断两个向量共起点以后的位置关系。
具体而言，如果(oldsymbol a)在(oldsymbol b)左侧，那么(oldsymbol{a} imes oldsymbol{b}<0)
IL double area(Point a,Point b,Point c){return cross(b-a,c-a);}
2.4 极角
利用C++函数库中的(mathtt{atan2(y,x)})可以帮助我们计算(arctan frac{y}{x})的值，其值域为((-pi,pi])，即当点位于三四象限时其返回一个负值。
IL double polar_angle(Vector A){return atan2(A.y,A.x);}
2.5 旋转
[egin{bmatrix}x yend{bmatrix} imesegin{bmatrix}cos heta & sin heta\ -sin heta &cos hetaend{bmatrix}
]
IL Vector rotate(Vector A,double rad){
return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}
三、直线&线段
3.1 直线的表示&参数
我们常用点向式（直线的参数方程）来表示一个直线：
[l:a+oldsymbol{v}t
]
为了方便，我们提供两种构造直线的方式（传入两点，或者传入一点和一个向量），并根据需要记录(a,b,oldsymbol{v},r)(极角)。
struct Line{
Point a,b;
Vector v;
double r;
Line(){}
Line(const Point &A,const Point &B,int op=0){
if(!op){a=A,b=B;v=b-a;r=polar_angle(v);}
else a=A,v=B;b=a+v;r=polar_angle(v);
}
Point point(double p){return a+v*p;}
};
3.2 点与直线（线段）的位置关系
如果一个点(p)满足(vec{ap} imes oldsymbol{v}=0)，我们就可以判定(p)在直线上。
如果我们需要判定点在线段上，还需要额外满足(vec{pa}cdot vec{pb}<0)
bool line_inc(const Point &p) const{
return sign(cross(p-a,v))==0;
}
bool seg_inc(const Point &p) const{
return sign(cross(a-p,b-p))==0&&sign(dot(a-p,b-p))<=0;
}
3.3 点在直线上的投影点(垂足)
用2.2中提到的方法计算，记得考虑要使用(oldsymbol{v})对应的单位向量。
IL Point point_line_proj(Point p,Line l){
return l.a+l.v*(dot(l.v,p-l.a)/dot(l.v,l.v));
}
3.4 点到直线（线段）的距离
我们采取面积/底边的方式计算点到直线的距离（高）。
如果要计算点到线段的距离，需要先判断高是否可以落在线段上，如果不在，应该取返回点与最近端点的距离。
IL double point_to_line(Point p,Line l){
return fabs(cross(l.v,p-l.a))/length(l.v);
}
IL double point_to_seg(Point p,Seg s){
if(s.a==s.b)return length(p-s.a);
Vector v1=s.b-s.a,v2=p-s.a,v3=p-s.b;
if(sign(dot(v1,v2))<0)return length(v2);
if(sign(dot(v1,v3))>0)return length(v3);
return point_to_line(p,Line(s.a,s.b));
}
3.5 直线（线段）间的位置关系
我们利用面积比来得到直线与直线的交点。
值得注意的是，为了防止被除数为零，我们需要预先判断两直线是否平行。
IL Point line_line_inter(Line x,Line y){//记得先判平行
Vector U=x.a-y.a;
double t=cross(y.v,U)/cross(x.v,y.v);//注意方向
return x.a+x.v*t;
}
如果两线段相交，则两线段必然相互跨立对方。我们通过跨立试验来判断，两线段是否相交（直线与线段是否相交类似）。我们同样需要预先判断两直线是否平行。
IL bool seg_seg_inter(Seg x,Seg y){
Point a1=x.a,a2=x.b,b1=y.a,b2=y.b;
double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1),
c3=cross(b2-b1,a1-b1),c4=cross(b2-b1,a2-b1);
return sign(c1)*sign(c2)<=0&&sign(c3)*sign(c4)<=0;
}
IL bool line_seg_inter(Line x,Seg y){
Point a1=x.a,a2=x.b,b1=y.a,b2=y.b;
double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1);
return sign(c1)*sign(c2)<=0;
}
四、圆
4.1 圆的表示&参数
我们可以使用圆心和半径表示一个圆。
struct Circle{
Point o;double r;
Circle(){}
Circle(const Point &O,double R):o(O),r(R){}
}
4.2 点与圆的位置关系
计算点与圆心的距离然后判断即可。
bool inc(const Point &p){
return dcmp(r,dis(o,p))>=0;
}
4.3 求三角形的外接圆
以任意两边中垂线的交点为圆心，交点到一点处的距离为半径作圆。
IL Line get_line(Point a,Point b){return Line((a+b)/2,rotate(b-a,pi/2),1);};
IL Circle get_circle(Point a,Point b,Point c){
Line u=get_line(a,b),v=get_line(a,c);
Point o=line_line_inter(u,v);
return Circle(o,dis(o,a));
}
4.4 最小圆覆盖
如果点(p)不在点集(S)的最小覆盖圆内，则(p)一定在({S}cup{p}) 的最小覆盖圆上。
利用这个定理，我们可以在(O(n))时间复杂度内求出一个点集的最小覆盖圆。
Circle min_circle(Point p[],int n){
random_shuffle(p,p+n);
Circle c=Circle(p[0],0);
for(int i=1;i<n;i++)
if(!c.inc(p[i])){
c=Circle(p[i],0);
for(int j=0;j<i;j++)
if(!c.inc(p[j])){
c=Circle((p[i]+p[j])/2,dis(p[i],p[j])/2);
for(int k=0;k<j;k++)
if(!c.inc(p[k]))
c=get_circle(p[i],p[j],p[k]);
}
}
return c;
}
4.5 伸缩变换
对于椭圆问题，我们可以对整个坐标系进行伸缩变换，再进行计算。
五、普通多边形
5.1 多边形的面积
我们用三角剖分求一个一般多边形的面积。
double poly_area(Point p[],int n){
double res=0;
for(int i=1;i<n-1;i++)
res+=cross(p[i]-p[0],p[i+1]-p[0]);
return res/2;
}
5.2 点与多边形的位置关系
我们采用射线法：考虑以该点为端点引出一条射线，如果这条射线与多边形有奇数个交点，则该点在多边形内部，否则该点在多边形外部，这被称为奇偶规则。
六、凸多边形
6.1 Andrew算法求凸包
我们采用Andrew算法可以(O(n))地计算一个点集的凸包。
int top,st[N];
bool used[N];
void Andrew(Point p[],int cnt){
sort(p+1,p+1+cnt);
for(int i=1;i<=cnt;i++){
while(top>1&&area(p[st[top-1]],p[st[top]],p[i])<=0)used[st[top--]]=0;
st[++top]=i;used[i]=1;
}
used[1]=0;
for(int i=cnt;i;i--){
if(used[i])continue;
while(top>1&&area(p[st[top-1]],p[st[top]],p[i])<=0)top--;
st[++top]=i;
}
}
6.2 点与凸多边形的位置关系
只要一个点在所有边的左侧，点就在凸包内，这可以遍历一遍实现或者用二分实现。
6.3 旋转卡壳
旋转卡壳可以用于求凸包的直径。
double rotating_caliper(Point p[]){
double res=0;
for(int i=1,a=3;i<top;i++){
Point d=p[st[i]],e=p[st[i+1]];
while(dcmp(area(d,e,p[st[a]]),area(d,e,p[st[a+1]]))<0)a=a%(top-1)+1;
res=max(res,max(dis(d,p[st[a]]),dis(e,p[st[a]])));
}
return res;
}
至于求平面最近点对，这不是个计算几何问题，但也放在这里，这个问题除了放在下面的解法外，还有分治的解法，可以去OI-wiki上查看。
struct cmp_y{
bool operator()(const Point &a, const Point &b)const{return a.y<b.y;}
};
multiset<Point,cmp_y>s;
double GetMin(){
double res=1e16;
for(int i=1,l=1;i<=n;i++){
while(l<i&&dcmp(p[i].x-p[l].x,res)>=0)s.erase(s.find(p[l++]));
for(auto it=s.lower_bound(Point(0,p[i].y-res));
it!=s.end()&&dcmp(it->y-p[i].y,res)<0;it++)
if(dcmp(dis(*it,p[i]),res)<0)res=dis(*it,p[i]);
s.insert(p[i]);
}
return res;
}
另外，平面最近最远点对问题有一个随机化的贪心算法：将整个坐标系随机旋转一个角度，再贪心计算。
6.4 半平面交
计算每条直线左侧的平面的交集。
bool cmp(const Line& x,const Line& y){
if(sign(x.r-y.r)==0)return area(x.a,x.b,y.b)<0;
return x.r<y.r;
}
IL bool on_right(Line a,Line b,Line c){
Point o=line_line_inter(b,c);
return sign(area(a.a,a.b,o))<=0;
}
double half_plane_inter(Line l[],int cnt){
sort(l+1,l+1+cnt,cmp);
int L=0,R=-1;
vector<int>q(cnt);
for(int i=1;i<=cnt;i++){
if(i>1&&sign(l[i].r-l[i-1].r)==0)continue;
while(L+1<=R&&on_right(l[i],l[q[R-1]],l[q[R]]))R--;
while(L+1<=R&&on_right(l[i],l[q[L]],l[q[L+1]]))L++;
q[++R]=i;
}
while(L+1<=R&&on_right(l[q[L]],l[q[R-1]],l[q[R]]))R--;
while(L+1<=R&&on_right(l[q[R]],l[q[L]],l[q[L+1]]))L++;
q[++R]=q[L];
vector<Point>ans(R);
int k=0;
for(int i=L;i<R;i++)
ans[++k]=line_line_inter(l[q[i]],l[q[i+1]]);
double res=0;
for(int i=2;i+1<=k;i++)res+=area(ans[1],ans[i],ans[i+1]);
return res/2;
}
七、面积并的计算
7.1 扫描线
计算图形的面积并时，一种常用的思路是运用扫描线——将一维取出排序，然后挨个统计对应的另一维的答案。
7.2 自适应辛普森法
总体思路是根据辛普森公式不断地近似计算函数，直到得到精度合适的值。
辛普森公式：
[int_{a}^{b} f(x) dx approx frac{b-a}{6}[f(a)+4f(frac{a+b}{2})+f(b)]
]
代码模板如下：
IL double f(double x){
}
IL double simpson(double a,double b){
return ((b-a)*(f(a)+f(b)+4*f((a+b)/2)))/6;
}
double divide(double L,double R,double ans){
double mid=(L+R)/2,l=simpson(L,mid),r=simpson(mid,R);
if(dcmp(l+r,ans)==0)return ans;
return divide(L,mid,l)+divide(mid,R,r);
}
如果本题的(f (x))计算较慢，需要根据题目需要将已经计算过的函数值作为参数传下去，避免重复计算。
如果精度不好保证，还需要将精度不断衰减保证效率。
附 :【模板】计算几何全家桶
#include<bits/stdc++.h>
using namespace std;
#define IL inline
const int N=1e5+5;
const double pi=acos(-1.0);
const double eps=1e-8;
IL int sign(double x){
if(fabs(x)<eps)return 0;
return (x<0)?-1:1;
}
IL int dcmp(double x,double y){return sign(x-y);}
struct Point{
double x,y;
Point(double X=0,double Y=0):x(X),y(Y){}
inline void in(){scanf("%lf%lf",&x,&y);}
inline void out(){printf("%.2lf %.2lf
",x,y);}
Point operator + (const Point &r) const{return Point(x+r.x,y+r.y);}
Point operator - (const Point &r) const{return Point(x-r.x,y-r.y);}
Point operator * (double p) const{return Point(x*p,y*p);}
Point operator / (double p) const{return Point(x/p,y/p);}
bool operator < (const Point &r) const{return x<r.x||(x==r.x&&y<r.y);}
bool operator == (const Point &r) const{return sign(x-r.x)==0&&sign(y-r.y)==0;}
};
typedef Point Vector;
IL double dot(Vector A,Vector B){return A.x*B.x+A.y*B.y;}
IL double cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}
IL double length(Vector A){return sqrt(dot(A,A));}
IL double dis(Point A,Point B){return length(A-B);}
IL double angle(Vector A,Vector B){return acos(dot(A,B)/length(A)/length(B));}
IL double project(Point a,Point b,Point c){return dot(b-a,c-a)/length(b-a);}
IL double polar_angle(Vector A){return atan2(A.y,A.x);}
IL double area(Point a,Point b,Point c){return cross(b-a,c-a);}
IL Vector rotate(Vector A,double rad){
return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));
}
struct Line{
Point a,b;
Vector v;
double r;
Line(){}
Line(const Point &A,const Point &B,int op=0){
if(!op){a=A,b=B;v=b-a;r=polar_angle(v);}
else a=A,v=B;b=a+v;r=polar_angle(v);
}
Point point(double p){return a+v*p;}
bool line_inc(const Point &p) const{return sign(cross(p-a,v))==0;}
bool seg_inc(const Point &p) const{
return sign(cross(a-p,b-p))==0&&sign(dot(a-p,b-p))<=0;
}
};
typedef Line Seg;
IL Point point_line_proj(Point p,Line l){
return l.a+l.v*(dot(l.v,p-l.a)/dot(l.v,l.v));
}
IL double point_to_line(Point p,Line l){
return fabs(cross(l.v,p-l.a))/length(l.v);
}
IL double point_to_seg(Point p,Seg s){
if(s.a==s.b)return length(p-s.a);
Vector v1=s.b-s.a,v2=p-s.a,v3=p-s.b;
if(sign(dot(v1,v2))<0)return length(v2);
if(sign(dot(v1,v3))>0)return length(v3);
return point_to_line(p,Line(s.a,s.b));
}
IL Point line_line_inter(Line x,Line y){//记得先判平行
Vector U=x.a-y.a;
double t=cross(y.v,U)/cross(x.v,y.v);//注意方向
return x.a+x.v*t;
}
IL bool seg_seg_inter(Seg x,Seg y){
Point a1=x.a,a2=x.b,b1=y.a,b2=y.b;
double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1),
c3=cross(b2-b1,a1-b1),c4=cross(b2-b1,a2-b1);
return sign(c1)*sign(c2)<=0&&sign(c3)*sign(c4)<=0;
}
IL bool line_seg_inter(Line x,Seg y){
Point a1=x.a,a2=x.b,b1=y.a,b2=y.b;
double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1);
return sign(c1)*sign(c2)<=0;
}
struct Circle{
Point o;double r;
Circle(){}
Circle(const Point &O,double R):o(O),r(R){}
bool inc(const Point &p){
return dcmp(r,dis(o,p))>=0;
}
};
IL Line get_line(Point a,Point b){return Line((a+b)/2,rotate(b-a,pi/2),1);};
IL Circle get_circle(Point a,Point b,Point c){
Line u=get_line(a,b),v=get_line(a,c);
Point o=line_line_inter(u,v);
return Circle(o,dis(o,a));
}
Circle min_circle(Point p[],int n){
random_shuffle(p,p+n);
Circle c=Circle(p[0],0);
for(int i=1;i<n;i++)
if(!c.inc(p[i])){
c=Circle(p[i],0);
for(int j=0;j<i;j++)
if(!c.inc(p[j])){
c=Circle((p[i]+p[j])/2,dis(p[i],p[j])/2);
for(int k=0;k<j;k++)
if(!c.inc(p[k]))
c=get_circle(p[i],p[j],p[k]);
}
}
return c;
}
double poly_area(Point p[],int n){
double res=0;
for(int i=1;i<n-1;i++)
res+=cross(p[i]-p[0],p[i+1]-p[0]);
return res/2;
}
int top,st[N];
void Andrew(Point p[],int cnt){
vector<bool>used(cnt);
sort(p+1,p+1+cnt);
for(int i=1;i<=cnt;i++){
while(top>1&&area(p[st[top-1]],p[st[top]],p[i])<=0)used[st[top--]]=0;
st[++top]=i;used[i]=1;
}
used[1]=0;
for(int i=cnt;i;i--){
if(used[i])continue;
while(top>1&&area(p[st[top-1]],p[st[top]],p[i])<=0)top--;
st[++top]=i;
}
}
double rotating_caliper(Point p[]){
double res=0;
for(int i=1,a=3;i<top;i++){
Point d=p[st[i]],e=p[st[i+1]];
while(dcmp(area(d,e,p[st[a]]),area(d,e,p[st[a+1]]))<0)a=a%(top-1)+1;
res=max(res,max(dis(d,p[st[a]]),dis(e,p[st[a]])));
}
return res;
}
bool cmp(const Line& x,const Line& y){
if(sign(x.r-y.r)==0)return area(x.a,x.b,y.b)<0;
return x.r<y.r;
}
IL bool on_right(Line a,Line b,Line c){
Point o=line_line_inter(b,c);
return sign(area(a.a,a.b,o))<=0;
}
double half_plane_inter(Line l[],int cnt){
sort(l+1,l+1+cnt,cmp);
int L=0,R=-1;
vector<int>q(cnt);
for(int i=1;i<=cnt;i++){
if(i>1&&sign(l[i].r-l[i-1].r)==0)continue;
while(L+1<=R&&on_right(l[i],l[q[R-1]],l[q[R]]))R--;
while(L+1<=R&&on_right(l[i],l[q[L]],l[q[L+1]]))L++;
q[++R]=i;
}
while(L+1<=R&&on_right(l[q[L]],l[q[R-1]],l[q[R]]))R--;
while(L+1<=R&&on_right(l[q[R]],l[q[L]],l[q[L+1]]))L++;
q[++R]=q[L];
vector<Point>ans(R);
int k=0;
for(int i=L;i<R;i++)
ans[++k]=line_line_inter(l[q[i]],l[q[i+1]]);
double res=0;
for(int i=2;i+1<=k;i++)res+=area(ans[1],ans[i],ans[i+1]);
return res/2;
}