###點的定義app
//考慮偏差的加法
double add(double a, double b)
{
if (abs(a + b) < EPS*(abs(a) + abs(b)))return 0;
return a + b;
}
struct P
{
double x, y;
P() {}
P(double x,double y):x(x),y(y){}
P operator+(P p) {
return P(add(x,p.x), add(y , p.y));
}
P operator-(P p) {
return P(add(x, -p.x), add(y, -p.y));
}
P operator*(double d) {
return P(x*d, y*d);
}
double dot(P p) {//內積，點積
return add(x*p.x, y*p.y);
}
double det(P p) {//外積，叉乘
return add(x*p.y, -y * p.x);
}
};
//判斷q點是否在線段p1-p2上
bool on_seg(P p1, P p2, P q) {
return (p1 - q).det(p2 - q) == 0 && (p1 - q).dot(p2 - q) <= 0;
}
//計算直線p1-p2與直線q1-q2的交點
P intersection(P p1, P p2, P q1, P q2) {
return p1 + (p2 - p1)*((q2 - q1).det(q1 - p1) / (q2 - q1).det(p2 - p1));
}
//判斷線段相交的代碼，此處是判斷n條線段相交的程序
for (int i = 0; i<n; i++) {
g[i][i] = 1;
for (int j = 0; j<i; j++) {
//判斷木棍i和木棍j是否有公共點
if ((p[i] - q[i]).det(p[j] - q[j]) == 0) {
//平行時
g[i][j] = g[j][i] = on_seg(p[i], q[i], p[j])
|| on_seg(p[i], q[i], q[j])
|| on_seg(p[j], q[j], p[i])
|| on_seg(p[j], q[j], q[i]);
}
else {
//不平行時
P r = intersection(p[i], q[i], p[j], q[j]);
g[i][j] = g[j][i] = on_seg(p[i], q[i], r) && on_seg(p[j], q[j], r);
}
}
}
###隊友提供spa
#include<stdio.h>
#include<math.h>
#include<string.h>
#define eps 1e-8
#define pi 3.141592653589793
using namespace std;
//二維點類
struct Point
{
double x,y;
Point(double a=0,double b=0){x=a;y=b;}
};
typedef Point Vector;
//二維直線類，通常方程ax+by+c=0
struct Line
{
double a,b,c,angle;
Point p1,p2;
Line(Point s,Point e)
{
a=s.y-e.y;
b=e.x-s.x;
c=s.x*e.y-e.x*s.y;
angle=atan2(e.y-s.y,e.x-s.x);
p1=s;p2=e;
}
Line(){}
};
//二維線段類
struct Segment
{
Point s,e;
Segment(Point a,Point b){s=a;e=b;}
Segment(double x1,double y1,double x2,double y2)
{
s=Point(x1,y1);
e=Point(x2,y2);
}
Segment(){}
};
//向量的加減及數乘
Vector operator + (Point a,Point b)
{
return Vector(a.x+b.x,a.y+b.y);
}
Vector operator - (Point a,Point b)
{
return Vector(a.x-b.x,a.y-b.y);
}
Vector operator * (Point a,double k)
{
return Vector(a.x*k,a.y*k);
}
Vector operator / (Point a,double k)
{
return Vector(a.x/k,a.y/k);
}
//求向量的模（長度）
double len(Vector a)
{
return sqrt(a.x*a.x+a.y*a.y);
}
//獲得sp-op和ep-op的叉積
//>0時ep在opsp的逆時針方向，<0時順時針，=0時共線
double Cross(Point &sp, Point &ep, Point &op)
{
return (sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y);
}
//兩向量求叉積，求三角形面積須要除以2
double Cross(Vector a,Vector b)
{
return a.x*b.y-b.x*a.y;
}
//兩向量求點積
double Dot(Vector a,Vector b)
{
return a.x*b.x+a.y*b.y;
}
//求最大最小值
double max(double a,double b)
{
if (a<b) return b; else return a;
}
double min(double a,double b)
{
if (a>b) return b; else return a;
}
//採用eps的精度判斷大/小於零
int epssgn(double x)
{
if (fabs(x)<eps) return 0;
else return x<0?-1:1;
}
//求兩點之間的直線距離
double dis(Point a,Point b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
//判斷兩直線是否平行
int Parallel(Line l1,Line l2)
{
return (fabs(l1.a*l2.b-l2.a*l1.b)<eps);
}
//判斷兩直線是否相等
int LineEqual(Line l1,Line l2)
{
if (!Parallel(l1,l2)) return 0;
else return (fabs(l1.a*l2.c-l2.a*l1.c)<eps && fabs(l1.b*l2.c-l2.b*l1.c)<eps);
}
//求點C到線段AB的距離
double PointToSegDist(Point A,Point B,Point C)
{
if (dis(A,B)<eps) return dis(B,C);
if (epssgn(Dot(B-A,C-A))<0) return dis(A,C);
if (epssgn(Dot(A-B,C-B))<0) return dis(B,C);
return fabs(Cross(B-A,C-A))/dis(B,A);
}
//求線段兩端AB到另外一線段CD的距離
double TwoSegMinDist(Point A,Point B,Point C,Point D)
{
return min(min(PointToSegDist(A,B,C),PointToSegDist(A,B,D)),
min(PointToSegDist(C,D,A),PointToSegDist(C,D,B)));
}
Point SymPoint(Point p,Line l) //求二維平面上點p關於直線p1p2的對稱點
{
Point result;
double a=l.p2.x-l.p1.x;
double b=l.p2.y-l.p1.y;
double t=((p.x-l.p1.x)*a+(p.y-l.p1.y)*b)/(a*a+b*b);
result.x=2*l.p1.x+2*a*t-p.x;
result.y=2*l.p1.y+2*b*t-p.y;
return result;
}
//判斷線段s1e1與s2e2是否相交（含端點）
//不含端點的話將下面的<=改爲<
int IsSegmentIntersect(Point s1, Point e1, Point s2, Point e2)
{
if( min(s1.x,e1.x)<= max(s2.x,e2.x) &&
min(s1.y,e1.y)<= max(s2.y,e2.y) &&
min(s2.x,e2.x)<= max(s1.x,e1.x) &&
min(s2.y,e2.y)<= max(s1.y,e1.y) &&
Cross(s2,e2,s1)*Cross(s2,e2,e1)<=0 &&
Cross(s1,e1,s2)*Cross(s1,e1,e2)<=0)
return 1;
return 0;
}
//知道直線上兩點p1p2，判斷直線與線段se是否相交，含頂點
int IsLineIntersectSegment(Point p1,Point p2,Point s,Point e)
{
if (Cross(p1,p2,s)*Cross(p1,p2,e)>eps) return 0;
else return 1;
}
int IsLineIntersectSegment(Line l1,Point s,Point e)
{
if (Cross(l1.p1,l1.p2,s)*Cross(l1.p1,l1.p2,e)>eps) return 0;
else return 1;
}
//求兩條直線l1和l2的交點
Point GetIntersect(Line l1, Line l2)
{
Point res;
res.x=(l1.b*l2.c-l2.b*l1.c)/(l1.a*l2.b-l2.a*l1.b);
res.y=(l1.c*l2.a-l2.c*l1.a)/(l1.a*l2.b-l2.a*l1.b);
return res;
}
//求質量均勻分佈的多邊形的重心
Point BaryCenter(Point *p,int n)
{
Point res(0,0);
double s=0.0,t;
int i;
for (i=1;i<n-1;i++)
{
t=Cross(p[i]-p[0],p[i+1]-p[0])/2; //分割成三角形，算面積
s+=t;
res.x+=(p[0].x+p[i].x+p[i+1].x)*t; //將面積做爲重量放在三角形的重心上
res.y+=(p[0].y+p[i].y+p[i+1].y)*t; //質量*座標，三角形重心所需的除以3放在後面
}
res.x/=(3*s);
res.y/=(3*s);
return res;
}
//求多邊形面積，要求點集按逆時針順序
double ConvexPolygonArea(Point *p,int n)
{
int i;
double area=0;
for (i=1;i<n-1;i++) area+=Cross(p[i]-p[0],p[i+1]-p[0]);
return area/2;
}
//求以原點爲圓心，過a、b，半徑爲r的扇形面積
double SectorArea(Point a,Point b,double r)
{
double theta=atan2(a.y,a.x)-atan2(b.y,b.x);
while (theta<=0) theta+=2*pi;
while (theta>=2*pi) theta-=2*pi;
theta=min(theta,2*pi-theta);
return r*r*theta/2;
}
//求以o爲圓心r爲半徑的圓與線段（直線）ab的交點
//返回的ret是交點，num是交點數量
//把判斷t範圍的兩個if去掉能夠計算圓與直線的交點
void CircleCrossLine(Point a,Point b,Point o,double r,Point ret[],int &num)
{
double x0=o.x,y0=o.y;
double x1=a.x,y1=a.y;
double x2=b.x,y2=b.y;
double dx=x2-x1,dy=y2-y1;
double A=dx*dx+dy*dy;
double B=2*dx*(x1-x0)+2*dy*(y1-y0);
double C=(x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)-r*r;
double delta=B*B-4*A*C;
num=0;
if (epssgn(delta)>=0)
{
double t1=(-B-sqrt(fabs(delta)))/(2*A);
double t2=(-B+sqrt(fabs(delta)))/(2*A);
if (epssgn(t1-1.0)<=0 && epssgn(t1)>=0)
ret[num++]=Point(x1+t1*dx,y1+t1*dy);
if (epssgn(t2-1.0)<=0 && epssgn(t2)>=0)
ret[num++]=Point(x1+t2*dx,y1+t2*dy);
}
}
//求以原點爲圓心半徑爲r的園，與原點、a、b組成的三角形的重疊面積
double Calc(Point a,Point b,double r)
{
Point p[2];
int num=0;
int ina=(epssgn(len(a)-r)<0);
int inb=(epssgn(len(b)-r)<0);
if (ina)
{
if (inb) return fabs(Cross(a,b))/2.0; //情形1：三角形徹底在圓內：直接求三角形面積
else //情形2：一個點在圓內，一個點在圓外：分割成一個徹底在圓內的三角形和一個扇形
{
CircleCrossLine(a,b,Point(0,0),r,p,num);
return SectorArea(b,p[0],r)+fabs(Cross(a,p[0]))/2.0;
}
}
else
{
CircleCrossLine(a,b,Point(0,0),r,p,num);
//情形2：一個點在圓內，一個點在圓外：分割成一個徹底在圓內的三角形和一個扇形
if (inb) return SectorArea(p[0],a,r)+fabs(Cross(p[0],b))/2.0;
else
{
//情形4：兩個點都在圓外，可是兩點的連線與圓有兩個焦點：分割成一個徹底在圓內的三角形和兩個扇形
if (num==2) return SectorArea(a,p[0],r)+SectorArea(p[1],b,r)+fabs(Cross(p[0],p[1]))/2.0;
//情形3：兩個點都在圓外，而且相交處爲一個扇形：直接求扇形面積
else return SectorArea(a,b,r);
}
}
}
//求質量均勻分佈的三角形的重心
Point TriangleMassCenter(Point a,Point b,Point c)
{
return (a+b+c)/3.0;
}
//求三角形的外心
Point TriangleCircumCenter(Point a,Point b,Point c)
{
Point res;
double a1=atan2(b.y-a.y,b.x-a.x)+pi/2;
double a2=atan2(c.y-b.y,c.x-b.x)+pi/2;
double ax=(a.x+b.x)/2;
double ay=(a.y+b.y)/2;
double bx=(c.x+b.x)/2;
double by=(c.y+b.y)/2;
double r1=(sin(a2)*(ax-bx)+cos(a2)*(by-ay))/(sin(a1)*cos(a2)-sin(a2)*cos(a1));
return Point(ax+r1*cos(a1),ay+r1*sin(a1));
}
//求三角形的垂心
Point TriangleOrthoCenter(Point a,Point b,Point c)
{
return TriangleMassCenter(a,b,c)*3.0-TriangleCircumCenter(a,b,c)*2.0;
}
//求三角形的心裏
Point TriangleInnerCenter(Point a,Point b,Point c)
{
Point res;
double la=len(b-c);
double lb=len(a-c);
double lc=len(a-b);
res.x=(la*a.x+lb*b.x+lc*c.x)/(la+lb+lc);
res.y=(la*a.y+lb*b.y+lc*c.y)/(la+lb+lc);
}
/*------Graham求凸包----*/
int cmp(Point &a,Point &b)
{
if (a.y==b.y) return (a.x<b.x);
return (a.y<b.y);
}
//對全部的點進行一次排序
void QSort(Point p[],int l,int r)
{
int i=l,j=r;
Point mid=p[(l+r)/2],swap;
while (i<=j)
{
while (cmp(p[i],mid)) i++;
while (cmp(mid,p[j])) j--;
if (i<=j)
{
swap=p[i];
p[i]=p[j];
p[j]=swap;
i++;j--;
}
}
if (i<r) QSort(p,i,r);
if (l<j) QSort(p,l,j);
}
//n爲原圖的點數，p[]爲原圖的點（0~n-1），top爲凸包點的數量（0~top-1），res[]爲凸包點的下標，。
int Graham(Point p[],int n,int res[])
{
int i,len,top;
top=1;
QSort(p,0,n-1);
for (i=0;i<3;i++) res[i]=i;
for (i=2;i<n;i++)
{
while (top && epssgn(Cross(p[i],p[res[top]],p[res[top-1]]))>=0) top--;
res[++top]=i;
}
len=top;
res[++top]=n-2;
for (i=n-3;i>=0;i--)
{
while (top!=len && epssgn(Cross(p[i], p[res[top]], p[res[top-1]]))>=0) top--;
res[++top]=i;
}
return top;
}
//判斷點x是否在凸包p[res[1~chnum]]中，返回1則爲在內部或邊界上
int PointInConvexHull(Point p[],int res[],int chnum,Point x)
{
//找一個凸包內的點
Point g=(p[res[0]]+p[res[chnum/3]]+p[res[2*chnum/3]])/3.0;
int l=0,r=chnum,mid;
//二分凸包
while (l+1<r)
{
mid=(l+r)/2;
if (epssgn(Cross(p[res[l]]-g,p[res[mid]]-g))>0)
{
if (epssgn(Cross(p[res[l]]-g,x-g))>=0
&& epssgn(Cross(p[res[mid]]-g,x-g))<0) r=mid;
else l=mid;
}
else
{
if (epssgn(Cross(p[res[l]]-g,x-g))<0
&& epssgn(Cross(p[res[mid]]-g,x-g))>=0) l=mid;
else r=mid;
}
}
r%=chnum;
if (epssgn(Cross(p[res[r]]-x,p[res[l]]-x))==-1) return 1; else return 0;
}
//旋轉卡殼求凸包p[res[1~chnum]]的直徑，對踵點數量appnum，存儲在app[][2]中
double Diameter(Point p[],int res[],int chnum,int app[][2],int &appnum) //AntiPodal Point
{
int i,j;
double ret=0,nowlen;
res[chnum]=res[0];
j=1;
appnum=0;
for (i=0;i<chnum;i++)
{
while (Cross(p[res[i]]-p[res[i+1]],p[res[j+1]]-p[res[i+1]])
<Cross(p[res[i]]-p[res[i+1]],p[res[j]]-p[res[i+1]]))
{
j++;
j%=chnum;
}
app[appnum][0]=res[i];
app[appnum][1]=res[j];
appnum++;
nowlen=dis(p[res[i]],p[res[j]]);
if (nowlen>ret) ret=nowlen;
nowlen=dis(p[res[i+1]],p[res[j+1]]);
if (nowlen>ret) ret=nowlen;
}
return ret;
}
//旋轉卡殼求凸包p[res[1~chnum]]內最大面積三角形
double ConvexHullMaxTriangleArea(Point p[],int res[],int chnum)
{
int i,j,k;
double area=0,tmp;
k=2;j=1;
res[chnum]=res[0];
for (i=0;i<chnum;i++)
{
//卡住i,j，若向前旋轉k能令面積更大的話就一直轉一下去
while (fabs(Cross(p[res[j]]-p[res[i]],p[res[(k+1)%chnum]]-p[res[i]]))
>fabs(Cross(p[res[j]]-p[res[i]],p[res[k]] -p[res[i]])))
k=(k+1)%chnum;
tmp=fabs(Cross(p[res[j]]-p[res[i]],p[res[k]]-p[res[i]]));
if (tmp>area) area=tmp;
//卡住i,k，若向前旋轉j能令面積更大的話就一直轉一下去
while (fabs(Cross(p[res[(j+1)%chnum]]-p[res[i]],p[res[k]]-p[res[i]]))
>fabs(Cross(p[res[j]]- p[res[i]],p[res[k]]-p[res[i]])))
j=(j+1)%chnum;
tmp=fabs(Cross(p[res[j]]-p[res[i]],p[res[k]]-p[res[i]]));
if (tmp>area) area=tmp;
//j,k轉到最大位置，向前轉i(i++)
}
return area/2;
}
//求兩凸包p、q間的最小距離，注意調用的時候要交換參數位置調用兩次
//卡住p（p步進），旋轉q，求能獲得的最小距離
double TwoConvexHullMinDist(Point P[],Point Q[],int n,int m)
{
int YMinP=0,YMaxQ=0,i;
double tmp,ans=999999999;
for (i=0;i<n;i++) if (P[i].y<P[YMinP].y) YMinP=i;
for (i=0;i<m;i++) if (Q[i].y>Q[YMaxQ].y) YMaxQ=i;
P[n]=P[0];
Q[m]=Q[0];
for (i=0;i<n;i++)
{
//注意，tmp保存的是>比較的結果
while (tmp=Cross(Q[YMaxQ+1]-P[YMinP+1],P[YMinP]-P[YMinP+1])
>Cross(Q[YMaxQ] -P[YMinP+1],P[YMinP]-P[YMinP+1]))
YMaxQ=(YMaxQ+1)%m;
if (tmp<0) ans=min(ans,PointToSegDist(P[YMinP],P[YMinP+1],Q[YMaxQ]));
else ans=min(ans,TwoSegMinDist (P[YMinP],P[YMinP+1],Q[YMaxQ],Q[YMaxQ+1]));
YMinP=(YMinP+1)%n;
}
return ans;
}
/*------Graham求凸包----*/
/*-----半平面交O(nlogn)模板-----*/
int cmp(const Line & l1,const Line & l2)
{
int d=epssgn(l1.angle-l2.angle);
if (!d) return (epssgn(Cross(l2.p1-l1.p1,l2.p2-l1.p1))>0); //極角相同時，將更靠半平面裏面的放在前面
return d<0;
}
void QSort(Line L[],int l,int r)
{
int i=l,j=r;
Line swap,mid=L[(l+r)/2];
while (i<=j)
{
while (cmp(L[i],mid)) i++;
while (cmp(mid,L[j])) j--;
if (i<=j)
{
swap=L[i];
L[i]=L[j];
L[j]=swap;
i++;j--;
}
}
if (i<r) QSort(L,i,r);
if (l<j) QSort(L,l,j);
}
//判斷l1與l2的交點是否在半平面hpl外
int IntersectionOutOfHalfPlane(Line &hpl,Line &l1,Line &l2)
{
Point p=GetIntersect(l1,l2);
return (epssgn(Cross(hpl.p1-p,hpl.p2-p))<0);
}
//半平面hpl向內推動dis長度
Line HalfPlaneMoveIn(Line &hpl,double &dis)
{
double dx=hpl.p1.x-hpl.p2.x;
double dy=hpl.p1.y-hpl.p2.y;
double ll=len(hpl.p1-hpl.p2);
Point pa=Point(dis*dy/ll+hpl.p1.x,hpl.p1.y-dis*dx/ll);
Point pb=Point(dis*dy/ll+hpl.p2.x,hpl.p2.y-dis*dx/ll);
return Line(pa,pb);
}
//求n個半平面l的半平面交，獲得的交點儲存在p中，交點數目返回到pn
//能夠將一個多邊形每一條邊當作半平面，求出來的交就是多邊形的核，要求pn>=3
void HalfPlaneIntersect(Line l[],int n,Point p[],int &pn)
{
int i,j;
int dq[MaxN],top,bot;
//排序是在知足全部半平面A*x+B*y+C>0或（<,<=,>=）,
//也就是全部半平面的符號均相同的狀況下對極角進行排序。
QSort(l,0,n-1);
//極角相同時，只保留最靠裏面的那條
for (i=j=0;i<n;i++) if (epssgn(l[i].angle-l[j].angle)>0) l[++j]=l[i];
n=j+1;
dq[0]=0; //雙端隊列
dq[1]=1;
top=1; //頂部和底部
bot=0;
for (i=2;i<n;i++)
{
//當棧頂的兩條直線交點在當前半平面外部時，彈棧
while (top>bot && IntersectionOutOfHalfPlane(l[i],l[dq[top]],l[dq[top-1]])) top--;
//因爲求的是一個凸多邊形，因此當半平面轉過接近一圈時，某個半平面知足上一個while的條件後，
//它又會影響到底部的兩條直線，當底部的兩條直線的交點，在當前的半平面外部時，底部彈棧
while (top>bot && IntersectionOutOfHalfPlane(l[i],l[dq[bot]],l[dq[bot+1]])) bot++;
dq[++top]=i; //當前半平面入棧
}
//當最頂部的兩條直線的交點不在最底部的半平面內時，頂部的那個半平面是多餘的，頂部彈棧
while (top>bot && IntersectionOutOfHalfPlane(l[dq[bot]],l[dq[top]],l[dq[top-1]])) top--;
//當最底部的兩條直線的交點不在最頂部的半平面內時，底部的那個半平面是多餘的，底部彈棧
while (top>bot && IntersectionOutOfHalfPlane(l[dq[top]],l[dq[bot]],l[dq[bot+1]])) bot++;
dq[++top]=dq[bot]; //將最底部的半平面放到最頂部來，方便下面求頂點
for (pn=0,i=bot;i<top;i++,pn++) p[pn]=GetIntersect(l[dq[i+1]],l[dq[i]]);
}
/*-----半平面交O(nlogn)模板-----*/
/*-----半平面交O(n^2)模板-----*/
//p是如今切出來的半平面的點，pn是點的數量，須要按順時針或者逆時針排序
//新的半平面直線爲ax+by+c=0
void HalfPlaneCut(Point p[],int &pn,double a,double b,double c)
{
int cnt=0,i;
Point tp[MaxN]; //temp_p
//如今交出來的點在新的半平面內，保留
for (i=1;i<=pn;i++) if (epssgn(a*p[i].x+b*p[i].y+c)>=0) tp[++cnt]=p[i];
else //不然若是其先後的點在半平面內，從新切割
{
if (epssgn(a*p[i-1].x+b*p[i-1].y+c)>0)
tp[++cnt]=GetIntersect(Line(p[i-1],p[i]),Line(a,b,c));
if (epssgn(a*p[i+1].x+b*p[i+1].y+c)>0)
tp[++cnt]=GetIntersect(Line(p[i],p[i+1]),Line(a,b,c));
}
pn=cnt;
for (i=1;i<=cnt;i++) p[i]=tp[i];
p[0]=p[cnt];
p[cnt+1]=p[1];
}
/*-----半平面交O(n^2)模板-----*/