计算几何
模板
cpp
namespace Geometry{
const int N=1e5+10;
const double PI=acos(-1);
const double eps=1e-12;
const double INF=1e18;
int sign(double x){//符号函数：精度为eps
if(fabs(x)<=eps) return 0;//x为0
if(x<0) return -1;//x为负数
return 1;//x为正数
}
struct Point{//点坐标
double x,y;
int id;
Point(int x=0,int y=0):x(x),y(y){}
Point operator-(Point p){return Point(x-p.x,y-p.y);}
Point operator+(Point p){return Point(x+p.x,y+p.y);}
Point operator*(double k){return Point(x*k,y*k);}
Point operator/(double k){return Point(x/k,y/k);}
double dis(){return x*x+y*y;}
int dis(int op){return x*x+y*y;}
Point norm(){return *this/length();}
Point vector(){return Point(sign(this->x),sign(this->y));}
double length(){return sqrt(dis());}
bool operator<(Point p)const{return x==p.x?y<p.y:x<p.x;}
bool operator==(Point p)const{return sign(x-p.x)==0&&sign(y-p.y)==0;}
bool operator!=(Point p)const{return sign(x-p.x)==0&&sign(y-p.y)==0;}
friend ostream& operator<<(ostream &os, const Point &p){
os<<"("<<p.x<<" "<<p.y<<")";
return os;
}
};
typedef Point Vector;
int cmp(double x,double y){//比较函数：精度为eps
if(fabs(x-y)<=eps) return 0;//x与y相等
if(x<y) return -1;//x小于y
return 1;//x大于y
}
double dot(Point a,Point b){//向量点积
return a.x*b.x+a.y*b.y;
}
double cross(Point a,Point b){//向量叉积
return a.x*b.y-b.x*a.y;
}
double get_angle(Point a,Point b){//求A->B的有向角
return acos(dot(a,b)/a.length()/b.length());
}
double area(Point a,Point b,Point c){//A为顶点，向量AB与向量AC的叉积，即三角形ABC面积的2倍（有向）
return cross(b-a,c-a);
}
Point rotate(Point a,double angle){
return Point(a.x*cos(angle)+a.y*sin(angle),-a.x*sin(angle)+a.y*cos(angle));
}
bool on_segment(Point p,Point a,Point b){//点p是否在线段ab上
return sign(cross(p-a,p-b))==0&&sign(dot(p-a,p-b))<=0;
}
Point get_line_intersection(Point p,Vector v,Point q,Vector w){//两直线的交点
//使用前提，直线必须有交点
//cross(v, w) == 0则两直线平行或者重合
if(!sign(cross(v,w))) return {INF,INF};
Vector u=p-q;
double t=cross(w,u)/cross(v,w);
Point o=p+v*t;
auto on_segment=[&](Point p,Point a,Point b){
return sign(dot(p-a,p-b))<=0;
};
//若判断线段/射线是否相交 还需判断交点是否在线段/射线上
if(!on_segment(o,p,p+v)||!on_segment(o,q,q+w))return {INF,INF};
return o;
}
double distance_to_line(Point p,Point a,Point b){//点p到直线ab的距离
Vector v1=b-a,v2=p-a;
return fabs(cross(v1,v2)/v1.length());
}
double distance_to_segment(Point p,Point a,Point b){//点p到线段ab的距离
if(a==b) return (p-a).length();
Vector v1=b-a,v2=p-a,v3=p-b;
if(sign(dot(v1,v2))<0) return v2.length();
if(sign(dot(v1,v3))>0) return v3.length();
return distance_to_line(p,a,b);
}
Point get_line_projection(Point p,Point a,Point b){//点p在直线ab上的投影
Vector v=b-a;
return a+v*(dot(v,p-a)/dot(v,v));
}
// bool segment_intersection(Point a1,Point a2,Point b1,Point b2){// 判断两个线段a1a2与b1b2是否相交
// //在同一直线上但不重合的两线段会错判为相交
// double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1);
// double c3=cross(b2-b1,a2-b1),c4=cross(b2-b1,a1-b1);
// return sign(c1)*sign(c2)<=0&&sign(c3)*sign(c4)<=0;
// }
Point segment_intersection(Point a1,Point a2,Point b1,Point b2){// 判断两个线段a1a2与b1b2是否相交
double c1=cross(a2-a1,b1-a1),c2=cross(a2-a1,b2-a1);
double c3=cross(b2-b1,a2-b1),c4=cross(b2-b1,a1-b1);
if(!(sign(c1)*sign(c2)<=0&&sign(c3)*sign(c4)<=0))return {INF,INF};
Point da=a2-a1,db=b2-b1;
double d=cross(da,db);
if(fabs(d)<=eps)return {INF,INF};
double t=cross(b1-a1,db)/d;
double u=cross(b1-a1,da)/d;
if(t>=-eps&&t<=1+eps&&u>=-eps&&u<=1+eps){
return a1+da*t;
}
return {INF,INF};
}
double polygon_area(vector<Point> &p){//求多边形面积(不一定是凸多边形)
double s=0;
for(int i=1;i+1<p.size();i++)
s+=cross(p[i]-p[0],p[i+1]-p[i]);
return s/2;
}
bool in_line(Point p1,Point p2,Point p3){ //判断三点是否共线
double res=(p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y);
if(fabs(res)<=eps)return true;//共线
return false;//不共线
}
bool is_line(vector<Point> &p){
for(int i=1;i<p.size();i++){
if(!in_line(p[0],p[1],p[i]))return false;
}
return true;
}
int stk[N],top;
bool used[N];
vector<Point> andrew(vector<Point> &q){// 求凸包
sort(q.begin(),q.end());
int n=q.size();
top=0;
for(int i=0;i<n;i++)used[i]=false;
for(int i=0;i<n;i++){
// 如果共线合法，改为<
while(top>=2&&sign(area(q[stk[top-1]],q[stk[top]],q[i]))<=0){
if(sign(area(q[stk[top-1]],q[stk[top]],q[i])<0))used[stk[top--]]=false;
else top--;
}
stk[++top]=i;
used[i]=true;
}
used[0]=false;
for(int i=n-1;i>=0;i--){
if(used[i])continue;
while(top>=2&&sign(area(q[stk[top-1]],q[stk[top]],q[i]))<=0)top--;
stk[++top]=i;
}
vector<Point> hull;
for(int i=1;i<top;i++)hull.push_back(q[stk[i]]);
return hull;
// 凸包的点为逆时针
}
// -------------------求凸包直径--------------------------------
int rotating_calipers(vector<Point> &hull){//旋转卡壳求凸包直径平方
if(is_line(hull))return (hull[0]-hull.back()).dis();
int res=0;
for(int i=0,j=2;i<hull.size();i++){
auto d=hull[i],e=hull[(i+1)%hull.size()];
while(area(d,e,hull[j])<area(d,e,hull[(j+1)%hull.size()]))j=(j+1)%hull.size();
res=max({res,(d-hull[j]).dis(0),(e-hull[j]).dis(0)});
}
return res;
}
//-------求最小矩形覆盖-------------------
double get_project(Point a, Point b, Point c){
return dot(b - a, c - a) / (b - a).length();
}
vector<Point> get_min_rectangle(vector<Point> hull,double &min_area){//旋转卡壳求最小矩形覆盖
if(is_line(hull)){
min_area=0;
return {hull[0],hull[0],hull.back(),hull.back()};
}
// andrew判断需要改成<=
vector<Point> ans(4);
int n=hull.size();
hull.push_back(hull[0]);
for(int i=0,a=2,b=1,c=2;i<n;i++){
auto d=hull[i],e=hull[i+1];
while(cmp(area(d,e,hull[a]),area(d,e,hull[a+1]))<0)a=(a+1)%n;
while(cmp(get_project(d,e,hull[b]),get_project(d,e,hull[b+1]))<0)b=(b+1)%n;
if(!i)c=a;
while(cmp(get_project(d,e,hull[c]),get_project(d,e,hull[c+1]))>0)c=(c+1)%n;
auto x=hull[a],y=hull[b],z=hull[c];
auto h=area(d,e,x)/(e-d).length();
auto w=dot(y-z,e-d)/(e-d).length();
if(h*w<min_area){
min_area=h*w;
ans[0] = d + (e - d).norm() * get_project(d, e, y);
ans[3] = d + (e - d).norm() * get_project(d, e, z);
auto u = (rotate(e - d, -PI / 2)).norm();
ans[1] = ans[0] + u * h;
ans[2] = ans[3] + u * h;
}
}
for(int i=0;i<4;i++){
if(!sign(ans[i].x))ans[i].x=0;
if(!sign(ans[i].y))ans[i].y=0;
}
return ans;
}
struct Line{//存储凸包上每条线段的起点和终点
Point st,ed;
Line(Point _st=Point(),Point _ed=Point()):st(_st),ed(_ed){}
};
double line_get_angle(const Line& a){
return atan2(a.ed.y-a.st.y,a.ed.x-a.st.x);
}
bool line_cmp(const Line& a, const Line& b){
double A=line_get_angle(a),B=line_get_angle(b);
if (!cmp(A,B)) return area(a.st,a.ed,b.ed)<0;
return A<B;
}
Point get_line_intersection(Line a,Line b){
return get_line_intersection(a.st,a.ed-a.st,b.st,b.ed-b.st);
}
bool on_right(Line& a, Line& b, Line& c){//判断bc的交点是否在a的右侧
auto o = get_line_intersection(b,c);
return sign(area(a.st,a.ed,o))<=0;
//需要保留经过半平面交的所有直线的时候，将该处等号去掉
}
int dq[N];// 双端队列
double half_plane_intersection(vector<Line> &line){
sort(line.begin(),line.end(),line_cmp);//对每条直线按照极坐标进行排序
int hh=0,tt=-1;//定义双端队列对头和队尾元素
for(int i=0;i<line.size();i++)//遍历每一条直线
{
if (i&&!cmp(line_get_angle(line[i]),line_get_angle(line[i-1]))) continue;//当两直线极坐标角度相同
while(hh+1<=tt&&on_right(line[i],line[dq[tt-1]],line[dq[tt]])) tt--;
//如果line[l[tt-1]]和line[l[tt]]的交点在line[i]的右侧,删除队尾元素
while(hh+1<=tt&&on_right(line[i],line[dq[hh]],line[dq[hh+1]])) hh++;
//如果line[l[hh]]和line[l[hh+1]]的交点在line[i]的右侧,删除队头元素
dq[++tt]=i;//将该直线加入队尾
}
//判断一下第一条直线与最后一条直线
while (hh+1<=tt&&on_right(line[dq[hh]],line[dq[tt-1]],line[dq[tt]])) tt--;//队尾更新队头
while (hh+1<=tt&&on_right(line[dq[tt]],line[dq[hh]],line[dq[hh+1]])) hh++;//队头更新队尾
dq[++tt]=dq[hh];//将第一条直线再次加入双端队列末尾
//计算半平面交交点
vector<Point> ans;
for (int i=hh;i<tt;i++){ //对于半平面交上的直线两两之间求出交点
ans.push_back(get_line_intersection(line[dq[i]],line[dq[i+1]]));
}
//计算半平面交面积
return polygon_area(ans);
}
// --------------------求三角形并(n^3logn)---------------------
double line_area(double a, int side,vector<vector<Point>> &tr){
vector<Point> q;
for (int i = 0; i < tr.size(); i ++ ){
auto t = tr[i];
if (cmp(t[0].x, a) > 0 || cmp(t[2].x, a) < 0) continue;
if (!cmp(t[0].x, a) && !cmp(t[1].x, a)){
if (side) q.emplace_back(t[0].y, t[1].y);
}
else if (!cmp(t[2].x, a) && !cmp(t[1].x, a)){
if (!side) q.emplace_back(t[2].y, t[1].y);
}
else{
static double d[3];
int u = 0;
for (int j = 0; j < 3; j ++ ){
auto o = get_line_intersection(t[j], t[(j + 1) % 3] - t[j], {a, -INF}, {0,INF * 2});
if (cmp(o.x, INF))d[u++] = o.y;
}
if (u){
sort(d, d + u);
q.emplace_back(d[0], d[u - 1]);
}
}
}
if(!q.size())return 0;
for (int i = 0; i < q.size(); i ++ ){
if (q[i].x > q[i].y)swap(q[i].x, q[i].y);
}
sort(q.begin(),q.end());
double res = 0, st = q[0].x, ed = q[0].y;
for (int i = 1; i < q.size(); i ++ ){
if (q[i].x <= ed) ed = max(ed, q[i].y);
else{
res += ed - st;
st = q[i].x, ed = q[i].y;
}
}
res += ed - st;
return res;
}
double range_area(double a, double b,vector<vector<Point>> &tr){
return (line_area(a, 1, tr) + line_area(b, 0, tr)) * (b - a) / 2;
}
double triangle_union_area(vector<vector<Point>> &tr){
vector<double> xs;
for(int i=0;i<tr.size();i++){
for(int j=0;j<3;j++){
xs.push_back(tr[i][j].x);
}
sort(tr[i].begin(),tr[i].begin()+3);
}
for(int i=0;i<tr.size();i++){
for(int j=i+1;j<tr.size();j++){
for(int x=0;x<3;x++){
for(int y=0;y<3;y++){
auto o=get_line_intersection(tr[i][x],tr[i][(x+1)%3]-tr[i][x],tr[j][y],tr[j][(y+1)%3]-tr[j][y]);
if(cmp(o.x,INF))xs.push_back(o.x);
}
}
}
}
sort(xs.begin(),xs.end());
double res=0;
for(int i=0;i+1<xs.size();i++){
if(cmp(xs[i],xs[i+1]))res+=range_area(xs[i],xs[i+1],tr);
}
return res;
}
struct Circle{
Point p;
double r;
};
pair<Point,Point> get_line(Point a,Point b){
//找到与向量ab垂直的向量的起点与终点
return {(a+b)/2,rotate(b-a,PI/2)};
}
Circle get_circle(Point a,Point b,Point c){
//根据不共线的三点求出圆的圆心坐标和半径
auto u=get_line(a,b),v=get_line(a,c);//求出垂直于ab和ac的向量
auto get_line_intersection=[&](Point p,Vector v,Point q,Vector w){
Vector u=p-q;
double t=cross(w,u)/cross(v,w);
return p+v*t;
};
auto p=get_line_intersection(u.first,u.second,v.first,v.second);//两向量交点即为所求圆的圆心
return {p,(p-a).length()};
}
//----------------求最小圆覆盖,期望O(n)---------------------
Circle smallest_circle(vector<Point> &po){
random_shuffle(po.begin(),po.end());
Circle c={po[0],0};//最开始时直接将第一个点看作圆心，半径为0
for(int i=1;i<po.size();i++){
if(cmp(c.r,(c.p-po[i]).length())<0){
//如果第i个点不在前i-1个点的最小圆覆盖内，则进入下一重循环
c={po[i],0};//将第i个点当作圆心，半径为0，重新确定圆
for(int j=0;j<i;j++){
if(cmp(c.r,(c.p-po[j]).length())<0){
c={(po[i]+po[j])/2,(po[i]-po[j]).length()/2};
for(int k=0;k<j;k++){
if(cmp(c.r,(c.p-po[k]).length())<0){//确定了最小圆覆盖上的三个点，直接求出圆
c=get_circle(po[i],po[j],po[k]);
}
}
}
}
}
}
return c;
}
// //------------求圆与三角形的面积并---------------------
// Point q[N],r;
// double R;
// Vector norm(Vector a){//法向量
// return a/get_length(a);
// }
// double get_circle_line_intersection(Point a, Point b, Point& pa, Point& pb)
// {
// auto e = get_line_intersection(a, b - a,r, rotate(b - a, PI / 2));
// auto mind = get_dist(r, e);
// if (!on_segment(e, a, b)) mind = min(get_dist(r, a), get_dist(r, b));
// if (cmp(R, mind) <= 0) return mind;
// auto len = sqrt(R * R - get_dist(r, e) * get_dist(r, e));
// pa = e + norm(a - b) * len;
// pb = e + norm(b - a) * len;
// return mind;
// }
// double get_sector(Point a, Point b)
// {
// auto angle = acos((dot(a, b)) / get_length(a) / get_length(b));
// if (sign(cross(a, b)) < 0) angle = -angle;
// return R * R * angle / 2;
// }
// double get_circle_triangle_area(Point a, Point b)
// {
// auto da = get_dist(r, a), db = get_dist(r, b);
// if (cmp(R, da) >= 0 && cmp(R, db) >= 0) return cross(a, b) / 2;
// if (!sign(cross(a, b))) return 0;
// Point pa, pb;
// auto mind = get_circle_line_intersection(a, b, pa, pb);
// if (cmp(R, mind) <= 0) return get_sector(a, b);
// if (cmp(R, da) >= 0) return cross(a, pb) / 2 + get_sector(pb, b);
// if (cmp(R, db) >= 0) return get_sector(a, pa) + cross(pa, b) / 2;
// return get_sector(a, pa) + cross(pa, pb) / 2 + get_sector(pb, b);
// }
// double work(int n)
// {
// double res = 0;
// for (int i = 0; i < n; i ++ )
// res += get_circle_triangle_area(q[i], q[(i + 1) % n]);
// return fabs(res);
// }
}基础知识
1. 前置知识点
(1) pi = acos(-1);
(2) 余弦定理 c^2 = a^2 + b^2 - 2abcos(t)
2. 浮点数的比较
const double eps = 1e-8;
int sign(double x) // 符号函数
{
if (fabs(x) < eps) return 0;
if (x < 0) return -1;
return 1;
}
int cmp(double x, double y) // 比较函数
{
if (fabs(x - y) < eps) return 0;
if (x < y) return -1;
return 1;
}
3. 向量
3.1 向量的加减法和数乘运算
3.2 内积（点积） A·B = |A||B|cos(C)
(1) 几何意义：向量A在向量B上的投影与B的长度的乘积。
(2) 代码实现
double dot(Point a, Point b)
{
return a.x * b.x + a.y * b.y;
}
3.3 外积（叉积） AxB = |A||B|sin(C)
(1) 几何意义：向量A与B张成的平行四边形的有向面积。B在A的逆时针方向为正。
(2) 代码实现
double cross(Point a, Point b)
{
return a.x * b.y - b.x * a.y;
}
3.4 常用函数
3.4.1 取模
double get_length(Point a)
{
return sqrt(dot(a, a));
}
3.4.2 计算向量夹角
double get_angle(Point a, Point b)
{
return acos(dot(a, b) / get_length(a) / get_length(b));
}
3.4.3 计算两个向量构成的平行四边形有向面积
double area(Point a, Point b, Point c)
{
return cross(b - a, c - a);
}
3.4.5 向量A顺时针旋转C的角度：
Point rotate(Point a, double angle)
{
return Point(a.x * cos(angle) + a.y * sin(angle), -a.x * sin(angle) + a.y * cos(angle));
}
4. 点与线
4.1 直线定理
(1) 一般式 ax + by + c = 0
(2) 点向式 p0 + vt
(3) 斜截式 y = kx + b
4.2 常用操作
(1) 判断点在直线上 A x B = 0
(2) 两直线相交
// cross(v, w) == 0则两直线平行或者重合
Point get_line_intersection(Point p, Vector v, Point q, vector w)
{
vector u = p - q;
double t = cross(w, u) / cross(v, w);
return p + v * t;
}
(3) 点到直线的距离
double distance_to_line(Point p, Point a, Point b)
{
vector v1 = b - a, v2 = p - a;
return fabs(cross(v1, v2) / get_length(v1));
}
(4) 点到线段的距离
double distance_to_segment(Point p, Point a, Point b)
{
if (a == b) return get_length(p - a);
Vector v1 = b - a, v2 = p - a, v3 = p - b;
if (sign(dot(v1, v2)) < 0) return get_length(v2);
if (sign(dot(v1, v3)) > 0) return get_length(v3);
return distance_to_line(p, a, b);
}
(5) 点在直线上的投影
double get_line_projection(Point p, Point a, Point b)
{
Vector v = b - a;
return a + v * (dot(v, p - a) / dot(v, v));
}
(6) 点是否在线段上
bool on_segment(Point p, Point a, Point b)
{
return sign(cross(p - a, p - b)) == 0 && sign(dot(p - a, p - b)) <= 0;
}
(7) 判断两线段是否相交
bool segment_intersection(Point a1, Point a2, Point b1, Point b2)
{
double c1 = cross(a2 - a1, b1 - a1), c2 = cross(a2 - a1, b2 - a1);
double c3 = cross(b2 - b1, a2 - b1), c4 = cross(b2 - b1, a1 - b1);
return sign(c1) * sign(c2) <= 0 && sign(c3) * sign(c4) <= 0;
}
5. 多边形
5.1 三角形
5.1.1 面积
(1) 叉积
(2) 海伦公式
p = (a + b + c) / 2;
S = sqrt(p(p - a) * (p - b) * (p - c));
5.1.2 三角形四心
(1) 外心，外接圆圆心
三边中垂线交点。到三角形三个顶点的距离相等
(2) 内心，内切圆圆心
角平分线交点，到三边距离相等
(3) 垂心
三条垂线交点
(4) 重心
三条中线交点（到三角形三顶点距离的平方和最小的点，三角形内到三边距离之积最大的点）
5.2 普通多边形
通常按逆时针存储所有点
5.2.1 定义
(1) 多边形
由在同一平面且不再同一直线上的多条线段首尾顺次连接且不相交所组成的图形叫多边形
(2) 简单多边形
简单多边形是除相邻边外其它边不相交的多边形
(3) 凸多边形
过多边形的任意一边做一条直线，如果其他各个顶点都在这条直线的同侧，则把这个多边形叫做凸多边形
任意凸多边形外角和均为360°
任意凸多边形内角和为(n?2)180°
5.2.2 常用函数
(1) 求多边形面积（不一定是凸多边形）
我们可以从第一个顶点除法把凸多边形分成n ? 2个三角形，然后把面积加起来。
double polygon_area(Point p[], int n)
{
double s = 0;
for (int i = 1; i + 1 < n; i ++ )
s += cross(p[i] - p[0], p[i + 1] - p[i]);
return s / 2;
}
(2) 判断点是否在多边形内（不一定是凸多边形）
a. 射线法，从该点任意做一条和所有边都不平行的射线。交点个数为偶数，则在多边形外，为奇数，则在多边形内。
b. 转角法
(3) 判断点是否在凸多边形内
只需判断点是否在所有边的左边（逆时针存储多边形）。
5.3 皮克定理
皮克定理是指一个计算点阵中顶点在格点上的多边形面积公式该公式可以表示为:
S = a + b/2 - 1
其中a表示多边形内部的点数，b表示多边形边界上的点数，S表示多边形的面积。
6. 圆
(1) 圆与直线交点
(2) 两圆交点
(3) 点到圆的切线
(4) 两圆公切线
(5) 两圆相交面积凸包
Andrew算法
- x，y从小到大排序
- 正向做一遍，反向做一遍，用栈维护
- 如果栈顶向量与当前点围成的面积为负，则弹出
- 凸包的点为逆时针
cpp
int stk[N],top;
bool used[N];
vector<Point> andrew(vector<Point> &q){
sort(q.begin(),q.end());
int n=q.size();
top=0;
for(int i=0;i<n;i++)used[i]=false;
for(int i=0;i<n;i++){
while(top>=2&&sign(area(q[stk[top-1]],q[stk[top]],q[i]))<=0){
if(sign(area(q[stk[top-1]],q[stk[top]],q[i])<0))used[stk[top--]]=false;
else top--;
}
stk[++top]=i;
used[i]=true;
}
used[0]=false;
for(int i=n-1;i>=0;i--){
if(used[i])continue;
while(top>=2&&sign(area(q[stk[top-1]],q[stk[top]],q[i]))<=0)top--;
stk[++top]=i;
}
vector<Point> hull;
for(int i=1;i<top;i++)hull.push_back(q[stk[i]]);
return hull;
}极角排序
atan2法(double时使用)
象限顺序3 4 1 2
cpp
double theta(auto p){
return atan2(p.y,p.x);
}
void psort(vector<Point> &p,Point c=Point(0,0)){
sort(p.begin(),p.end(),[&](auto a,auto b){
return cmp(theta(a-c),theta(b-c))<0;
});
}叉乘法1(整数时使用)
象限顺序1 2 3 4
cpp
int quadrand(Point p){
return (p.y<0)<<1|(p.x<0)^(p.y<0);
}
void psort(vector<Point> &p,Point c=Point(0,0)){
sort(p.begin(),p.end(),[&](auto a,auto b){
return quadrand(a-c)<quadrand(b-c)||quadrand(a-c)==quadrand(b-c)&&cross(a-c,b-c)>0;
});
}叉乘法2
象限顺序2 3 4 1,如果两个点共线，那么较近的点更小
cpp
sort(tmp.begin(),tmp.end(),[&](auto x,auto y){
int c=cross(x-p[t],y-p[t]);
if(c==0)return dot(x-p[t],x-p[t])<dot(y-p[t],y-p[t]);
if(((x-p[t]).x<0)!=((y-p[t]).x<0))return (x-p[t]).x<(y-p[t]).x;
return c>0;
});平面最近点对
cpp
Point tmp[N];
void Merge(int l,int r,vector<Point> &a){
int mid=l+r>>1;
int i=l,j=mid+1;
int k=l;
while(i<=mid&&j<=r){
if(a[i].y<a[j].y)tmp[k++]=a[i++];
else tmp[k++]=a[j++];
}
while(i<=mid)tmp[k++]=a[i++];
while(j<=r)tmp[k++]=a[j++];
for(int i=l;i<=r;i++)a[i]=tmp[i];
}
double NearestPointDis(int l,int r,vector<Point> &a){
if(l>=r)return LLONG_MAX;
if(l+1==r){
if(a[l].y>a[r].y)swap(a[l],a[r]);
return (a[l]-a[r]).dis();
}
int mid=l+r>>1;
int t=a[mid].x;
double d=min(NearestPointDis(l,mid,a),NearestPointDis(mid+1,r,a));
Merge(l,r,a);
vector<int> pos;
for(int i=l;i<=r;i++){
if(1ll*(a[i].x-t)*(a[i].x-t)<d)pos.push_back(i);
}
for(int i=0;i<pos.size();i++){
for(int j=i+1;j<pos.size()&&1ll*(a[pos[i]].y-a[pos[j]].y)*(a[pos[i]].y-a[pos[j]].y)<d;j++){
d=min(d,(a[pos[i]]-a[pos[j]]).dis());
}
}
return d;
}平面最近最远点对随机法
cpp
pair<double,double> calc(int arg,vector<Point> &p){
int n=p.size();
for(int i=0;i<n;i++){
double x=p[i].x,y=p[i].y;
p[i].x=x*cos(arg)-y*sin(arg);
p[i].y=x*sin(arg)+y*cos(arg);
}
sort(p.begin(),p.end());
static const int S=10,T=10;
double mn=LLONG_MAX,mx=-LLONG_MAX;
for(int i=0;i<n;i++){
for(int j=i+1;j<n&&j<=i+S;j++)mn=min(mn,(p[i]-p[j]).dis());
for(int j=n-1;j>i&&j>=n-T-1;j--)mx=max(mx,(p[i]-p[j]).dis());
}
return {mn,mx};
}
pair<double,double> get(int T,vector<Point> &p){
double mn=LLONG_MAX,mx=-LLONG_MAX;
for(int i=0;i<T;i++){
auto [a,b]=calc(i*PI/T,p);
mn=min(mn,a),mx=max(mx,b);
}
return {mn,mx};
}旋转卡壳求凸包直径
cpp
int rotating_calipers(vector<Point> &hull){//旋转卡壳求凸包直径平方
if(is_line(hull))return (hull[0]-hull.back()).dis();
int res=0;
for(int i=0,j=2;i<hull.size();i++){
auto d=hull[i],e=hull[(i+1)%hull.size()];
while(area(d,e,hull[j])<area(d,e,hull[(j+1)%hull.size()]))j=(j+1)%hull.size();
res=max({res,(d-hull[j]).dis(0),(e-hull[j]).dis(0)});
}
return res;
}旋转卡壳求最小矩形覆盖
cpp
double get_project(Point a, Point b, Point c){
return dot(b - a, c - a) / (b - a).length();
}
vector<Point> get_min_rectangle(vector<Point> hull,double &min_area){
// 旋转卡壳求最小矩形覆盖
// andrew判断需要改成<=
// ans[4]四个顶点逆序
if(is_line(hull)){
min_area=0;
return {hull[0],hull[0],hull.back(),hull.back()};
}
vector<Point> ans(4);
int n=hull.size();
hull.push_back(hull[0]);
for(int i=0,a=2,b=1,c=2;i<n;i++){
auto d=hull[i],e=hull[i+1];
while(cmp(area(d,e,hull[a]),area(d,e,hull[a+1]))<0)a=(a+1)%n;
while(cmp(get_project(d,e,hull[b]),get_project(d,e,hull[b+1]))<0)b=(b+1)%n;
if(!i)c=a;
while(cmp(get_project(d,e,hull[c]),get_project(d,e,hull[c+1]))>0)c=(c+1)%n;
auto x=hull[a],y=hull[b],z=hull[c];
auto h=area(d,e,x)/(e-d).length();
auto w=dot(y-z,e-d)/(e-d).length();
if(h*w<min_area){
min_area=h*w;
ans[0] = d + (e - d).norm() * get_project(d, e, y);
ans[3] = d + (e - d).norm() * get_project(d, e, z);
auto u = (rotate(e - d, -PI / 2)).norm();
ans[1] = ans[0] + u * h;
ans[2] = ans[3] + u * h;
}
}
for(int i=0;i<4;i++){
if(!sign(ans[i].x))ans[i].x=0;
if(!sign(ans[i].y))ans[i].y=0;
}
return ans;
}k远点对
- 运行k次如下操作
- 找到凸包直径的对应两个点，记录这两点到其他点的距离，然后删去
- 小根堆维护k大值
半平面交(求凸多边形的交集)
cpp
struct Line{//存储凸包上每条线段的起点和终点
Point st,ed;
Line(Point _st=Point(),Point _ed=Point()):st(_st),ed(_ed){}
};
double line_get_angle(const Line& a){
return atan2(a.ed.y-a.st.y,a.ed.x-a.st.x);
}
bool line_cmp(const Line& a, const Line& b){
double A=line_get_angle(a),B=line_get_angle(b);
if (!cmp(A,B)) return area(a.st,a.ed,b.ed)<0;
return A<B;
}
Point get_line_intersection(Line a,Line b){
return get_line_intersection(a.st,a.ed-a.st,b.st,b.ed-b.st);
}
bool on_right(Line& a, Line& b, Line& c){//判断bc的交点是否在a的右侧
auto o = get_line_intersection(b,c);
return sign(area(a.st,a.ed,o))<=0;
//需要保留经过半平面交的所有直线的时候，将该处等号去掉
}
int dq[N];// 双端队列
double half_plane_intersection(vector<Line> &line){
sort(line.begin(),line.end(),line_cmp);//对每条直线按照极坐标进行排序
int hh=0,tt=-1;//定义双端队列对头和队尾元素
for(int i=0;i<line.size();i++)//遍历每一条直线
{
if (i&&!cmp(line_get_angle(line[i]),line_get_angle(line[i-1]))) continue;//当两直线极坐标角度相同
while(hh+1<=tt&&on_right(line[i],line[dq[tt-1]],line[dq[tt]])) tt--;
//如果line[l[tt-1]]和line[l[tt]]的交点在line[i]的右侧,删除队尾元素
while(hh+1<=tt&&on_right(line[i],line[dq[hh]],line[dq[hh+1]])) hh++;
//如果line[l[hh]]和line[l[hh+1]]的交点在line[i]的右侧,删除队头元素
dq[++tt]=i;//将该直线加入队尾
}
//判断一下第一条直线与最后一条直线
while (hh+1<=tt&&on_right(line[dq[hh]],line[dq[tt-1]],line[dq[tt]])) tt--;//队尾更新队头
while (hh+1<=tt&&on_right(line[dq[tt]],line[dq[hh]],line[dq[hh+1]])) hh++;//队头更新队尾
dq[++tt]=dq[hh];//将第一条直线再次加入双端队列末尾
//计算半平面交交点
vector<Point> ans;
for (int i=hh;i<tt;i++){ //对于半平面交上的直线两两之间求出交点
ans.push_back(get_line_intersection(line[dq[i]],line[dq[i+1]]));
}
//计算半平面交面积
return polygon_area(ans);
}最小覆盖圆
性质1：最小覆盖圆是唯一的 性质2：若P不在S的最小覆盖圆内部，则P在{P}∪S的最小覆盖圆边上
三维计算几何基础
1. 三维向量表示(x, y, z)
2. 向量加减法、数乘运算，与二维相同
3. 模长 |A| = sqrt(x * x + y * y + z * z)
4. 点积
(1) 几何意义：A·B = |A| * |B| * cos(C)
(2) 代数求解：(x1, y1, z1) · (x2, y2, z2) = (x1x2, y1y2, z1z2);
5. 叉积
(1) 几何意义：AxB = |A| * |B| * sin(C)，方向：右手定则
(2) 代数求解：AxB = (y1z2 - z1y2, z1x2 - x1z2, x1y2 - x2y1)
6. 如何求平面法向量
任取平面上两个不共线的向量A、B：AxB
7. 判断点D是否在平面里
任取平面上两个不共线的向量A、B：先求法向量C = AxB，然后求平面上任意一点到D的向量E与C的点积，判断点积是否为0。
8. 求点D到平面的距离
任取平面上两个不共线的向量A、B：先求法向量C = AxB。然后求平面上任意一点到D的向量E在C上的投影长度即可。即：E·C / |C|
9. 多面体欧拉定理
顶点数 - 棱长数 + 表面数 = 2
10. 三维凸包三维凸包
cpp
namespace Geometry_3D{
const double pi=acos(-1);
const double eps=1e-12;
double rand_eps(){
return ((double)rand()/RAND_MAX-0.5)*eps;
}
struct Point_3D{
double x,y,z;
Point_3D(double x=0,double y=0,double z=0):x(x),y(y),z(z){}
void shake(){
x+=rand_eps(),y+=rand_eps(),z+=rand_eps();
}
double len(){
return sqrt(x*x+y*y+z*z);
}
};
typedef Point_3D Vector_3D;
Vector_3D operator+(const Vector_3D &A, const Vector_3D &B) {
return Vector_3D(A.x + B.x, A.y + B.y, A.z + B.z);
}
Vector_3D operator-(const Point_3D &A, const Point_3D &B) {
return Vector_3D(A.x - B.x, A.y - B.y, A.z - B.z);
}
Vector_3D operator*(const Vector_3D &A, double p) {
return Vector_3D(A.x * p, A.y * p, A.z * p);
}
Vector_3D operator/(const Vector_3D &A, double p) {
return Vector_3D(A.x / p, A.y / p, A.z / p);
}
double Dot_3D(const Vector_3D &A, const Vector_3D &B) {
return A.x * B.x + A.y * B.y + A.z * B.z;
}
Vector_3D Cross_3D(const Vector_3D &A, const Vector_3D &B) {
return Vector_3D(A.y * B.z - A.z * B.y, A.z * B.x - A.x * B.z, A.x * B.y - A.y * B.x);
}
int cmp(double x,double y){//比较函数：精度为eps
if(fabs(x-y)<eps) return 0;//x与y相等
if(x<y) return -1;//x小于y
return 1;//x大于y
}
Point_3D q[N];
struct Plane{
int v[3];
Vector_3D Normal(){
return Cross_3D(q[v[1]] - q[v[0]], q[v[2]] - q[v[0]]);
}
double area(){
return Normal().len()/2;
}
bool above(int i){
return Dot_3D(q[i]-q[v[0]],Normal())>=0;
}
};
Plane plane[N],np[N];
int m=0;
bool g[N][N];
void get_convex_3d(int n){
plane[m++]={0,1,2};
plane[m++]={2,1,0};
for(int i=3;i<n;i++){
int cnt=0;
for(int j=0;j<m;j++){
bool t=plane[j].above(i);
if(!t)np[cnt++]=plane[j];
for(int k=0;k<3;k++){
g[plane[j].v[k]][plane[j].v[(k+1)%3]]=t;
}
}
for(int j=0;j<m;j++){
for(int k=0;k<3;k++){
int a=plane[j].v[k],b=plane[j].v[(k+1)%3];
if(g[a][b]&&!g[b][a]){
np[cnt++]={a,b,i};
}
}
}
m=cnt;
for(int j=0;j<m;j++)plane[j]=np[j];
}
}
}
using namespace Geometry_3D;
void solve(){
int n;
cin>>n;
for(int i=0;i<n;i++){
cin>>q[i].x>>q[i].y>>q[i].z;
q[i].shake();
}
get_convex_3d(n);
double res=0;
for(int i=0;i<m;i++){
res+=plane[i].area();
}
cout<<setprecision(3)<<fixed<<res<<"\n";
}三角形面积并
cpp
Point get_line_intersection(Point p,Vector v,Point q,Vector w){//两直线的交点
//使用前提，直线必须有交点
//cross(v, w) == 0则两直线平行或者重合
if(!sign(cross(v,w))) return {INF,INF};
Vector u=p-q;
double t=cross(w,u)/cross(v,w);
Point o=p+v*t;
auto on_segment=[&](Point p,Point a,Point b){
return sign(dot(p-a,p-b))<=0;
};
//若判断线段/射线是否相交 还需判断交点是否在线段/射线上
if(!on_segment(o,p,p+v)||!on_segment(o,q,q+w))return {INF,INF};
return o;
}
//-------------三角形面积并(n^3logn)--------------------------
double line_area(double a, int side,vector<vector<Point>> &tr){
vector<Point> q;
for (int i = 0; i < tr.size(); i ++ ){
auto t = tr[i];
if (cmp(t[0].x, a) > 0 || cmp(t[2].x, a) < 0) continue;
if (!cmp(t[0].x, a) && !cmp(t[1].x, a)){
if (side) q.emplace_back(t[0].y, t[1].y);
}
else if (!cmp(t[2].x, a) && !cmp(t[1].x, a)){
if (!side) q.emplace_back(t[2].y, t[1].y);
}
else{
static double d[3];
int u = 0;
for (int j = 0; j < 3; j ++ ){
auto o = get_line_intersection(t[j], t[(j + 1) % 3] - t[j], {a, -INF}, {0,INF * 2});
if (cmp(o.x, INF))d[u++] = o.y;
}
if (u){
sort(d, d + u);
q.emplace_back(d[0], d[u - 1]);
}
}
}
if(!q.size())return 0;
for (int i = 0; i < q.size(); i ++ ){
if (q[i].x > q[i].y)swap(q[i].x, q[i].y);
}
sort(q.begin(),q.end());
double res = 0, st = q[0].x, ed = q[0].y;
for (int i = 1; i < q.size(); i ++ ){
if (q[i].x <= ed) ed = max(ed, q[i].y);
else{
res += ed - st;
st = q[i].x, ed = q[i].y;
}
}
res += ed - st;
return res;
}
double range_area(double a, double b,vector<vector<Point>> &tr){
return (line_area(a, 1, tr) + line_area(b, 0, tr)) * (b - a) / 2;
}
double triangle_union_area(vector<vector<Point>> &tr){
vector<double> xs;
for(int i=0;i<tr.size();i++){
for(int j=0;j<3;j++){
xs.push_back(tr[i][j].x);
}
sort(tr[i].begin(),tr[i].begin()+3);
}
for(int i=0;i<tr.size();i++){
for(int j=i+1;j<tr.size();j++){
for(int x=0;x<3;x++){
for(int y=0;y<3;y++){
auto o=get_line_intersection(tr[i][x],tr[i][(x+1)%3]-tr[i][x],tr[j][y],tr[j][(y+1)%3]-tr[j][y]);
if(cmp(o.x,INF))xs.push_back(o.x);
}
}
}
}
sort(xs.begin(),xs.end());
double res=0;
for(int i=0;i+1<xs.size();i++){
if(cmp(xs[i],xs[i+1]))res+=range_area(xs[i],xs[i+1],tr);
}
return res;
}最小圆覆盖
cpp
struct Circle{
Point p;
double r;
};
pair<Point,Point> get_line(Point a,Point b){
//找到与向量ab垂直的向量的起点与终点
return {(a+b)/2,rotate(b-a,PI/2)};
}
Circle get_circle(Point a,Point b,Point c){
//根据不共线的三点求出圆的圆心坐标和半径
auto u=get_line(a,b),v=get_line(a,c);//求出垂直于ab和ac的向量
auto get_line_intersection=[&](Point p,Vector v,Point q,Vector w){
Vector u=p-q;
double t=cross(w,u)/cross(v,w);
return p+v*t;
};
auto p=get_line_intersection(u.first,u.second,v.first,v.second);//两向量交点即为所求圆的圆心
return {p,(p-a).length()};
}
//----------------求最小圆覆盖,期望O(n)---------------------
Circle smallest_circle(vector<Point> &po){
random_shuffle(po.begin(),po.end());
Circle c={po[0],0};//最开始时直接将第一个点看作圆心，半径为0
for(int i=1;i<po.size();i++){
if(cmp(c.r,(c.p-po[i]).length())<0){
//如果第i个点不在前i-1个点的最小圆覆盖内，则进入下一重循环
c={po[i],0};//将第i个点当作圆心，半径为0，重新确定圆
for(int j=0;j<i;j++){
if(cmp(c.r,(c.p-po[j]).length())<0){
c={(po[i]+po[j])/2,(po[i]-po[j]).length()/2};
for(int k=0;k<j;k++){
if(cmp(c.r,(c.p-po[k]).length())<0){//确定了最小圆覆盖上的三个点，直接求出圆
c=get_circle(po[i],po[j],po[k]);
}
}
}
}
}
}
return c;
}自适应辛普森积分
cpp
const double eps=1e-12;
double a,b,c,d,L,R;
double f(double x){
return (c*x+d)/(a*x+b);
}
double simpson(double l,double r){
double mid=(l+r)/2;
return (r-l)*(f(l)+4*f(mid)+f(r))/6;
}
double asr(double l,double r,double s){
double mid=(l+r)/2;
double left=simpson(l,mid),right=simpson(mid,r);
if(fabs(left+right-s)<eps)return left+right;
return asr(l,mid,left)+asr(mid,r,right);
}
void solve(){
cin>>a>>b>>c>>d>>L>>R;
cout<<setprecision(6)<<fixed<<asr(L,R,simpson(L,R))<<"\n";
}自适应辛普森积分求圆的面积并
cpp
const double eps=1e-12;
int n;
struct Point{
double x,y;
bool operator<(const Point &t)const{
return x!=t.x?x<t.x:y<t.y;
}
};
struct Circle{
Point p;
double R;
} c[N];
Point q[N];
int dcmp(double x,double y){
if(fabs(x-y)<eps)return 0;
if(x<y)return -1;
return 1;
}
double f(double x){
int cnt=0;
for(int i=0;i<n;i++){
double X=fabs(x-c[i].p.x),R=c[i].R;
if(dcmp(X,R)<0){
double Y=sqrt(R*R-X*X);
q[cnt++]={c[i].p.y-Y,c[i].p.y+Y};
}
}
if(!cnt)return 0;
sort(q,q+cnt);
double res=0,st=q[0].x,ed=q[0].y;
for(int i=1;i<cnt;i++){
if(q[i].x<=ed)ed=max(ed,q[i].y);
else{
res+=ed-st;
st=q[i].x,ed=q[i].y;
}
}
return res+ed-st;
}
double simpson(double l,double r,double sl,double sr){
double mid=(l+r)/2;
return (r-l)*(sl+4*f(mid)+sr)/6;
}
double asr(double l,double r,double s,double sl,double sr){
double mid=(l+r)/2;
double S=f(mid);
double left=simpson(l,mid,sl,S),right=simpson(mid,r,S,sr);
if(fabs(left+right-s)<eps)return left+right;
return asr(l,mid,left,sl,S)+asr(mid,r,right,S,sr);
}
double dis(Point a,Point b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
void init(){ //预处理，删除被包含的圆
int cnt = 0;
for (int i = 0, j; i < n; ++i) {
for (j = i+1; j < n; ++j) {
if(c[i].R-c[j].R + dis(c[i].p, c[j].p) < eps){
break;
}
}
if(j >= n) c[cnt++] = c[i];
}
n = cnt;
}
vector<Point> segs;
void merge(){
sort(segs.begin(),segs.end());
vector<Point> res;
double st=segs[0].x,ed=segs[0].y;
for(int i=1;i<segs.size();i++){
if(dcmp(segs[i].x,ed)<=0)ed=max(ed,segs[i].y);
else{
res.push_back({st,ed});
st=segs[i].x,ed=segs[i].y;
}
}
res.push_back({st,ed});
segs.swap(res);
}
void solve(){
cin>>n;
for(int i=0;i<n;i++)cin>>c[i].p.x>>c[i].p.y>>c[i].R,segs.push_back({c[i].p.x - c[i].R, c[i].p.x + c[i].R});
sort(c,c+n,[&](Circle a,Circle b){
return a.R<b.R;
});
init();
merge();
double res=0;
for(int i=0;i<segs.size();i++){
auto [l,r]=segs[i];
res+=asr(l,r,f((l+r)/2),f(l),f(r));
}
cout<<setprecision(3)<<fixed<<res<<"\n";
}计算积分(发散收敛)
试计算积分
保留至小数点后位。若积分发散，请输出 。
cpp
const double eps=1e-12;
double a;
double f(double x){
return pow(x,a/x-x);
}
double simpson(double l,double r,double sl,double sr){
double mid=(l+r)/2;
return (r-l)*(sl+4*f(mid)+sr)/6;
}
double asr(double l,double r,double s,double sl,double sr){
double mid=(l+r)/2;
double S=f(mid);
double left=simpson(l,mid,sl,S),right=simpson(mid,r,S,sr);
if(fabs(left+right-s)<eps)return left+right;
return asr(l,mid,left,sl,S)+asr(mid,r,right,S,sr);
}
void solve(){
cin>>a;
if(a<0)return cout<<"orz",void();
cout<<setprecision(5)<<fixed<<asr(eps,20,f((eps+20)/2),f(eps),f(20));
}