# 高三進階程式設計期末報告-最近點對解法研究及複雜度探討

## 題目
給平片上相異$N$點$(x_1,y_1),(x_2,y_2),...,(x_N,y_n)$
請計算最近點對的距離的平方
也就是找出$i,j$，滿足$1≤i<j≤N$，且$(x_i-x_j)^2+(y_i-y_j)^2$為最小

輸入說明：
第一行包含一個正整數$N$，代表平面上的點數。
接下來$N$行，每行包含兩整樹$x_i,y_i$代表第$i$個點的座標
$1≤N≤2×10^5$
$-10^9≤x_i,y_i≤10^9$
$(x_i,y_i)≠(x_j,y_j)≠j$，不存在兩點座標相同

輸出說明：
輸出正整數，代表最近點對的距離的平方

範例輸入
4
1 3
7 5
-4 9
2 -5

範例輸出
40

## 算法一 暴力枚舉
### 想法
共n個點
枚舉$n(n-1)/2$組點對的最短距離，取最小值就是答案
### 複雜度$O(n^2)$
因為要枚舉$n(n-1)/2$組，故為$O(n^2)$
## 算法二 掃描線算法
有點像枚舉算法
枚舉前先將座標以x排序
邊向左枚舉並維護當前最小距離(d)
且當計算的兩點的x座標距離大於d就break換下一個點枚舉(若x座標距離大於最小值則不可能為新的最近點對)
可以想像成向右擴展d的範圍內找點枚舉，且當d變動擴展範圍也會跟著變動

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define x first
#define y second

int dis(pair<int,int> a, pair<int,int> b){
    return ((a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y));
}
signed main(){
    int n;
    while(cin >> n){
        vector<pair<int,int> > p;
        for(int i = 0; i < n; i++){
            int x,y;
            cin >> x >> y;
            p.push_back(make_pair(x,y));
        }
        sort(p.begin(),p.end());
        int d = 9223372036854775807;
        for(int i = 0; i < n; i++){
            for(int j = i+1; j < n; j++){
                if((p[i].x - p[j].x)*(p[i].x - p[j].x) > d)break;
                d = min(d, dis(p[i],p[j]));
            }
        }
        cout << d << '\n';
    }
}

### 時間複雜度 $O(nlogn) 最差O(n^2)$
有用到sort所以只少有$O(nlogn)$
最差情況$O(n^2)$，全部點都在同一條平行y軸的直線上(點的x都一樣)時就無法break掉任何點對
不過一般情況下還是蠻快的

## 算法三 掃描線算法(加上set對y座標做二分搜)
### 想法
做法和上一個掃描線作法很像，只是sort完後不是直接枚舉
而是將向左邊和上下延伸d的範圍內的點都放進set(以y座標比較)
左界維護方式:每次for迴圈都檢查並erase
上下界維護方式:用二分搜
然後二分搜找到區域內的上下界後再將此區域的點進行枚舉

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define x first
#define y second
int dis(pair<int,int> a, pair<int,int> b){
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
signed main(){ios::sync_with_stdio(0);cin.tie(0);
    int n;
    while(cin >> n){
        vector<pair<int,int> > p;
        for(int i = 0; i < n; i++){
            int x, y;cin >> x >> y;
            p.push_back(make_pair(x,y));
        }
        sort(p.begin(),p.end());
        set<pair<int,int> > st;
        st.insert(make_pair(p[0].y,p[0].x));
        int d = 5e18;
        for(int i = 1, l = 0; i < n;i++){
            while(l < i && (p[i].x - p[l].x)*(p[i].x - p[l].x) > d){
                st.erase(make_pair(p[l].y,p[l].x));
                l++;
            }
            auto itdown = st.lower_bound(make_pair(p[i].y - d,0));
            auto itup   = st.upper_bound(make_pair(p[i].y + d,0));
            for(auto it = itdown; it != itup; it++){
                d = min(d, dis(p[i],make_pair(it->y,it->x)));//因為存進set時是y放前面所以it->y才是x
            }
            st.insert(make_pair(p[i].y,p[i].x));
        }
        cout << d << '\n';
    }
}

### 時間複雜度 $O(nlogn)$
使用set後最差情況也被優化成$O(nlogn)$了
但這個算法和下面的要介紹的分治法比起來常數太大了
如果題目設計得比較好，是有可能被卡掉的
## 算法四 分治法
可以將平面分成兩半，以簡化問題
最近點對可能會出現在三個情況
1左半部
2右半部
3兩點橫跨左右半部
所以這題的關鍵就是3要怎麼算

我們可以用類似掃描線的方法處理狀況三
1和2情況會有一組最小距離(d)
儲存從中間向左右展開d的空間內的所有點
接者再對這些點的y排序，並對這些點枚舉，並且只要超過y座標的差超過d就break

### 實作

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define x first
#define y second

int n;
vector<pair<int,int> > p;
vector<pair<int ,int> > tmp;

bool cmp(pair<int,int> a, pair<int,int> b){
    return a.y < b.y;
}
int dis(pair<int,int> a, pair<int,int> b){
    return ((a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y));
}
#define m (l + r) / 2
int cqd(int l, int r){
    if(r - l == 1)return 9223372036854775807;
    int d = min(cqd(l,m),cqd(m,r));
    tmp.clear();
    for(int i = l; i < r; i++){
        if((p[i].x - p[m].x)*(p[i].x - p[m].x) < d){
            tmp.push_back(p[i]);
        }
    }

    sort(tmp.begin(),tmp.end(),cmp);
    for(int i = 0; i < tmp.size(); i++){
        for(int j = i + 1; j < tmp.size(); j++){
            if((tmp[i].y - tmp[j].y)*(tmp[i].y - tmp[j].y)) > d){
                break;

            }
            d = min(d, dis(tmp[i],tmp[j]));
        }
    }
    return d;

}

