这个模板的缺点:有很多模板都是非常C语言化的,虽然可读性比较差但是应该性能非常感人。鉴于ACM可以开O2所以方向用STL就好了。但是有的地方不好,比如半平面交的排序,对同一个向量多次判断倾斜角,其实预处理的时候要是需要,就把倾斜角初始化就好了。
基础
//不要输出-0.0之类的数
const double eps=1e-8;
const double pi=acos(-1.0);
//判断浮点数的符号
inline int cmp(double x){
return (fabs(x)<eps)?0:((x>0.0)?1:-1);
}
inline double sqr(double x){
return x*x;
}
多边形
点
struct Point {
double x,y;
Point() {};
Point(const double x,const double y):x(x),y(y) {};
friend Point operator+(const Point &a,const Point &b) {
return Point(a.x+b.x,a.y+b.y);
}
friend Point operator-(const Point &a,const Point &b) {
return Point(a.x-b.x,a.y-b.y);
}
friend Point operator*(const Point &p,const double k) {
return Point(p.x*k,p.y*k);
}
friend Point operator*(const double k,const Point &p) {
return Point(p.x*k,p.y*k);
}
friend Point operator/(const Point &p,const double k) {
return Point(p.x/k,p.y/k);
}
friend bool operator==(const Point &a,const Point &b) {
return cmp(a.x-b.x)==0&&cmp(a.y-b.y)==0;
}
Point rotate(double A) {
//向量绕原点旋转A弧度
return Point(x*cos(A)-y*sin(A),x*sin(A)+y*cos(A));
}
double norm() {
return sqrt(sqr(x)+sqr(y));
}
};
double det(const Point &a,const Point &b) {
return a.x*b.y-a.y*b.x;
}
double dot(const Point &a,const Point &b) {
return a.x*b.x+a.y*b.y;
}
double dist(const Point &a,const Point &b) {
return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
}
线段
struct Line {
Point a,b;
Line() {};
Line(const Point &a,const Point &b):a(a),b(b) {};
Line move_dist(const double &d) {
//向法向平移d单位长度
//单位法向量n,从a指向b
Point n=b-a;
n=n/n.norm();
//左旋90度
n=n.rotate(pi/2.0);
return Line(a+n*d,b+n*d);
}
};
double dist_point_to_line(const Point &p,const Line &l) {
Point a=l.a,b=l.b;
//当a与b可以重合时,这里要加上下面的语句
/*if(a==b)
return a.dist(p);*/
if(cmp(dot(p-a,b-a))<0)
return dist(p,a);
if(cmp(dot(p-b,a-b))<0)
return dist(p,b);
return fabs(det(a-p,b-p)/dist(a,b));
}
Point point_project_on_line(const Point &p,const Line &l) {
Point a=l.a,b=l.b;
double r=dot(b-a,p-a)/dot(b-a,b-a);
return a+(b-a)*r;
}
bool point_on_line(const Point &p,const Line &l) {
Point a=l.a,b=l.b;
//这里的line是线段
//第一个cmp意思是叉积等于0,意味着直线穿过该点
//第二个cmp的<=意思是点在线段内(含端点),当改为<为点在线段内(不含端点)
return cmp(det(p-a,b-a))==0&&cmp(dot(p-a,p-b))<=0;
}
bool parallel(const Line &tl,const Line &l) {
Point a=tl.a,b=tl.b;
//叉积等于0,意味着向量平行
return !cmp(det(a-b,l.a-l.b));
}
bool intersect(const Line &tl,const Line &l,Point &p) {
Point a=tl.a,b=tl.b;
//判断直线是否相交,相交则求出交点(不需要交点可以直接return)
if(parallel(tl,l))
return false;
double s1=det(a-l.a,l.b-l.a);
double s2=det(b-l.a,l.b-l.a);
p=(b*s1-a*s2)/(s1-s2);
return true;
}
多边形
const int MAXN=105;
struct Polygon {
int n;
Point a[MAXN];
Polygon() {};
double perimeter() {
double sum=0.0;
a[n]=a[0];
for(int i=0; i<n; i++)
sum+=(a[i+1]-a[i]).norm();
return sum;
}
double area() {
double sum=0.0;
a[n]=a[0];
for(int i=0; i<n; i++)
sum+=det(a[i+1],a[i]);
return sum/2.0;
}
Point masscenter(){
Point ans(0.0,0.0);
//在这里,当多边形面积为0,返回的是原点
if(cmp(area())==0)
return ans;
a[n]=a[0];
for(int i=0;i<n;i++)
ans=ans+(a[i]+a[i+1])*det(a[i+1],a[i]);
return ans/area()/6.0;
}
//下面两个只有格点多边形能用
int border_point_num(){
int num=0;
a[n]=a[0];
for(int i=0;i<n;i++)
num+=__gcd(abs(int(a[i+1].x-a[i].x)),abs(int(a[i+1].y-a[i].y)));
return num;
}
int inside_point_num(){
return (int)area()+1-border_point_num()/2;
}
};
int point_in_polygon(Point &p,Polygon &po) {
Point *a=po.a;
int n=po.n;
int num=0,d1,d2,k;
a[n]=a[0];
for(int i=0; i<n; i++) {
if(point_on_line(p,Line(a[i],a[i+1])))
return 2;
k=cmp(det(a[i+1]-a[i],p-a[i]));
d1=cmp(a[i].y-p.y);
d2=cmp(a[i+1].y-p.y);
if(k>0&&d1<=0&&d2>0)
num++;
if(k<0&&d2<=0&&d1>0)
num--;
}
return num!=0;
}
凸包
struct Polygon_Convex {
vector<Point> P;
Polygon_Convex(int Size=0) {
P.resize(Size);
}
Polygon to_polygon() {
//注意多边形的最大点数要够
Polygon p;
p.n=P.size();
for(int i=0; i<p.n; i++) {
p.a[i]=P[i];
}
return p;
}
double diameter(int &First,int &Second){
//旋转卡壳求直径,O(n)
vector<Point> &p=P;
int n=P.size();
double maxd=0.0;
if(n==1){
First=Second=0;
return maxd;
}
#define next(i) ((i+1)%n)
for(int i=0,j=1;i<n;++i){
while(cmp(det(p[next(i)]-p[i],p[j]-p[i])-det(p[next(i)]-p[i],p[next(j)]-p[i]))<0)
j=next(j);
double d=dist(p[i],p[j]);
if(d>maxd){
maxd=d;
First=i;
Second=j;
}
d=dist(p[next(i)],p[next(j)]);
if(d>maxd){
maxd=d;
First=i;
Second=j;
}
}
#undef next(i)
return maxd;
}
};
bool comp_less(const Point&a,const Point &b) {
//水平序
return cmp(a.x-b.x)<0||cmp(a.x-b.x)==0&&cmp(a.y-b.y)<0;
}
Polygon_Convex convex_hull(vector<Point> a) {
Polygon_Convex res(2*a.size()+5);
sort(a.begin(),a.end(),comp_less);
a.erase(unique(a.begin(),a.end()),a.end());
int m=0;
for(int i=0; i<a.size(); ++i) {
while(m>1&&cmp(det(res.P[m-1]-res.P[m-2],a[i]-res.P[m-2]))<=0)
--m;
res.P[m++]=a[i];
}
int k=m;
for(int i=int(a.size())-2; i>=0; --i) {
while(m>k&&cmp(det(res.P[m-1]-res.P[m-2],a[i]-res.P[m-2])<=0))
--m;
res.P[m++]=a[i];
}
//当只有一个点时,凸包保留一个点,否则结尾和开头重复了
res.P.resize(m-(a.size()>1));
return res;
}
int point_in_polygon_convex(const Point &p,const Polygon_Convex &pc) {
//0在外部,1在内部,2在边界上
//包括边界
int n=pc.P.size();
const vector<Point> &P=pc.P;
//找一个内部点
Point g=(P[0]+P[n/3]+P[2*n/3])/3.0;
int l=0,r=n;
while(l+1<r){
int mid=(l+r)>>1;
if(cmp(det(P[l]-g,P[mid]-g))>0){
if(cmp(det(P[l]-g,p-g))>=0&&cmp(det(P[mid]-g,p-g))<0)
r=mid;
else
l=mid;
}
else{
if(cmp(det(P[l]-g,p-g))<0&&cmp(det(P[mid]-g,p-g))>=0)
l=mid;
else
r=mid;
}
}
r%=n;
int z=cmp(det(P[r]-p,P[l]-p));
//z==0在边界上,三点共线
//z==1在凸包外
//z==-1在凸包内
return (z+2)%3;
}
半平面
强烈建议加上四周围的inf边,这样传入的v既不会是空的,也不会有相反向量叉积做除法(转180度之前会遇到inf边)。
小心三线共点返回的半平面,有点东西。(好像可以被反向平行检测掉)
vector<Halfplane> Hp;
Hp.push_back(Halfplane(Point(-inf,-inf),Point(inf,-inf)));
Hp.push_back(Halfplane(Point(inf,-inf),Point(inf,inf)));
Hp.push_back(Halfplane(Point(inf,inf),Point(-inf,inf)));
Hp.push_back(Halfplane(Point(-inf,inf),Point(-inf,-inf)));
struct Halfplane {
//向量first->second的左侧
Point first,second;
Halfplane() {};
Halfplane(Point p1,Point p2):first(p1),second(p2) {};
};
inline int satisfy(Point a,Halfplane p) {
return cmp(det(a-p.first,p.second-p.first))<=0;
}
Point intersect_point(const Halfplane &a,const Halfplane &b) {
double k=det(b.first-b.second,a.first-b.second);
double t=det(b.first-b.second,a.second-b.second);
//反向向量已经被半平面交制裁了,这个函数调用之前就要先保证不平行,否则后果自负
//把边界也放进来了,所以反向向量至少转180度会遇到边界
k=k/(k-t);
return a.first+(a.second-a.first)*k;
}
inline bool compare(const Halfplane &a,const Halfplane &b) {
int res=cmp((a.second-a.first).angle()-(b.second-b.first).angle());
return res==0?satisfy(a.first,b):res<0;
}
inline bool parallel(const Halfplane &a,const Halfplane &b){
Point pa=a.second-a.first;
Point pb=b.second-b.first;
return !cmp(det(pa,pb));
}
//半平面交,O(nlogn)
vector<Point> halfplane_intersection(vector<Halfplane> v) {
//半平面把边界放进来,不可能是空的
if(v.empty())
return vector<Point>();
sort(v.begin(),v.end(),compare);
deque<Halfplane> q;
deque<Point> ans;
q.push_back(v[0]);
int vs=v.size();
for(int i=1; i<vs; ++i) {
if(cmp((v[i].second-v[i].first).angle()-(v[i-1].second-v[i-1].first).angle())==0)
continue;
while(!ans.empty()&&!satisfy(ans.back(),v[i])) {
ans.pop_back();
q.pop_back();
}
while(!ans.empty()&&!satisfy(ans.front(),v[i])) {
ans.pop_front();
q.pop_front();
}
//前面已经去掉平行了,再平行就是反平行
if(parallel(q.back(),v[i]))
return vector<Point>();
ans.push_back(intersect_point(q.back(),v[i]));
q.push_back(v[i]);
}
while(!ans.empty()&&!satisfy(ans.back(),q.front())) {
ans.pop_back();
q.pop_back();
}
while(!ans.empty()&&!satisfy(ans.front(),q.back())) {
ans.pop_front();
q.pop_front();
}
ans.push_back(intersect_point(q.back(),q.front()));
return vector<Point>(ans.begin(),ans.end());
}