寻找一个点，它到一条线上所有给定点的距离之和为 K

问题描述

给定平面直角坐标系上 $n$ 个点 $(x_i,y_i)$，现在要在坐标系上找到一个点 $(x,y)$，使得它到直线 $ax + by + c = 0$ 上所有给定点 $(x_i,y_i)$ 的距离之和为 $K$。

解题思路

寻找一个点到一条直线上点的距离，在坐标系中比较麻烦。我们可以将直线平移并旋转，使得直线变成 $y=0$ 的形式。平移和旋转的操作是线性的，可以直接用矩阵乘法实现。

对于世界坐标系上的点 $P(x,y)$，做这样一次平移和旋转变换后，变成了坐标系上的点 $P'(x',y')$：

$$ \begin{pmatrix} x' \ 1 \end{pmatrix}

\begin{pmatrix} 1 & 0 & -x \ 0 & 1 & -y \ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} \cos\theta & \sin\theta & 0 \ -\sin\theta & \cos\theta & 0 \ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x \ y \ 1 \end{pmatrix} $$

经过这样的变换，直线 $ax+by+c=0$ 变成了 $y=0$，原来的点 $(x_i,y_i)$ 变成了点 $(x_i',y_i')$。此时，点 $(x,y)$ 到直线 $y=0$ 上的距离为 $y$，点 $(x_i',y_i')$ 到直线上的距离也是 $y_i'$。因此，我们只需要求出点 $(x,y)$ 到点 $(x_i',y_i')$ 的距离之和即可。

假设点 $A(x_1,y_1)$ 与点 $B(x_2,y_2)$ 的距离为 $L_{AB}$，则有以下公式：

$$ L_{AB} = \sqrt{(x_2-x_1)^2 + (y_2-y_1)^2} $$

我们可以用这个公式来求点 $(x,y)$ 到点 $(x_i',y_i')$ 的距离。具体来说，我们可以把公式展开，得到：

$$ L_{P'P_i} = \sqrt{(x_i'-x')^2 + y_i'^2} = \sqrt{(x_i-x+\cos\theta y_i-\sin\theta)^2 + (\sin\theta x_i+\cos\theta y_i)^2} $$

于是，题目就转化成了在平面直角坐标系上寻找一个点 $(x,y)$，使得它到 $n$ 个给定点的距离之和等于 $K$。可以用数学方法解决这个问题，但是时间复杂度比较高。

还有一种比较高效的做法，就是先假设点 $(x,y)$ 位于直线的左侧，再二分它到直线的距离 $y$。这样可以将时间复杂度降为 $O(n \log n)$。

代码实现

class Point:
    def __init__(self, x: float, y: float):
        self.x = x
        self.y = y

    def __sub__(self, other: "Point") -> "Vector":
        return Vector(self.x - other.x, self.y - other.y)

    def dist(self, other: "Point") -> float:
        return float((self - other).abs())

    def transform(self, matrix: "Matrix") -> "Point":
        x = self.x * matrix[0][0] + self.y * matrix[1][0] + matrix[2][0]
        y = self.x * matrix[0][1] + self.y * matrix[1][1] + matrix[2][1]
        return Point(x, y)

class Vector:
    def __init__(self, x: float, y: float):
        self.x = x
        self.y = y

    def __add__(self, other: "Vector") -> "Vector":
        return Vector(self.x + other.x, self.y + other.y)

    def __sub__(self, other: "Vector") -> "Vector":
        return Vector(self.x - other.x, self.y - other.y)

    def __mul__(self, other: Union[int, float, "Vector"]) -> Union[float, "Vector"]:
        if isinstance(other, (int, float)):
            return Vector(self.x * other, self.y * other)
        elif isinstance(other, Vector):
            return self.x * other.y - self.y * other.x
        else:
            raise TypeError()

    def abs(self) -> float:
        return math.sqrt(self.x ** 2 + self.y ** 2)

class Matrix:
    def __init__(self, data: List[List[float]]):
        self.data = data

    def __mul__(self, other: "Matrix") -> "Matrix":
        m = len(self.data)  # size of A
        n = len(other.data[0])  # size of B[0]

        AB = [[0] * n for i in range(m)]
        for i in range(m):
            for j in range(n):
                for k in range(len(other.data)):
                    AB[i][j] += self.data[i][k] * other.data[k][j]

        return Matrix(AB)

    @staticmethod
    def identity(n: int) -> "Matrix":
        I = [[float(i == j) for j in range(n)] for i in range(n)]
        return Matrix(I)

    @staticmethod
    def rotate(theta: float) -> "Matrix":
        R = [[math.cos(theta), math.sin(theta), 0],
            [-math.sin(theta), math.cos(theta), 0],
            [0, 0, 1]]
        return Matrix(R)

    @staticmethod
    def translate(dx: float, dy: float) -> "Matrix":
        T = [[1, 0, dx],
            [0, 1, dy],
            [0, 0, 1]]
        return Matrix(T)

def solve(points: List[Point], a: float, b: float, c: float, K: float) -> Optional[Point]:
    eps = 1e-7
    n = len(points)
    theta = math.atan2(b, a)
    matrix = Matrix.identity(3) * Matrix.rotate(-theta) * Matrix.translate(-c / math.sqrt(a ** 2 + b ** 2), 0)
    transformed_points = [pt.transform(matrix) for pt in points]
    miny = min(pt.y for pt in transformed_points)
    maxy = max(pt.y for pt in transformed_points)
    dist_sum = sum(pt.dist(Point(0, pt.y)) for pt in transformed_points)

    if K > dist_sum:
        return None

    l, r = miny, maxy
    while r - l > eps:
        mid = (l + r) / 2
        d = sum(abs(pt.y - mid) for pt in transformed_points)
        if d <= K:
            l = mid
        else:
            r = mid

    ans = Point(0, l)
    return ans.transform(Matrix.translate(c / math.sqrt(a ** 2 + b ** 2), 0) * Matrix.rotate(theta))

参考资料

• 平面几何学入门及其在编程竞赛中的应用