signed main(){
    while(cin >> n){
        p.clear();
        tmp.clear();
        int x, y;
        for(int i = 0; i < n; i++){
            cin >> x >> y;
            p.push_back(make_pair(x,y));
        }
        sort(p.begin(),p.end());
        cout << cqd(0,n) << '\n';
    }

}

### 複雜度 $O(nlog^2n)$
分治本身複雜度是$O(nlogn)$
每次還要sort $O(logn)$
而第三的狀況的處理是常數
故此算法複雜度為 $O(nlog^2n)$

:::info
:::spoiler 狀況3為常數的證明
在蒐集完左右展開的所有點，並對y排序後，開始枚舉從下向上枚舉。
枚舉中有一個很重要的動作，y座標的差超過d後就break
所以事實上需要枚舉的點就只有點$P$向上、左、右展開d的範圍

這個範圍是$d*2d$的矩形。先假設P的x在正中心。
將矩形分成兩個$d*d$的正方形，在這個正方形內所有點都要保持距離大於d的性質(因為來自同一邊遞迴)
再將這個正方形切成四個$d/2 * d/2$的正方形，每個小正方形內都只能有1個點(因為小正方形對角線長d/sqrt(2)，所以小正方形內任意兩點都小於d)
再算上p自己，最多就只能有7個點

那如果P不再中心呢
那就可能會發生小正方形內有兩個點的狀況，但就只有兩個小正方形有機會，就是覆蓋中線的那兩個小正方形。
故最多就只能有9個點
:::
## 算法五 隨機算法
### 算法
1.將全部點打亂
2.將d(當前最短距離)初始為$P_1\&P_2$的距離
3.已邊長為d/2畫網格
4.將一個一個點加入網格，並檢查新加入的點能不能更新最近點對，若成功更新最近點對，則回到第3步，重新來(前面加過的點只要重新放入網格就好，不用重新檢查了)

檢查方法如下
檢查當前網格(新加入點的所在網格)向上下左右個展開兩格的範圍內能不能更新最近點對(如圖)
只要超過這個5*5範圍外的點就不可能能更新最近點對

### 實作細節
1.利用$c++$除法取整的特性，直接將座標除以$d/2$，不用真的開空間存
2.在檢查$5*5$的空間時，可以用$unorder\_map$或$unorder\_set$的$find$功能
:::warning
$underder$的有用到$hash$，所以如果裡面有包$pair$，$C++$內建的$hash$會爛掉，要自己寫一個，或自己定義
覺得麻煩可以將點設成1e9*x+y(要注意題目最大點的範圍)
用這個方法要將全部點平移到第一象限
:::
### 實作

#include<bits/stdc++.h>
using namespace std;
#define x first
#define y second
#define int long long
int n;
vector<pair<int,int> > p;
struct pair_hash {
    template <class T1, class T2>
    std::size_t operator () (const std::pair<T1, T2>& p) const {
        auto h1 = std::hash<T1>{}(p.first);
        auto h2 = std::hash<T2>{}(p.second);
        return h1 ^ h2;
    }
};

int dx[25] = {-2,-2,-2,-2,-2,-1,-1,-1,-1,-1,0,0,0,0,0,1,1,1,1,1,2,2,2,2,2};
int dy[25] = {-2,-1,0,1,2,-2,-1,0,1,2,-2,-1,0,1,2,-2,-1,0,1,2,-2,-1,0,1,2};

int dis(pair<int,int> a, pair<int,int> b){
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

signed main(){
    while(cin >> n){
        p.clear();
        for(int i = 0; i < n; i++){
            int x, y;cin >> x >> y;
            p.push_back(make_pair(x, y));
        }
        shuffle(p.begin(),p.end(),mt19937(random_device()()));
        int d = dis(p[0],p[1]);
        long double gd = sqrt(d) / 2;

        unordered_map<pair<int,int>, pair<int,int> ,pair_hash> st;
        st.insert(make_pair(make_pair(p[0].x/gd,p[0].y/gd) ,p[0]));
        st.insert(make_pair(make_pair(p[1].x/gd,p[1].y/gd) ,p[1]));
        for(int i = 2; i < n; i++){
            bool flag = 0;
            int nx = p[i].x/gd;
            int ny = p[i].y/gd;
            for(int j = 0; j < 25; j++){
                auto it = st.find(make_pair(nx+dx[j],ny+dy[j]));
                if(it != st.end()){
                    if(dis(p[i],it->y) < d){
                        d = dis(p[i],it->y);
                        gd = sqrt(d) / 2;
                        flag = 1;
                    }
                }
            }
            if(flag){
                st.clear();
                for(int j = 0; j <= i; j++){
                    st.insert(make_pair(make_pair(p[j].x/gd,p[j].y/gd) ,p[j]));
                }
            }
            else{
                st.insert(make_pair(make_pair(p[i].x/gd,p[i].y/gd) ,p[i]));
            }


        }
        cout << d << '\n';
    }
}

### 複雜度 期望$O(n)$
這個算法會影響時間複雜度的就是更新最近點對，而什麼時候更新是隨機的，所以只能算期望複雜度
再加入第i個點時更新的機率為$\frac{2}{i}$(只有一組點為最近點對，且兩點順序可以交換)
且當第i個成功組成新的最近點對，就要花$O(i)$的時間更新網格
故每新增一個點的期望複雜度為$O(i)*\frac{2}{i}$為常數，故檢查n個點的總複雜度期望為$O(n)$
若運氣非常不好，點shuffle成每次都要更新，就會變為$O(n^2)$