UOJ#277 【清华集训2016】定向越野

Posted by EtaoinWu on 周四 07 九月 2017

调了 1+day?

卡精度 ?

5-10-60-90-95-100?

或称最慢提交 ?

心累

#pragma GCC optimize(3)
#include "bits/stdc++.h"
//#define SHOW_CONSOLE
//#include <graphics.h>
using namespace std;
namespace io
{
#define USE_FEAD_IO
#ifdef USE_FEAD_IO
#define BUF_SIZE 5000010
  char buf[BUF_SIZE]; int cur = BUF_SIZE; FILE * in = stdin;
  inline char gnc()
  {
    if (cur == BUF_SIZE) { fread(buf, BUF_SIZE, 1, in); cur = 0; }
    return buf[cur++];
  }
#else
  inline char gnc()
  {
    return (char)getchar();
  }
#endif
  template<typename T>
  inline void gi(T &dx)
  {
    dx = 0;
    int yc = gnc();
    bool nega = false;
    while (yc < '0' || yc > '9') { nega = (yc == '-' ? true : nega); yc = gnc(); }
    while (yc >= '0' && yc <= '9') { dx = (T)(dx * 10 + yc - '0'); yc = gnc(); }
    if (nega) { dx = -dx; }
  }
  void gc(char &a)
  {
    do a = gnc(); while (!isgraph(a));
  }
  void gpc(char &a)
  {
    do a = gnc(); while (!isprint(a));
  }
  void gss(char * jc)
  {
    *jc = gnc();
    while (!isgraph(*jc)) * jc = gnc();
    while (isgraph(*jc)) * ++jc = gnc();
    *jc = 0;
  }
#if __cplusplus >= 201103l || (defined(_MSVC_LANG) && _MSVC_LANG >= 201103l)
  template<typename T, typename ...Args> inline void gi(T &a, Args &...args)
  { gi(a); gi(args...); }
#else
  template<typename t1, typename t2> inline void gi(t1 &a, t2 &b) { gi(a); gi(b); }
  template<typename t1, typename t2, typename t3> inline void gi(t1 &a, t2 &b, t3 &c) { gi(a); gi(b); gi(c); }
  template<typename t1, typename t2, typename t3, typename t4> inline void gi(t1 &a, t2 &b, t3 &c, t4 &d) { gi(a); gi(b); gi(c); gi(d); }
  template<typename t1, typename t2, typename t3, typename t4, typename t5> inline void gi(t1 &a, t2 &b, t3 &c, t4 &d, t5 &e) { gi(a); gi(b); gi(c); gi(d); gi(e); }
#endif
}
using namespace io;

