看完闫老师的视频写的。
向量方法确实简便。

const double eps = 1e-9;
vector<vector<double>> getCenters(vector<int> p1, vector<int> p2, int r) {
    double x1 = p1[0], y1 = p1[1];
    double x2 = p2[0], y2 = p2[1];
    double xc = (x1 + x2) / 2;
    double yc = (y1 + y2) / 2;
    double d = sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
    if (d / 2 >  r + eps) {
        return {};
    }

    vector<double> v = {(x2 - x1) / d, (y2 - y1) / d};
    /*
        cos, -sin  --> 0, -1  -->  0, 1 
        sin,  cos      1,  0      -1, 0

        (x, y) -> (y, -x) and (-y, x)
    */
    double h = sqrt(r * r - (d / 2) * (d / 2));

    vector<vector<double>> res;
    res.push_back({xc + v[1] * h, yc - v[0] * h});
    res.push_back({xc - v[1] * h, yc + v[0] * h});
    return res;
}

int numPoints(vector<vector<int>>& points, int r) {
    int n = points.size();
    int res = 1;
    for (int i = 0; i < n; i++) {
        for (int j = i + 1; j < n; j++) {
            vector<vector<double>> centers = getCenters(points[i], points[j], r);
            for (auto center : centers) {
                int cnt = 0;
                double cx = center[0];
                double cy = center[1];
                for (int k = 0; k < n; k++) {
                    double x = points[k][0];
                    double y = points[k][1];
                    if ((cx - x) * (cx - x) + (cy - y) * (cy - y) < r * r + eps) {
                        cnt++;
                    }
                }
                res = max(res, cnt);
            }
        }
    }
    return res;
}