二维几何
点和向量化为坐标Coord进行运算,使用stl中的complex实现。
复数相乘的几何意义为长度相乘,极角相加。
用直线上的一点p和方向向量v表示一条经过p的直线,直线上的所有点q满足q=p+t*v,其中t是参数;当限制t≥0时,该参数方程表示射线;限制0≤t≤1时,该参数方程表示线段。
此外,如果已知线段端点a1和a2,可以通过Line(a1,a2-a1)来得到对应的参数形式。
Morley定理:三角形每个内角的三等分线相交成等边三角形。
欧拉定理:平面图的点数V、边数E和面数F满足V+F-E=2。
const double EPS=1e-9,PI=acos(-1);
typedef complex<double> Coord;
#define X real()
#define Y imag()
struct Line
{
Coord p,v;
Line(Coord p=Coord(),Coord v=Coord()):p(p),v(v) {}
Coord point(double t)
{
return p+v*t;
}
};
struct Circle
{
Coord c;
double r;
Circle(Coord c=Coord(),double r=0):c(c),r(r) {}
Coord point(double t)//t为参数,幅角
{
return c+polar(r,t);
}
};
/*
Coord(double x=0,double y=0);//构造函数
double real(Coord a);//a的实部(复平面的横坐标),也可写作a.real()
double imag(Coord a);//a的虚部(复平面的纵坐标),也可写作a.imag()
double abs(Coord a);//向量a的模长,或是点a到原点的距离
double norm(Coord a);//abs的平方,比abs快,但是要注意浮点数精度溢出
double arg(Coord a);//a的幅角,与atan2(a.real(),a.imag())等价
Coord polar(double r,double t);//极坐标生成方式,r为幅值,t为幅角
//运算符重载+、-、*、/(以及对应的赋值运算,但是赋值运算不能写在表达式中,详见参考地址)、<<、>>(输出括号形式的坐标)
*/
int sgn(double d)
{
return (d>EPS)-(d<-EPS);
}
bool operator!=(const Coord &A,const Coord &B)//不等运算符,涉及到浮点数比较要重写
{
return sgn(A.X-B.X)||sgn(A.Y-B.Y);
}
bool operator==(const Coord &A,const Coord &B)
{
return !(A!=B);
}
bool cmpCoord(const Coord &A,const Coord &B)//复数没有小于运算,只能这样定义一个比较函数
{
return sgn(A.X-B.X)?
A.X<B.X:
A.Y+EPS<B.Y;
}
bool cmpLine(const Line &A,const Line &B)//按极角排序,求凸包中使用
{
return arg(A.v)<arg(B.v);
}
double Dot(Coord A,Coord B)
{
return A.X*B.X+A.Y*B.Y;
}
double Cross(Coord A,Coord B)
{
return A.X*B.Y-B.X*A.Y;
}
double Angle(Coord A,Coord B)
{
return acos(Dot(A,B)/abs(A)/abs(B));
}
double Area2(Coord A,Coord B,Coord C)//三角形ABC有向面积的两倍
{
return Cross(B-A,C-A);
}
Coord Rotate(Coord A,double rad)//向量A逆时针旋转rad弧度
{
return A*polar(1.0,rad);
}
Coord Normal(Coord A)//A的法向量,把A逆时针旋转九十度并长度化为1
{
double L=abs(A);
return Coord(-A.Y/L,A.X/L);
}
bool onLeft(Coord P,Line L)//p是否在有向直线L左侧,不含线上
{
return Cross(L.v,P-L.p)>0;
}
double DistanceToLine(Coord P,Line L)//点到直线距离(有向)
{
return Cross(L.v,P-L.p)/abs(L.v);
}
double DistanceToLine(Coord P,Coord A,Coord B)
{
return DistanceToLine(P,Line(A,B-A));
}
double DistanceToSegment(Coord P,Coord A,Coord B)//点到线段的距离(无向)
{
if(A==B)return abs(P-A);
Coord v1=B-A,v2=P-A,v3=P-B;
if(sgn(Dot(v1,v2))<0)return abs(v2);
if(sgn(Dot(v1,v3))>0)return abs(v3);
return fabs(DistanceToLine(P,Line(A,B-A)));
}
Coord getLineProjection(Coord P,Line L)//点在直线上的投影
{
return L.point(Dot(L.v,P-L.p)/norm(L.v));
}
Coord getLineProjection(Coord P,Coord A,Coord B)
{
return getLineProjection(P,Line(A,B-A));
}
Coord getSymmetry(Coord P,Coord O)//P关于O的对称点
{
return O+O-P;
}
Coord getSymmetry(Coord P,Line L)//P关于L的对称点
{
return getSymmetry(P,getLineProjection(P,L));
}
Coord getLineIntersection(Line L1,Line L2)//直线交点,须确保两直线相交
{
return L1.point(Cross(L2.v,L1.p-L2.p)/Cross(L1.v,L2.v));
}
Coord getLineIntersection(Coord A1,Coord A2,Coord B1,Coord B2)
{
return getLineIntersection(Line(A1,A2-A1),Line(B1,B2-B1));
}
bool SegmentProperIntersection(Coord A1,Coord A2,Coord B1,Coord B2)//线段相交判定,交点不在一条线段的端点
{
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 sgn(C1)*sgn(C2)<0&&sgn(C3)*sgn(C4)<0;
}
bool onSegment(Coord P,Coord A1,Coord A2)//判断点是否在线段上,不包含端点
{
return sgn(Dot(A1-P,A2-P))<0&&!sgn(Cross(A1-P,A2-P));
}
double PolygonArea(const vector<Coord> &p)//计算多边形的有向面积,凸多边形即为面积
{
double area=0;
for(int i=1,n=p.size()-1; i<n; ++i)
area+=Area2(p[0],p[i],p[i+1]);
return area/2;
}
int inPolygon(Coord p,const vector<Coord> &poly)//点在多边形内的判定,转角法,正值为内部,0为外部,-1在边界上
{
int ans=0;
for(int i=0,k,d1,d2,n=poly.size(); i!=n; ++i)
{
if(onSegment(p,poly[i],poly[(i+1)%n]))return -1;//在边界上
k=sgn(Cross(poly[(i+1)%n]-poly[i],p-poly[i]));
d1=sgn(poly[i].Y-p.Y);
d2=sgn(poly[(i+1)%n].Y-p.Y);
if(k>0&&d1<=0&&d2>0)++ans;
if(k<0&&d2<=0&&d1>0)++ans;
}
return ans;
}
int ConvexHull(vector<Coord> p,vector<Coord> &sol)//获得凸包;不希望凸包的边上有输入点,把两个<=改成<
{
sort(p.begin(),p.end(),cmpCoord);//先比横坐标再比纵坐标
for(int i=0; i!=p.size(); ++i)
{
while(sol.size()>1&&
Area2(sol[sol.size()-2],sol[sol.size()-1],p[i])<=0)
sol.pop_back();
sol.push_back(p[i]);
}
for(int i=sol.size()-2,k=sol.size(); i>=0; --i)
{
while(sol.size()>k&&
Area2(sol[sol.size()-2],sol[sol.size()-1],p[i])<=0)
sol.pop_back();
sol.push_back(p[i]);
}
if(p.size()>1)sol.pop_back();
return sol.size();
}
vector<Coord> cutPolygon(const vector<Coord> &poly,Coord A,Coord B)//用有向直线A->B切割多边形poly, 返回“左侧”。 如果退化,可能会返回一个单点或者线段,复杂度O(n^2)
{
vector<Coord> newpoly;
for(int i=0,n=poly.size(); i!=n; ++i)
{
Coord C=poly[i],D=poly[(i+1)%n];
if(sgn(Cross(B-A,C-A))>=0)newpoly.push_back(C);
if(!sgn(Cross(B-A, C-D)))
{
Coord ip=getLineIntersection(Line(A,B-A),Line(C,D-C));
if(onSegment(ip,C,D))newpoly.push_back(ip);
}
}
return newpoly;
}
vector<Coord> getHalfPlaneIntersection(vector<Line> &L)//半平面交
{
sort(L.begin(),L.end(),cmpLine);//按极角排序
int first,last;//双端队列的第一个元素和最后一个元素
vector<Coord> p(L.size(),Coord()); //p[i]为q[i]和q[i+1]的交点
vector<Line> q(L.size(),Line());//双端队列
q[first=last=0]=L[0]; //队列初始化为只有一个半平面L[0]
for(int i=0,n=L.size(); i!=n; ++i)
{
while(first<last&&!onLeft(p[last-1],L[i]))
--last;
while(first<last&&!onLeft(p[first],L[i]))
++first;
q[++last]=L[i];
if(!sgn(Cross(q[last].v,q[last-1].v)))
{
--last;
if(onLeft(L[i].p,q[last]))
q[last]=L[i];
}
if(first<last)
p[last-1]=getLineIntersection(q[last-1], q[last]);
}
while(first<last&&!onLeft(p[last-1],q[first]))
--last;//删除无用平面
if(last-first<=1)return vector<Coord>();//空集
p[last]=getLineIntersection(q[last],q[first]);
return vector<Coord>(p.begin()+first,p.begin()+last+1);//从deque复制到输出中
}
int getLineCircleIntersection(Line L,Circle C,vector<Coord> &sol)
{
double a=L.v.X,
b=L.p.X-C.c.X,
c=L.v.Y,
d=L.p.Y-C.c.Y,
e=a*a+c*c,
f=2*(a*b+c*d),
g=b*b+d*d-C.r*C.r,
delta=f*f-4*e*g;
if(sgn(delta)<0)return 0;
if(!sgn(delta))
return sol.push_back(L.point(-f/(2*e))),1;
sol.push_back(L.point((-f-sqrt(delta))/(2*e)));
sol.push_back(L.point((-f+sqrt(delta))/(2*e)));
return 2;
}
int getCircleIntersection(Circle C1,Circle C2,vector<Coord> &sol)
{
double d=abs(C1.c-C2.c);
if(!sgn(d))
return sgn(C1.r-C2.r)?0:-1;//重合返回-1
if(sgn(C1.r+C2.r-d)<0||sgn(fabs(C1.r-C2.r)-d)>0)//外离或内含
return 0;
double a=arg(C2.c-C1.c),
da=acos((C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d));
Coord p1=C1.point(a-da),p2=C1.point(a+da);
sol.push_back(p1);
if(p1==p2)return 1;//相切
return sol.push_back(p2),2;
}
Line getTangent(Coord C,Coord P)//圆心C,圆上一点P处切线
{
return Line(P,Normal(C-P));
}
int getTangents(Coord p,Circle C,vector<Coord> &sol)//点到圆的切点,返回个数
{
Coord u=p-C.c;
double d=abs(u);
if(d<C.r)return 0;//点在圆内
if(!sgn(d-C.r))//点在圆上
return sol.push_back(p),1;
double base=arg(u),ang=acos(C.r/d);
sol.push_back(C.point(base+ang));
sol.push_back(C.point(base-ang));
return 2;
}
int getTangents(Circle A,Circle B,vector<Coord> &a,vector<Coord> &b)//公共切线的切点
{
int cnt=0;
if(A.r<B.r)
swap(A,B),swap(a,b);//有时需标记交换
double d=abs(A.c-B.c),
rdiff=A.r-B.r,
rsum=A.r+B.r;
if(sgn(d-rdiff)<0)return 0;//内含
double base=arg(B.c-A.c);
if(!sgn(d)&&!sgn(rdiff))return -1;//重合,无穷多条切线
if(!sgn(d-rdiff))//内切,外公切线
{
a.push_back(A.point(base));
b.push_back(B.point(base));
return 1;
}
//有外公切线的情形
double ang=acos(rdiff/d);
a.push_back(A.point(base+ang));
b.push_back(B.point(base+ang));
a.push_back(A.point(base-ang));
b.push_back(B.point(base-ang));
cnt+=2;
if(!sgn(d-rsum))
{
a.push_back(A.point(base));
b.push_back(B.point(base+PI));
++cnt;
}
else if(sgn(d-rsum)>0)
{
double ang_in=acos(rsum/d);
a.push_back(A.point(base+ang_in));
b.push_back(B.point(base+ang_in+PI));
a.push_back(A.point(base-ang_in));
b.push_back(B.point(base-ang_in+PI));
cnt+=2;
}
return cnt;
}
double AreaCircleWithTriangle(Circle C,Coord A,Coord B)//C和三角形OAB的相交面积,如果三角形顶点不在O上则把圆和三角形同时平移,直到有一个顶点在O上
{
int sg=sgn(Cross(A,B));
if(!sg||A==C.c||B==C.c)return 0;
double OA=abs(A-C.c),OB=abs(B-C.c),angle=Angle(A,B),
d=DistanceToLine(Coord(),A,B);
if(sgn(OA-C.r)<=0&&sgn(OB-C.r)<=0)
return Cross(A,B)/2;
if(sgn(OA-C.r)>=0&&sgn(OB-C.r)>=0&&sgn(d-C.r)>=0)
return sg*C.r*C.r*angle/2;
if(sgn(OA-C.r)>=0&&sgn(OB-C.r)>=0&&sgn(d-C.r)<0)
{
Coord prj=getLineProjection(Coord(),A,B);
if(!onSegment(prj,A,B))return sg*C.r*C.r*angle/2;
vector<Coord> p;
Line L=Line(A,B-A);
getLineCircleIntersection(L,C,p);
double s1=C.r*C.r*angle/2,
s2=C.r*C.r*Angle(p[0],p[1])/2;
s2-=fabs(Cross(p[0],p[1])/2);
s1=s1-s2;
return sg*s1;
}
if(sgn(OB-C.r)<0)swap(A,B);
Line L=Line(A,B-A);
vector<Coord> inter;
getLineCircleIntersection(L,C,inter);
Coord inter_point=inter[!onSegment(inter[0],A,B)];
double s=fabs(Cross(inter_point,A)/2);
s+=C.r*C.r*Angle(inter_point,B)/2;
return s*sg;
}
double AreaCircleWithPolygon(Circle C,const vector<Coord> &p)
{
double ans=0;
for(int i=0; i<p.size(); ++i)
ans+=AreaCircleWithTriangle(C,p[i],p[(i+1)%p.size()]);
return fabs(ans);
}
Coord getGravityCenter(const vector<Coord> &p)//多边形重心
{
Coord a(0,0);
double am=0,mj;
for(int i=0; i<p.size(); ++i)
{
mj=Cross(p[i],p[(i+1)%p.size()]);
a+=mj*(p[i]+p[(i+1)%p.size()]);
am+=mj;
}
return a/am/3.0;
}