計算幾何是一個偏冷門的板塊（基本只有 SCOI 考/gg），但是基本的了解還是要有的。

平面幾何基礎

向量與點

由于點和向量的運算和性質等的相似性，我們將其一并定義。
定義向量與點

注意到二者有類似的形式。事實上，我們經常將二者視為同一種東西，因此經常直接將點代做向量進行向量運算。例如點 \(A\)，我們可以直接將其視作 \(\vec A\) 帶入運算。
向量的基本運算分為向量與實數運算和向量間運算。
其中，向量間運算主要分為叉積和點積兩種。

叉積的幾何意義為 \(\vec a\) 與 \(\vec b\) 之間夾成的三角形的面積的兩倍，點積的幾何意義是一個向量在另一個向量上的投影長度。
若 \(\vec b\) 在 \(\vec a\) 的逆時針方向，叉積為正，反之為負，共線為 0。這一個性質相對關鍵，比較常用。
而如果兩個向量的夾角為鈍角，則點積為負，反之為正，恰好垂直則為 0。

在旋轉之前，我們先介紹極角。極角是指 \(x\) 軸正半軸與向量逆時針方向夾成的角。這個還是比較好理解的。

向量的旋轉可以通過和角公式推導。設旋轉角度為 \(\theta\)，\(\vec a\) 的長度為 \(r\)，極角為 \(\varphi\)，有

這就是向量的旋轉公式，還比較好推。

直線

由于可能存在豎直的線，因此直線與線段一般都不使用斜率來存，而是用線上的兩個點來確定線的形態。
同時，由于計算幾何的特殊性，我們是要求直線與線段是有方向性的，其方向的規定與向量類似。
定義一條直線被點 \(s,t\) 確定，方向是 \(s\to t\)。

點 \(B\) 到直線的距離 \(L\) 可以通過常見的面積法得到。由于叉積的幾何意義，我們直接有

其中 \(dis_{u,v}\) 表示兩點間的距離。

