平面内有n个点，如何快速求出距离最近的点对？【分治法】【屈婉玲教材】

题目如标题，很容易想到劈成两半，每一半各求最近的，最后merge

难点在f(n)如何估计，假设两个子问题的点集分别是left和right，两部分的最近距离分别是d1和d2，所以有可能的最小距离δ=min(d1,d2)。如果left中的每一个点都要和right中的每一个点求距离，那么f(n)=n^2，分治算法是不会简化的。幸运的是如果left中的每一个点最多只需要和right中的最多6个点比较【使用下图中的鸽巢原理证明】

所以还有T(n)=O(nlog(n))，代码如下：

import random
import math
import sys

def random_point():
    x,y = random.uniform(0,1000),random.uniform(0,1000)
    return (x,y)

def distance(a,b):
    return math.sqrt((a[0]-b[0])**2+(a[1]-b[1])**2)

def brute_force(points):
    res,a,b = sys.maxsize,None,None

    l = len(points)
    for i in range(1,l):
        for j in range(i):
            dis = distance(points[i],points[j])
            if (dis < res):
                res,a,b = dis,points[i],points[j]
    return res,a,b

def divide_algo(points):
    def divide(points):
        l = len(points)
        if (l <= 1):
            return sys.maxsize,None,None
        elif (l == 2):
            return distance(points[0],points[1]),points[0],points[1]
        mid = l >> 1
        d1,a1,b1 = divide(points[:mid])
        d2,a2,b2 = divide(points[mid:])
        d,a,b = (d1,a1,b1) if d1<d2 else (d2,a2,b2)
        line = points[mid][0]
        left,right = [],[]
        for x,y in points:
            if (line-d<=x<line):
                left.append((x,y))
            elif (line<=x<=line+d): #bug1: line<=x的情况应该划分给右边
                right.append((x,y))
        sorted_right = sorted(right, key=lambda pts: (pts[1],pts[0]))
        right_len = len(sorted_right)
        for lpx,lpy in left:
            target,ll,rr = lpy-d,0,right_len-1
            while (ll <= rr):
                mm = ll + ((rr-ll)>>1)
                if (sorted_right[mm][1] >= target):
                    rr = mm - 1
                else:
                    ll = mm + 1
            ##sorted_right[ll] is the first y >= lpy-d
            for i in range(7):
                if ll+i >= right_len:
                    break
                dis = distance((lpx,lpy),sorted_right[ll+i])
                if (dis<d):
                    d,a,b = dis,(lpx,lpy),sorted_right[ll+i]
        return d,a,b
    points.sort()
    return divide(points)

if __name__ == '__main__':
    points = [random_point() for _ in range(600)]
    #points = [(863.3684416788369, 625.7796544031397), (245.20062985148627, 620.1477779298687), (242.46118979625464, 177.28076233614954)]
    print(points)
    print('brute_force={}'.format(brute_force(points)))
    print('divide_algo={}'.format(divide_algo(points)))

 最后致敬已经在天堂的屈老师