const double eps = 1e-3;
const double pi = acos(-1);
inline double fuck(double x)
{
  while (x < 0)x += pi;
  while (x >= pi)x -= pi;
  return x;
}
struct Point
{
  double x, y;
  Point() :x(0), y(0) {}
  Point(double x, double y) :x(x), y(y) {}
};
typedef Point Vector;
inline Vector operator +(const Vector &A, const Vector &B) { return Vector(A.x + B.x, A.y + B.y); }
inline Vector operator -(const Point &A, const Point &B) { return Vector(A.x - B.x, A.y - B.y); }
inline Vector operator *(const Vector &A, double p) { return Vector(A.x*p, A.y*p); }
inline Vector operator /(const Vector &A, double p) { return Vector(A.x / p, A.y / p); }
inline ostream &operator<<(ostream&a, const Vector &p)
{
  return a << "(" << p.x << "," << p.y << ")";
}
//  考虑 eps
inline int sgn(double x)
{
  if (fabs(x) < eps) return 0; else return x < 0 ? -1 : 1;
}
inline bool operator<(const Point&a, const Point&b)
{
  return sgn(a.x - b.x) < 0 || (sgn(a.x - b.x) == 0 && sgn(a.y - b.y) < 0);
}
inline bool operator==(const Point&a, const Point&b)
{
  return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0;
}
// 向量极角 
inline double angle(const Vector &v) { return atan2(v.y, v.x); }
// 点积 
inline double dot(const Vector & A, const Vector & B) { return A.x*B.x + A.y*B.y; }
inline double length(const Vector & A) { return sqrt(dot(A, A)); }
inline double angle(const Vector & A, const Vector & B) { return acos(dot(A, B) / length(A) / length(B)); }
// 叉积 
inline double cross(const Vector & A, const Vector & B) { return A.x*B.y - A.y*B.x; }
inline double area2(const Point &A, const Point & B, const Point & C) { return cross(B - A, C - A); }
// 向量旋转 
inline Vector rotate(const Vector & A, double rad)
{
  return Vector(A.x*cos(rad) - A.y*sin(rad), A.x*sin(rad) + A.y*cos(rad));
}
// 计算向量的单位法线 , 左转 90 度以后把长度归一化 , 调用前确保 A 不是零向量 
inline Vector norm(const Vector & A)
{
  double L = length(A);
  return Vector(-A.y / L, A.x / L);
}
// 直线交点 
inline Point getIntersection(const Point & P, const Vector & v, const Point & Q, const Vector & w)
{
  Vector u = P - Q;
  double t = cross(w, u) / cross(v, w);
  return P + v*t;
}
// 点到直线的距离 
inline double distanceToLine(const Point & P, const Point & A, const Point & B)
{
  Vector v1 = B - A, v2 = P - A;
  return fabs(cross(v1, v2)) / length(v1);
}
// 点到线段的距离 
inline double distanceToSegment(const Point & P, const Point & A, const Point & B)
{
  if (A == B) return length(P - A);
  Vector v1 = B - A, v2 = P - A, v3 = P - B;
  if (sgn(dot(v1, v2)) < 0) return length(v2);
  else if (sgn(dot(v1, v3)) < 0) return length(v3);
  else return fabs(cross(v1, v2)) / length(v1);
}
// 线段相交判定 ( 不允许端点处相交 )
inline bool segmentIntersectionp(const Point & a1, const Point & a2, const Point & b1, const Point & b2)
{
  double c1 = cross(a2 - a1, b1 - a1), c2 = cross(a2 - a1, b2 - a1),
    c3 = cross(b2 - b1, a1 - b1), c4 = cross(b2 - b1, a2 - b1);
  return sgn(c1)*sgn(c2) < 0 && sgn(c3)*sgn(c4) < 0;
}
// 线段相交判定 ( 允许端点处相交 )
inline bool on(const Point & p, const Point & a1, const Point & a2)
{
  return sgn(cross(a1 - p, a2 - p)) == 0 && sgn(dot(a1 - p, a2 - p)) < 0;
}
// 线段相交判定 ( 允许端点处相交 )
bool segmentIntersection(const Point & a1, const Point & a2, const Point & b1, const Point & b2)
{
  double c1 = cross(a2 - a1, b1 - a1), c2 = cross(a2 - a1, b2 - a1),
    c3 = cross(b2 - b1, a1 - b1), c4 = cross(b2 - b1, a2 - b1);
  if (sgn(c1)*sgn(c2) < 0 && sgn(c3)*sgn(c4) < 0) return true;
  if (sgn(c1) == 0 && on(b1, a1, a2)) return true;
  if (sgn(c2) == 0 && on(b2, a1, a2)) return true;
  if (sgn(c3) == 0 && on(a1, b1, b2)) return true;
  if (sgn(c4) == 0 && on(a2, b1, b2)) return true;
  return false;
}
// 多边形的有向面积 
double PolygonArea(Point *p, int n)
{
  double area = 0;
  for (int i = 1; i < n; i++) {
    area += cross(p[i] - p[0], p[i + 1] - p[0]);
  }
  return area / 2;
}
struct Circle
{
  Point c;
  double r;
  Circle() {}
  Circle(const Point & c, double r) :c(c), r(r) {}
  Point point(double a) const
  {
    return Point(c.x + r*cos(a), c.y + r*sin(a));
  }
};
struct Line
{
  Point p;
  Vector v;
  Line(const Point & p, const Vector & v) :p(p), v(v) {}
  Point point(double t) const
  {
    return p + v*t;
  }
  Line move(double d) const
  {// 平移 
    return Line(p + norm(v)*d, v);
  }
};
inline bool on(const Point & p, const Line &a)
{
  return on(p, a.p, a.p + a.v);
}
inline bool on(const Point & p, const Circle &c)
{
  return !(sgn(length(p - c.c) - c.r));
}
inline double length(const Line & a)
{
  return length(a.v);
}
// 弧度转角度 
inline double lineAngleDegree(const Vector &v)
{//0<=x<180
  double ang = angle(v) * 180.0 / pi;
  while (sgn(ang)<0) ang += 360.0;
  while (sgn(ang - 180) >= 0) ang -= 180.0;
  return ang;
}
// 构造直线 
inline Line getLine(const Point & p1, const Point & p2)
{
  return Line(p1, p2 - p1);
}
// 构造直线 
inline Line getLine(double x1, double y1, double x2, double y2)
{
  return getLine(Point(x1, y1), Point(x2, y2));
}
// 求直线交点 
inline Point getIntersection(const Line & a, const Line & b)
{
  return getIntersection(a.p, a.v, b.p, b.v);
}
// 直线与圆的交点 
int getIntersection(const Line &L, const Circle & C, double &t1, double &t2, vector<Point>& sol)
{
  double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;
  double e = a*a + c*c, f = 2 * (a*b + c*d), cir = b*b + d*d - C.r*C.r;
  double delta = f*f - 4 * e*cir;
  if (sgn(delta)<0) return 0;
  if (sgn(delta) == 0) {
    t1 = t2 = -f / (2 * e); sol.push_back(L.point(t1));
    return 1;
  }
  t1 = (-f - sqrt(delta)) / (2 * e); sol.push_back(L.point(t1));
  t2 = (-f + sqrt(delta)) / (2 * e); sol.push_back(L.point(t2));
  return 2;
}
// 直线与圆求交 
int getIntersection(const Line &L, const Circle & C)
{
  double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;
  double e = a*a + c*c, f = 2 * (a*b + c*d), cir = b*b + d*d - C.r*C.r;
  double delta = f*f - 4 * e*cir;
  if (sgn(delta)<0) return 0;
  double t1, t2;
  if (sgn(delta) == 0) {
    t1 = t2 = -f / (2 * e);
    if (t1 < 0 || t1 > 1) return 0;
    return 1;
  }
  int ans = 2;
  t1 = (-f - sqrt(delta)) / (2 * e);
  t2 = (-f + sqrt(delta)) / (2 * e);
  if (t1 < 0 || t1 > 1) --ans;
  if (t2 < 0 || t2 > 1) --ans;
  return ans;
}
// 两圆相交的代码 
int getIntersection(const Circle & C1, const Circle & C2, vector<Point>& sol)
{
  double d = length(C1.c - C2.c);
  if (sgn(d) == 0) {
    if (sgn(C1.r - C2.r) == 0) return -1;// 两圆重合 
    return 0;
  }
  if (sgn(C1.r + C2.r - d) < 0) return 0;
  if (sgn(fabs(C1.r - C2.r) - d) > 0) return 0;
  double a = angle(C2.c - C1.c);
  double da = acos((C1.r*C1.r + d*d - C2.r*C2.r) / (2 * C1.r*d));
  Point p1 = C1.point(a - da), p2 = C1.point(a + da);

  sol.push_back(p1);
  if (p1 == p2) return 1;
  sol.push_back(p2);
  return 2;
}
// 过点 p 到圆 C 的切线 。v[i] 是第 i 条切线的向量 。 返回切线条数 
int getTangents(const Point & p, const Circle &C, vector<Vector>& v)
{
  Vector u = C.c - p;
  double dist = length(u);
  if (dist < C.r) return 0;
  else if (sgn(dist - C.r) == 0) {
    v.push_back(rotate(u, pi / 2));
    return 1;
  } else {
    double ang = asin(C.r / dist);
    v.push_back(rotate(u, -ang));
    v.push_back(rotate(u, +ang));
    return 2;
  }
}
//  两圆的公切线 。
//  函数返回切线的条数 。-1 表示无穷条切线 。
// a[i] 和 b[i] 分别是第 i 条切线在圆 A 和圆 B 上的切点 
int getTangents(Circle A, Circle B, Point* a, Point* b)
{
  int cnt = 0;
  if (A.r<B.r) {
    swap(A, B); swap(a, b);
  }
  double d2 = (A.c.x - B.c.x)*(A.c.x - B.c.x) + (A.c.y - B.c.y)*(A.c.y - B.c.y);
  double rdiff = A.r - B.r;
  double rsum = A.r + B.r;
  if (d2 < rdiff * rdiff) return 0;// 内含 

  double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
  if (d2 == 0 && A.r == B.r) return -1;// 无限多条切线 
  if (d2 == rdiff * rdiff) {// 内切 ,1 条切线 
    a[cnt] = A.point(base); b[cnt] = B.point(base); cnt++;
    return 1;
  }
  // 有外共切线 
  double ang = acos((A.r - B.r) / sqrt(d2));
  a[cnt] = A.point(base + ang); b[cnt] = B.point(base + ang); cnt++;
  a[cnt] = A.point(base - ang); b[cnt] = B.point(base - ang); cnt++;
  if (d2 == rsum*rsum) {// 一条内公切线 
    a[cnt] = A.point(base); b[cnt] = B.point(pi + base); cnt++;
  } else if (d2>rsum*rsum) {// 两条内公切线 
    ang = acos(rsum / sqrt(d2));
    a[cnt] = A.point(base + ang); b[cnt] = B.point(pi + base + ang); cnt++;
    a[cnt] = A.point(base - ang); b[cnt] = B.point(pi + base - ang); cnt++;
  }
  return cnt;
}
// 求三角形外接圆 
Circle circumscribedCircle(const Point & p1, const Point & p2, const Point & p3)
{
  double Bx  = p2.x - p1.x, By = p2.y - p1.y;
  double Cx  = p3.x - p1.x, Cy = p3.y - p1.y;
  double D  = 2 * (Bx*Cy - By*Cx);
  double cx  = (Cy*(Bx*Bx + By*By) - By*(Cx*Cx + Cy*Cy)) / D + p1.x;
  double cy  = (Bx*(Cx*Cx + Cy*Cy) - Cx*(Bx*Bx + By*By)) / D + p1.y;
  Point p = Point(cx, cy);
  return Circle(p, length(p1 - p));
}
// 求三角形内接圆 
Circle inscribedCircle(const Point & p1, const Point & p2, const Point & p3)
{
  double a  = length(p2 - p3);
  double b  = length(p3 - p1);
  double c  = length(p1 - p2);
  Point p  = (p1*a + p2*b + p3*c) / (a + b + c);
  return Circle(p, distanceToLine(p, p1, p2));
}
const int maxn = 666;
extern Point f[];
struct graph
{
  static const int vv = maxn * maxn * 4;
  static const int ee = vv * 6;
  int beg[vv], nxt[ee], to[ee];
  double len[ee];
  int e;
  void aa(int x, int y, double z)
  {
    ++e; nxt[e] = beg[x]; to[e] = y; len[e] = z; beg[x] = e;
  }
  void ae(int x, int y, double z)
  {
    if (x == y) return;
    aa(x, y, z); aa(y, x, z);
  }
} g;
Point f[graph::vv];
int n;
Circle cir[maxn];
Point zza[maxn], zzb[maxn];
Point s, t;
int cc;
map<Point, int> tra;
vector<Vector> zv;
int tr(const Point & x)
{
  if (tra[x]) return tra[x];
  ++cc;
  f[cc] = x;
  return tra[x] = cc;
}
bool ok(const Line &a)
{
  for (int i = 1; i <= n; ++i) {
    if (getIntersection(a, cir[i]) == 2) {
      return false;
    }
  }
  return true;
}
vector<int> v;
vector<double> cac;
vector<int> sb;
bool sbcmp(int x, int y)
{
  return cac[x] < cac[y];
}
double dis[graph::vv];
bool vis[graph::vv];
struct info
{
  int x;
  double d;
  bool operator<(const info &a) const
  {
    return d > a.d;
  }
};
void spfa(int s, int t)
{
  fill(dis, dis + graph::vv, 1e20);
  dis[s] = 0;
  //memset(vis, false, sizeof vis);
  priority_queue<info> q;
  q.push(info({ s, 0 }));
  while (!q.empty()) {
    int x = q.top().x; q.pop();
    if (x == t) return;
    if (vis[x]) continue;
    vis[x] = true;
    for (int i = g.beg[x]; i; i = g.nxt[i]) {
      int y = g.to[i];
      double w = g.len[i];
      if (dis[y] > dis[x] + w) {
        dis[y] = dis[x] + w;
        q.push(info({ y,dis[y] }));
      }
    }
  }
  return;
}
int main()
{
  gi(s.x, s.y, t.x, t.y);
  gi(n);
  for (int i = 1; i <= n; ++i) {
    gi(cir[i].c.x, cir[i].c.y, cir[i].r);
  }
  int sss = tr(s), ttt = tr(t);
  {
    if (ok(getLine(s, t))) {
      g.ae(sss, ttt, length(t - s));
    }
  }
  for (int i = 1; i <= n; ++i) {
    getTangents(s, cir[i], zv);
    for (int j = 0; j < zv.size(); ++j) {
      Line l(s, zv[j]);
      double f1, f2;
      vector<Point> sol;
      /*assert*/(getIntersection(l, cir[i], f1, f2, sol) == 1);
      if (ok(Line(s, zv[j] * f1))) {
        Point ss = s + zv[j] * f1;
        g.ae(tr(s), tr(ss), length(zv[j] * f1));
      }
    }
    zv.clear();
    getTangents(t, cir[i], zv);
    for (int j = 0; j < zv.size(); ++j) {
      Line l(t, zv[j]);
      double f1, f2;
      vector<Point> sol;
      /*assert*/(getIntersection(l, cir[i], f1, f2, sol) == 1);
      if (ok(Line(t, zv[j] * f1))) {
        Point ss = t + zv[j] * f1;
        g.ae(tr(t), tr(ss), length(zv[j] * f1));
      }
    }
    zv.clear();
  }
  for (int i = 1; i <= n; ++i) {
    for (int j = i + 1; j <= n; ++j) {
      int x = getTangents(cir[i], cir[j], zza, zzb);
      // assert x = 2 or 3 or 4
      /*assert*/(x >= 2);
      for (int k = 0; k < x; ++k) {
        int aa = tr(zza[k]), bb = tr(zzb[k]);
        Line t = getLine(zza[k], zzb[k]);
        if (!ok(t)) continue;
        g.ae(aa, bb, length(t));
      }
    }
  }
  for (int k = 1; k <= n; ++k) {
    v.clear(); cac.clear(); sb.clear();
    for (int i = 1; i <= cc; ++i) {
      if (on(f[i], cir[k])) {
        v.push_back(i);
        cac.push_back(angle(f[i] - cir[k].c));
      }
    }
    for (int i = 0; i < v.size(); ++i) sb.push_back(i);
    sort(sb.begin(), sb.end(), sbcmp);
    for (int ii = 1; ii < v.size(); ++ii) {
      int x = sb[ii - 1], y = sb[ii];
      g.ae(tr(f[v[x]]), tr(f[v[y]]), cir[k].r * fuck(cac[y] - cac[x]));
    }
    {
      int x = sb[v.size() - 1], y = sb[0];
      g.ae(tr(f[v[x]]), tr(f[v[y]]), cir[k].r * fuck(cac[y] - cac[x]));
    }
  }
  spfa(sss, ttt);
  printf("%.1lf\n", dis[ttt]);
  return 0;
}