直線交點由相似法得到。假設第一條直線的起終點分別為 \(A,B\)，第二條的為 \(C,D\)，設原點為 \(O'\)，交點為 \(O\)，則

其中分子上是 \(\triangle ADC\) 面積的兩倍，分母是四邊形 \(ACBD\) 的面積。這個分數相當于是求出了 \(\frac{AO'}{AB}\) 的值。

一個點 \(A\) 向一條直線作垂線的垂足，實際上就是投影向量，可以通過點積快速得到。
設垂足為 \(O\)，有

線段

基本與直線類似。需要注意的是對于線段而言，一個點到線段的距離與直線是不同的。
如果一個點不在線段“對應”的平面內，那么距離就不是直線，而是到最近的端點的距離。所謂“對應”就是指是否這個點與兩個端點的兩條連線與線段本身的夾角都小于等于九十度。這個東西可以用點積的性質快速判。

多邊形

一般我們存儲多邊形的方式都是通過點集。存儲的相鄰兩個點之間存在連邊，最后一個點與第一個點連邊，所有邊組成一個封閉圖形。
簡潔起見，我們定義第 \(i\) 個點在多邊形上的下一個點為 \(nxt (i)\)，這主要是為了方便最后一個點與第一個點間連邊的表示。
在代碼實現中，由于使用了 vector 存儲，因此點的編號是從 0 開始的。而在算式中，我們為了方便從 1 開始編號。一般將多邊形的點積設為 \(P\)，多邊形上的第 \(i\) 個點表示為 \(P_i\)，一共有 \(n\) 個點。

我們在這里約定一般的多邊形都是點逆時針旋轉的方向存儲的。

多邊形的周長是好算的，直接將兩個相鄰點的距離加起來即可。

多邊形的面積則是一個相對困難的點。由于沒有保證多邊形是一個凸多邊形，因此我們不能采用類似將多邊形分成若干個三角形的方式快速計算。我們希望能有一個簡潔的公式能夠快速計算。事實上，我們有

確實比較簡潔，但證明并不顯然。對于凸多邊形，這個公式還是相對顯然的。討論原點在多邊形內部與外部兩種情況，可以通過畫圖比較明顯的得出來。
而對于凹多邊形，情況就變得復雜了起來。畫個圖可以發現，由于叉積自帶符號的特性，這個式子容斥的意味變得比較深。

其中紅色的是負，藍色的是正。在感性上可能能理解，在理性上感覺不嚴謹。

叉積性質

\(1\over2\) 叉積即帶符號的三角形面積，多邊形每條邊對應的面積為正/負，稱之為 \(+1\) 邊 / \(-1\) 邊。

對于邊 \(AB\)（\(A=A_i,B=A_{(i+1)\bmod n}\)），引射線 \(OA\)，看 \(AB\) 相對于 \(A\) 之后的射線是順時針還是逆時針轉（轉動 \(<180^\circ\)），順時針就是 \(-1\)，逆時針就是 \(+1\)。

射線性質

若 \(O\) 在多邊形上，就施以輕微擾動，讓它去多邊形內/外。以 \(O\) 為起點，引一條射線 \(r\)，要求 \(r\) 不經過多邊形的頂點（自然也不與多邊形的邊重合），稱這樣的射線為“可用射線”，此時 \(r\) 交的邊都“跨過”\(r\)。

由于“叉積性質”，跨過 \(r\) 的邊必定是 \(+1\)/\(-1\) 交替出現。

若點 \(O\) 在多邊形外，由反證法可得若有邊跨過 \(r\)，則第一條跨過 \(r\)（交點離 \(O\) 最近）的邊必定是 \(-1\) 邊（反證法：保證是第一條，還要形成閉合回路，那么必定會繞成順時針，與假設矛盾）。同理，\(O\) 在多邊形內時，第一條跨過 \(r\) 的邊必定是 \(+1\) 邊。

貢獻系數

要求對平面內每個點 \(B\)，若其在多邊形內，則貢獻系數為 \(1\)，否則（多邊形上/多邊形外）貢獻系數為 \(0\)。首先忽略多邊形上的點、不在可用射線上的點。

從 \(O\) 向 \(B\) 引射線，稱射線在 \(B\) 之后的部分為 \(r'\)。跨過 \(r'\) 的邊的 \(\pm1\) 之和即該點的貢獻系數。

\(r'\) 可看作“\(B\) 逃逸的過程”，每有一條邊跨過則必定會轉換“是否在多邊形內”的狀態，否則不轉換。因為最終必定在多邊形外，所以若 \(B\) 在多邊形外則恰有偶數條邊跨過 \(r'\)，否則恰有奇數條邊跨過 \(r'\)。

于是多邊形外的點 \(B\) 貢獻系數為 \(0\)。又由于“射線性質”，無論 \(O\) 在何處，多邊形內的點 \(B\) 對應的 \(r'\) 都會首先被 \(+1\) 邊跨過，故貢獻系數是 \(+1\)。

由于我對容斥的理解很差，因此這里給出一個同學的證明qwq。

凸包

求凸包

最先就是求點集的凸包。個人習慣是用 Andrew 算法。大體思路是將凸包分為兩部分，分為上凸殼與下凸殼，然后拼起來。這個還是比較符合直覺的。

以求上凸包為例。先直接將點按橫坐標排序。然后直接從左到右掃，如果新的點與棧頂的上一個點組成的直線在棧頂的前兩個點形成的向量的逆時針方向，那棧頂就可以彈掉了。
下凸包同理，直接反著做一遍即可。

半平面交

然后是重頭戲，求半平面交。半平面交是指給定一些平面然后要我們求這些平面的交。實際操作的時候一般是通過給出直線來給定半平面然后求這些半平面的交。

半平面是指一條直線及其某一邊的平面。這個時候之前提到的“直線的方向性”的作用就體現于此。具體而言，一般欽定直線左邊的平面被稱作半平面。更準確地說，是直線的“逆時針方向”的平面。求半平面交的過程之所以比較繁復，就體現在定義的隱性，容易把自己轉暈。

然后我們要引入一個叫極角排序的東西。上面在向量的內容中已經粗略介紹了極角的概念。極角是指向量與 \(x\) 正半軸逆時針方向所成的夾角。我們可以通過 c++ 內置的 atan2(y,x) 函數（注意這個函數是縱坐標在前）輕松求出極角的大小。不過需要注意的是，如果極角大小超過了 \(\pi\)，那么這個函數的值會通過負數來表示，也就是變成了從 \(x\) 軸正半軸逆時針旋轉一個大小為負的角度。這個東西顯然也是符合定義的，也比較符合直覺，但是這樣會對之后的極角排序帶來一定的影響。

極角排序，顧名思義，以極角的大小作為關鍵字來為向量排序。這個排序的效果是令所有的線與 \(x\) 軸正半軸的角度依次增大。（注意直線的方向性）求半平面交的過程其實比較樸素，就是將直線一條一條地放在平面上，然后圍成一個圖形。極角排序就可以規整這個“放直線”的順序，像我們畫出凸殼一樣逆時針旋轉著來放這些直線。

也就是說，我們實際上只需要寫一個 sort，然后重定義一下直線排序的關鍵字，然后里面就是判斷直線的極角，小的在前大的在后。
但我們還要考慮一下直線平行（這里的平行指代極角相等而不是斜率相等，還是注意方向性）的情況。我們發現由于是求半平面交，因此如果兩條直線平行，那么更右邊的那條直線就沒用了。因此我們讓更右邊的直線排在后面，排完序后去掉即可。至于如何判斷哪條在右邊，可以先不管，看代碼的時候依照代碼畫一下圖即可。粗略地講就是用叉積判一下逆時針還是順時針。

然后就是“一條一條放直線”的過程。我們在這個過程中實際上干兩件事：

• 算一下當前的直線與前一條直線的交點。（由于我們極角排序了，因此最后形成的圖形的拐點一定是由我們放的相鄰兩條直線交成的）
• 刪去一定沒有貢獻的直線。（這里說的沒有貢獻是指其上的任意一個點都一定不在最后的半平面交內部）

一條線什么時候會完全沒用呢？

如圖。由于極角排序，因此我們插入直線的順序是黑到紫到藍。顯然插入了藍色線以后紫色的線就一定沒有任何作用了，刪掉了也是一樣的。我們發現，還是由于極角排序規定了插入直線的方向性，因此只需要很方便的判定上兩條直線的交點是否在當前直線的右邊，也就是順時針方向即可。如果是，那么就可以直接將上一條線去掉。例如圖中黑色紫色的交點在藍色線右邊，那紫色線就可以直接去掉了。

還有一種情況是我們已經畫了大半了，這個時候插入的線碰到了最前面插入的直線。我們可以用類似于上面的方法判定前面的是否完全沒用。可以之后去看代碼，不再贅述。
上述操作在實現的時候就是用雙端隊列即可。

最后的問題是，我們一直在用后插入的線去判斷之前插入的線是否沒用。最后我們要做的事情就是用之前插入的線去判斷后插入的線是否沒用。

然后就可以愉快的將所有的拐點合成一個多邊形。順帶一提，如果半平面交是一個封閉圖形，那么一定是一個凸殼。

旋轉卡殼

旋轉卡殼求凸包直徑的算法，也就是求凸包內最遠的兩點。

這個算法感性理解很簡單，就是對于每一個點“最遠”的線段，然后全部取 \(\max\) 即可。這個比較有單調性，隨著選的點在轉，線段一起轉即可。
但是嚴謹證明實際上需要很多條件，我其實也扯不清楚。這里推薦一個博客：http://www.rzrgm.cn/Sinktank/p/18308988 ，說的很清楚，圖也很棒。

圓

三點確定一圓

可以用求中垂線的一次函數然后求交點的方法。但是顯然誤差巨大，同時沒有斜率的問題要單獨分討（豎直的線）。

于是我們要拋棄初中方法，擁抱高中數學，掏出圓的方程，然后解出來很大一坨，看起來復雜度是常數，實際慢的離譜。于是我們可以先把算式中的一些重復東西先預處理，重復的浮點數運算是很慢的。

最小圓覆蓋

實際上是個超級大暴力，然后分析復雜度期望是線性的。這種看起來很劣的東西是用類似隨機化的東西保證的期望復雜度。
其思想就是：

• 以每個點為起點開始枚舉，畫個以其為圓心半徑為 0 的圓。
• 如果某個點不在畫的這個圓內，就把圓擴大，直到有三個點已經確定了圓的形態。

顯然這種方法的正確性是對的，因為一旦沒有了圓其就再找一個更大的圓。如果運氣差其將枚舉每一種可能的圓。
于是其用什么保證復雜度呢？我們將點集的順序隨機打亂。可以證明復雜度期望會收斂成為線性的。但反正我個人不是很喜歡這種很不符合直覺的算法。可以去看一下代碼感受其有多暴力。

code

幾乎所有操作都有注釋，可以放心食用。同時定義了一些便于書寫的函數，如 get(s,t) 就是求 \(s\to t\) 的向量。
同時特別提醒一下叉積是用“^”符號替代。

#include<bits/stdc++.h>
using namespace std;

#define ld double
const ld eps=1e-6,pi=acos(-1.0),inf=2e9+7;
mt19937 rnd(time(0));

int dcmp(ld x){return x<-eps?-1:(x>eps);}
ld Abs(ld x){return dcmp(x)*x;}

namespace bs{//基礎的線段向量直線
	struct pt{//點或向量 
		ld x,y;
		void rd(){cin>>x>>y;}
		void put(){cout<<x<<' '<<y<<'\n';}
	};
	bool operator < (const pt a,const pt b){return (a.x!=b.x)?a.x<b.x:a.y<b.y;} //以 x 為第一關鍵字排序，用于求凸包 
	
	pt operator * (const pt a,const ld b){return (pt){a.x*b,a.y*b};}//向量的常數運算 
	pt operator / (const pt a,const ld b){return (pt){a.x/b,a.y/b};}
	
	pt operator + (const pt a,const pt b){return (pt){a.x+b.x,a.y+b.y};}//向量間運算 
	pt operator - (const pt a,const pt b){return (pt){a.x-b.x,a.y-b.y};}
	ld operator * (const pt a,const pt b){return a.x*b.x+a.y*b.y;}//點積 
	ld operator ^ (const pt a,const pt b){return a.x*b.y-a.y*b.x;}//叉積 
	bool operator == (const pt a,const pt b){return (dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0);}//相等 
		
	ld len(pt a){return sqrt(a.x*a.x+a.y*a.y);}//模長
	ld ang(pt a){return atan2(a.y,a.x);}       //與 x 軸的夾角 
	ld dis(pt a,pt b){return len(a-b);}        //兩點間距離 
	int side(pt a,pt b){return dcmp(b^a);}     //判斷向量 a 在 b 的哪邊，1 逆，-1 順，0 共線
	pt get(pt a,pt b){return b-a;}             //得到兩點間的向量，b 為尾 
	pt rotate(pt a,ld theta){                  //向量轉動一個特定角度 
		ld c=cos(theta),s=sin(theta);
		return (pt){a.x*c-a.y*s,a.x*s+a.y*c};
	}
	
	struct seg{//線段 
		pt s,t;
	};
	seg get_seg(pt a,pt b){return (seg){a,b};}
	ld len(seg a){return dis(a.s,a.t);}  //線段長度 
	ld dis(seg a,pt b){//點到線段距離 
		pt s=a.s,t=a.t;
		if(dcmp(get(s,b)*get(s,t))<0||dcmp(get(t,b)*get(t,s))<0)return min(dis(s,b),dis(t,b));
		return (get(s,b)^get(t,b))/len(a);
	} 
	pt close_pt(seg a,pt b){//線段上最近的點 
		pt s=a.s,t=a.t;
		if(dcmp(get(s,b)*get(s,t))<0||dcmp(get(t,b)*get(t,s))<0)return dis(s,b)<dis(t,b)?s:t;
		pt v=get(a.s,a.t); 
		return (v*(((b-a.s)*v)/(v*v)))+a.s; 
	}
	struct line{//直線 
		pt s,t; //半平面統一在逆時針方向（左邊） 
	};
	line get_line(pt a,pt b){return (line){a,b};}
	int side(line a,pt b){return side(b-a.s,a.t-a.s);}  //判斷點 b 在直線 a 的哪邊  注意直線是有方向的 
	pt cross_line(line a,line b){   //直線交點 
		pt A=a.s,B=a.t,C=b.s,D=b.t;
		return A+(get(A,B)*((get(A,C)^get(A,D))/(get(A,B)^get(C,D))));
	}
	pt ver_line(line a,pt b){  //一個點向直線做垂線的垂足 
		pt v=get(a.s,a.t);  //直線的方向向量
		return (v*(((b-a.s)*v)/(v*v)))+a.s; 
	}
	pt sym_line(line a,pt b){  //一個點關于直線的對稱點 
		return (ver_line(a,b)*2)-b;
	}
	bool operator < (const line a,const line b){  //極角排序的關鍵字 注意到這里的排序并非嚴格的按照極角從小到大的順序排列，只是相對位置不變（因為極角為負的會在前） 
		ld a1=ang(get(a.s,a.t)),a2=ang(get(b.s,b.t));//注意在 sort 的關鍵字中不能出現等號 
		if(dcmp(a2-a1)!=0)return dcmp(a2-a1)>0;   //按極角本身大小排序 
		return dcmp(get(a.s,b.t)^get(a.s,a.t))>0; //如果 a 在 b 的逆時針方向，那么 b 就沒有用就排在前面然后在去重的時候去掉 
	}
}
using namespace bs;

namespace poly{//多邊形
	struct polygon{
		vector <pt> pts;//原則上逆時針排序 
		void rd(int n){for(int i=1;i<=n;i++){pt tmp;tmp.rd();pts.push_back(tmp);}}
		void put(){cout<<"poly_put\n";for(pt tmp:pts)tmp.put();}
		
		pt operator [](int x){return pts[x];}
		int nxt(int x){return x<pts.size()-1?x+1:0;}
		
		ld C(){//多邊形周長 
			ld ans=0;
			for(int i=0;i<pts.size();i++)ans+=dis(pts[i],pts[nxt(i)]);
			return ans;
		}
		ld S(){//多邊形面積 
			ld ans=0;
			for(int i=0;i<pts.size();i++)ans+=pts[i]^pts[nxt(i)];
			return Abs(ans)/2;
		}
	};
}
using namespace poly;

namespace tb{//凸包 
	void andrew(polygon &p){  //求點集凸包 
		sort(p.pts.begin(),p.pts.end());
		vector <pt> q(p.pts.size()<<1);
		int top=-1,siz=p.pts.size();q[++top]=p[0];
		for(int i=1;i<siz;i++){
			while(top&&side(get(q[top-1],q[top]),get(q[top-1],p[i]))>=0)top--;
			q[++top]=p[i];
		}
		int now=top;
		for(int i=siz-2;i>=0;i--){
			while(top>now&&side(get(q[top-1],q[top]),get(q[top-1],p[i]))>=0)top--;
			q[++top]=p[i];
		}
		if(siz>1)top--;
		p.pts.clear();for(int i=0;i<=top;i++)p.pts.push_back(q[i]);
	}
	ld kk(polygon p){  //旋轉卡殼求凸包直徑
		ld ans=0;int beg=0,siz=p.pts.size(),tmp=0;
		for(int i=0;i<siz;i++){  //找到對于 0 最遠的線段 
			if(dis(get_seg(p[i],p[p.nxt(i)]),p[0])>tmp)
			beg=i,tmp=dis(get_seg(p[i],p[p.nxt(i)]),p[0]);
		}
		for(int i=0;i<siz;i++){
			while(dis(get_seg(p[beg],p[p.nxt(beg)]),p[i])<dis(get_seg(p[p.nxt(beg)],p[p.nxt(p.nxt(beg))]),p[i]))beg=p.nxt(beg);
			ans=max({ans,dis(p[beg],p[i]),dis(p[p.nxt(beg)],p[i])});
		}
		return ans;
	}
	polygon si(vector <line> Q){  //求半平面交 
		int siz=Q.size(),len=0;sort(Q.begin(),Q.end());//極角排序 
		vector <line> q(siz+1),lst(siz+1);vector <pt> p(siz+1);//p 表示交點坐標 
		for(int i=0;i<siz;i++){   //去重 
			if(i&&dcmp(ang(get(Q[i].s,Q[i].t))-ang(get(Q[i-1].s,Q[i-1].t)))==0)continue;
			q[++len]=Q[i];
		}
		int l=1,r=0;              //雙端隊列 
		for(int i=1;i<=len;i++){
			while(l<r&&dcmp(get(p[r],q[i].t)^get(p[r],q[i].s))>0)--r;
			while(l<r&&dcmp(get(p[l+1],q[i].t)^get(p[l+1],q[i].s))>0)++l;
			lst[++r]=q[i];if(r>l)p[r]=cross_line(lst[r],lst[r-1]);
		}
		while(l<r&&dcmp(get(p[r],lst[l].t)^get(p[r],lst[l].s))>0)--r;  //由于前面只有后插入的來彈出前插入的，于是這里要用之前插入的去彈出后插入的 
//		while(l<r&&dcmp(get(p[l+1],lst[r].t)^get(p[l+1],lst[r].s))>0)++l;
		
		p[r+1]=cross_line(lst[l],lst[r]),++r;  //不知道為什么用 r++ 會錯，感覺判定有點奇怪（實際上是因為加的時候仍然是 p_r 而不是 p_r+1 被賦值） 
		polygon P;for(int i=l+1;i<=r;i++)P.pts.push_back(p[i]);
		return P;
	}
}
using namespace tb;

namespace circle{//圓 
	struct cir{  
		pt o;ld r;
		bool in(pt a){return dcmp(r-dis(a,o))>=0;}
	};
	cir get_cir(pt o,ld r){return (cir){o,r};}  //直接得到圓 
	cir get_cir(pt u,pt v,pt w){  //用圓的解析式三點確定一圓，精度相對高 
	    ld A,B,C,D;
	    ld tmp1=v.x-w.x,tmp2=v.x*w.y,tmp3=w.x*v.y,tmp4=(u.x*u.x+u.y*u.y),tmp5=(v.x*v.x+v.y*v.y),tmp6=(w.x*w.x+w.y*w.y),tmp7=v.y-w.y;
	    A=u.x*(tmp7)-u.y*(tmp1)+tmp2-tmp3;
	    B=tmp4*(-tmp7)+tmp5*(u.y-w.y)+tmp6*(v.y-u.y);
	    C=tmp4*(tmp1)+tmp5*(w.x-u.x)+tmp6*(u.x-v.x);
	    D=tmp4*(tmp3-tmp2)+tmp5*(u.x*w.y-w.x*u.y)+tmp6*(v.x*u.y-u.x*v.y);
	    return get_cir((pt){-B/(2.0*A),-C/(2.0*A)},sqrt((B*B+C*C-4*A*D)/(4*A*A)));
	}
	cir cir_cov(vector <pt> p){  //最小圓覆蓋 
		shuffle(p.begin(),p.end(),rnd);int siz=p.size(); //注意要隨機打亂才能保證復雜度 
		cir c=get_cir((pt){0,0},0);
		for(int i=0;i<siz;i++){
			if(!c.in(p[i])){
				c=get_cir(p[i],0);
				for(int j=0;j<i;j++){
					if(!c.in(p[j])){
						c=get_cir((p[i]+p[j])/2.0,dis(p[i],p[j])/2.0);
						for(int k=0;k<j;k++)if(!c.in(p[k]))c=get_cir(p[i],p[j],p[k]);
					}
				}
			}
		}
		return c;
	}
}
using namespace circle;

const int N=2e6+7;  //記得對答案特判 -0 的情況避免被創飛 

signed main(){
	ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
	
	return 0;
}