CTU Open 2016 C Cable Connection gym101505
题意:
n个点,坐标 ∈ [ 1 , 10000 ] \in [1,10000] ∈[1,10000],问不穿过(原文:run between)两个点 (包括原点)的直线被坐标轴截得的最短距离是多少。
题解
分析:
题目意思就是求 第一象限中包括原点的凸包的切线被坐标轴截得的线段的最小值。
这个线段有两个性质
- 最优线段必过凸包上一点。
-
过 ( x 0 , y 0 ) (x_0,y_0) (x0,y0)点的直线被坐标轴截得的距离 l l l 为:(其中k为斜率)
l = ( y 0 − k x 0 ) 2 + ( x 0 − y 0 k ) 2 l=(y_0-kx_0)^2+(x_0-\frac{y_0}{k})^2 l=(y0−kx0)2+(x0−ky0)2
对上式求导,并令 k 0 = y 0 x 0 k_0=\frac{y_0}{x_0} k0=x0y0得到
l ′ = k 0 2 − k 4 k 0 − k + k 0 l'=k_0^2-\frac{k^4}{k_0}-k+k_0 l′=k02−k0k4−k+k0
k = k 3 0 k=\sqrt[3]k_0 k=3k0为一个极值。上面的公式比赛的时候不会解,可以考虑基本不等式:该结论的初等数学证明论文。
实现
这里有三种实现。最后一个超级短,但是不会证明。
-
解析几何法。
对于凸包上的每个点,用上面的求导结论O(1)求出最值。不过这个最值可能取不到,因为取到最值时的直线可能穿过凸包。对于这种情况,注意到函数的单调性,我们直接取其相邻两点和它连成的两条直线中的更优的解即可。
实现时,我们直接把合法的线段极小值与过相邻两点的线段长度取个min即可。
/**
* CTU Open 2016
* Problem Solution: Cable
* @author Mimino
*/
#include <cstdio>
#include <cmath>
#include <vector>
#include <algorithm>
using namespace std;
const int MAXN = 1000000;
const double EPS = 1e-9;
int N;
struct Vector {
double x, y;
Vector(double _x=0, double _y=0): x(_x), y(_y) {}
double size() const {
return sqrt(x*x+y*y);
}
bool operator<(const Vector & v) const {
return x < v.x || (x == v.x && y < v.y);
}
} ps[MAXN];
Vector operator-(const Vector & v1, const Vector & v2) {
return Vector(v1.x - v2.x, v1.y - v2.y);
}
double cross_product(const Vector & v1, const Vector & v2) {
return v1.x * v2.y - v1.y * v2.x;
}
double direction(const Vector & segment1, const Vector & segment2, const Vector & point) {
return cross_product(point - segment1, segment2 - segment1);
}
int main() {
while (scanf("%d", &N) == 1) {
for (int i = 0; i < N; ++i)
scanf("%lf%lf", &ps[i].x, &ps[i].y);
// convex hull
sort(ps, ps+N);
vector<Vector> hull;
for (int i = 0; i < N; ++i) {
while (hull.size() >= 2 && direction(hull[hull.size()-2], hull[hull.size()-1], ps[i]) <= 0)
hull.pop_back();
hull.push_back(ps[i]);
}
double result = 1e12;
// between two points
for (int i = 0; i+1 < hull.size(); ++i) {
double dx = hull[i+1].x - hull[i].x,
dy = hull[i+1].y - hull[i].y;
if (dx > 0 && dy < 0) {
Vector p1(0, hull[i].y - dy / dx * hull[i].x),
p2(hull[i+1].x - dx / dy * hull[i+1].y, 0);
result = min(result, (p1-p2).size());
}
}
// on a single point
for (int i = 0; i < hull.size(); ++i) {
double xx = pow(hull[i].x * hull[i].y * hull[i].y, 1.0/3.0),
yy = pow(hull[i].x * hull[i].x * hull[i].y, 1.0/3.0);
Vector p1(0, hull[i].y+yy), p2(hull[i].x+xx, 0);
if ((!i || direction(p1, p2, hull[i-1])+EPS >= 0) && (i+1 == hull.size() || direction(p1, p2, hull[i+1])+EPS >= 0))
result = min(result, (p1-p2).size());
}
printf("%.3lf\n", result);
}
return 0;
}
-
暴力分类的三分法
枚举凸包上的点的思路和上面一样,只是求极值时用三分(先要证明是一个单峰曲线,所以还是要先求导,否则只是凭感觉瞎搞)。用三分时需要一个区间,所以我们对每个点,算出它和前后两个点的连线的斜率作为三分的左右边界。这一步要分6类:
k f < 0 , k g < 0 k_f<0 ,k_g<0 kf<0,kg<0 k f > 0 , k g < 0 k_f>0 ,k_g<0 kf>0,kg<0 k f > 0 , k g > 0 k_f>0 ,k_g>0 kf>0,kg>0 k f > 0 , k g > 0 , g ⃗ ⋅ f ⃗ > = 0 k_f>0 ,k_g>0,\vec{g}\cdot \vec f>=0 kf>0,kg>0,g⋅f>=0 k f = 0 , k g = 0 k_f=0 ,k_g=0 kf=0,kg=0 k f = − ∞ , k g = − ∞ k_f=-\infty ,k_g=-\infty kf=−∞,kg=−∞ 把大于0的k变为0 或者 − ∞ -\infty −∞, 注意两个k都大于0还要判一下向量夹角,如果是锐角就要让k取 ( − ∞ , 0 ) (-\infty,0) (−∞,0)
极角排序的时候还要按照长度排序,(注意凸包板子下标,我又忘了orz)(我照着数据分了一晚上类orz)
死亡样例:
2 10000 10000 1 1 2 1 1 1 10000
#include<cstdio> #include<algorithm> #include<string> #include<iostream> #include<map> #include<vector> #include<cmath> using namespace std; typedef long double db; inline int abs(int x) { return x>0 ? x : -x; } const int maxn = 1e6 + 5; const db eps = 1e-10; int n, q, a; long double ans[maxn]; bool eq(double a, double b) { return abs(a - b) <= eps; } struct V { double x, y; V(double a = 0.0, double b = 0.0) :x(a), y(b) {} void sc() { scanf("%lf%lf", &x, &y); } double L() { return x * x + y * y; } double operator |(V const &o)const { return x * o.y - o.x * y; } double operator *(V const &o)const { return x * o.x + o.y * y; } bool operator <(V const &o)const { if (eq(x, o.x))return y < o.y; return x < o.x; } V operator -(V const &o)const { return V(x - o.x, y - o.y); } void pr() { printf("%lf %lf\n", x, y); } }st[maxn]; bool cmpr(V const &a, V const &b) { V v1 = a - st[0], v2 = b - st[0]; if (abs(v1 | v2) < eps)return v1.L() < v2.L(); return (v1 | v2) < -eps; } vector<V> ch; void getCH() { sort(st + 1, st + n+1, cmpr); ch.push_back(st[0]); for(int i=1;i<=n;i++) { while (ch.size() > 1 && (st[i] - ch.back() | ch.back() - ch[ch.size() - 2]) < eps)ch.pop_back(); ch.push_back(st[i]); } } int in(int x) { return x ? x - 1 : ch.size() - 1; } int out(int x) { return x + 1 == (int)ch.size() ? 0 : x + 1; } db f(db k, db x0, db y0) { return (y0 - k * x0)*(y0 - k * x0) + (x0 - y0 / k)*(x0 - y0 / k); } double ternary_search(double l, double r,db x0,db y0) { double eps = 1e-10; //set the error limit here while (r - l > eps) { double m1 = l + (r - l) / 3; double m2 = r - (r - l) / 3; double f1 = f(m1,x0,y0); //evaluates the function at m1 double f2 = f(m2,x0,y0); //evaluates the function at m2 if (f1 > f2) l = m1; else r = m2; } return f(l,x0,y0); //return the minimum of f(x) in [l, r] } map<string, char> mp; int main() { //freopen("cable.in", "r", stdin); //freopen("cable.out", "w", stdout); while (cin >> n) { ch.clear(); for (int i = 1; i <= n; i++) { st[i].sc(); } st[0] = { 0,0 }; getCH(); int p = 1; //for (auto t : ch)t.pr(); double ans = 1e18; if (n == 1) { ans = ternary_search(-1e60, -1e-60, st[1].x,st[1].y); printf("%.3lf\n", sqrt(ans)); continue; } while (p != 0) { V v0 = ch[in(p)]- ch[p],v1=ch[out(p)]-ch[p]; //all x are different,except most left ones,k0>k1 db k1 = (v1.x == 0 ? -1e60 : v1.y / v1.x); db k0 = (v0.x == 0 ? -1e60 : v0.y / v0.x); if ((v0 * v1) >= 0 && k1 >= 0 && k0 >= 0)k1 = -1e60, k0 = 0;// int fi2 = 1; if (k0 < -eps|| k1 < -eps) { if (k1 == -1e60&&k0 == -1e60)continue; if (k0 >= -eps)k0 = 0; if (fabs(k1) < eps)k1 = 0; if (k1 == 0 && k0 == 0)continue; if (k1 >eps)k1 = -1e60; ans=min(ans,ternary_search(k1+1e-60, k0-1e-60, ch[p].x, ch[p].y)); } p = out(p); } printf("%.3lf\n", sqrt(ans)); } }
-
优雅的三分法
猜想:凸包的切线被坐标轴截得的线段的长度是关于其截距的单峰函数。
于是我们直接对答案关于截距 x x x进三分得到如下ac代码
#include <cstdio> #include <cmath> #include <vector> #include <algorithm> #include<iostream> #include<set> using namespace std; const int MAXN = 1000001; const double EPS = 1e-7; int N,i; int x[MAXN], y[MAXN]; double f(double x0) { double y0 = 0; for (i = 1; i <= N; i++) { y0 = max(y0, (x0*y[i]) / (x0 - x[i])); } return y0 * y0 + x0 * x0; } int main() { while (~scanf("%d",&N)) { int mx = 0; for (i = 1; i <= N; i++) {scanf("%d%d",&x[i],&y[i]);mx=max(mx,x[i]); } double l = mx, r = 30000; while (r - l > EPS) { double m1 = l + (r - l) / 3; double m2 = r - (r - l) / 3; if (f(m1) < f(m2))r = m2; else l = m1; } printf("%.3lf\n", sqrt(f(l))); } return 0; }