## 最近點對問題 Closest Pair of Points
reference: APCS進階程式設計課程 課本 ==6-2 分⽽治之==

---
緣起:
> 有位同學問了，
> 老師，為什麼我寫了時間 ==複雜度較低== 的程式，但 APCS 跑出來的 ==時間== 卻更久
> :star: 感謝 `20803 江旻達` 同學 2024/06 提供的 [非常好問題](https://zh.wikipedia.org/zh-tw/早上好中国，现在我有冰淇淋#:~:text=激情9，因為-,非常好電影,-，動作非常好) 與兩個不同方法的 `程式碼` :star:

[TOC]
### 題目
10000 個二維座標點 $(0 \le x, y \le 40000)$ 找最近的兩個點的 ==距離==

### 簡述策略
首先這題分而治之的大方向是:

// 1. 從 x 座標分割成左右兩個半邊
mid = (left + right) / 2;

// 2. 取得兩半邊各自的最近點對 "距離" 後
min_left  = get_min_distance(arr, left, mid);
min_right = get_min_distance(arr, mid, right);

// 3. 針對交錯的點對計算最短距離 (一個點在左，一個點在右的 pair)
// 3.1 根據推論(請查看課本)，發現一個好特性，
// 只需要抓分割線左右邊 - 距離 d 的點來判斷就好了
d = min(min_left, min_right);

// 3.2 找交錯 pair 的最小距離
// 在 X 座標於 [mid - d, mid + d] 這個區間中的點集中起來
// 找最小距離

他寫的兩個版本（如同課本介紹的順序）分別是
1. 將中間的所有點兩兩比較 $O(n^2)$
2. 將中間的所有點 沿著 y 軸排序 $O(n \log n)$，再從上到下檢查：每個點，與後七個點的距離$O(n)$

:::spoiler 1.cpp

// old method
#include <iostream>
#include<cmath>
#include<algorithm>
#include<iomanip>
#include <chrono>
using namespace std;
struct P{
    double x; double y;
};
bool cmpx(P a,P b){
    return a.x<b.x;
}
bool cmpy(P a,P b){
    return a.y<b.y;
}
double d(P a,P b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
double closestpair(P arr[],int h,int l){
    int n=h-l+1;
    if(n==1)return(-1);
    if(n==2)return(d(arr[l],arr[h]));
    if(n==3)return(min(d(arr[l],arr[h]),min(d(arr[l+1],arr[h]),d(arr[l],arr[l+1]))));

    int mid=(h+l)/2;
    double dl=closestpair(arr,mid,l);
    double dr=closestpair(arr,h,mid+1);
    double ans=min(dl,dr);

    int midl=mid;
    int midh=mid;
    int midx=arr[mid].x;
    for(int i=mid;i>=l;i--){
        if(arr[i].x<midx-ans) break;
        midl=i;
    }
    for(int i=mid;i<=h;i++){
        if(arr[i].x>midx+ans) break;
        midh=i;
    }
    for(int i=midl;i<=midh;i++){
        for(int j=i+1;j<=midh;j++){
            ans=min(ans,d(arr[i],arr[j]));
        }
    }
    // int k=0;
    // P arry[n];
    // for(int i=l;i<=h;i++){
    //     if(arr[i].x>=arr[mid].x-ans&&arr[i].x<=arr[mid].x+ans){
    //         arry[k].x=arr[i].x;
    //         arry[k].y=arr[i].y; 
    //         k++;
    //     }  
    // }
    // sort(arry,arry+k,cmpy);
    // for(int i=0;i<k;i++){
    //     for(int j=0;j<7;j++){
    //         if(i+j+1<k)ans=min(ans,d(arry[i],arry[i+j+1]));
    //     }
    // }
    // k=0;
    return(ans);
}
int main(){
    auto start = std::chrono::high_resolution_clock::now();

    P arr[10001];
    int N; 
    while(cin>>N){
        if(N==0) break;
        else{
            for(int i=0;i<N;i++){
                cin>>arr[i].x>>arr[i].y;
            }
        }
        sort(arr,arr+N,cmpx);
        double d=closestpair(arr,N-1,0);
        if(d>10000||d==-1){
            cout<<"INFINITY "<<endl;
        }
        //else cout << setiosflags (ios::fixed) << setprecision (4)<< floor (d * 10000)/10000<< endl;
        else cout<<fixed<<setprecision(4)<<d<<endl;
    }
    auto end = std::chrono::high_resolution_clock::now();

    // 計算時間差
    std::chrono::duration<double> elapsed = end - start;

    // 顯示經過的時間（秒）
    std::cout << "Elapsed time: " << elapsed.count() << " seconds\n";
}

:::

:::spoiler 2.cpp

// new method
#include <iostream>
#include<cmath>
#include<algorithm>
#include<iomanip>
#include <chrono>
using namespace std;
struct P{
    double x; double y;
};
bool cmpx(P a,P b){
    return a.x<b.x;
}
bool cmpy(P a,P b){
    return a.y<b.y;
}
double d(P a,P b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
double closestpair(P arr[],int h,int l){
    int n=h-l+1;
    if(n==1)return(-1);
    if(n==2)return(d(arr[l],arr[h]));
    if(n==3)return(min(d(arr[l],arr[h]),min(d(arr[l+1],arr[h]),d(arr[l],arr[l+1]))));

    int mid=(h+l)/2;
    double dl=closestpair(arr,mid,l);
    double dr=closestpair(arr,h,mid+1);
    double ans=min(dl,dr);

    // int midl=mid;
    // int midh=mid;
    // int midx=arr[mid].x;
    // for(int i=mid;i>=l;i--){
    //     if(arr[i].x<midx-ans) break;
    //     midl=i;
    // }
    // for(int i=mid;i<=h;i++){
    //     if(arr[i].x>midx+ans) break;
    //     midh=i;
    // }
    // for(int i=midl;i<=midh;i++){
    //     for(int j=i+1;j<=midh;j++){
    //         ans=min(ans,d(arr[i],arr[j]));
    //     }
    // }
    int k=0;
    P arry[n];
    for(int i=l;i<=h;i++){
        if(arr[i].x>=arr[mid].x-ans&&arr[i].x<=arr[mid].x+ans){
            arry[k].x=arr[i].x;
            arry[k].y=arr[i].y; 
            k++;
        }  
    }
    sort(arry,arry+k,cmpy);
    for(int i=0;i<k;i++){
        for(int j=0;j<7;j++){
            if(i+j+1<k)ans=min(ans,d(arry[i],arry[i+j+1]));
        }
    }
    k=0;
    return(ans);
}
int main(){
    auto start = std::chrono::high_resolution_clock::now();

    P arr[10001];
    int N; bool con=true;
    while(cin>>N){
        if(N==0) break;
        else{
            for(int i=0;i<N;i++){
                cin>>arr[i].x>>arr[i].y;
                if(arr[i].x>40000||arr[i].y>40000||arr[i].x<0||arr[i].y<0)con =false;
            }
        }
        sort(arr,arr+N,cmpx);
        double d=closestpair(arr,N-1,0);
        if(d>10000||d==-1||con==false){
            cout<<"INFINITY "<<endl;
            con=true;
        }
        //else cout << setiosflags (ios::fixed) << setprecision (4)<< floor (d * 10000)/10000<< endl;
        else cout<<fixed<<setprecision(4)<<d<<endl;
    }
    auto end = std::chrono::high_resolution_clock::now();

    // 計算時間差
    std::chrono::duration<double> elapsed = end - start;

    // 顯示經過的時間（秒）
    std::cout << "Elapsed time: " << elapsed.count() << " seconds\n";
}

:::
p.s. 我在他的 main 頭尾加上了計時的程式碼 `auto start = ...`

---
理論上來說，方法 2 會比較好
但送到 APCScalss 後，反而 ==方法 1 的執行時間更短==，多次測試都得到一樣的結果

這就是問題:
:::danger
為什麼時間複雜度低的好方法，反而需要更長的時間執行?
:::

---
### 時間複雜度的意義
$O(n^2)$ 是否永遠都比 $O(n \log n)$ 還慢？
1. $10\ n \log n > 0.1\ n^2,$ ` if n < 200`
演算法為了降低時間複雜度，通常會多做一些事情，減緩最耗時的步驟，
當 n 不夠大時，這些耗時的部分就導致時間過久了
:::info
以這題來說，n=10000 已經夠大了，不是這個問題
:::
2. 特例：==平均狀況== 與 ==最壞狀況==
我們平常(或之前)討論的是 ==最壞狀況== 的時間複雜度，最久會需要多長的時間
但這個最壞狀況不見得常常發生，比如:
quicksort 的最壞狀況就是：每次選擇的切割點都是最大/最小值，導致算法會變成 $O(n^2)$

---
來到這個最近點對問題，可以先思考第二種方法是想解決什麼問題。
他想解決的問題是:
> 第一種方法會導致 $O(n^2)$ 的搜尋次數 (搜尋範圍是 `[mid-d, mid+d]` 中的點。
> 所以當範圍中的點 $n$ 很少，時間並不會太久；當範圍中的點 $n$ 很多(==最壞狀況==)，搜尋時間才會過久
:::info
那我就試著生成一個中間點很多的測資，然後分別用這兩種程式來跑就可以看出差異了 :+1:
:::

### 生成特例測資
其實生測資的程式，我也是請 AI 生的 :+1:
:::spoiler gen.cpp

/*
幫我寫一個 CPP 程式，生成最近點對問題的特例測資
10000 個點
0 <= xy <= 40000

我需要的特例測資是，面對以下算法會過久的測資:

distance = min(left_min, right_min);
min_i <- min i that point[i].x >= point[mid].x - distance
max_i <- max i that ...
for(i = min_i; i<=max_i; i++) {
  for(j=i+1; j<=max_i; j++) {
    distance = min(distance(new_distance_from(point[i], point[j]));
  }
}
*/

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstdlib>

using namespace std;

int main() {
    ofstream fout("closest_pair_hardcase.txt");

    const int N = 10000;
    const int x_start = 20000;
    const int x_range = 8; // 所有 x 在 [20000, 20008] 之間
    const int y_step = 4;   // 保持距離不太近
    const int y_start = 0;

    fout << N << "\n";
    for (int i = 0; i < N; ++i) {
        int x = x_start + (rand() % x_range); // 產生接近的 x 值
        int y = y_start + i * y_step;         // 線性增加 y 值，保證距離夠大
        fout << x << " " << y << "\n";
    }

    fout.close();
    cout << "Test case generated: closest_pair_hardcase.txt" << endl;

    return 0;
}

:::
他做的事情就是，將點沿著 y 軸平分，所以大約間隔 4 個單位 $(40000 / 10000)$
> 以確保左右半邊找到的解 (`min_left`, `min_right`) 足夠大

接著將這些點的 x 座標均勻分布在 `[mid-4. mid+4]` 的範圍
> 以確保所有點都在中間線附近 (距離 4 以內)

### 測試兩份程式碼
1. 使用上述程式碼生成測資 `closest_pair_hardcase.txt`
2. 分別使用兩份 code 來執行，
但為了避免貼上測資太慢導致的速度差異，會使用 ==輸出入轉向== 這個技術來輸入測資
(細節在後)
3. $O(n^2)$ 的方法一，明顯跑比較久!
:::success
結案 :+1:，特例的最壞狀況，確實會導致 `方法一` 速度太慢
:::

#### 細節: ==輸出入轉向==
> 使用 `codeblocks` 或 `devc++` 等軟體，執行程式後都會生成一個 `.exe` 檔案，
> 這個檔案就是你寫的程式的 ==執行檔==，先假設為 `code.exe` 好了

> 你可以在 命令提示自元 (CMD) 找到這個 `code.exe` 並執行他，後續的流程就跟平常用 `codeblock` 跳出的黑框框一樣

> 但我們要做到的是輸入轉向，不是從鍵盤輸入，而是從 `.txt` 檔案輸入

> 方法是 `code.exe < closest_pair_hardcase.txt`
> 就能將文件的內容當作輸入，放到 code.exe 中

<div style="height: 500px;"></div>