# 2020資訊之芽—計算幾何(Computational Geometry)
> [name=peienwu]
> [time=Mon, Aug 23, 2021 12:06 PM]
###### tags: `2020資訊之芽`
###### tags: `2021暑期筆記`
###### tags: `計算幾何`
計算幾何：「以數學方式精確與高效的計算出幾何結構與其屬性的嚴謹應用」。我們可以延伸解決n維的幾何問題，但這篇討論的範圍都是二維平面上的幾何。
## 課程內容
### 座標與向量
* 長度、角度、座標、向量
* 最常以座標、向量表示
* Define x first, y second
* 內積Dot：$A \cdot B = A_xB_x+A_yB_y$
* 外積Cross：$A \times B = A_xB_y-A_yB_x$
* 運算子重載：加減乘除、取純量（$abs()$）
### 有向面積
* 用外積算面積（有正負：有向面積）
* 逆時針為正、順時針為負
* 多邊形面積：任選平面上一點A，將所有點與A連線
* 透過順、逆加總有向面積（p0,p1...pn,其中p0=pn）
* AREA = $\frac{1}{2}\sum_{i=0}^{N-1}\vec{P_i}\times\vec{P_{i+1}}$
### 線段相交
* 線段不平行：$P_1$ 與$P_2$ 會在線段$P_3$、$P_4$異側（方向函數）
* 線段平行：檢驗是否共線、並確認某一個點是否在線段上
### 誤差分析：EPS
* 使用二進位儲存：必產生誤差
* 精度：float $10^{-7}$, double $10^{-16}$, long double $10^{-19}$
* 誤差容忍值$eps$，將 $x\pm eps$視為 $x$
* 重載運算子：==（視為相等）,>,<，加上誤差範圍
* Eps大小：多落在 $10^{-6}$到 $10^{-12}$ 之間
* 誤差：加減法，絕對誤差相加；乘除法，相對誤差相加
* 下界：數字範圍為V內時，$eps$ 至少要VK 乘上資料型態本身誤差
* 上界：題目一般會給
* 避免誤差大法：非到最後關頭，否則都用整數運算！
### 三角函數
* 泰勒展開式逼近，時間並非O(1)
* $atan2(y,x) = \theta$，回傳值域$(-pi,pi]$
* 回傳long double ：使用$atan2l(y,x)$
* 常數大，不建議使用
### IEEE 754
* Sign, Exponent, Mantissa
* 正負號、指數部分、小數部分
### 極角排序
* 給定很多點，依照與某特定點（原點）的角度進行排序
* Sort by cross，依照內積排序
* 題目：平面上n個點，問一條直線最多通過幾個點
### 凸包
* 多邊形：簡單多邊形（邊不相交）、凸多邊形（內角都≤180）、凹多邊形（有內角>180）
* 能包住所有點的凸多邊形
* 凸包求法：Monotone Chain（二維平面）、DC（三維）
* **Monotone Chain**
* 將所有點按照(x,y)排序
* 把下凸包、上凸包「圍」出來
* 合併下凸包、上凸包
* 開一個vector（功能為stack）紀錄當前下半凸包
* 檢查新加入的點會讓哪些點不再是凸包上的點
### 模擬退火(SA)
* 尋找空間中近似最優解
* 一個隨機算法
* 例題：給你平面上 N 個點，請你找出一個點，使得這個點連到這 N 個點的距離總和最短
## 計算幾何函式模板
計算幾何最重要的莫過於座標上的點，實作方式可以用 $std::pair$ 或是自己定義一個類別，將點的資訊以及相關的運算式定義出來。我們總共需要重載點的大於、等於、小於的運算子，以及加法減法、向量外積內積等，同時還有很多功能是可以繼續定義下去，例如向量乘上一個定值，可以繼續加入類別中。
### 點的模板及運算子重載
```cpp=
struct pt{
int x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(abs(x-b.x)<=eps && abs(y-b.y)<=eps)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
pt operator*(int k) {return {x * k, y * k};} //向量乘法（係數積）
pt operator/(int k) {return {x / k, y / k};} //向量除法
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
```
以上是預設點的x座標y座標都是整數的情況，如果要改成使用自定義型別，可以改用樣板(Template)自定義資料型別，根據題目的要求，使用整數或是浮點數進行運算。
```cpp=+
template <typename T>
class pt{
T x,y;
//下方運算子重載與上方相同
}
pt<int> p[N]; //點座標宣告為整數
pt<double> pp[N]; //點的不同資料型別宣告
```
### 方向函數
針對 $\vec a$ 以及 $\vec b$ 外積的結果，可以知道兩者之間相對的方向。如果 $\vec a$ 和 $\vec b$ 共線，則回傳0，$\vec a$ 轉向 $\vec b$ 如果是順時針則回傳1，其餘回傳-1。
```cpp=+
int dir(pt a, pt b, pt o) { //方向函數
int cross = (a - o) ^ (b - o);
if(fabs(cross) <= eps) return 0;
else if(cross > 0) return 1;
else return -1;
}
```
### 點與線段關係
以下函式可以判斷點o是否在 $\overline{AB}$ 上，首先利用外積是否為0判斷是 $\overline{OA}$ 與 $\overline{OB}$ 是否平行；接著以內積判斷是否在線段中，而非線段的兩側（平行的條件下內積只可能是1或是-1）。
```cpp=+
bool onseg(pt a, pt b, pt o){ //o是否在ab線段上
int cross = (a - o) ^ (b - o); //是否平行
int dot = (a - o) * (b - o); //是否在線段中
return (cross == 0)&&(dot <= 0);
}
```
### 線段相交
以下函式為給定四個點$A,B,C,D$，判斷 $\overline{AB}$ 是否相交於 $\overline{CD}$。首先是特例的判斷，線段的其中一端點在另一線段上，利用上方點與線段關係的函式完成這個判斷。
特例判斷完成之後，我們需要用到上方方向函式判斷線段兩端點是否在另一條線段的異側，即以下的關係式：
$$dir(A,B,C)\times dir(A,B,D) < 0$$
我們要檢查兩條線段，其相乘結果必須皆為負數，表示處於線段的異側！最後是平行線的判斷，如果兩線平行且相交，表示兩線共線，這可以在特例時就被判斷出來。
```cpp=+
bool Intersection(pt a, pt b, pt c, pt d){ //線段ab是否與cd相交
if(onseg(a,b,c)||onseg(a,b,d))return true; //點c、d是否洽在線段ab上
if(onseg(c,d,a)||onseg(c,d,b))return true; //點a、b是否洽在線段cd上
if(dir(a,b,c)*dir(a,b,d)==-1 && dir(c,d,a)*dir(c,d,b)==-1)
return true; //對於線段兩端點看另外兩端點必須方向相反
return false;
}
```
### 極角排序
在進行全點對線段共線問題的判斷時，使用極角排序通常會比單純暴力枚舉更快速。極角也就是極座標中每一個跟原點的夾角。如果兩個點位在左半平與右半平面，則先將其判斷出來，如此才能確定起始的角度為何。
如果位在同一個左右半平面，則透過外積的方式比較兩個向量的先後順序。以下程式碼是從座標平面270度的地方開始逆時針掃一圈依序經過的點。
```cpp=+
bool cmp(pt a, pt b){
bool f1 = a < pt{0,0};
bool f2 = b < pt{0,0};
if(f1 != f2)return f1 < f2;
return (a ^ b) > 0;
//逆時針將點進行極角排序，從270度開始逆時針
}
sort(p.begin(),p.end(),cmp); //以id為原點進行極角排序
```
### 凸包函數（使用Monotone Chain）
首先先將所有點依照x座標進行排序，之後用掃描線由左而右的將符合要求的點推入維護的單調容器中，維護下凸包，接著利用reverse()函數將所有點逆序，也就是x座標由大到小讓掃描線由右而左掃過一遍，將上凸包也圍起來。時間複雜度為$O(n\log n)$。
```cpp=+
bool check(pt a, pt b, pt o){
pt aa = a - o;
pt bb = b - o;
return (aa ^ bb) >= 0;
}
vector<pt> convex_hull(){
vector<pt> hull;
sort(p.begin(),p.end(),cmp); //首先對x進行排序
for(auto i : p){ //依序走訪，如果遇到外積<0則不在凸包上
while(hull.size() > 1 && check(i,hull[hull.size()-1],hull[hull.size()-2])){
hull.pop_back();
}
hull.push_back(i); //在凸包hull的每一點都符合外積小於0
}
int down_hull = hull.size();
hull.pop_back(); //x最大的點會在凸包上，不用做兩次先pop一次
reverse(p.begin(),p.end()); //將所有點逆序之後做一次上面的凸包
for(auto i: p){
while(hull.size() > down_hull && check(i,hull[hull.size()-1],hull[hull.size()-2])){
hull.pop_back();
}
hull.push_back(i);
}
return hull; //起點會經過兩次，剛好來算有向面積
}
```
### 旋轉卡尺
旋轉卡尺可以被應用在尋找最遠點對、面積最大三角形等問題。利用兩條平行的線中間夾著凸包，繞一圈的過程中更新需要求的數值。實作上來說就是使用兩個指針，分別指向旋轉卡尺的平行線所在的兩個點，依照旋轉的方向進行增減的動作！
```cpp=+
bool check2(pt a,pt b,pt c,pt d){
int aa = abs((a - c)^(b - c));
int bb = abs((a - d)^(b - d));
return aa < bb;
}
void solve(){
int ans = 0,d = h,sz = hull.size();
rep(i,0,sz-1){
while(check2(hull[i],hull[(i+1)%sz],hull[d],hull[(d+1)%sz]))
d = (d+1)%sz;
ans = max(ans,(hull[i]-hull[d]).dis());
ans = max(ans,(hull[(i+1)%sz]-hull[d]).dis());
}
}
```
## 基本數學知識
計算幾何圍繞著幾個主軸，向量運算、內積、外積，利用它們進行角度、共線與否、距離等等的判斷。之前有稍微接觸過向量，不過內積與外積是第一次碰到的主題。
因為目前絕大部分的討論都是在二維平面上進行，因此以下都是以二維平面為前提所進行的討論！
### 內積（點積）
內積跟 $\cos\theta$ 有關，因此主要可以幫助我們判斷線段是否垂直（等於零）、以及共線時點位於正向或是反向的判斷（等於正負一）。兩個向量 $\vec u$ 以及 $\vec v$ 的內積可以寫成以下關係式：
$$\vec u\cdot \vec v = |\vec u||\vec v|\cos\theta = u_1v_1 + u_2v_2$$
兩個向量做內積的正負號會跟餘弦函數的正負變化相同（如下圖），值域為 $[-|\vec u||\vec v|,|\vec u||\vec v|]$，並且在 $\theta = 0,\pi$ 有最大最小值，$\theta = \frac{\pi}{2},\frac{3\pi}{2}$ 的值為零。
### 外積（叉積）
外積跟 $\sin\theta$ 有關，主要可以判斷兩向量方向關係（順逆時針旋轉）、是否平行、比較角度大小等。外積的應用十分廣泛，找凸包以及旋轉卡尺都會用到外積判斷兩個向量角度關係。兩個向量 $\vec u$ 以及 $\vec v$ 的外積可以寫成以下關係式：
$$\vec u\times \vec v = |\vec u||\vec v|\sin\theta = \begin{bmatrix}u_1&u_2\\v_1&v_2\end{bmatrix}=u_1v_2 - u_2v_1$$
用更簡單的方式理解外積 $\vec u\times \vec v$ ，其正負值可以想像成 $\vec u$ 轉向 $\vec v$ 所經的劣弧順逆時鐘方向。順時針為正、逆時針為負。
關於外積的數值變化，與正弦函數的變化是一樣的（下圖），其值域跟內積一樣，不過最大最小值發生在 $\theta = \frac{\pi}{2},\frac{3\pi}{2}$，並在 $\theta = 0,\pi$ 時兩向量叉積為零。
### 面積
#### 測量師公式（行列式）
這是一個從給定多邊形的座標推得面積的公式，寫成很多個三角形有向面積的總和。以下圖來說，$\triangle FBC$、$\triangle FCD$ 、$\triangle FDE$ 的有向面積皆大於零，而 $\triangle FAB$、$\triangle FEA$ 都會因為有向面積是負的（逆時針旋轉）而被扣除掉，運算的總和即是多邊形 $ABCDE$ 的面積！
一般化的公式，多邊形上總共有 $N$ 個點，令第 $N+1$ 個點為第1個點（為了要繞一圈計算面積），多邊形面積為：
$$AREA = \frac{1}{2}\sum_{i=1}^{N}\vec{P_i}\times\vec{P_{i+1}}$$
#### 三角形外積面積公式
三角形面積有非常多算法，不過利用外積的公式還是第一次聽到。以下是公式推導過程：
先從高中三角函數的三角形公式開始：
$$\begin{split}\triangle ABC &= \frac{1}{2}\overline{AB}\,\overline{AC}\cdot\sin A
\\&=\frac{1}{2}\sqrt{\overline{AB}^2\,\overline{AC}^2\,(1-\cos^2A)}
\\&=\frac{1}{2}\sqrt{\overline{AB}^2\,\overline{AC}^2-(\overline{AB}\cdot\overline{AC})^2}
\\&=\frac{1}{2}\sqrt{(x_1^2+y_1^2)(x_2^2+y_2^2)-(x_1\,x_2+y_1\,y_2)^2}
\\&=\frac{1}{2}\sqrt{(x_1^2\,y_2^2)+(x_2^2\,y_1^2)-2x_1\,y_2\,x_2\,y_1}
\\&=\frac{1}{2}\sqrt{[(x_1\,y_2)-(x_2\,y_1)]^2}
\\&=\frac{1}{2}|\overrightarrow{AB}\times \overrightarrow{AC}|\end{split}$$
將三角形其中一點對另外兩點的向量做外積，除以2即為三角形面積。這個公式會在旋轉卡尺的地方使用到！
#### 平行四邊形面積
根據上面三角形面積公式的推導，可以相對應得知道兩向量所夾平行四邊形面積公式：
$$AREA = |\overrightarrow{AB}\times \overrightarrow{AC}|$$
## 上機作業
- [x] 向量加法
- [x] 等長線段對
- [x] 向左轉向右轉
- [x] 線段相交
- [x] 最小凸多邊形
- [x] TIOJ 1178 Convex Hull
- [x] 來吧，遊戲開始了
- [x] 遊戲：最終回
- [x] TIOJ 1205 直角三角形
- [x] TIOJ 1105 H.PS3
- [x] ZJ b288: 夏季大三角
- [x] TIOJ 1500 Clean up on aisle 3
- [x] ZJ a871: Museum Area
- [x] TIOJ 1280 領土 (Territory)
- [x] TIOJ 1678 剪多邊形（molding）
- [x] ZJ d269: 11579 - Triangle Trouble
### 向量加法
[題目連結](https://neoj.sprout.tw/problem/398/)
[Submission](https://neoj.sprout.tw/challenge/178462/)
> 題目敘述：
> 給你n個數字（0≤i<1，小數點精度到末九位），想知道到底有多少組 $(i,j,k)$ 滿足 $v_i+v_j=v_k$，其中 $i,j,k$ 可以重複。
這題其實跟計算幾何沒什麼關係，直接用unordered_map去做（有點像two sum，不過下面的code好像也不用開到multi），簡單！不過我在浮點數的地方吃了一些WA，最後算了直接改用字串處理這個惱人的東西XD
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define int long long
#define double long double
#define rep(i,l,r) for(int i=l;i<=r;i++)
#define rep2(i,l,r) for(int i=l;i<r;i++)
#define pii pair<int,int>
#define x first
#define y second
#define eps (1e-9)
#define INF 1e10
#define N 2001
#define ll long long
#define ld long double
#define int long long
using namespace std;
int n;
signed main(){
cin>>n;
vector<int> vec(n);
rep2(i,0,n){
string s;cin>>s;
int num = 0,times = 1000000000;
for(int i=2;i<=10;i++){
num += (s[i]-'0')*times;
times/=10;
}
vec[i] = num;
}
unordered_multimap<int,int> mp;
for(int i=0;i<n;i++){
mp.insert({vec[i],i});
}
int ans = 0;
for(int i=0;i<n;i++){
for(int j=0;j<n;j++){
ans += mp.count(vec[i] + vec[j]);
}
}
cout<<ans<<endl;
}
```
### 等長線段對
[題目連結](https://neoj.sprout.tw/problem/399/)
[Submission](https://neoj.sprout.tw/challenge/178471/)
> 題目敘述：
> 給定平面上很多個點，求出有幾對線段等長（輸入有重複的點）。
既然n≤500，那就直接枚舉吧，沒啥特別難度。
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a,i<=b;i++)
#define rep2(i,a,b) for(int i=a;i<b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 1003
#define eps 1e-9
#define x first
#define y second
using namespace std;
int n;
pii p[N];
int dist(pii a,pii b){
int x = a.x-b.x,y = a.y-b.y;
return x*x+y*y;
}
signed main(){
Orz;
cin>>n;
rep2(i,0,n)cin>>p[i].x>>p[i].y;
map<int,int>mp;
for(int i=0;i<n;i++){
for(int j=i+1;j<n;j++){
int dis = dist(p[i],p[j]);
mp[dis]+=1;
}
}
int ans = 0;
for(auto i:mp){
int cnt = i.second;
ans+=((cnt*(cnt-1))/2);
}
cout<<ans<<endl;
}
```
### 向左轉向右轉
[題目連結](https://neoj.sprout.tw/problem/400/)
[Submission](https://neoj.sprout.tw/challenge/178524/)
> 題目敘述
> 給你平面上n個點，依序走訪每一個點，試問走訪過程中共執行幾次的左轉、右轉以及迴轉。
很特別，計算幾何讓電腦可以處理平常我們所看到的平面圖形，可以利用向量內積、外積等方式判斷方向。這一題最重要的就是**方向函數**。傳入3個點$(A,B,O)$，方向函數會會回傳$\stackrel\longrightarrow{OA}\times \stackrel\longrightarrow{OB}$ 的正負數值。
下圖為外積$\stackrel\longrightarrow{OA}\times \stackrel\longrightarrow{OB}$ 的結果，當 $\sin\theta$的結果為負，也就是下圖的情況，從B走到A就需要往左邊走；反之亦然。
至於如何判斷當兩個向量的方向呈現一直線時，也就是外積回傳的值為0時（$\sin\theta = 0$），應該是同向還是異向呢？這時候就需要搭配向量內積（這我想了很久），因為內積公式是$A\cdot B = |A||B|\cos\theta$，將兩個向量內積之後就可以很明確的判斷到底是朝原本的方向走，還是反方向的行走！
:::info
**內積、外積公式**
有一點數學，不過蠻有趣的。可以利用$\sin$與$\cos$達到計算角度的目的，利用兩者不同的值域，互相搭配，就可以更輕鬆的進行判斷！注意到外積的正負就代表著A到B是順時針或是逆時鐘。
$$A\cdot B = |A||B|\cos\theta = A_xB_x+A_yB_y\\A\times B = |A||B|\sin\theta = A_xB_y-A_yB_x$$
<br>
**方向函數**
當我們要判斷方向的時候，會利用正弦函數，逆時針正、順時針為負進行判斷！
```cpp=
int dir(pt a, pt b, pt o) {
int cross = (a - o) ^ (b - o);
if(fabs(cross) <= eps) return 0;
else if(cross > 0) return 1;
else return -1;
}
```
注意到此時在判斷是否為平行的時候（cross==0），使用到$fabs()$這個函數，目的是為了避免誤差而導致判斷錯誤，因此需要進行誤差的處理（其實不用也沒差啦，只是這樣嚴謹一點）
:::
以下是AC Code：
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 1003
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
int x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(x-b.x<=eps && y-b.y<=eps)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
vector<pt> a;
int dir(pt a, pt b, pt o) {
int cross = (a - o) ^ (b - o);
if(fabs(cross) <= eps) return 0;
else if(cross > 0) return 1;
else return -1;
}
int n,t;
signed main(){
Orz;
cin>>n;
a.resize(n+2);
rep(i,1,n)cin>>a[i].x>>a[i].y;
int right = 0,left = 0,turn = 0;
pt pre = a[1],from = a[2];
for(int i=3;i<=n;i++){
int ori = dir(a[i],from,pre);
if(ori == 1)right+=1;
else if(ori == -1)left+=1;
else if(ori == 0 && ((a[i]-from)*(from-pre))<0)turn+=1;
pre = from;from = a[i];
}
cout<<left<<" "<<right<<" "<<turn<<endl;
}
```
### 線段相交
[題目連結](https://neoj.sprout.tw/problem/401/)
[Submission](https://neoj.sprout.tw/challenge/178537/)
線段相交 = ~~線段香蕉~~，自動選字永遠都是香蕉，有點煩XDD
如何判斷兩線段是否相交？首先需要一個函數可以判斷點是否在一個線段上，如此一來就可以判斷端點在另一條線段上的特殊情況。以下程式碼為判斷點$P_o$ 是否在 $\overline{P_aP_b}$ 上。利用向量外積可以判斷兩線段是否平行，而使用內積公式可以判斷$P_o$是否在線段中，而非線段的兩側！
```cpp=
bool onseg(pt a, pt b, pt o){ //o是否在ab線段上
int cross = (a - o) ^ (b - o); //是否平行
int dot = (a - o) * (b - o); //是否在線段中
return (cross == 0)&&(dot <= 0);
}
```
說明：由點$P_o$指向a和b的向量必須呈現180度角（也就是異向），才可確保在ab線段中（跟a,b重合也算是跟ab線段相交）。
接下來是主要的部分，首先先確認4個端點是否恰好在另外一條線段上，判斷完之後就是處理一般相交的情況。若線段 $\overline{P_1P_2}$ 與 $\overline{P_3P_4}$ 相交，則點 $P_1$ 與點 $P_2$ 會在線段$\overline{P_3P_4}$ 的異側。用方向函數表示：$dir(a,b,c)\times dir(a,b,d)<0$。確認完兩個線段之後即完成線段相交的判斷！
```cpp=
bool Intersection(pt a, pt b, pt c, pt d){ //線段ab是否與cd相交
if(onseg(a,b,c)||onseg(a,b,d))return true; //點c、d是否洽在線段ab上
if(onseg(c,d,a)||onseg(c,d,b))return true; //點a、b是否洽在線段cd上
if(dir(a,b,c)*dir(a,b,d)==-1 && dir(c,d,a)*dir(c,d,b)==-1)
return true; //對於線段兩端點看另外兩端點必須方向相反
return false;
}
```
由下圖可以得到上面的結論，當兩線段相交時，方向函數得到的值（用外積，也就是下圖 $\theta_1$ 以及 $\theta_2$）的方向），會呈現一正一負，從兩個相反的方向看同一條線段得出來的結論！
AC Code:
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 1003
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
int x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(x-b.x<=eps && y-b.y<=eps)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
vector<pt> point;
int dir(pt a, pt b, pt o) { //方向函數
int cross = (a - o) ^ (b - o);
if(fabs(cross) <= eps) return 0;
else if(cross > 0) return 1;
else return -1;
}
bool onseg(pt a, pt b, pt o){ //o是否在ab線段上
int cross = (a - o) ^ (b - o); //是否平行
int dot = (a - o) * (b - o); //是否在線段中
return (cross == 0)&&(dot <= 0);
}
bool Intersection(pt a, pt b, pt c, pt d){ //線段ab是否與cd相交
if(onseg(a,b,c)||onseg(a,b,d))return true; //點c、d是否洽在線段ab上
if(onseg(c,d,a)||onseg(c,d,b))return true; //點a、b是否洽在線段cd上
if(dir(a,b,c)*dir(a,b,d)==-1 && dir(c,d,a)*dir(c,d,b)==-1)
return true; //對於線段兩端點看另外兩端點必須方向相反
return false;
}
int n,t;
void solve(){
point.assign(4,{0,0});
rep(i,0,3)cin>>point[i].x>>point[i].y;
if(Intersection(point[0],point[1],point[2],point[3])){
cout<<"Yes"<<endl;
}
else cout<<"No"<<endl;
}
signed main(){
Orz;
int t;cin>>t;
while(t--){
solve();
}
}
```
### TIOJ 1178 Convex Hull
[題目連結](https://tioj.ck.tp.edu.tw/problems/1178)
[Submission](https://tioj.ck.tp.edu.tw/submissions/262532)
> 題目敘述
給定n個二維平面的點，找出位在凸包上的所有點的個數
最小凸多邊形 = 凸包，要找出能包住所有點的最小凸多邊形，簡稱凸包。聽說最好寫的凸包演算法是：Andrew's Monotone Chain，翻成中文叫做Andrew's 單調鍊？有一點單調+鍊的味道。下圖是我用照片合成起來的GIF，大致模擬出使用Andrew's Monotone Chain 找凸包的方法。
:::info
**Andrew's Monotone Chain**
這個演算法的時間複雜度是 $O(n\log n)$，空間複雜度 $O(n)$，資料說它可以解決了凸包有重疊的點、共線的點、退化成線段和點的情況。它的名字叫做「單調鍊」，要維護一個有點像單調隊列的東西，對於在容器中第 $i$ 個位置的點都滿足 $\stackrel\longrightarrow{P_i P_{i+1}}\times \stackrel\longrightarrow{P_{i+1} P_{i+2}} > 0$ ，如果有點做外積後的結果小於等於0，則它會被pop掉（這是依照上圖逆時針完成凸包的描述，如果方向相反則會變號）。
<br>
以下是此演算法的執行步驟：
1. 先把所有的點按照 $(x,y)$ 排序
2. 將下凸包圍出來，有點像維護單調隊列，對所有新加入的點i計算點i-2、i-1與i之間的外積，如果不符合情況代表圍不到新加入的點，需要將舊的點pop出來
3. 將原本已經排序好的點逆序
3. 再把上凸包由x座標大到小圍出來，將上下合併就是凸包了（必須注意起終點被push的次數問題）！
:::
一般會用一個vector儲存在凸包上面的點（不包含在邊上的點，只有位於轉折點的點），在頭尾的部分（x座標最大與最小）需要特別處理，讓每一個點最多近到vector一次。
:::warning
**實作細節**
以下是確認是否需要將vector中元素pop出來的關鍵，對向量$\stackrel\longrightarrow{OA}\times \stackrel\longrightarrow{OB}$ 做外積的結果，必須排除外積結果為0的情況，如果將0也納入，會造成一個點被push進去很多次，在數量和計算上出現問題。
等於的情況代表這些點會共線，因此如果題目要求要輸出所有凸包上的點，就必須在等號發生的情況也一併紀錄下去。對一個極端的多邊形來說（所有點都線線時），必須要處理去重的情況（用set完成），才不會有一個點被重複push進去的問題。
```cpp=
bool check(pt a,pt b,pt o){
int cross = (a - o)^(b - o);
return cross >= 0; //這裡很關鍵，別吃WA
}
```
除此之外，上凸包在範圍限制上是需要注意的。假設x座標最大的點i，當在圍上凸包的過程中i是不可以被pop出去的，因此vector的大小必須大於下凸包的大小。
凸包使用第i-1跟第i個點的向量去看第i到第i+1個點的向量，決定一個點要不要被推入vector中。當我們逆序從x座標最大的點往前看時，要確保每一輪結束之後在i點後都必須要有至少一個點，設定hull.size() > down_hull的原因是防止在下凸包的點被圍上凸包的過程更新到。
```cpp=+
int down_hull = hull.size(); //圍上凸包的程式碼片段
for(auto i: p){
while(hull.size() > down_hull
&& check(i,hull[hull.size()-1],hull[hull.size()-2])){
hull.pop_back();
}
hull.push_back(i);
}
```
:::
以下是AC Code:
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 100001
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
int x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(x-b.x<=eps && y-b.y<=eps)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
bool cmp(pt a, pt b){
if(a.x == b.x)return a.y < b.y;
return a.x < b.x;
}
vector<pt> p;
bool check(pt a,pt b,pt o){
int cross = (a - o)^(b - o);
return cross >= 0; //這裡很關鍵，別吃WA
}
int n,t;
vector<pt> convex_hull(){
vector<pt> hull;
sort(p.begin(),p.end(),cmp); //首先對x進行排序
for(auto i : p){ //依序走訪，如果遇到外積<0則不在凸包上
while(hull.size()>=2 && check(i,hull[hull.size()-1],hull[hull.size()-2])){
hull.pop_back();
}
hull.push_back(i); //在凸包hull的每一點都符合外積小於0
}
int down_hull = hull.size();
hull.pop_back(); //x最大的點會在凸包上，不用做兩次先pop一次
reverse(p.begin(),p.end()); //將所有點逆序之後做一次上面的凸包
for(auto i: p){
while(hull.size() > down_hull && check(i,hull[hull.size()-1],hull[hull.size()-2])){
hull.pop_back();
}
hull.push_back(i);
}
return hull; //起點會經過兩次，剛好來算有向面積
}
signed main(){
Orz;
cin>>n;
p.assign(n,{0,0});
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
vector<pt> hull = convex_hull();
cout<<hull.size()-1<<endl;
}
```
### 最小凸多邊形
[題目連結](https://neoj.sprout.tw/problem/402/)
[Submission](https://neoj.sprout.tw/challenge/178589/)
> 題目敘述
找出二維平面上n個點的凸包所圍出來的面積為何？
跟上一題類似，在找到全部在凸包上面的點後，就可以利用**有向面積**把凸包面積算出來，有一個公式可以計算多邊形面積，利用外積得到正負值，轉一圈後得到面積！對於多邊形的頂點 $P_0,P_1,...,P_{n-1},P_n=P_0$ 的面積如下：
$$Area = \frac{1}{2}\sum_{i=0}^{n-1}\stackrel\longrightarrow{P_i}\times \stackrel\longrightarrow{P_{i+1}}$$
其中最後一個點會回到起點，形成一個封閉的迴路。
以下是AC Code：
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 100001
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
int x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(x-b.x<=eps && y-b.y<=eps)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
bool cmp(pt a, pt b){
if(a.x == b.x)return a.y < b.y;
return a.x < b.x;
}
vector<pt> p;
bool check(pt a,pt b,pt o){
int cross = (a - o)^(b - o);
return cross >= 0;
}
int n,t;
vector<pt> convex_hull(){
vector<pt> hull;
sort(p.begin(),p.end(),cmp); //首先對x進行排序
for(auto i : p){ //依序走訪，如果遇到外積<0則不在凸包上
while(hull.size() > 1 && check(i,hull[hull.size()-1],hull[hull.size()-2])){
hull.pop_back();
}
hull.push_back(i); //在凸包hull的每一點都符合外積小於0
}
int down_hull = hull.size();
hull.pop_back(); //x最大的點會在凸包上，不用做兩次先pop一次
reverse(p.begin(),p.end()); //將所有點逆序之後做一次上面的凸包
for(auto i: p){
while(hull.size() > down_hull && check(i,hull[hull.size()-1],hull[hull.size()-2])){
hull.pop_back();
}
hull.push_back(i);
}
return hull; //起點會經過兩次，剛好來算有向面積
}
signed main(){
Orz;
cin>>t;
while(t--){
cin>>n;
p.assign(n,{0,0});
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
vector<pt> hull = convex_hull();
int area = 0,len = hull.size();
for(int i=0;i<len-1;i++)area += (hull[i]^hull[i+1]);
cout<<fixed<<setprecision(1)<<((ld)area/2)<<endl;
}
}
```
### 來吧，遊戲開始了。
[題目連結](https://neoj.sprout.tw/problem/790/)
[Submission](https://neoj.sprout.tw/challenge/178691/)
[GGB模擬](https://www.geogebra.org/graphing/h4fxdquw)
> 題目敘述
給你二維平面上n個點(n≤2400)，每一個點座標皆不相同，求出總共可以圍出多少個三角形？
這是NEOJ上的加分題，好像是一個題組吧，反正總共有三題，這是第一題。如果 $O(n^3)$ 的枚舉，複雜度會爆炸（量級約$10^{10}$），根據電神的說法，這一題要用極角排序以及雙指標找到共線，接著就可以利用排列組合把因為共線而不能形成三角形的組合扣掉，就是答案了。
這一題的核心概念是找共線，具體來說的作法是枚舉每一個點的同時，以它為原點對其他的點進行排序，如果遇到有相同的極角座標表示這些點共線，同時利用陣列cnt[x]統計共線點數為x的線段總共有幾條。
以下的GIF就是大致上程式執行的樣子。因為一條長度為x的線段會因為枚舉x次的關係，在最後扣掉的情況會重複x次因此需要除掉。
:::info
**共線與三角形**
一般情況下（任三點不共線），總共可以形成 $C^n_3$ 個三角形，如果有一條m個點共線的情況下（其他點不共線），則可以形成的三角形數量就必須扣除共線限制的情況，變成 $C^n_3-C^m_3$ 個三角形。
:::
時間複雜度為：枚舉每一個點 $O(n)$，極角排序 $O(n\log n)$，總時間複雜度 $O(n^2\log n)$
:::warning
**實作小細節**
**1. 維護共線連續區間**
我們要想辦法讓有共線的點們所在位置是一個連續的位置。三個點共線可能為在對角線的象限中，也就是點差了180度，如此一來就沒辦法讓共線的點為在連續的區間。為了達到這個目的，我們將所有位於下半平面的點都移到上半平面（在上半平面找到有相同 $\tan\theta$ 值的位置），接著就能利用雙指針找極角座標排序後有相同極角的區間之最大值！
**2. 特例判斷**
如果有一點y座標為0但x座標為負，要將其移到x軸正向的地方，不能把這種情況涵蓋為一般情況，否則原本在x軸正向的點會被移到x軸負向，沒有達到預期的效果。
:::
以下是AC Code：
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 100001
#define all(x) x.begin(),x.end()
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
int x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(x-b.x == 0 && y-b.y == 0)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
vector<pt> p;
vector<int> cnt;
int n,ans = 0;
bool cmp(pt a, pt b){
bool f1 = a < pt{0,0};
bool f2 = b < pt{0,0};
if(f1 != f2)return f1 < f2;
return (a ^ b) > 0;
//逆時針將點進行極角排序，從270度開始逆時針
}
//用cnt[i]統計區間長度為i的線段數量
void solve(pt id){
vector<pt> pp;
for(auto i : p){ //以id為原點
pt cur = i-id;
if(cur == pt{0,0})continue;
if(cur.y < 0){cur.x = -cur.x;cur.y = -cur.y;}
if(cur.x < 0 && cur.y==0){cur.x = -cur.x;}
pp.push_back(cur);
}
sort(all(pp),cmp); //將id當作原點進行排序
int p1 = 0,p2 = 0,len = pp.size(); //雙指針找共線區間
while(p1 < n-1){ //最大化區間
while(p2+1 < len && (pp[p1]^pp[p2+1]) == 0)p2++;
cnt[p2-p1+2]+=1;
p1 = p2+1;
}
}
signed main(){
Orz;
cin>>n;
p.assign(n,{0,0});
cnt.resize(n+1,0);
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
rep(i,0,n-1)solve(p[i]);
int ans = (n*(n-1)*(n-2))/6;
rep(i,3,n)ans-=(cnt[i]*(i-1)*(i-2))/6;
cout<<ans<<endl;
}
```
### 遊戲：最終回
[題目連結](https://neoj.sprout.tw/problem/792/)
[Submission](https://neoj.sprout.tw/challenge/178786/)
> 題目敘述
共有n個二維平面上的格子點，這些點會形成簡單多邊形。試求或在簡單多邊形內部的格線總長（包括垂直與水平格線）。
這邊有一個不嚴謹的推導方式，不過他是正確的。令多邊形內部格線長度為S，多邊形的邊落在的格線長度為T，多邊形面積T，則有以下關係式：
$$S = 2A-\frac{T}{2}$$
詳細的公式推導可以可以參閱下圖，平行四邊形（斜線部分）內部**垂直**的格線長度為： 大矩形 $(x_1+x_2)(y_1+y_2)$ 扣掉左右上下共四個三角形兩兩拼成一個矩形 $x_1y_1$ 以及 $x_2y_2$，還有左上右下兩個正方形 $2x_2y_1$，整理之後會發現其實跟面積是一樣的。對於垂直部分也是類似的情況。
好像隱約發現到面積與格線長度有十分密切的關係，算出面積，把在格線上的邊進行特判扣掉，就可以得到格線長度。
這一題我想了很久，一直看不出來關係式到底長怎樣，直到大神提點才發現原來有這樣的關係，我反應好遲鈍:cry:
以下是AC Code：
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 100001
#define all(x) x.begin(),x.end()
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
int x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(x-b.x == 0 && y-b.y == 0)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
vector<pt> p;
vector<int> cnt;
int n,edge,ans;
int solve(){
int area = 0;
rep(i,0,n-1){
area += (p[i]^p[i+1]);
if(p[i].y == p[i+1].y)edge += abs(p[i].x-p[i+1].x);
if(p[i].x == p[i+1].x)edge += abs(p[i].y-p[i+1].y);
}
area = abs(area);
return area;
}
signed main(){
Orz;
cout<<fixed<<setprecision(1);
while(cin>>n){
p.assign(n+1,{0,0});
edge = 0;
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
p[n] = p[0];
ans = solve();
cout<<ans-((ld)edge/2)<<endl;
}
}
```
### TIOJ 1205 直角三角形
[題目連結](https://tioj.ck.tp.edu.tw/problems/1205)
[Submission](https://tioj.ck.tp.edu.tw/submissions/262842)
> 題目敘述
給你N（N≤1500）個座標平面上的點，請問總共可形成多少個直角三角形呢？
從極角排序後的第一個點開始逆時針進行雙指針的枚舉。這邊使用到一個很特別的手法，對於共線的情況我們先透過預處理的方式將共線的點合併起來，並用cnt[x]陣列紀錄第x個點是由幾個點所合併起來的，如此一來，在進行計算的時候就不會有共線要分別處理的問題（不需擔心是不是可以跟之前的點形成直角三角形，因為相同斜率的點已經被合併剩下一個），直接將數量相乘就可以知道直角三角形的數量！
時間複雜度：枚舉所有點 $O(n)\times$ 進行極角排序$O(n\log n)$ 以及雙指標$O(n)$，總時間複雜度為 $O(n^2\log n)$。
:::success
**實作小細節**
雙指針進行枚舉的過程中，很有可能會指標指向的索引值會超出範圍。解決的方法有兩種：
1. 超出了即代表繞了一圈，只需要對索引值取餘數即可。
2. 除了取餘數的方法之外，其實也可以直接在點集後面將所有點再推入一次，讓角度從360延伸成720度，就不會有超出範圍的問題！
:::
以下是AC Code：
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 100001
#define all(x) x.begin(),x.end()
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
int x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(x-b.x == 0 && y-b.y == 0)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
vector<pt> p,temp,pp;
vector<int> cnt;
int n,ans = 0;
bool cmp(pt a, pt b){
bool f1 = a < pt{0,0};
bool f2 = b < pt{0,0};
if(f1 != f2)return f1 < f2;
return (a ^ b) > 0;
//逆時針將點進行極角排序，從270度開始逆時針
}
//O(n)枚舉每個點當直角情況
int solve(pt id){
pp.clear();cnt.clear();temp.clear();
for(pt i : p){
pt cur = i - id;
if(cur == pt{0,0})continue;
temp.push_back(cur);
}
sort(all(temp),cmp); //以id為原點進行極角排序
pp.push_back(temp[0]); //pp每一角度只存至多一個點
cnt.push_back(1); //考慮每個點共線情況
int len = temp.size();
rep(i,1,len-1){
int cross = temp[i]^temp[i-1],dot = temp[i]*temp[i-1];
if(cross == 0 && dot >= 0)cnt[cnt.size()-1] += 1; //共線數量+=1
else {pp.push_back(temp[i]);cnt.push_back(1);} //非共線設定數量為1
}
len = pp.size(); //考慮橫跨一周的情況
rep(i,0,len-1){ //雙指針i,p1可能會超過一圈
pp.push_back(pp[i]); //將點再繞一圈
cnt.push_back(cnt[i]);
}
int ans = 0,p1 = 0;
rep(i, 0, len-1){
while(p1 < i+len && (pp[i]^pp[p1]) >= 0 && (pp[i]*pp[p1]) > 0)p1 += 1;
//夾銳角的情況要p1+=1
if((pp[i]^pp[p1]) > 0 && (pp[i]*pp[p1]) == 0)ans += cnt[i]*cnt[p1];
//正向的直角三角形，若共線則兩者數量相乘
}
return ans;
}
signed main(){
Orz;
while(cin>>n){
if(n == 0)break;
p.assign(n,{0,0});
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
int ans = 0;
rep(i,0,n-1){
ans += solve(p[i]);
}
cout<<ans<<endl;
}
}
```
### TIOJ 1105 H.PS3
[題目連結](https://tioj.ck.tp.edu.tw/problems/1105)
[Submission $O(n^2)$](https://tioj.ck.tp.edu.tw/submissions/262930)
[Submission $O(n\log n)$](https://tioj.ck.tp.edu.tw/submissions/262947)
> 題目敘述
給你平面上N個點（N≤3000），請求出最遠點對的索引值（小的在前、大的在後）
我做了一份[最近點對：不同複雜度之解決方式](/hVplrqxCRdGiMkn2lwXGAA)的筆記，共有四種方法可以解決那個問題，這一題要求的是最遠點對，作法與最近點對其實差蠻遠的。由上幾題知道凸包的求法，因為凸包是可以圍住所有點的多邊形，因此最遠點對也應該在凸包上，而且所在的位置會為在凸包的兩側上（如果不落在凸包上，一定可以把點向兩側延伸到凸包上，且移動過後的點對距離一定比原始的點對距離大）。
找完凸包之後，可以用旋轉卡尺的方式尋找最遠點對。想像兩條平行線中間夾著凸包，逆時鐘旋轉繞行凸包一圈，過程不斷更新最遠點對的距離。在實作上兩條平行線可以被想像成由 $P_i$ 指向 $P_{i+1}$ 的向量，透過外積三角形面積公式決定卡尺該如何移動。
$$AREA = |\overrightarrow{AB}\times \overrightarrow{AC}|$$
以下圖為例，我們要找 $\overline{HM}$ 為底可以形成的最大三角形面積的頂點，因為在同底的情況下面積就代表點與邊的垂直距離，最大的垂直距離意味著這條底邊可以垂直延伸的最遠距離。因為凸包必定是凸多邊形，因此三角形的面積會呈現單峰函數，因此只需要從下一個三角形面積的大小，決定雙指針中比較快的指標的移動情況。
如果仔細來看，以下圖為例，當前較快的指標指向的位置是 $D$ 點，考慮一條與與 $\overline{HM}$ 平行的直線，若下一個點 $J$ 在平行線段的另外一側，則將指標移往 $J$ 點。可能會有一個疑問，如果比較下圖的線段長度，會發現到 $\overline{DH}$ 的長度比經過 $J$ 點的兩條線段都還要長，那為何還要更新至 $J$ 點？舉這個例子不太好，不過可以想像當旋轉卡尺轉到以 $\overline{FH}$ 為底的時候，會將最遠點對的距離更新成 $\overline{HD}$ 的長度。如果今天 $H$ 的左側又多加了一個新點 $P$，則最遠點對會變成 $\overline{PD}$ 的距離。
簡單來說，最遠點對一定會發生對角的凸包點上面，即使現在以 $\overline{HM}$ 為底最遠點並非 $J$ 而是 $D$ ，但在旋轉卡尺旋轉到 $\overline{FH}$ 時就能將距離更新成 $\overline{HD}$ 的距離。
:::success
**實作小細節**
這一題有點麻煩，因為他要輸出的是最遠點對的索引值，而不是最遠點對之間的距離。在尋找凸包的過程中，會對所有點進行排序，因此原有的索引值順序會被打亂，需要在一開始輸入的時後就好好維護每一個座標的索引值。
:::
以下是AC Code：
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 100001
#define all(x) x.begin(),x.end()
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
int x,y,ind;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(x-b.x == 0 && y-b.y == 0)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
int dis() {return x*x + y*y;}
};
vector<pt> p,hull;
pt pt_ans;
int n,h;
bool check(pt a, pt b, pt o){
pt aa = a - o;
pt bb = b - o;
return (aa ^ bb) >= 0;
}
bool check2(pt a,pt b,pt c,pt d){
int aa = abs((a - c)^(b - c));
int bb = abs((a - d)^(b - d));
return aa < bb;
}
bool cmp(pt a, pt b){
if(a == b)return a.ind < b.ind;
if(a.x == b.x)return a.y < b.y;
return a.x < b.x;
}
void convex_hull(){
stable_sort(all(p),cmp);
rep(i,0,n-2)if(p[i] == p[i+1])p[i+1].ind = p[i].ind;
hull.clear();
for(auto i : p){
while(hull.size() > 1 && check(i,hull[hull.size()-1],hull[hull.size()-2]))
hull.pop_back();
hull.push_back(i);
}
int sz = hull.size();
h = hull.size()-1;
hull.pop_back();
reverse(all(p));
for(auto i : p){
while(hull.size() > sz && check(i,hull[hull.size()-1],hull[hull.size()-2]))
hull.pop_back();
hull.push_back(i);
}
hull.pop_back();
}
void solve(){
int ans = 0,d = h,sz = hull.size();
rep(i,0,sz-1){
while(check2(hull[i],hull[(i+1)%sz],hull[d],hull[(d+1)%sz]))
d = (d+1)%sz;
if(ans < (hull[i]-hull[d]).dis()){
ans = (hull[i]-hull[d]).dis();
int a = hull[i].ind,b = hull[d].ind;if(a > b)swap(a,b);
pt_ans = {a,b};
}
else if(ans == (hull[i]-hull[d]).dis()){
int a = hull[i].ind,b = hull[d].ind;if(a > b)swap(a,b);
if(pt_ans > (pt){a,b})pt_ans = {a,b};
}
if(ans < (hull[(i+1)%sz]-hull[d]).dis()){
ans = (hull[(i+1)%sz]-hull[d]).dis();
int a = hull[(i+1)%sz].ind,b = hull[d].ind;if(a > b)swap(a,b);
pt_ans = {a,b};
}
else if(ans == (hull[(i+1)%sz]-hull[d]).dis()){
int a = hull[(i+1)%sz].ind,b = hull[d].ind;if(a > b)swap(a,b);
if(pt_ans > (pt){a,b})pt_ans = {a,b};
}
}
}
signed main(){
Orz;
while(cin>>n){
if(n == 0)break;
pt_ans = (pt){0,0};
p.resize(n,{0,0});
rep(i,0,n-1){
cin>>p[i].x>>p[i].y;
p[i].ind = i;
}
convex_hull();
solve();
cout<<pt_ans.x<<" "<<pt_ans.y<<endl;
}
}
/*
5
9 1
1 5
1 2
9 9
5 1
*/
```
### ZJ b288: 夏季大三角
[題目連結](https://zerojudge.tw/ShowProblem?problemid=b288)
[解題報告](https://zerojudge.tw/ShowThread?postid=26741&reply=0)
> 題目敘述
請輸出在N個二維平面的座標，挑選3顆出來成組成三角形的最大面積
比較一下兩個複雜度的作法，第一個是使用 $O(n^3)$ 枚舉所有的點並計算面積，所需要的時間是0.4sec，而且需要特別注意不能使用到海龍公式計算面積，否則有很大的機會會超時。
以下作法是先進行 $O(n\log n)$ 找尋凸包，因為面積最大的三角形必定三個點都在凸包上，因此用 $n^2$ 的時間進行枚舉，旋轉卡尺（類似最遠點對的作法）找面積最大的第三個點，就能在總時間複雜度 $O(n^2)$ 完成！（會再更少，因為只要枚舉凸包上的點）
以下是AC Code：
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 100001
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
ld x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(abs(x-b.x)<=eps && abs(y-b.y)<=eps)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
ld operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
ld operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
bool cmp(pt a, pt b){
if(a.x == b.x)return a.y < b.y;
return a.x < b.x;
}
vector<pt> p,hull;
int n,t,h;
ld ans;
bool check(pt a,pt b,pt o){
int cross = (a - o)^(b - o);
return cross >= 0; //這裡很關鍵，別吃WA
}
bool check2(pt a,pt b,pt c,pt d){
ld aa = (a - c)^(b - c);
ld bb = (a - d)^(b - d);
return aa < bb;
}
ld area(pt a,pt b){
return abs(a^b)/2;
}
void convex_hull(){
hull.clear();
sort(p.begin(),p.end(),cmp); //首先對x進行排序
for(auto i : p){ //依序走訪，如果遇到外積<0則不在凸包上
while(hull.size()>=2 && check(i,hull[hull.size()-1],hull[hull.size()-2])){
hull.pop_back();
}
hull.push_back(i); //在凸包hull的每一點都符合外積小於0
}
int down_hull = hull.size();
h = down_hull-1;
hull.pop_back(); //x最大的點會在凸包上，不用做兩次先pop一次
reverse(p.begin(),p.end()); //將所有點逆序之後做一次上面的凸包
for(auto i: p){
while(hull.size() > down_hull && check(i,hull[hull.size()-1],hull[hull.size()-2])){
hull.pop_back();
}
hull.push_back(i);
}
hull.pop_back();
}
void solve(){
int d,sz = hull.size();
rep(i,0,sz-1){
d = i+1;
rep(j,i+1,sz-1){
while(check2(hull[i],hull[(j)%sz],hull[d],hull[(d+1)%sz]))
d = (d+1)%sz;
ans = max(ans,area((hull[d]-hull[i]),(hull[d]-hull[j])));
}
}
}
signed main(){
Orz;
cin>>n;
p.assign(n,{0,0});
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
convex_hull();
ans = 0;
solve();
cout<<fixed<<setprecision(6);
cout<<ans<<endl;
}
```
### TIOJ 1500 Clean up on aisle 3
[題目連結](https://tioj.ck.tp.edu.tw/problems/1500)
[Submission](https://tioj.ck.tp.edu.tw/submissions/262966)
> 題目敘述
平面上n個點找最近點對的距離
最近點對真的有超多種作法的，枚舉、掃描線、分治、隨機都可以做！這邊有[一篇筆記](/hVplrqxCRdGiMkn2lwXGAA)比較各種時間複雜度的最近點對作法，這邊不多做贅述！
以下程式碼是掃描線演算法，最差情況下的時間複雜度是 $O(n^2)$，因為需要排序，所以下限為 $\Omega(n\log n)$！
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 50005
#define all(x) x.begin(),x.end()
#define INF 5e18
#define eps 1e-9
#define x first
#define y second
using namespace std;
int n;
pii p[N];
ld dis(pii a, pii b){
ld x = a.x-b.x, y = a.y-b.y;
return sqrt(x*x + y*y);
}
signed main(){
Orz;
cout<<fixed<<setprecision(6);
while(cin>>n){
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
sort(p,p+n);
ld d = INF;
rep(i,0,n-1){
rep(j,i+1,n-1){
if(p[j].x > p[i].x + d)break;
d = min(d, dis(p[i],p[j]));
}
}
cout<<d<<endl;
}
}
```
### TIOJ 1280 領土 (Territory)
[題目連結](https://tioj.ck.tp.edu.tw/problems/1280)
[Submission](https://tioj.ck.tp.edu.tw/submissions/262848)
> 題目敘述
一個國家有 n 個安全哨，每一個都有座標 $(x,y)$ ，代表在座標軸上的位置。輸出該國安全哨所能圍出的最大領土。
n個點所能圍成的最大面積，其實等價於凸包的面積。與前幾題的**最小凸多邊形**是一模一樣的題目！
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 100001
#define all(x) x.begin(),x.end()
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
int x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(x-b.x == 0 && y-b.y == 0)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
vector<pt> p,temp,pp;
vector<int> cnt;
int n,ans = 0;
bool cmp(pt a, pt b){
if(a.x == b.x)return a.y < b.y;
return a.x < b.x;
}
bool check(pt a,pt b,pt o){
pt aa = a - o;
pt bb = b - o;
return (aa^bb) >= 0;
}
vector<pt> solve(){
sort(all(p),cmp);
vector<pt> h;
for(pt i : p){
while(h.size()>=2 && check(i,h[h.size()-1],h[h.size()-2]))
h.pop_back();
h.push_back(i);
}
int sz = h.size();
h.pop_back();
reverse(all(p));
for(auto i : p){
while(h.size()>sz && check(i,h[h.size()-1],h[h.size()-2]))
h.pop_back();
h.push_back(i);
}
return h;
}
signed main(){
Orz;
cin>>n;
p.resize(n,{0,0});
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
vector<pt> hull = solve();
int area = 0,sz = hull.size();
rep(i,0,sz-2){
area += (hull[i]^hull[i+1]);
}
cout<<((area%2)?(area/2)+1:(area/2))<<endl;
}
```
### ZJ a871: Museum Area
[題目連結](https://zerojudge.tw/ShowProblem?problemid=a871)
> 題目敘述
n個點圍成的多邊形，求面積
水題，直接套行列式公式即可算出答案！
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pdd pair<double,double>
#define int long long
#define ld long double
#define N 15
#define x first
#define y second
using namespace std;
int n;
pdd p[N];
ld check(pdd a, pdd b){
return a.x*b.y - a.y*b.x;
}
signed main(){
Orz;
while(cin>>n){
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
p[n] = p[0];
ld area = 0.0;
rep(i,0,n-1)area += check(p[i],p[i+1]);
ld ans = (ld)area/2;
if(ans<0)ans = -ans;
cout<<fixed<<setprecision(2);
cout<<ans<<endl;
}
}
```
### TIOJ 1678 剪多邊形（molding）
[題目連結TIOJ](https://tioj.ck.tp.edu.tw/problems/1678)
[TIOJ Submission](https://tioj.ck.tp.edu.tw/submissions/262849)
[題目連結ZJ](https://zerojudge.tw/ShowProblem?problemid=d546)
> 題目敘述
間單來說是求出多邊形面積以及凸包面積的差，詳細可以點上面題目連結。
題目說多邊形需要才剪下的面積，我們就算凸包面積以及多邊形面積，兩者的差去除上題目給的色塊面積即是答案！
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define int long long
#define ll long long
#define ld long double
#define N 100001
#define all(x) x.begin(),x.end()
#define eps 1e-9
#define x first
#define y second
using namespace std;
struct pt{
int x,y;
bool operator < (pt b){
if(x == b.x)return y < b.y;
return x < b.x;
}
bool operator > (pt b){
if(x == b.x)return y > b.y;
return x > b.x;
}
bool operator == (pt b){
if(x-b.x == 0 && y-b.y == 0)return true;
return false;
}
pt operator+(pt b) {return {x + b.x, y + b.y};} //向量相加
pt operator-(pt b) {return {x - b.x, y - b.y};} //向量相減
int operator^(pt b) {return x * b.y - y * b.x;} //向量外積cross
int operator*(pt b) {return x * b.x + y * b.y;} //向量內積dot
};
vector<pt> p,temp,pp;
vector<int> cnt;
int n,a,ans = 0;
bool cmp(pt a, pt b){
if(a.x == b.x)return a.y < b.y;
return a.x < b.x;
}
bool check(pt a,pt b,pt o){
pt aa = a - o;
pt bb = b - o;
return (aa^bb) >= 0;
}
vector<pt> solve(){
sort(all(p),cmp);
vector<pt> h;
for(pt i : p){
while(h.size()>=2 && check(i,h[h.size()-1],h[h.size()-2]))
h.pop_back();
h.push_back(i);
}
int sz = h.size();
h.pop_back();
reverse(all(p));
for(auto i : p){
while(h.size()>sz && check(i,h[h.size()-1],h[h.size()-2]))
h.pop_back();
h.push_back(i);
}
return h;
}
int Polygon(){
vector<pt> temp(p);
temp.push_back(temp[0]);
int area = 0;
rep(i,0,n-1)area += (temp[i]^temp[i+1]);
return abs(area);
}
signed main(){
Orz;
cin>>n>>a;
p.resize(n,{0,0});
rep(i,0,n-1)cin>>p[i].x>>p[i].y;
ld small = (ld)Polygon()/2;
vector<pt> hull = solve();
int area = 0,sz = hull.size();
rep(i,0,sz-2)area += (hull[i]^hull[i+1]);
ld big = (ld)area/2;
int ans = ceil((big-small)/a);
cout<<ans<<endl;
}
```
### ZJ d269: 11579 - Triangle Trouble
[題目連結](https://zerojudge.tw/ShowProblem?problemid=d269)
> 題目敘述
有一個三角形工廠有一個很大的問題。給你一些邊的邊長，想辦法找出用這些邊長圍出最大的三角形。
根據海龍公式，三角形面積：
$$\triangle ABC = \sqrt{s\cdot(s-a)\cdot(s-b)\cdot(s-c)}$$
可以利用貪婪法，將所有邊長由大到小進行排序，每一次拿最大的三個邊長進行枚舉，即可算出最大的三角形面積。不難理解，當換上一個比較大的邊，算出來的s也會比較大，跟邊相減的值也會比較大，總面積自然較大（好啦，這是非常不嚴謹的證明XD）
在想題過程中，我有思考到，如果周長一樣的情況下，到底何種面積的三角形面積會比較大？答案是正三角形！
:::info
**三角形周長固定下面積的比較**
根據海龍公式：
$$s = \frac{1}{2}(a+b+c)$$
想要比較在周長固定下三角形的面積，可以用算幾不等式比較，因為 $s$ 是定值，所以可以列出以下式子：
$$\frac{(s-a)+(s-b)+(s-c)}{3} ≥ \sqrt[3]{(s-a)(s-b)(s-c)}$$
等好成立時，$a=b=c$。因為$s = \frac{a+b+c}{2}$，因此：
$$(\frac{a}{2})^2 ≥ (s-a)(s-b)(s-c)$$
得到海龍公式
$$\triangle ABC = \sqrt{s\cdot(s-a)\cdot(s-b)\cdot(s-c)} ≤ \sqrt{\frac{3a}{2}\cdot\frac{a^3}{8}}=\frac{\sqrt{3}}{4}a^2$$
:::
以下是使用貪婪法的AC Code：
```cpp=
#include <bits/stdc++.h>
#define Orz ios::sync_with_stdio(0),cin.tie(0)
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define pii pair<int,int>
#define pdd pair<double,double>
#define ll long long
#define ld double
#define N 100001
#define all(x) x.begin(),x.end()
#define eps 1e-9
#define x first
#define y second
using namespace std;
int t,n;
vector<ld> p;
ld area(ld a ,ld b, ld c){
if(a > b + c)return -1;
ld p = (a+b+c)/2;
return p*(p-a)*(p-b)*(p-c);
}
signed main(){
Orz;
cin>>t;
while(t--){
cin>>n;
p.assign(n,0);
rep(i,0,n-1)cin>>p[i];
sort(all(p),greater<>());
ld ans = 0;
rep(i,0,n-3)
ans = max(ans,area(p[i],p[i+1],p[i+2]));
cout<<fixed<<setprecision(2);
cout<<sqrt(ans)<<endl;
}
}
```
### 1224 . 矩形覆蓋面積計算
[題目連結](https://tioj.ck.tp.edu.tw/problems/1224)
[Submission2:AC](https://tioj.ck.tp.edu.tw/submissions/261541)
> 題意：給你平面上n個矩形，請求出它們覆蓋的總表面積。
這一題所使用的技巧是<font color="#f00">掃描線</font>以及<font color="#f00">線段樹</font>，下圖中的水平藍色線即為掃描線，由y=0開始往上掃描，當遇到了矩形的邊，利用線段樹查詢區間內當前的矩形寬度，乘上兩掃描線的高度差即為面積。當然，掃描線也可以使用垂直方向的線段由左而右的掃描，實作細節是一樣的。
#### 線段樹維護
##### 方法一
我們可以定義線段樹$seg[cur]$為區間$[l,r]$中有被**矩形覆蓋的大小**有多大，也就是圖中當前掃描線對應到的區域的寬度。這樣子維護有一個問題，當我們直接用$seg[cur]$儲存答案，我們在修改的時候沒有辦法確切知道這段區間被覆蓋的情況。
下圖為一種模擬的情況，每一個區間的數字代表著非0的數字個數，也就是它的寬度。今天我們要對區間$[4,6]$加減值，將區間拆成$[4,4]$跟$[5,6]$，這時候區間$[3,4]$的數值是1，我們卻不知道到底是3還是4是有被覆蓋到的，必須要遞迴下去到葉節點才能得到完整的覆蓋情況，這時候每一次加減值的複雜就會提升到$O(n)$，因此不能以這種方式維護。
##### 方法二
有別於第一種方法對$seg[id]$進行維護，我們可以多開一個區間 $tag$ 來紀錄被矩形覆蓋的情況。下圖有3個矩形，其中的數字代表每一塊區域被覆蓋的情況，這邊使用了$tag$來紀錄（他是附在區間上的，不會像圖中一樣的方式呈現）。tag的數值為非負整數，紀錄當前區間有多少矩形覆蓋在上面，用$tag$來輔助維護$seg[id]$可以在$O(logn)$的時間進行修改與查詢。
以下程式碼是是 $tag$ 的轉移，當大的區間的tag值不為0，代表有一個矩形曾完整覆蓋這個區間，這時候可以直接回傳區間大小，否則即回傳左右節點的$sed[left],seg[right]$的數值。
這邊定義$seg[id]$為：「考慮 id 的子孫們（不含 id 本身）的所有 tag 值，假設這些子孫只有被tag值作用過，共有多少非0的數字」。
```cpp=
seg[cur].val = (seg[2*cur].tag?mid-l:seg[2*cur].val)
+(seg[2*cur+1].tag?r-mid:seg[2*cur+1].val);
```
#### 實作方法
##### 矩形維護
首先是維護矩形的方法。我們一個矩形總共要維護四個東西：矩形左界x1、矩形右界x2、矩形上下界的y座標（分上下兩條），這兩條邊是下界或是上界val。為什麼要水平方向要分兩條討論？是因為下界代表進入，當掃描線掃到這一條邊的時候表示我們要新增區間 $[x1,x2)$ 進入線段樹；反之如果掃到了上界，則表示離開這個矩形，在線段樹中扣掉區間 $[x1,x2)$。
```cpp=
struct Node{ //每一個矩陣分成上下兩條邊
int x1; //矩形左界x1
int x2; //矩形右界x2
int y; //矩形y座標（分上下兩邊）
int val; //val = ±1(進入代表1、離開代表-1)
}arr[2*N];
```
上下界我們利用val維護，當 $val=1$ 時表示是矩形的下界； $val=-1$ 則是矩形上界，這兩個搭配在一起剛好就可以用線段樹區間加值的方式進行操作！總共有 $n$ 個矩形，因此我們要掃描線總共掃描 $2n$ 條線段。
##### 線段樹
一樣對值域（這題是1000000）的4倍開了線段樹，同時維護一個非負整數 $tag$ 表示區間被覆蓋的情況。當每一次修改完成之後，我們可以直接取用根節點 $seg[1]$ 的數值表示寬度（非0的個數）！
```cpp=+
//seg[i]表示i的左右兩子樹的區間非0的個數
struct node{ //建立線段樹
int val; //維護非0個數
int tag; //使用tag紀錄區間被覆蓋次數
}seg[4*M];
```
接下來就是在程式執行的過程中將 $2n$ 條邊依照y座標進行排序 $O(nlogn)$，接著依序使用掃描線搭配線段樹的修改，計算矩形的面積。最後就是輸出加起來的答案。
:::warning
**Debug 小錯誤**
[Submission1:WA](https://tioj.ck.tp.edu.tw/submissions/261512)
可以看到有一筆測資過不了，95分QQQ
後來debug之後發現到，因為我是對每一個矩形先輸入下界之後才是上界，當我在排序的過程中，上界有可能有機會跑到下界之前，造成 $tag$ 被扣到負的情況，但在定義中可以清楚知道 $tag$ 是非負整數造成錯誤。因此只要把排序的過程改成 stable_sort() 即可！
```cpp=
stable_sort(arr,arr+(n<<1),cmp);
```
:::
最後終於是程式碼的部分，以下：
```cpp=
#include <bits/stdc++.h>
#define ios ios::sync_with_stdio(0),cin.tie(0);
#define N 100005
#define M 1000001
#define lld long long
using namespace std;
int n;
struct Node{ //每一個矩陣分成上下兩條邊
int x1; //矩形左界x1
int x2; //矩形右界x2
int y; //矩形y座標（分上下兩邊）
int val; //val = ±1(進入代表1、離開代表-1)
}arr[2*N];
//seg[i]表示i的左右兩子樹的區間非0的個數
struct node{ //建立線段樹
int val; //維護非0個數
int tag; //使用tag紀錄區間被覆蓋次數
}seg[4*M];
bool cmp(Node a, Node b){
return a.y<b.y;
}
//對區間[ql,qr)進行加值val
void modify(int cur,int l,int r,int ql,int qr,int val){
if(r <= l || ql >= r || qr <= l)return;
if(ql <= l && qr >= r){
seg[cur].tag += val;
return;
}
int mid = (l+r)/2;
modify(2*cur,l,mid,ql,qr,val);
modify(2*cur+1,mid,r,ql,qr,val);
//左右節點如有tag表示被完全覆蓋，直接加上區間大小，否則加上seg[左右子樹]
seg[cur].val = (seg[2*cur].tag?mid-l:seg[2*cur].val)
+(seg[2*cur+1].tag?r-mid:seg[2*cur+1].val);
}
int main(){
ios;
memset(arr,0,sizeof(arr));
memset(seg,0,sizeof(seg));
cin>>n; //依序輸入左右下上：x1,x2,y1,y2
for(int i=0;i<(n<<1);i+=2){
int x1,x2,y1,y2;cin>>x1>>x2>>y1>>y2;
arr[i] = (Node){x1,x2,y1,1}; //插入矩形下邊，帶入val = 1
arr[i+1] = (Node){x1,x2,y2,-1}; //上邊要val = -1
}
stable_sort(arr,arr+(n<<1),cmp); //依照y座標由小到大排序
int y0 = 0,val = 0; //有下而上的枚舉所有水平邊
lld ans = 0LL; //上一條y的座標，計算高，val為矩形結合起來的寬
for(int i=0;i<(n<<1);i++){ //枚舉2n條y的邊
ans += (lld)(arr[i].y-y0)*val; //計算面積（寬*高）
modify(1,0,M,arr[i].x1,arr[i].x2,arr[i].val);
y0 = arr[i].y;
val = seg[1].val; //修改後（下一輪）的矩陣寬度
}
cout<<ans<<'\n';
}
```
## CSES Geometry
- [x] Point Location Test
- [x] Line Segment Intersection
- [x] Polygon Area
- [x] Point in Polygon
- [x] Polygon Lattice Points
- [x] Minimum Euclidean Distance
- [x] Convex Hull
### Point Location Test
[CODE](https://cses.fi/paste/3a495a0f4bdec8af2df7b4/)
給你三個點 $P_1,P_2,P_3$ ，判斷出向量 $(P_1,P_2)$ 之於 $P_3$ 的方向為何？
關係共有相交、位於左側以及位於右側。
### Line Segment Intersection
[CODE](https://cses.fi/paste/854258b8f67b255c2df828/)
給你兩條線段，判斷他們是否相交。線段相交裸題。
### Polygon Area
[CODE](https://cses.fi/paste/49f2a020a47797b42df861/)
給你N個點，計算出此多邊形圍成的面積為何。帶入行列式公式可解。
### Point in Polygon
[CODE](https://cses.fi/paste/d0edfaf7380124852e075e/)
給你N個點構成的一個簡單多邊形，判斷一個點是否在多邊形內部。
參考資料：[greekforgreek](https://www.geeksforgeeks.org/how-to-check-if-a-given-point-lies-inside-a-polygon/)
凸多邊形的情況，我們可以利用內角和是360度轉一圈的的方法，判斷一點是否在多邊形內部。但本題是簡單多邊形，因此用這種方法是不可行的。
這張圖是作法，我們做出一個由x點出發向x軸正向的射線，計算這一條射線跟多邊形共有幾個交點。如果是奇數，表示在內部；反之則是在外部。交在多邊形上的點要特別注意，因為他可能會被統計到兩次，因此我們特別處理當點相交一點時，只計算那條邊的兩端點中，y座標比較大的那一點是否相交的情況。
最後，我們發現遇上如圖中g點的狀況直接就被處理好了！他會被算到零次，因為在看兩條邊加上去之後，又因為頂點是y座標比較大的，會被剪掉兩次。
### Polygon Lattice Points
[CODE](https://cses.fi/paste/7fa78ccd60f7877f2e07dc/)
給定頂點座標均是整點（或正方形格子點）的簡單多邊形，皮克定理說明了其面積 $A$ 和內部格點數目 $i$、邊上格點數目 $b$ 的關係：
$$A = i + \frac{b}{2} - 1$$
[PROOF](http://episte.math.ntu.edu.tw/articles/sm/sm_25_10_1/index.html)
### Minimum Euclidean Distance
[CODE](https://cses.fi/paste/a728ec171111ce792e091f/)
[最近點對問題](https://hackmd.io/@peienwu/closest_pair)。不過，既然是計算幾何，我們就用掃描線做最近點對。掃描線有兩個做法（可以參考那個連結），至於搭配set輔助步驟就是：
1. 將點輸入並且排序，X座標為主，Y座標為輔。
2. 使用set，並以Y座標為排序基準（pair的首項），以儲存第 $i$ 點的左方、水平距離小於等於d的點。
3. 右掃描線依序窮舉各點作為右端點。
(1) Erase與右端點水平距離大於d的點們（左掃描線右移）
(2) 用二分搜找出與第 $i$ 點垂直距離小於d的點，並嘗試更新
(3) 將第 $i$ 點加入set中。
### Convex Hull
[CODE](https://cses.fi/paste/ac97fa8545e467692e09aa/)
凸包裸題。
## Leetcode：Geography
Leetcode 總共有27題的tag是計算幾何的，[題單在這裡](https://leetcode.com/tag/geometry/)。有三題被鎖起來不能看，所以總共有24題。
### 149 Max Points on a Line
[題目連結](https://leetcode.com/problems/max-points-on-a-line)
```cpp=
class Solution {
public:
int gcd(int a,int b){
if(a > b)swap(a,b);
if(a == 0)return b;
return gcd(b % a, a);
}
int maxPoints(vector<vector<int>>& points) {
int n = points.size(),ans = 0;
for(int i = 0;i < n;i++){
map<pair<int,int>,int> mp;
int same = 0,hori = 0;
for(int j = i+1;j < n;j++){
int dx = points[i][0] - points[j][0];
int dy = points[i][1] - points[j][1];
if(dx == 0 && dy == 0){same++;continue;}
if(dx == 0 && dy != 0){hori++;continue;}
int g = gcd(abs(dx),abs(dy));
if(dy < 0 || (dy == 0 && dx < 0)){dx = -dx;dy = -dy;}
mp[{dx/g,dy/g}]++;
}
int sum = hori + same + 1;
for(auto it : mp)sum = max(sum,it.second + 1);
ans = max(ans,sum);
}
return ans;
}
};
```
### 223 Rectangle Area
[題目連結](https://leetcode.com/problems/rectangle-area)
```cpp=
class Solution {
public:
int computeArea(int ax1, int ay1, int ax2, int ay2, int bx1, int by1, int bx2, int by2) {
int x = max(min(ax2,bx2) - max(ax1,bx1),0);
int y = max(min(ay2,by2) - max(ay1,by1),0);
int ans = (ax2 - ax1)*(ay2 - ay1)+(bx2 - bx1)*(by2 - by1)-x * y;
return ans;
}
};
```
### 335 Self Crossing
[題目連結](https://leetcode.com/problems/self-crossing)
```cpp=
class Solution {
public:
bool isSelfCrossing(vector<int>& d) {
int n = d.size();
if(n <= 3)return false;
for(int i = 3;i < n;i++){
//第四條、第五條、第六條交第一條
if(d[i] >= d[i-2] && d[i-1] <= d[i-3])return true;
if(i >= 4 && d[i-1] == d[i-3] && d[i] + d[i-4] >= d[i-2])return true;
if(i >= 5 && d[i-2] >= d[i-4] && d[i] >= d[i-2] - d[i-4]
&& d[i-1] >= d[i-3]-d[i-5] && d[i-1] <= d[i-3])return true;
}
return false;
}
};
```
### 587 Erect the Fence
[題目連結](https://leetcode.com/problems/erect-the-fence)
```cpp=
#define pii pair<int,int>
#define ff first
#define ss second
class Solution {
public:
bool check(pii a,pii b,pii o){
pii aa = {a.ff - o.ff,a.ss - o.ss};
pii bb = {b.ff - o.ff,b.ss - o.ss};
int cross = aa.ff * bb.ss - aa.ss * bb.ff;
return cross > 0;
}
vector<vector<int>> outerTrees(vector<vector<int>>& point) {
vector<pii> h;
int n = point.size();
vector<pii> p(n);
for(int i = 0;i < n;i++)p[i] = {point[i][0],point[i][1]};
sort(p.begin(),p.end());
for(auto i : p){
while(h.size() > 1 && check(i,h[h.size()-1],h[h.size()-2]))
h.pop_back();
h.push_back(i);
}
int down = h.size();
h.pop_back();
reverse(p.begin(),p.end());
for(auto i : p){
while(h.size() > down && check(i,h[h.size()-1],h[h.size()-2]))
h.pop_back();
h.push_back(i);
}
set<pii> s;for(auto i : h)s.insert(i);
n = s.size();
vector<vector<int>> ans;ans.resize(n);
for(int i=0;i<n;i++)ans[i].resize(2);
int id = 0;
for(auto i : s){ans[id][0] = i.ff;ans[id][1] = i.ss;id++;}
return ans;
}
};
```
### 593 Valid Square
[題目連結](https://leetcode.com/problems/valid-square)
```cpp=
class Solution {
public:
int dis(vector<int>& p1, vector<int>& p2){
int x = p1[0] - p2[0],y = p1[1] - p2[1];
return x * x + y * y;
}
bool validSquare(vector<int>& p1, vector<int>& p2, vector<int>& p3, vector<int>& p4) {
map<int,int>mp; //長度、個數
mp[dis(p1,p2)]++;
mp[dis(p1,p3)]++;
mp[dis(p1,p4)]++;
mp[dis(p2,p3)]++;
mp[dis(p2,p4)]++;
mp[dis(p3,p4)]++;
return mp.size() == 2 && mp.begin()->second == 4;
}
};
```
### 812 Largest Triangle Area
[題目連結](https://leetcode.com/problems/largest-triangle-area)
```cpp=
class Solution {
public:
int cross(int x1,int y1,int x2,int y2){
return x1 * y2 - x2 * y1;
}
double area(int a,int b,int c,int d,int e,int f){
double sum = 0;
sum += cross(a,b,c,d);
sum += cross(c,d,e,f);
sum += cross(e,f,a,b);
return abs(sum / 2.0);
}
double largestTriangleArea(vector<vector<int>>& points) {
int n = points.size();
double ans = 0;
for(int i = 0;i < n;i++){
for(int j = i + 1;j < n;j++){
for(int p = j + 1;p < n;p++){
ans = max(ans,area(points[i][0],points[i][1]
,points[j][0],points[j][1]
,points[p][0],points[p][1]));
}
}
}
return ans;
}
};
```
### 836 Rectangle Overlap
[題目連結](https://leetcode.com/problems/rectangle-overlap)
```cpp=
class Solution {
public:
bool isRectangleOverlap(vector<int>& rec1, vector<int>& rec2) {
//最小的右端點 - 最大的左端點
int x = min(rec1[2],rec2[2]) - max(rec1[0],rec2[0]);
//最小的上端點 - 最大的下端點
int y = min(rec1[3],rec2[3]) - max(rec1[1],rec2[1]);
return x > 0 && y > 0;
}
};
```
### 858 Mirror Reflection
[題目連結](https://leetcode.com/problems/mirror-reflection)
```cpp=
class Solution {
public:
int gcd(int a,int b){
if(a == 0)return b;
return gcd(b % a,a);
}
int mirrorReflection(int p, int q) {
int len = p * q / gcd(q,p);
int a = (len / p) % 2,b = (len / q) % 2;
if(a == 0 && b == 1)return 0;
else if(a == 1 && b == 1)return 1;
else return 2;
}
};
```
### 478 Generate Random Point in a Circle
[題目連結](https://leetcode.com/problems/generate-random-point-in-a-circle)
直接Random半徑會出事（可能不夠亂，或是半徑太小），如果random面積之後算半徑才OK。
```cpp=
class Solution {
public:
double R,X,Y;
Solution(double radius, double x_center, double y_center) {
R = radius;X = x_center;Y = y_center;
srand(time(NULL));
}
vector<double> randPoint() {
double Area = rand() * R * R * M_PI / (RAND_MAX + 1.0);
double r = sqrt(Area / M_PI);
double theta = 2.0 * M_PI * rand() / (RAND_MAX + 1.0);
vector<double> ans;
ans.push_back(X + r * cos(theta));
ans.push_back(Y + r * sin(theta));
return ans;
}
};
/**
* Your Solution object will be instantiated and called as such:
* Solution* obj = new Solution(radius, x_center, y_center);
* vector<double> param_1 = obj->randPoint();
*/
```
### 883 Projection Area of 3D Shapes
[題目連結](https://leetcode.com/problems/projection-area-of-3d-shapes)
三種不同的投影對應到三種不同的角度看圖形。x-y的面積即為由上而下看有方格的個數。x-z是從前方看，因此對應到的是每一行的最大方塊個數。
```cpp=
class Solution {
public:
int projectionArea(vector<vector<int>>& grid) {
int n = grid.size(),ans = 0;
for(int i = 0;i < n;i++){
int maxR = 0,maxC = 0;
for(int j = 0;j < n;j++){
if(grid[i][j] > 0)ans++; //由上往下看
maxR = max(maxR,grid[i][j]); //由側邊看
maxC = max(maxC,grid[j][i]); //由前面看
}
ans += maxR + maxC;
}
return ans;
}
};
```
### 892 Surface Area of 3D Shapes
[題目連結](https://leetcode.com/problems/surface-area-of-3d-shapes)
```cpp=
class Solution {
public:
int surfaceArea(vector<vector<int>>& grid) {
int ans = 0,n = grid.size();
for(int i = 0;i < n;i++){
for(int j = 0;j < n;j++){
if(grid[i][j] > 0)ans += 4 * grid[i][j] + 2;
if(i < n - 1){
ans -= min(grid[i][j],grid[i+1][j])*2;
}
if(j < n - 1){
ans -= min(grid[i][j],grid[i][j+1])*2;
}
}
}
return ans;
}
};
```
### 939 Minimum Area Rectangle
[題目連結](https://leetcode.com/problems/minimum-area-rectangle)
```cpp=
class Solution {
public:
int minAreaRect(vector<vector<int>>& points) {
vector<vector<int>>p = points;
int n = points.size();
set<pair<int,int>>s;
for(int i=0;i<n;i++)s.insert({p[i][0],p[i][1]});
int ans = INT_MAX;
for(int i = 0;i < n;i++){
for(int j = i+1;j < n;j++){
if(s.find({p[i][0],p[j][1]})!=s.end()
&& s.find({p[j][0],p[i][1]})!=s.end()){
if(p[i][0] == p[j][0] || p[i][1] == p[j][1])continue;
ans = min(ans,abs(p[i][0]-p[j][0])*abs(p[i][1]-p[j][1]));
}
}
}
if(ans == INT_MAX)return 0;
else return ans;
}
};
```
### 963 Minimum Area Rectangle II
[題目連結](https://leetcode.com/problems/minimum-area-rectangle-ii)
```cpp=
class Solution {
public:
double minAreaFreeRect(vector<vector<int>>& p){
set<pair<int,int>> s;
int n = p.size();
for(int i = 0;i < n;i++)s.insert({p[i][0],p[i][1]});
double ans = 1e9;
for(int i = 0;i < n;i++){
for(int j = i+1;j < n;j++){
int x1 = p[j][0] - p[i][0];
int y1 = p[j][1] - p[i][1];
for(int k = j+1;k < n;k++){
int x2 = p[k][0] - p[i][0];
int y2 = p[k][1] - p[i][1];
if(x1 * x2 + y1 * y2 != 0)continue;
int nx = p[k][0] + x1,ny = p[k][1] + y1;
if(s.find({nx,ny}) != s.end()){
ans = min(ans,sqrt(x1*x1+y1*y1) * sqrt(x2*x2+y2*y2));
}
}
}
}
if(ans == 1e9)return 0;
return ans;
}
};
```
### 973 K Closest Points to Origin
[題目連結](https://leetcode.com/problems/k-closest-points-to-origin)
```cpp=
class Solution {
public:
vector<vector<int>> kClosest(vector<vector<int>>& points, int k) {
int n = points.size();
multimap<int,int> mp;
for(int i = 0;i < n;i++){
int x = points[i][0];
int y = points[i][1];
mp.insert({x * x + y * y,i});
}
auto it = mp.begin();
vector<vector<int>> ans;
for(int i=0;i<k;i++){
int id = it->second;
ans.push_back(points[id]);
it++;
}
return ans;
}
};
```
### 1030 Matrix Cells in Distance Order
[題目連結](https://leetcode.com/problems/matrix-cells-in-distance-order)
```cpp=
class Solution {
public:
vector<vector<int>> allCellsDistOrder(int rows, int cols, int rCenter, int cCenter) {
int n = rows,m = cols;
multimap<int,pair<int,int>> mp;
for(int i = 0;i < n;i++){
for(int j = 0;j < m;j++){
int dis = abs(i - rCenter) + abs(j - cCenter);
mp.insert({dis,{i,j}});
}
}
vector<vector<int>> ans;ans.resize(n*m);
for(int i=0;i<n*m;i++)ans[i].resize(2,0);
int id = 0;
for(auto it : mp){
ans[id][0] = (it.second.first);
ans[id][1] = (it.second.second);
id++;
}
return ans;
}
};
```
### 1037 Valid Boomerang
[題目連結](https://leetcode.com/problems/valid-boomerang)
```cpp=
class Solution {
public:
bool isBoomerang(vector<vector<int>>& points) {
int x1 = points[1][0]-points[0][0],y1 = points[1][1]-points[0][1];
int x2 = points[2][0]-points[1][0],y2 = points[2][1]-points[1][1];
int cross = x1 * y2 - x2 * y1;
if(cross == 0)return false;
return true;
}
};
```
### 1232 Check If It Is a Straight Line
[題目連結](https://leetcode.com/problems/check-if-it-is-a-straight-line)
```cpp=
class Solution {
public:
bool checkStraightLine(vector<vector<int>>& coordinates) {
vector<vector<int>> p = coordinates;
int n = p.size();
for(int i = 0;i < n-2;i++){
int x1 = p[i+1][0] - p[i][0],y1 = p[i+1][1] - p[i][1];
int x2 = p[i+2][0] - p[i+1][0],y2 = p[i+2][1] - p[i+1][1];
if(x1 * y2 != x2 * y1)return false;
}
return true;
}
};
```
### 1266 Minimum Time Visiting All Points
[題目連結](https://leetcode.com/problems/minimum-time-visiting-all-points)
```cpp=
class Solution {
public:
int minTimeToVisitAllPoints(vector<vector<int>>& points) {
int n = points.size(),ans = 0;
vector<vector<int>> p = points;
for(int i = 0;i < n-1;i++){
ans += max(abs(p[i+1][0] - p[i][0]),abs(p[i+1][1] - p[i][1]));
}
return ans;
}
};
```
### 1401 Circle and Rectangle Overlapping
[題目連結](https://leetcode.com/problems/circle-and-rectangle-overlapping)
```cpp=
class Solution {
public:
bool checkOverlap(int radius, int x_center, int y_center, int x1, int y1, int x2, int y2) {
int x = clamp(x_center,x1,x2) - x_center;
int y = clamp(y_center,y1,y2) - y_center;
return x*x + y*y <= radius * radius;
}
};
```
### 1453 Maximum Number of Darts Inside of a Circular Dartboard
[題目連結](https://leetcode.com/problems/maximum-number-of-darts-inside-of-a-circular-dartboard)
```cpp=
#define ld long double
#define pdd pair<ld,ld>
#define x first
#define y second
#define exp 1e-6
class Solution {
public:
double dis(pdd a,pdd b){
ld x = a.x - b.x,y = a.y - b.y;
return sqrt(x * x + y * y);
}
pair<pdd,pdd> get_center(pdd a,pdd b,ld R){
pdd mid = {(a.x + b.x) / 2,(a.y + b.y) / 2};
ld theta = atan2(a.y - b.y, b.x - a.x);
ld tmp = dis(a,b) / 2, d = sqrt(R * R - tmp * tmp);
pair<pdd,pdd> ans;
ans.x = {mid.x - d * sin(theta),mid.y - d * cos(theta)};
ans.y = {mid.x + d * sin(theta),mid.y + d * cos(theta)};
return ans;
}
int numPoints(vector<vector<int>>& point, int R) {
int n = point.size(),ans = 1;
pdd p[n];for(int i=0;i<n;i++){p[i] = {point[i][0],point[i][1]};}
for(int i = 0;i < n;i++){
for(int j = i+1;j < n;j++){
if(dis(p[i],p[j]) - 2.0 * R >= exp)continue;
pair<pdd,pdd> cur = get_center(p[i],p[j],R);
int cnt = 0;
for(int k = 0;k < n;k++)
if(dis(p[k], cur.x) - R<= exp)cnt ++;
ans = max(ans, cnt);cnt = 0;
for(int k = 0;k < n;k++)
if(dis(p[k], cur.y) - R <= exp)cnt ++;
ans = max(ans, cnt);
}
}
return ans;
}
};
```
### 1515 Best Position for a Service Centre
[題目連結](https://leetcode.com/problems/best-position-for-a-service-centre)
這一題是模擬退火，非常酷，之前沒有寫過的東西。
```cpp=
#define eps 1e-6
#define ld long double
#define pdd pair<ld,ld>
#define x first
#define y second
class Solution {
public:
pdd p[105];int n;
ld dis_all(pdd mid){
ld sum;
for(int i = 0;i < n;i++){
ld x = p[i].x - mid.x,y = p[i].y - mid.y;
sum += sqrt(x * x + y * y);
}
return sum;
}
double getMinDistSum(vector<vector<int>>& pos) {
n = pos.size();
for(int i = 0;i < n;i++)p[i] = {pos[i][0],pos[i][1]};
pdd cur = p[0];ld mid_dis = dis_all(p[0]);
ld test_size = 100;
ld dx[4] = {1.0,0.0,-1.0,0.0},dy[4] = {0.0,1.0,0.0,-1.0};
bool flag = 0;
while(test_size > eps){
flag = 0;
for(int i = 0;i < 4;i++){
pdd newp = cur;
newp.x += dx[i] * test_size;
newp.y += dy[i] * test_size;
ld new_dis = dis_all(newp);
if(new_dis < mid_dis){
mid_dis = new_dis;
cur = newp;
flag = 1;
break;
}
}
if(flag == 0)test_size /= 2.0;
}
return mid_dis;
}
};
```
換另外一種迭代方式
```cpp=
#define eps 1e-8
#define ld long double
#define pdd pair<ld,ld>
#define x first
#define y second
class Solution {
public:
pdd p[105];int n;
ld dis_all(pdd mid){
ld sum;
for(int i = 0;i < n;i++){
ld x = p[i].x - mid.x,y = p[i].y - mid.y;
sum += sqrt(x * x + y * y);
}
return sum;
}
double getMinDistSum(vector<vector<int>>& pos) {
srand(time(NULL));
n = pos.size();
for(int i = 0;i < n;i++)p[i] = {pos[i][0],pos[i][1]};
pdd cur = p[0];ld mid_dis = dis_all(p[0]);
ld test_size = 150.0;
while(test_size > eps){
pdd newp = cur;
int temp = rand();
newp.x += cos(temp) * test_size;
newp.y += sin(temp) * test_size;
ld new_dis = dis_all(newp);
if(new_dis < mid_dis){
mid_dis = new_dis;
cur = newp;
}
else test_size *= 0.99;
}
return mid_dis;
}
};
```
### 1610 Maximum Number of Visible Points
[題目連結](https://leetcode.com/problems/maximum-number-of-visible-points)
```cpp=
#define ld long double
const int N = 1e5+5;
class Solution {
public:
int visiblePoints(vector<vector<int>>& points, int angle, vector<int>& loc) {
vector<ld> ang;
int overlap = 0,n = points.size(),ans = 0;
for(auto p : points){
if(p[0] == loc[0] && p[1] == loc[1])overlap++;
else ang.push_back(atan2l(p[1]-loc[1],p[0]-loc[0]) * 180 / (ld)M_PI);
}
int sz = ang.size();
for(int i = 0;i < sz;i++)ang.push_back(ang[i] + 360);
sort(ang.begin(),ang.end());
sz = ang.size();
for(int i = 0,it2 = 0;i < sz;i++){
while(it2 < sz && ang[it2] - ang[i] <= angle)it2++;
ans = max(ans,it2 - i);
}
return ans + overlap;
}
};
```
### 1828 Queries on Number of Points Inside a Circle
[題目連結](https://leetcode.com/problems/queries-on-number-of-points-inside-a-circle)
```cpp=
class Solution {
public:
vector<int> countPoints(vector<vector<int>>& points, vector<vector<int>>& queries) {
int n = queries.size(),m = points.size();
vector<int> ans;ans.resize(n);
for(int i = 0;i < n;i++){
int rx = queries[i][0],ry = queries[i][1],sum = 0;
for(int j = 0;j < m;j++){
int x = points[j][0] - rx,y = points[j][1] - ry;
if(x * x + y * y <= queries[i][2] * queries[i][2])sum++;
}
ans[i] = sum;
}
return ans;
}
};
```
### 2101. Detonate the Maximum Bombs
[題目連結](https://leetcode.com/problems/detonate-the-maximum-bombs/)
把圖形考慮成一張圖，符合條件的就增加一條有向邊，接著對每一個點做一次DFS即可。時間 $O(n^2)$。
```cpp=
#define pii pair<int,int>
#define ld long double
class Solution {
public:
vector<int> E[105];
ld dis(pii a,pii b){
int x = a.first - b.first,y = a.second - b.second;
return sqrt((long long)x * x + (long long)y * y);
}
int sum = 1;
bool vis[105];
void dfs(int now){
vis[now] = 1;
for(auto i : E[now]){
if(vis[i])continue;
sum++;vis[i] = 1;
dfs(i);
}
}
int maximumDetonation(vector<vector<int>>& bomb) {
int n = bomb.size(),ans = 0;
for(int i = 0;i < n;i++){
for(int j = 0;j < n;j++){
pii a = {bomb[i][0],bomb[i][1]},b = {bomb[j][0],bomb[j][1]};
int r1 = bomb[i][2],r2 = bomb[j][2];
ld d = dis(a,b);
if(r1 >= d)E[i].push_back(j);
if(r2 >= d)E[j].push_back(i);
}
}
for(int i = 0;i < n;i++){
fill(vis,vis+105,0);
sum = 1;
dfs(i);
ans = max(ans,sum);
}
return ans;
}
};
```
## 延伸問題
- [ ] 給定平面上N個點，問一條直線最多能穿越幾個點
- [ ] 線段相交座標
- [ ] 給你一個線段及一個圓，判斷圓跟線段的最短距離
- [ ] 給你一個三角形（其中一頂點為圓心）和圓形，求出兩者的交集面積
- [ ] 兩圓交點
- [ ] 圓跟多邊形交集面積
- [ ] 最小包覆圓
- [ ] 半平面交
- [ ] Bentley–Ottmann Algorithm
- [ ] Voronoi Diagram
- [ ] Delaunay Triangulation
## 心得
計算幾何，顧名思義就是在電腦完成幾何的運算，要怎麼把平面的東西轉化成電腦看得懂的東西就是計算幾何在做的事情。常常我們覺得很容易判斷的事情，例如判斷線段是否相交，我們可以利用肉眼直輕易判斷出來，因為我們有強大的空間感幫助我們進行判斷，但換作是電腦就必須用一些數學的技巧，對於不同的情況做各自的判斷，才能讓電腦正確回答兩條線段的相交情形。
除此之外，在寫題過程中，使用到ggb進行輔助，讓我可以對程式的執行過程有更是覺化的概念，也幫助我在解題時能更理解解題的策略！上面一題三角形個數的判斷，就使用了ggb判定將點搬移的所有情況。利用它我抓到了當點的y座標為零時並沒進行好特殊情況的判斷，這也是一個視覺化之後的好處！
有一題沒有做的是模擬退火的題目實作，要求圓與三角形的交集面積，感覺超級複雜，以後有時間來慢慢實作！
Sign in with Wallet
Connect another wallet