hdu 5027 Help!（计算几何 三分求极值）

hdu 5027 Help!

根据题意，在陆地上的路线解分为两部分：①从起点直线到达；②无法直线到达。救到人后返回陆地的时间可以先求出来。
可以求出起点到凸多边形的相切点，在两点间的边可直线到达。找出使总时间最小的入水点，可以用三分搜索求极值（列出公式求导后可发现是一个单峰函数）。
考虑入水点无法直线到达的情况，原理与上类似

计算几何的实现还是相当蛋疼的。有几个比较烦人的地方：

1、找切点。对于每个点 p 求出 sz[p] = ( vector(p, p+1) ^
vector(p, point(Xo, Yo)） )>=0，若sz[p] != sz[p+1]，则说明点 p+1 是切点。而且，若起点就在多边形的边上，则切点为其所处的边的两端点。
2、水点无法直线到达的情况下，应先到达切点，再沿着边寻找入水点
3、对于求出的切点 a，b。最好改变其位置使 a->b 之间的点属于①，b->a间的点属于②。

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;
const int MAXN = 50010;
const double eps = 1e-4, inf = 1e100;
double Ts, Vr, Vs, rres, aans, len;
int n, sz[MAXN];
inline int cmp(double a, double b)
{
return ((a-b)>eps) - ((a-b)<-eps);
}
struct point
{
double x, y;
point() {}
point(double x, double y)
{
this->x = x;
this->y = y;
}
double operator ^ (const point & a)
{
return x*a.y-y*a.x;
}
double operator * (const point & a)
{
return x*a.x + y*a.y;
}
double dis2()
{
return (x*x + y*y);
}
double dis()
{
return sqrt(x*x + y*y);
}
point operator - (const point &a )
{
return point(x-a.x, y-a.y);
}
point operator + (const point &a )
{
return point(x+a.x, y+a.y);
}
point operator / (const double &a )
{
return point(x/a, y/a);
}
point operator * (const double a )
{
return point(x*a, y*a);
}
bool operator == (const point &a) const
{
return cmp(x,a.x)==0 && cmp(y,a.y) == 0;
}
point trunc(double r)
{
double l=dis();
if (!cmp(l,0))return *this;
return point(x*r/l,y*r/l);
}
} da[MAXN], OO, PP, frm;
double persToedge(point o, point p1, point p2)
{
if (p1 == p2)
{
return (o-p1).dis();
}
point l1 = o-p1, l2 = o-p2, l3 = p1-p2;
if (cmp(l1*l3, 0)==1)
{
return l1.dis();
}
if (cmp(l2*l3, 0)==-1)
{
return l2.dis();
}
return fabs(l1^l2)/(l3.dis());
}
double dis_Point_to_Line(point o, point a, point b)
{
return (fabs((a-o)^(b-o)))/((a-b).dis());
}
void circle_X_line(point o, double r, point a, point b, int &nc, point *res)
{
int k ;
double dis = dis_Point_to_Line(o, a, b);
k = cmp(r, dis);
if (k==-1)
{
nc = 0;
return;
}
point xx = a+((b-a)*(((b-a)*(o-a))/ ((b-a).dis2()) ));
if (k == 0)
{
nc = 1;
res[0]=xx;
return;
}
dis = r*r - dis*dis;
nc = 2;
res[0] = xx-((b-a).trunc(dis));
res[1] = xx+((b-a).trunc(dis));
if (res[0].x > res[1].x) swap(res[0], res[1]);
}
inline bool isPointOnSeg(point o, point a, point b)
{
if (a.x > b.x) swap(a, b);
return cmp(b.x, o.x)>=0 && cmp(o.x, a.x)>=0;
}
double getNum(point x, int cs = -1)
{
double tp = (PP-x).dis()/Vs;
if (cs > 0)
return tp += ((len + (x-frm).dis()) / Vr);
return (OO-x).dis()/Vr + tp;
}
double tripleSearch(point bg, point ed, int cs=-1)
{
point ml, mr;
while ((bg-ed).dis() > eps)
{
ml = (bg*2+ed)/3.0;
mr = (bg+ed*2)/3.0;
if (getNum(ml, cs) < getNum(mr, cs)) ed = mr;
else bg = ml;
}
return getNum(bg, cs);
}

void cal(int bg, int ed, int cs)
{
int nc, i = bg;
point ap[2], p1, p2;
while (true)
{
nc = 0;
p1 = da[(i+n)%n];
p2 = da[(i+1+n)%n];
if (cs == 1) frm = p1;
else frm = p2;
circle_X_line(PP, (Ts-aans)*Vs, p1, p2, nc, ap);
if (nc == 1)
{
if (isPointOnSeg(ap[0], p1, p2))
rres = min(rres, (len + (ap[0]-p1).dis())/Vr + (PP-ap[0]).dis()/Vs);
}
else if (nc == 2)
{
point e1, e2;
if (cmp(ap[0].x, max(p1.x, p2.x))>0 ||
cmp(ap[1].x, min(p1.x, p2.x))<0) continue;
if (isPointOnSeg(p1, ap[0], ap[1])) e1 = p1;
else e1 = ap[0];
if (isPointOnSeg(p2, ap[0], ap[1])) e2 = p2;
else e2 = ap[1];
rres = min(rres, tripleSearch(e1, e2, 1));
}
len += (p1-p2).dis();
if (i == ed) break;
i = (i+cs+n)%n;
}
}
void solve()
{
aans = inf;
rres = inf;
int bg, ed, nc, ssa[2], hc = 0;
point ap[2], p1, p2, ssp[2];
if (cmp( (da[1]-da[0]) ^ (da[2]-da[1]), 0) == -1)
reverse(da, da+n);
for (int i = 0; i< n; ++i)
aans = min(aans, persToedge(PP, da[i], da[(i+1)%n])/Vs*2);
if (aans > Ts)
{
puts("-1");
return;
}
for (int i = 0; i< n; ++i)
sz[i] = cmp( (da[(i+1)%n]-da[i]) ^ (OO-da[i]), 0)>=0;
for (int i = 0, j = 0; i< n; ++i)
if (sz[i] != sz[(i+1)%n])
ssa[j] = (i+1)%n, ssp[j] = da[ssa[j]], ++j, hc = 1;
if (!hc)    //
{
// rres = 0;
for (int i = 0; i< n; ++i)
if (cmp(dis_Point_to_Line(OO, da[i], da[(i+1)%n]), 0) == 0)
{
ssa[0] = i;
ssa[1] = (i+1)%n;
ssp[0] = da[i];
ssp[1] = da[(i+1)%n];
break;
}
}
else
{
if (sz[ssa[0]] != 0) swap(ssa[0], ssa[1]), swap(ssp[0], ssp[1]);
bg = ssa[0], ed = ssa[1] < ssa[0]?ssa[1]+n:ssa[1];
for (int i = bg; i<ed; ++i)
{
nc = 0;
p1 = da[i%n];
p2 = da[(i+1)%n];
if (p1.x > p2.x) swap(p1, p2);
circle_X_line(PP, (Ts-aans)*Vs, p1, p2, nc, ap);
if (nc == 0) continue;
else if (nc == 1)
{
if (isPointOnSeg(ap[0], p1, p2))
rres = min(rres, (ap[0]-OO).dis()/Vr+(ap[0]-PP).dis()/Vs);
}
else
{
point e1, e2;
if (cmp(ap[0].x, p2.x)>0 || cmp(ap[1].x, p1.x) < 0) continue;
if (isPointOnSeg(p1, ap[0], ap[1])) e1 = p1;
else e1 = ap[0];
if (isPointOnSeg(p2, ap[0], ap[1])) e2 = p2;
else e2 = ap[1];
rres = min(rres, tripleSearch(e1, e2));
}
}
}
len = (ssp[1]-OO).dis();
cal(ssa[1], (ssa[0]-1+n)%n, 1);
len = (ssp[0]-OO).dis();
cal((ssa[0]-1+n)%n, ssa[1], -1);
if (rres == inf) puts("-1");
else printf("%.2lf\n", rres + aans);
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
//freopen("out.txt", "w", stdout);
#endif // ONLINE_JUDGE
int t;
scanf("%d", &t);
while (t--)
{
scanf("%lf%lf%lf",  &Ts, &Vr, &Vs);
scanf("%lf%lf%lf%lf", &OO.x, &OO.y, &PP.x, &PP.y);
scanf("%d", &n);
for (int i = 0; i< n; ++i)
scanf("%lf%lf", &da[i].x, &da[i].y);
solve();
}
return 0;
}