【學習筆記】計算幾何全家桶
本來是不想碼的，但總是忘記一些基本操作，還是記下來比較好。
一：【准備工作】
#define LD double
#define Vector Point
#define Re register int
const LD eps=1e-8;//據說：出題的大學生們基本上都是用的這個值
inline int dcmp(LD a){return a<-eps?-1:(a>eps?1:0);}//處理精度
inline LD Abs(LD a){return a*dcmp(a);}//取絕對值
struct Point{
LD x,y;Point(LD X=0,LD Y=0){x=X,y=Y;}
inline void in(){scanf("%lf%lf",&x,&y);}
inline void out(){printf("%.2lf %.2lf\n",x,y);}
};
二：【向量】
1.【模長】
對於 \(\vec{a}=(x,y),\) \(|\vec{a}|=\sqrt{x^2+y^2}\) \(=\sqrt{|\vec{a}|^{2}}\) \(=\sqrt{\vec{a} \cdot \vec{a}}\) 。
inline LD Len(Vector a){return sqrt(Dot(a,a));}//【模長】
2.【向量加減】
對於 \(\vec{a}=(x_{1},y_{1}),\) \(\vec{b}=(x_{2},y_{2}),\) \(\vec{a}+\vec{b}=(x_{1}+x_{2},y_{1}+y_{2})\) 。
對於 \(\vec{a}=(x_{1},y_{1}),\) \(\vec{b}=(x_{2},y_{2}),\) \(\vec{a}-\vec{b}=(x_{1}-x_{2},y_{1}-y_{2})\) 。
inline Vector operator+(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
inline Vector operator-(Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
3.【向量數乘（放縮）】
對於 \(\vec{a}=(x,y),\) \(\lambda \vec{a}=(\lambda x,\lambda y)\) 。
除法也可以理解為數乘：\(\frac{\vec{a}}{\lambda}=\frac{1}{\lambda}\vec{a}=(\frac{1}{\lambda} x,\frac{1}{\lambda} y)\) 。
inline Vector operator*(Vector a,LD b){return Vector(a.x*b,a.y*b);}
4.【點積（內積）（數量積）】
\(\vec{a} \cdot \vec{b}=|\vec{a}||\vec{b}| \cos \theta\) \((\theta=\langle\vec{a}, \vec{b}\rangle)\) 。
對於 \(\vec{a}=(x_{1}, y_{1}), \vec{b}=(x_{2}, y_{2}),\) \(\vec{a} \cdot \vec{b}=x_{1} x_{2}+y_{1} y_{2}\) 。
夾角 \(\theta\) 與點積大小的關系：
\((1).\) 若 \(\theta=0^{\circ},\) \(\vec{a} \cdot \vec{b}=|\vec{a}||\vec{b}|\) 。
\((2).\) 若 \(\theta=180^{\circ},\) \(\vec{a} \cdot \vec{b}=-|\vec{a}||\vec{b}|\) 。
\((3).\) 若 \(\theta < 90^{\circ},\) \(\vec{a} \cdot \vec{b}>0\) 。
\((4).\) 若 \(\theta=90^{\circ},\) \(\vec{a} \cdot \vec{b}=0\) 。
\((5).\) 若 \(\theta > 90^{\circ},\) \(\vec{a} \cdot \vec{b}<0\) 。
inline LD Dot(Vector a,Vector b){return a.x*b.x+a.y*b.y;}//【點積】
5.【叉積（外積）（向量積）】
\(\vec{a} \times \vec{b}=|\vec{a}||\vec{b}| \cos \theta\) \((\theta=\langle\vec{a}, \vec{b}\rangle)\) 。
對於 \(\vec{a}=(x_{1}, y_{1}), \vec{b}=(x_{2}, y_{2}),\) \(\vec{a} \times \vec{b}=x_{1} y_{2}-y_{1} x_{2}\) 。
向量位置與叉積大小的關系：
\((1).\) 若 \(\vec{a} \| \vec{b},\) \(\vec{a} \times \vec{b}=0\) 。
\((2).\) 若 \(\vec{a}\) 在 \(\vec{b}\) 右側，\(\vec{a} \times \vec{b}>0\) 。
\((3).\) 若 \(\vec{a}\) 在 \(\vec{b}\) 左側，\(\vec{a} \times \vec{b}<0\) 。
inline LD Cro(Vector a,Vector b){return a.x*b.y-a.y*b.x;}//【叉積】
三：【點、向量的位置變換】
1.【點、向量的旋轉】
\((1).\) 對於點 \(P=(x,y)\) 或向量 \(\vec{a}=(x,y)\)，將其順時針旋轉 \(\theta\) 角度（點：關於原點，向量：關於起點）： \(\begin{vmatrix}x&y\end{vmatrix}\) \(\times\) \(\begin{vmatrix}cos \theta & -sin \theta\\ sin \theta & cos \theta \end{vmatrix}\) \(=\) \(\begin{vmatrix}xcos \theta +ysin \theta &-xsin \theta + ycos \theta \end{vmatrix}\) 。
inline Point turn_P(Point a,LD theta){//【點A\向量A順時針旋轉theta(弧度)】
LD x=a.x*cos(theta)+a.y*sin(theta);
LD y=-a.x*sin(theta)+a.y*cos(theta);
return Point(x,y);
}
\((2).\) 將點 \(A(x,y)\) 繞點 \(B(x_0,y_0)\) 順時針旋轉 \(\theta\) 角度：\(\begin{vmatrix}(x\!-\!x_0)cos \theta +(y\!-\!y_0)sin \theta + x_0 &-(x\!-\!x_0)sin \theta + (y\!-\!y_0)cos \theta + y_0 \end{vmatrix}\) 。
inline Point turn_PP(Point a,Point b,LD theta){//【將點A繞點B順時針旋轉theta(弧度)】
LD x=(a.x-b.x)*cos(theta)+(a.y-b.y)*sin(theta)+b.x;
LD y=-(a.x-b.x)*sin(theta)+(a.y-b.y)*cos(theta)+b.y;
return Point(x,y);
}
四：【圖形與圖形之間的關系】
1.【點與線段】
\((1).\) 判斷點 \(P\) 是否在線段 \(AB\) 上：
inline int pan_PL(Point p,Point a,Point b){//【判斷點P是否在線段AB上】
return !dcmp(Cro(p-a,b-a))&&dcmp(Dot(p-a,p-b))<=0;//做法一
// return !dcmp(Cro(p-a,b-a))&&dcmp(min(a.x,b.x)-p.x)<=0&&dcmp(p.x-max(a.x,b.x))<=0&&dcmp(min(a.y,b.y)-p.y)<=0&&dcmp(p.y-max(a.y,b.y))<=0;//做法二
//PA,AB共線且P在AB之間(其實也可以用len(p-a)+len(p-b)==len(a-b)判斷，但是精度損失較大)
}
\((2).\) 點 \(P\) 到線段 \(AB\) 的距離：
inline bool operator==(Point a,Point b){return !dcmp(a.x-b.x)&&!dcmp(a.y-b.y);}//兩點坐標重合則相等
inline LD dis_PL(Point p,Point a,Point b){//【點P到線段AB距離】
if(a==b)return Len(p-a);//AB重合
Vector x=p-a,y=p-b,z=b-a;
if(dcmp(Dot(x,z))<0)return Len(x);//P距離A更近
if(dcmp(Dot(y,z))>0)return Len(y);//P距離B更近
return Abs(Cro(x,z)/Len(z));//面積除以底邊長
}
2.【點與直線】
\((1).\) 判斷點 \(P\) 是否在直線 \(AB\) 上：
inline int pan_PL_(Point p,Point a,Point b){//【判斷點P是否在直線AB上】
return !dcmp(Cro(p-a,b-a));//PA,AB共線
}
\((2).\) 點 \(P\) 到直線 \(AB\) 的垂足 \(F\)：
inline Point FootPoint(Point p,Point a,Point b){//【點P到直線AB的垂足】
Vector x=p-a,y=p-b,z=b-a;
LD len1=Dot(x,z)/Len(z),len2=-1.0*Dot(y,z)/Len(z);//分別計算AP,BP在AB,BA上的投影
return a+z*(len1/(len1+len2));//點A加上向量AF
}
\((3).\) 點 \(P\) 關於直線 \(AB\) 的對稱點：
inline Point Symmetry_PL(Point p,Point a,Point b){//【點P關於直線AB的對稱點】
return p+(FootPoint(p,a,b)-p)*2;//將PF延長一倍即可
}
3.【線與線】
\((1).\) 兩直線 \(AB,CD\) 的交點 \(Q\)：
inline Point cross_LL(Point a,Point b,Point c,Point d){//【兩直線AB,CD的交點】
Vector x=b-a,y=d-c,z=a-c;
return a+x*(Cro(y,z)/Cro(x,y));//點A加上向量AF
}
\((2).\) 判斷直線 \(AB\) 與線段 \(CD\) 是否相交：
inline int pan_cross_L_L(Point a,Point b,Point c,Point d){//【判斷直線AB與線段CD是否相交】
return pan_PL(cross_LL(a,b,c,d),c,d);//直線AB與直線CD的交點在線段CD上
}
\((3).\) 判斷兩線段 \(AB,CD\) 是否相交：
inline int pan_cross_LL(Point a,Point b,Point c,Point d){//【判斷兩線段AB,CD是否相交】
LD c1=Cro(b-a,c-a),c2=Cro(b-a,d-a);
LD d1=Cro(d-c,a-c),d2=Cro(d-c,b-c);
return dcmp(c1)*dcmp(c2)<0&&dcmp(d1)*dcmp(d2)<0;//分別在兩側
}
4.【點與多邊形】
\((1).\) 判斷點 \(A\) 是否在任意多邊形 \(Poly\) 以內（射線法）：
inline int PIP(Point *P,Re n,Point a){//【射線法】判斷點A是否在任意多邊形Poly以內
Re cnt=0;LD tmp;
for(Re i=1;i<=n;++i){
Re j=i<n?i+1:1;
if(pan_PL(a,P[i],P[j]))return 2;//點在多邊形上
if(a.y>=min(P[i].y,P[j].y)&&a.y<max(P[i].y,P[j].y))//縱坐標在該線段兩端點之間
tmp=P[i].x+(a.y-P[i].y)/(P[j].y-P[i].y)*(P[j].x-P[i].x),cnt+=dcmp(tmp-a.x)>0;//交點在A右方
}
return cnt&1;//穿過奇數次則在多邊形以內
}
\((2).\) 判斷點 \(A\) 是否在凸多邊形 \(Poly\) 以內（二分法）：
inline int judge(Point a,Point L,Point R){//判斷AL是否在AR右邊
return dcmp(Cro(L-a,R-a))>0;//必須嚴格以內
}
inline int PIP_(Point *P,Re n,Point a){//【二分法】判斷點A是否在凸多邊形Poly以內
//點按逆時針給出
if(judge(P[1],a,P[2])||judge(P[1],P[n],a))return 0;//在P[1_2]或P[1_n]外
if(pan_PL(a,P[1],P[2])||pan_PL(a,P[1],P[n]))return 2;//在P[1_2]或P[1_n]上
Re l=2,r=n-1;
while(l<r){//二分找到一個位置pos使得P[1]_A在P[1_pos],P[1_(pos+1)]之間
Re mid=l+r+1>>1;
if(judge(P[1],P[mid],a))l=mid;
else r=mid-1;
}
if(judge(P[l],a,P[l+1]))return 0;//在P[pos_(pos+1)]外
if(pan_PL(a,P[l],P[l+1]))return 2;//在P[pos_(pos+1)]上
return 1;
}
5.【線與多邊形】
\((1).\) 判斷線段 \(AB\) 是否在任意多邊形 \(Poly\) 以內：不相交且兩端點 \(A,B\) 均在多邊形以內。
\((2).\) 判斷線段 \(AB\) 是否在凸多邊形 \(Poly\) 以內：兩端點 \(A,B\) 均在多邊形以內。
6.【多邊形與多邊形】
\((1).\) 判斷任意兩個多邊形是否相離：屬於不同多邊形的任意兩邊都不相交 且 一個多邊形上的任意點都不被另一個多邊形所包含。
inline int judge_PP(Point *A,Re n,Point *B,Re m){//【判斷多邊形A與多邊形B是否相離】
for(Re i1=1;i1<=n;++i1){
Re j1=i1<n?i1+1:1;
for(Re i2=1;i2<=m;++i2){
Re j2=i2<m?i2+1:1;
if(pan_cross_LL(A[i1],A[j1],B[i2],B[j2]))return 0;//兩線段相交
if(PIP(B,m,A[i1])||PIP(A,n,B[i2]))return 0;//點包含在內
}
}
return 1;
}
五：【圖形面積】
1.【任意多邊形面積】
inline LD PolyArea(Point *P,Re n){//【任意多邊形P的面積】
LD S=0;
for(Re i=1;i<=n;++i)S+=Cro(P[i],P[i<n?i+1:1]);
return S/2.0;
}
2.【圓的面積並】
自適應辛普森法亂搞。
- 【模板】 \(\text{CIRU - The area of the union of circles [SP8073]}\) \ 圓的面積並 \(\text{[Bzoj2178]}\)
【題解】 \(\text{Xing_Ling}\)
3.【三角形面積並】
自適應辛普森法亂搞。
或者掃描線？wo太菜了不會寫。
六：【凸包】
1.【求凸包】
\((1).\) \(\text{Graham}\) 掃描法
inline bool cmp1(Vector a,Vector b){return a.x==b.x?a.y<b.y:a.x<b.x;};//按坐標排序
inline int ConvexHull(Point *P,Re n,Point *cp){//【Graham掃描法】求凸包
sort(P+1,P+n+1,cmp1);
Re t=0;
for(Re i=1;i<=n;++i){//下凸包
while(t>1&&dcmp(Cro(cp[t]-cp[t-1],P[i]-cp[t-1]))<=0)--t;
cp[++t]=P[i];
}
Re St=t;
for(Re i=n-1;i>=1;--i){//上凸包
while(t>St&&dcmp(Cro(cp[t]-cp[t-1],P[i]-cp[t-1]))<=0)--t;
cp[++t]=P[i];
}
return --t;//要減一
}
2.【旋轉卡殼】
（這玩意兒有 \(2*2*2*3=24\) 種讀音）
Rd Ans=Len(cp[2]-cp[1]);
for(Re i=1,j=3;i<=cnt;++i){
while(dcmp(Cro(cp[i+1]-cp[i],cp[j]-cp[i])-Cro(cp[i+1]-cp[i],cp[j+1]-cp[i]))<0)j=j<cnt?j+1:1;//注意是<0，如果寫<=0的話可能會被兩個點的數據卡掉
Ans=max(Ans,max(Len(cp[j]-cp[i]),Len(cp[j]-cp[i+1])));//求最遠距離
}
3.【半平面交】
struct Line{
Point a,b;LD k;Line(Point A=Point(0,0),Point B=Point(0,0)){a=A,b=B,k=atan2(b.y-a.y,b.x-a.x);}
inline bool operator<(const Line &O)const{return dcmp(k-O.k)?dcmp(k-O.k)<0:judge(O.a,O.b,a);}//如果角度相等則取左邊的
}L[N],Q[N];
inline Point cross(Line L1,Line L2){return cross_LL(L1.a,L1.b,L2.a,L2.b);}//獲取直線L1,L2的交點
inline int judge(Line L,Point a){return dcmp(Cro(a-L.a,L.b-L.a))>0;}//判斷點a是否在直線L的右邊
inline int halfcut(Line *L,Re n,Point *P){//【半平面交】
sort(L+1,L+n+1);Re m=n;n=0;
for(Re i=1;i<=m;++i)if(i==1||dcmp(L[i].k-L[i-1].k))L[++n]=L[i];
Re h=1,t=0;
for(Re i=1;i<=n;++i){
while(h<t&&judge(L[i],cross(Q[t],Q[t-1])))--t;//當隊尾兩個直線交點不是在直線L[i]上或者左邊時就出隊
while(h<t&&judge(L[i],cross(Q[h],Q[h+1])))++h;//當隊頭兩個直線交點不是在直線L[i]上或者左邊時就出隊
Q[++t]=L[i];
}
while(h<t&&judge(Q[h],cross(Q[t],Q[t-1])))--t;
while(h<t&&judge(Q[t],cross(Q[h],Q[h+1])))++h;
n=0;
for(Re i=h;i<=t;++i)P[++n]=cross(Q[i],Q[i<t?i+1:h]);
return n;
}
4.【閔可夫斯基和】
Vector V1[N],V2[N];
inline int Mincowski(Point *P1,Re n,Point *P2,Re m,Vector *V){//【閔可夫斯基和】求兩個凸包{P1},{P2}的向量集合{V}={P1+P2}構成的凸包
for(Re i=1;i<=n;++i)V1[i]=P1[i<n?i+1:1]-P1[i];
for(Re i=1;i<=m;++i)V2[i]=P2[i<m?i+1:1]-P2[i];
Re t=0,i=1,j=1;V[++t]=P1[1]+P2[1];
while(i<=n&&j<=m)++t,V[t]=V[t-1]+(dcmp(Cro(V1[i],V2[j]))>0?V1[i++]:V2[j++]);
while(i<=n)++t,V[t]=V[t-1]+V1[i++];
while(j<=m)++t,V[t]=V[t-1]+V2[j++];
return t;
}
5.【動態凸包】
struct ConvexHull{
int op;set<Point>s;
inline int PIP(Point P){
IT it=s.lower_bound(Point(P.x,-inf));//找到第一個橫坐標大於P的點
if(it==s.end())return 0;
if(it->x==P.x)return (P.y-it->y)*op<=0;//比較縱坐標大小
if(it==s.begin())return 0;
IT j=it,k=it;--j;return Cro(P-*j,*k-*j)*op>=0;//看叉姬
}
inline int judge(IT it){
IT j=it,k=it;
if(j==s.begin())return 0;--j;
if(++k==s.end())return 0;
return Cro(*it-*j,*k-*j)*op>=0;//看叉姬
}
inline void insert(Point P){
if(PIP(P))return;//如果點P已經在凸殼上或凸包里就不插入了
IT tmp=s.lower_bound(Point(P.x,-inf));if(tmp!=s.end()&&tmp->x==P.x)s.erase(tmp);//特判橫坐標相等的點要去掉
s.insert(P);IT it=s.find(P),p=it;
if(p!=s.begin()){--p;while(judge(p)){IT tmp=p--;s.erase(tmp);}}
if((p=++it)!=s.end())while(judge(p)){IT tmp=p++;s.erase(tmp);}
}
}up,down;
int x,y,T,op;
int main(){
// freopen("123.txt","r",stdin);
in(T),up.op=1,down.op=-1;
while(T--){
in(op),P.In();
if(op<2)up.insert(P),down.insert(P);//插入點P
else puts((up.PIP(P)&&down.PIP(P))?"YES":"NO");//判斷點P是否在凸包內
}
}
七：【圓】
1.【三點確定一圓】
\((1).\) 暴力解方程：
設 \(x^2+y^2+Dx+Ey+F=0\)，圓心為 \(O\)，半徑為 \(r\)，帶入三點 \(A(x_1,y_1),B(x_2,y_2),C(x_3,y_3)\)，解得：
\(\begin{cases} D=\frac{(x_2^2+y_2^2-x_3^2-y_3^2)(y_1-y_2)-(x_1^2+y_1^2-x_2^2-y_2^2)(y_2-y_3)}{(x_1-x_2)(y_2-y_3)-(x_2-x_3)(y_1-y_2)} \\ E=\frac{x_1^2+y_1^2-x_2^2-y_2^2+D(x_1-x_2)}{y_2-y_1} \\ F=-(x_1^2+y_1^2+Dx_1+Ey_1) \\ O=(-\frac{D}{2},-\frac{E}{2}) \\ r=\frac{D^2+E^2-4F}{4} \end{cases}\)
#define S(a) ((a)*(a))
struct Circle{Point O;LD r;Circle(Point P,LD R=0){O=P,r=R;}};
inline Circle getCircle(Point A,Point B,Point C){//【三點確定一圓】暴力解方程
LD x1=A.x,y1=A.y,x2=B.x,y2=B.y,x3=C.x,y3=C.y;
LD D=((S(x2)+S(y2)-S(x3)-S(y3))*(y1-y2)-(S(x1)+S(y1)-S(x2)-S(y2))*(y2-y3))/((x1-x2)*(y2-y3)-(x2-x3)*(y1-y2));
LD E=(S(x1)+S(y1)-S(x2)-S(y2)+D*(x1-x2))/(y2-y1);
LD F=-(S(x1)+S(y1)+D*x1+E*y1);
return Circle(Point(-D/2.0,-E/2.0),sqrt((S(D)+S(E)-4.0*F)/4.0));
}
\((2).\) 向量求三角形垂心：
inline Circle getcircle(Point A,Point B,Point C){//【三點確定一圓】向量垂心法
Point P1=(A+B)*0.5,P2=(A+C)*0.5;
Point O=cross_LL(P1,P1+Normal(B-A),P2,P2+Normal(C-A));
return Circle(O,Len(A-O));
}
2.【最小覆蓋圓】
【定理】 如果點 \(p\) 不在集合 \(\{S\}\) 的最小覆蓋圓內，則 \(p\) 一定在 \(\{S\}\cup{p}\) 的最小覆蓋圓上。
inline int PIC(Circle C,Point a){return dcmp(Len(a-C.O)-C.r)<=0;}//判斷點A是否在圓C內
inline void Random(Point *P,Re n){for(Re i=1;i<=n;++i)swap(P[i],P[rand()%n+1]);}//隨機一個排列
inline Circle Min_Circle(Point *P,Re n){//【求點集P的最小覆蓋圓】
// random_shuffle(P+1,P+n+1);
Random(P,n);Circle C=Circle(P[1],0);
for(Re i=2;i<=n;++i)if(!PIC(C,P[i])){
C=Circle(P[i],0);
for(Re j=1;j<i;++j)if(!PIC(C,P[j])){
C.O=(P[i]+P[j])*0.5,C.r=Len(P[j]-C.O);
for(Re k=1;k<j;++k)if(!PIC(C,P[k]))C=getcircle(P[i],P[j],P[k]);
}
}
return C;
}
3.【三角剖分】
inline LD calc(Point A,Point B,Point O,LD R){//【三角剖分】
if(A==O||B==O)return 0;
Re op=dcmp(Cro(A-O,B-O))>0?1:-1;LD ans=0;
Vector x=A-O,y=B-O;
Re flag1=dcmp(Len(x)-R)>0,flag2=dcmp(Len(y)-R)>0;
if(!flag1&&!flag2)ans=Abs(Cro(A-O,B-O))/2.0;//兩個點都在里面
else if(flag1&&flag2){//兩個點都在外面
if(dcmp(dis_PL(O,A,B)-R)>=0)ans=R*R*Angle(x,y)/2.0;//完全包含了圓弧
else{//分三段處理 △+圓弧+△
if(dcmp(Cro(A-O,B-O))>0)swap(A,B);//把A換到左邊
Point F=FootPoint(O,A,B);LD lenx=Len(F-O),len=sqrt(R*R-lenx*lenx);
Vector z=turn_P(F-O,Pi/2.0)*(len/lenx);Point B_=F+z,A_=F-z;
ans=R*R*(Angle(A-O,A_-O)+Angle(B-O,B_-O))/2.0+Cro(B_-O,A_-O)/2.0;
}
}
else{//一個點在里面，一個點在外面
if(flag1)swap(A,B);//使A為里面的點，B為外面的點
Point F=FootPoint(O,A,B);LD lenx=Len(F-O),len=sqrt(R*R-lenx*lenx);
Vector z=turn_P(F-O,Pi/2.0)*(len/lenx);Point C=dcmp(Cro(A-O,B-O))>0?F-z:F+z;
ans=Abs(Cro(A-O,C-O))/2.0+R*R*Angle(C-O,B-O)/2.0;
}
return ans*op;
}
【計算幾何全家桶 Code】
#include<algorithm>
#include<cstdio>
#include<cmath>
/*一：【准備工作】*/
#define LD double
#define LL long long
#define Re register int
#define Vector Point
using namespace std;
const int N=262144+3;
const LD eps=1e-8,Pi=acos(-1.0);
inline int dcmp(LD a){return a<-eps?-1:(a>eps?1:0);}//處理精度
inline LD Abs(LD a){return a*dcmp(a);}//取絕對值
struct Point{
LD x,y;Point(LD X=0,LD Y=0){x=X,y=Y;}
inline void in(){scanf("%lf%lf",&x,&y);}
inline void out(){printf("%.2lf %.2lf\n",x,y);}
};
/*二：【向量】*/
inline LD Dot(Vector a,Vector b){return a.x*b.x+a.y*b.y;}//【點積】
inline LD Cro(Vector a,Vector b){return a.x*b.y-a.y*b.x;}//【叉積】
inline LD Len(Vector a){return sqrt(Dot(a,a));}//【模長】
inline LD Angle(Vector a,Vector b){return acos(Dot(a,b)/Len(a)/Len(b));}//【兩向量夾角】
inline Vector Normal(Vector a){return Vector(-a.y,a.x);}//【法向量】
inline Vector operator+(Vector a,Vector b){return Vector(a.x+b.x,a.y+b.y);}
inline Vector operator-(Vector a,Vector b){return Vector(a.x-b.x,a.y-b.y);}
inline Vector operator*(Vector a,LD b){return Vector(a.x*b,a.y*b);}
inline bool operator==(Point a,Point b){return !dcmp(a.x-b.x)&&!dcmp(a.y-b.y);}//兩點坐標重合則相等
/*三：【點、向量的位置變換】*/
/*1.【點、向量的旋轉】*/
inline Point turn_P(Point a,LD theta){//【點A\向量A順時針旋轉theta(弧度)】
LD x=a.x*cos(theta)+a.y*sin(theta);
LD y=-a.x*sin(theta)+a.y*cos(theta);
return Point(x,y);
}
inline Point turn_PP(Point a,Point b,LD theta){//【將點A繞點B順時針旋轉theta(弧度)】
LD x=(a.x-b.x)*cos(theta)+(a.y-b.y)*sin(theta)+b.x;
LD y=-(a.x-b.x)*sin(theta)+(a.y-b.y)*cos(theta)+b.y;
return Point(x,y);
}
/*四：【圖形與圖形之間的關系】*/
/*1.【點與線段】*/
inline int pan_PL(Point p,Point a,Point b){//【判斷點P是否在線段AB上】
return !dcmp(Cro(p-a,b-a))&&dcmp(Dot(p-a,p-b))<=0;//做法一
// return !dcmp(Cro(p-a,b-a))&&dcmp(min(a.x,b.x)-p.x)<=0&&dcmp(p.x-max(a.x,b.x))<=0&&dcmp(min(a.y,b.y)-p.y)<=0&&dcmp(p.y-max(a.y,b.y))<=0;//做法二
//PA,AB共線且P在AB之間(其實也可以用len(p-a)+len(p-b)==len(a-b)判斷，但是精度損失較大)
}
inline LD dis_PL(Point p,Point a,Point b){//【點P到線段AB距離】
if(a==b)return Len(p-a);//AB重合
Vector x=p-a,y=p-b,z=b-a;
if(dcmp(Dot(x,z))<0)return Len(x);//P距離A更近
if(dcmp(Dot(y,z))>0)return Len(y);//P距離B更近
return Abs(Cro(x,z)/Len(z));//面積除以底邊長
}
/*2.【點與直線】*/
inline int pan_PL_(Point p,Point a,Point b){//【判斷點P是否在直線AB上】
return !dcmp(Cro(p-a,b-a));//PA,AB共線
}
inline Point FootPoint(Point p,Point a,Point b){//【點P到直線AB的垂足】
Vector x=p-a,y=p-b,z=b-a;
LD len1=Dot(x,z)/Len(z),len2=-1.0*Dot(y,z)/Len(z);//分別計算AP,BP在AB,BA上的投影
return a+z*(len1/(len1+len2));//點A加上向量AF
}
inline Point Symmetry_PL(Point p,Point a,Point b){//【點P關於直線AB的對稱點】
return p+(FootPoint(p,a,b)-p)*2;//將PF延長一倍即可
}
/*3.【線與線】*/
inline Point cross_LL(Point a,Point b,Point c,Point d){//【兩直線AB,CD的交點】
Vector x=b-a,y=d-c,z=a-c;
return a+x*(Cro(y,z)/Cro(x,y));//點A加上向量AF
}
inline int pan_cross_L_L(Point a,Point b,Point c,Point d){//【判斷直線AB與線段CD是否相交】
return pan_PL(cross_LL(a,b,c,d),c,d);//直線AB與直線CD的交點在線段CD上
}
inline int pan_cross_LL(Point a,Point b,Point c,Point d){//【判斷兩線段AB,CD是否相交】
LD c1=Cro(b-a,c-a),c2=Cro(b-a,d-a);
LD d1=Cro(d-c,a-c),d2=Cro(d-c,b-c);
return dcmp(c1)*dcmp(c2)<0&&dcmp(d1)*dcmp(d2)<0;//分別在兩側
}
/*4.【點與多邊形】*/
inline int PIP(Point *P,Re n,Point a){//【射線法】判斷點A是否在任意多邊形Poly以內
Re cnt=0;LD tmp;
for(Re i=1;i<=n;++i){
Re j=i<n?i+1:1;
if(pan_PL(a,P[i],P[j]))return 2;//點在多邊形上
if(a.y>=min(P[i].y,P[j].y)&&a.y<max(P[i].y,P[j].y))//縱坐標在該線段兩端點之間
tmp=P[i].x+(a.y-P[i].y)/(P[j].y-P[i].y)*(P[j].x-P[i].x),cnt+=dcmp(tmp-a.x)>0;//交點在A右方
}
return cnt&1;//穿過奇數次則在多邊形以內
}
inline int judge(Point a,Point L,Point R){//判斷AL是否在AR右邊
return dcmp(Cro(L-a,R-a))>0;//必須嚴格以內
}
inline int PIP_(Point *P,Re n,Point a){//【二分法】判斷點A是否在凸多邊形Poly以內
//點按逆時針給出
if(judge(P[1],a,P[2])||judge(P[1],P[n],a))return 0;//在P[1_2]或P[1_n]外
if(pan_PL(a,P[1],P[2])||pan_PL(a,P[1],P[n]))return 2;//在P[1_2]或P[1_n]上
Re l=2,r=n-1;
while(l<r){//二分找到一個位置pos使得P[1]_A在P[1_pos],P[1_(pos+1)]之間
Re mid=l+r+1>>1;
if(judge(P[1],P[mid],a))l=mid;
else r=mid-1;
}
if(judge(P[l],a,P[l+1]))return 0;//在P[pos_(pos+1)]外
if(pan_PL(a,P[l],P[l+1]))return 2;//在P[pos_(pos+1)]上
return 1;
}
/*5.【線與多邊形】*/
/*6.【多邊形與多邊形】*/
inline int judge_PP(Point *A,Re n,Point *B,Re m){//【判斷多邊形A與多邊形B是否相離】
for(Re i1=1;i1<=n;++i1){
Re j1=i1<n?i1+1:1;
for(Re i2=1;i2<=m;++i2){
Re j2=i2<m?i2+1:1;
if(pan_cross_LL(A[i1],A[j1],B[i2],B[j2]))return 0;//兩線段相交
if(PIP(B,m,A[i1])||PIP(A,n,B[i2]))return 0;//點包含在內
}
}
return 1;
}
/*五：【圖形面積】*/
/*1.【任意多邊形面積】*/
inline LD PolyArea(Point *P,Re n){//【任意多邊形P的面積】
LD S=0;
for(Re i=1;i<=n;++i)S+=Cro(P[i],P[i<n?i+1:1]);
return S/2.0;
}
/*2.【圓的面積並】*/
/*3.【三角形面積並】*/
/*六：【凸包】*/
/*1.【求凸包】*/
inline bool cmp1(Vector a,Vector b){return a.x==b.x?a.y<b.y:a.x<b.x;};//按坐標排序
inline int ConvexHull(Point *P,Re n,Point *cp){//【Graham掃描法】求凸包
sort(P+1,P+n+1,cmp1);
Re t=0;
for(Re i=1;i<=n;++i){//下凸包
while(t>1&&dcmp(Cro(cp[t]-cp[t-1],P[i]-cp[t-1]))<=0)--t;
cp[++t]=P[i];
}
Re St=t;
for(Re i=n-1;i>=1;--i){//上凸包
while(t>St&&dcmp(Cro(cp[t]-cp[t-1],P[i]-cp[t-1]))<=0)--t;
cp[++t]=P[i];
}
return --t;//要減一
}
/*2.【旋轉卡殼】*/
/*3.【半平面交】*/
struct Line{
Point a,b;LD k;Line(Point A=Point(0,0),Point B=Point(0,0)){a=A,b=B,k=atan2(b.y-a.y,b.x-a.x);}
inline bool operator<(const Line &O)const{return dcmp(k-O.k)?dcmp(k-O.k)<0:judge(O.a,O.b,a);}//如果角度相等則取左邊的
}L[N],Q[N];
inline Point cross(Line L1,Line L2){return cross_LL(L1.a,L1.b,L2.a,L2.b);}//獲取直線L1,L2的交點
inline int judge(Line L,Point a){return dcmp(Cro(a-L.a,L.b-L.a))>0;}//判斷點a是否在直線L的右邊
inline int halfcut(Line *L,Re n,Point *P){//【半平面交】
sort(L+1,L+n+1);Re m=n;n=0;
for(Re i=1;i<=m;++i)if(i==1||dcmp(L[i].k-L[i-1].k))L[++n]=L[i];
Re h=1,t=0;
for(Re i=1;i<=n;++i){
while(h<t&&judge(L[i],cross(Q[t],Q[t-1])))--t;//當隊尾兩個直線交點不是在直線L[i]上或者左邊時就出隊
while(h<t&&judge(L[i],cross(Q[h],Q[h+1])))++h;//當隊頭兩個直線交點不是在直線L[i]上或者左邊時就出隊
Q[++t]=L[i];
}
while(h<t&&judge(Q[h],cross(Q[t],Q[t-1])))--t;
while(h<t&&judge(Q[t],cross(Q[h],Q[h+1])))++h;
n=0;
for(Re i=h;i<=t;++i)P[++n]=cross(Q[i],Q[i<t?i+1:h]);
return n;
}
/*4.【閔可夫斯基和】*/
Vector V1[N],V2[N];
inline int Mincowski(Point *P1,Re n,Point *P2,Re m,Vector *V){//【閔可夫斯基和】求兩個凸包{P1},{P2}的向量集合{V}={P1+P2}構成的凸包
for(Re i=1;i<=n;++i)V1[i]=P1[i<n?i+1:1]-P1[i];
for(Re i=1;i<=m;++i)V2[i]=P2[i<m?i+1:1]-P2[i];
Re t=0,i=1,j=1;V[++t]=P1[1]+P2[1];
while(i<=n&&j<=m)++t,V[t]=V[t-1]+(dcmp(Cro(V1[i],V2[j]))>0?V1[i++]:V2[j++]);
while(i<=n)++t,V[t]=V[t-1]+V1[i++];
while(j<=m)++t,V[t]=V[t-1]+V2[j++];
return t;
}
/*5.【動態凸包】*/
/*七：【圓】*/
/*1.【三點確定一圓】*/
#define S(a) ((a)*(a))
struct Circle{Point O;LD r;Circle(Point P,LD R=0){O=P,r=R;}};
inline Circle getCircle(Point A,Point B,Point C){//【三點確定一圓】暴力解方程
LD x1=A.x,y1=A.y,x2=B.x,y2=B.y,x3=C.x,y3=C.y;
LD D=((S(x2)+S(y2)-S(x3)-S(y3))*(y1-y2)-(S(x1)+S(y1)-S(x2)-S(y2))*(y2-y3))/((x1-x2)*(y2-y3)-(x2-x3)*(y1-y2));
LD E=(S(x1)+S(y1)-S(x2)-S(y2)+D*(x1-x2))/(y2-y1);
LD F=-(S(x1)+S(y1)+D*x1+E*y1);
return Circle(Point(-D/2.0,-E/2.0),sqrt((S(D)+S(E)-4.0*F)/4.0));
}
inline Circle getcircle(Point A,Point B,Point C){//【三點確定一圓】向量垂心法
Point P1=(A+B)*0.5,P2=(A+C)*0.5;
Point O=cross_LL(P1,P1+Normal(B-A),P2,P2+Normal(C-A));
return Circle(O,Len(A-O));
}
/*2.【最小覆蓋圓】*/
inline int PIC(Circle C,Point a){return dcmp(Len(a-C.O)-C.r)<=0;}//判斷點A是否在圓C內
inline void Random(Point *P,Re n){for(Re i=1;i<=n;++i)swap(P[i],P[rand()%n+1]);}//隨機一個排列
inline Circle Min_Circle(Point *P,Re n){//【求點集P的最小覆蓋圓】
// random_shuffle(P+1,P+n+1);
Random(P,n);Circle C=Circle(P[1],0);
for(Re i=2;i<=n;++i)if(!PIC(C,P[i])){
C=Circle(P[i],0);
for(Re j=1;j<i;++j)if(!PIC(C,P[j])){
C.O=(P[i]+P[j])*0.5,C.r=Len(P[j]-C.O);
for(Re k=1;k<j;++k)if(!PIC(C,P[k]))C=getcircle(P[i],P[j],P[k]);
}
}
return C;
}
/*3.【三角剖分】*/
inline LD calc(Point A,Point B,Point O,LD R){//【三角剖分】
if(A==O||B==O)return 0;
Re op=dcmp(Cro(A-O,B-O))>0?1:-1;LD ans=0;
Vector x=A-O,y=B-O;
Re flag1=dcmp(Len(x)-R)>0,flag2=dcmp(Len(y)-R)>0;
if(!flag1&&!flag2)ans=Abs(Cro(A-O,B-O))/2.0;//兩個點都在里面
else if(flag1&&flag2){//兩個點都在外面
if(dcmp(dis_PL(O,A,B)-R)>=0)ans=R*R*Angle(x,y)/2.0;//完全包含了圓弧
else{//分三段處理 △+圓弧+△
if(dcmp(Cro(A-O,B-O))>0)swap(A,B);//把A換到左邊
Point F=FootPoint(O,A,B);LD lenx=Len(F-O),len=sqrt(R*R-lenx*lenx);
Vector z=turn_P(F-O,Pi/2.0)*(len/lenx);Point B_=F+z,A_=F-z;
ans=R*R*(Angle(A-O,A_-O)+Angle(B-O,B_-O))/2.0+Cro(B_-O,A_-O)/2.0;
}
}
else{//一個點在里面，一個點在外面
if(flag1)swap(A,B);//使A為里面的點，B為外面的點
Point F=FootPoint(O,A,B);LD lenx=Len(F-O),len=sqrt(R*R-lenx*lenx);
Vector z=turn_P(F-O,Pi/2.0)*(len/lenx);Point C=dcmp(Cro(A-O,B-O))>0?F-z:F+z;
ans=Abs(Cro(A-O,C-O))/2.0+R*R*Angle(C-O,B-O)/2.0;
}
return ans*op;
}
int main(){}