最近点对数学分析Solution

算法原理

我们以(0, 0)为坐标原点，将整个坐标系旋转某个角度。

假设此时估价函数evaluationFunction()以x坐标为例，如果我们进行无限次操作，那么一定会有作为答案的点对在x轴上的投影重合。

那么根据数学直觉，如果旋转的角度合适，最近点对的估价值不会差距太大。所以我们进行T次操作，并且每次旋转后按照evaluationFuncion()的估价值将n个点排序，并且对于每一个点向后验证checkNum个点，只要参数设置合理，就有很大概率能找到最近点对。

checkIfCorrect(limit, dataSiz)

对于检验程序，我们只需要关注checkIfCorrect(limit, dataSiz)，其中limit是实验的组数，dataSiz是每组实验的数据规模（通过GetData()函数随机生成并赋值个point数组）。

同时，我们使用以验证过的期望时间复杂度为$O(nlog(n))$的分治法进行对拍，检验结果是否正确。

Hyper Parameters

本方法涉及的超参数以及更改的方法：

1. evaluationFunction()：估价函数，第59行代码中可进行修改，目前以$x\times y$作为衡量指标
2. T：180度的划分份数，在代码第90行被设置为常量，目前为3
3. checkNum：对于每个点向后检查的点数，在代码第90行被设置为常量，目前为50

已尝试实验结果：对于上述默认超参数，设置组数为10000000，每组数据及大小为1000，结果是数学分析法完全正确。

鼓励更多人来尝试！

Code

#include <iostream>
#include <vector>
#include <limits.h>
#include <cmath>
#include <random>
#include <iostream>
#include <limits.h>
#include <stdio.h>
#define MAXN 10000000
using namespace std;
class Point {
public:
    double x, y;
    double getDistance(const Point& node) const {
        return sqrt(pow(x-node.x, 2) + pow(y-node.y, 2));
    }
    friend double getDistance(const Point& a, const Point& b) { 
        return sqrt(pow(a.x-b.x, 2) + pow(a.y-b.y, 2));
    }
} point[MAXN];
void GetData(const int& siz) { // 数据生成器
    srand(time(0));
    for(int i=0; i<siz; i++) {
        point[i].x = (unsigned long long)rand()*rand()*rand(), 
      point[i].y = (unsigned long long)rand()*rand()*rand();
    }
}  
// ---------- 分治法 ----------
Point auxiliary[MAXN]; // 辅助数组
double DivideAndConquerSolution(const Point* arr, const int& l, const int& r) {
    if(r-l==1) return getDistance(arr[l], arr[r]);
    if(r-l==2) {
        double d1 = arr[l].getDistance(arr[l+1]), d2 = arr[l].getDistance(arr[l+2]), d3 = arr[l+1].getDistance(arr[l+2]);
        return std::min(d1, std::min(d2, d3));
    }
    int mid = (l+r) >> 1, a, b;
    double d1 = DivideAndConquerSolution(arr, l, mid), d2 = DivideAndConquerSolution(arr, mid+1, r);
    double d = std::min(d1, d2);
    int index = 0;
    for(int i=mid; (i>=l)&&(arr[mid].x-arr[i].x<d); i--) auxiliary[index] = arr[i], index++;
    for(int i=mid+1; (i<=r)&&(arr[i].x-arr[mid].x<d); i++) auxiliary[index] = arr[i], index++;
    std::sort(auxiliary, auxiliary+index, [&](const Point& a, const Point& b){
        return a.y < b.y;
    });
    for(int i=0; i<index; i++) {
        for(int j=i+1; j<index; j++) {
            if(auxiliary[j].y - auxiliary[i].y>=d) break;
            d = std::min(d, getDistance(auxiliary[i], auxiliary[j]));
        }
    }
	return d; 
}

// ---------- 数学法 ----------
struct NODE {
    double a[4];
    bool operator < (const NODE &t) const { // evaluationFunction()
        return a[0]*a[1] < t.a[0]*t.a[1];
    }
} p[MAXN]; 
double MathSolution(const Point* arr, const int& len, const int& T, const int& checkNum) { // 数学法
    double z, w, x, y, x_, y_, answer = __DBL_MAX__;
    int k = std::min(checkNum, len);
    for(int i=1; i<=T; i++) {
        double ang = M_PI/2/i;
        z = sin(ang), w = cos(ang);
        for(int i=0; i<len; i++){
            x = arr[i].x, y = arr[i].y;
            p[i] = NODE({x*w-y*z, x*z+y*w, x, y});
        }
        std::stable_sort(p,p+len);
        for(int i=0; i<len; i++)
            for(int j=1; j<=k&&i+j<len; j++){
                x=p[i].a[2], y=p[i].a[3];
                x_=p[i+j].a[2], y_=p[i+j].a[3];
                z=sqrt((x-x_)*(x-x_)+(y-y_)*(y-y_));
                answer=std::min(answer, z);
            }
    }
    return answer;
}
// ----------- 验证 ----------
void checkIfCorrect(const int& limit, const int& dataSiz) { // 验证程序
		#define showGap 10000
  	// showGap用于设定隔多少组显示一次当前进度，用于观测运行情况
    vector< pair<double, double> > dif;
    for(int i=1; i<=limit; i++) {
        GetData(dataSiz);
      	const int T = 3, checkNum = 50; // !!!!!!!
        double a = MathSolution(point, dataSiz, T, checkNum);
        sort(point, point+dataSiz, [&](const Point& a, const Point& b) { // lambda function, -std=c++11
            return a.x < b.x;
        });
        double b = DivideAndConquerSolution(point, 0, dataSiz-1);
        if(a!=b) dif.push_back((pair<int, int>){a, b});
        if(i%showGap==0) printf("当前已经进行了%d组\n",i);
    }
    printf("共有%d组不同:\n", (int)dif.size());
    for(int i=0; i<dif.size(); i++) printf("%d %.5lf %.5lf\n", i+1, dif[i].first, dif[i].second);
}
int main() {
	checkIfCorrect(1000000, 1000); // example
	return 0;
}