Codechef 13.9 TWOROADS 数学+计算几何+拉格朗日乘数法
题目大意：
给定平面上的$n$个点，求两条直线使得每个点到这两条直线距离的最小值的平方和最小。
数据范围$1\leq{n}\leq{100}$，坐标范围的绝对值$\leq{1000}$。
算法讨论：
考虑任意两条不重合直线，则整个平面会被这两条直线的两条垂直的角平分线分为四个区域，其中两个相对的区域距离其中一条直线比较近，另外两个相对的区域距离另一条直线比较近。
其中一条角平分线将平面上的点集分成了两个部分，考虑两部分分别形成的凸包的公切线必定是两部分各选出一个点的连线，并且这条直线等价于一开始的角平分线，因此我们可以枚举任意两个点的连线作为第一条角平分线，我们将所有点按照这个点在第一条角平分线上的投影长度排序，这样第二条角平分线对于平面起到的划分效果只有$O(n)$种，于是一共只有$O(n^3)$种状态，我们依次枚举取最优解即可。
关于如何求最优解，我们实际上需要考虑这样一个问题，给定平面上若干个点，求一条直线使得这条直线到所有点距离的平方和最小，这个问题可以参考http://www.zhihu.com/question/37942571我的回答。这个公式是支持$O(1)$时间内在点集里加点和删点的，于是我们枚举第一条角平分线，再枚举第二条角平分线位置的同时就能顺便维护两个区域的答案了。
时空复杂度：
时间复杂度$O(n^3)$，空间复杂度$O(n)$。
代码：
#include <cstdio>
#include <cstring>
#include <cctype>
#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long double ldb;
static const ldb eps = 1e-8;
int dcmp(const ldb &x) {
if (fabs(x) < eps)
return 0;
else {
return x > 0 ? 1 : -1;
}
}
#define N 110
struct Point {
ldb x, y;
Point() {}
Point(ldb _x, ldb _y) : x(_x), y(_y) {}
void read() {
cin >> x >> y;
}
friend Point operator - (const Point &a, const Point &b) {
return Point(b.x - a.x, b.y - a.y);
}
friend ldb dot(const Point &a, const Point &b) {
return a.x * b.x + a.y * b.y;
}
friend ldb cross(const Point &a, const Point &b) {
return a.x * b.y - a.y * b.x;
}
friend ldb dist(const Point &a, const Point &b) {
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
friend ldb dist2(const Point &a, const Point &b) {
return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}
}P[N];
struct Solver {
ldb x, y, x2, y2, xy;
int num;
void reset() {
x = y = x2 = y2 = xy = num = 0;
}
void add(Point p, int type) {
x += type * p.x;
y += type * p.y;
x2 += type * (p.x * p.x);
y2 += type * (p.y * p.y);
xy += type * (p.x * p.y);
num += type;
}
ldb solve() {
static ldb A, B, C, delta;
if (num == 0)
return 0;
A = x2 - x * x / num;
B = (x * y / num - xy) * 2;
C = y2 - y * y / num;
delta = 16 * (A + C) * (A + C) - 4 * (-4) * (B * B - 4 * A * C);
return (sqrt(delta) - 4 * (A + C)) / -8;
}
}ldown_rup, lup_rdown;
ldb calc() {
return ldown_rup.solve() + lup_rdown.solve();
}
struct Note {
Point p;
int type;
ldb dis;
Note() {}
friend bool operator < (const Note &a, const Note &b) {
return a.dis < b.dis;
}
}A[N];
int id;
int main() {
#ifndef ONLINE_JUDGE
freopen("tt.in", "r", stdin);
#endif
int n;
scanf("%d", &n);
int i, j, k;
for (i = 1; i <= n; ++i)
P[i].read();
ldb ans = 1e10, h, dis;
for (i = 1; i <= n; ++i) {
for (j = i + 1; j <= n; ++j) {
ldown_rup.reset();
lup_rdown.reset();
id = 0;
for (k = 1; k <= n; ++k) {
++id;
if (k == i) {
A[id].p = P[i];
A[id].type = 1;
A[id].dis = 0;
lup_rdown.add(P[i], 1);
}
else if (k == j) {
A[id].p = P[j];
A[id].type = -1;
A[id].dis = dist(P[i], P[j]);
ldown_rup.add(P[j], 1);
}
else {
h = fabs(cross(P[i] - P[j], P[i] - P[k])) / dist(P[i], P[j]);
dis = sqrt(dist2(P[i], P[k]) - h * h);
A[id].p = P[k];
A[id].type = dcmp(cross(P[i] - P[j], P[i] - P[k]));
A[id].dis = dis * dcmp(dot(P[i] - P[j], P[i] - P[k]));
if (A[id].type == 1)
lup_rdown.add(P[k], 1);
else
ldown_rup.add(P[k], 1);
}
}
sort(A + 1, A + id + 1);
ans = min(ans, calc());
for (k = 1; k <= id; ++k) {
if (A[k].type == 1) {
lup_rdown.add(A[k].p, -1);
ldown_rup.add(A[k].p, 1);
}
else {
ldown_rup.add(A[k].p, -1);
lup_rdown.add(A[k].p, 1);
}
ans = min(ans, calc());
}
}
}
printf("%.10lf\n", (double)ans / n);
return 0;
}
评论 (0)