fork download
copy 
  1. /**
  2.  * Solution to WF/2017/A/Airport.
  3.  * Find the longest line contained inside a polygon, allowing for lines touching
  4.  * vertices or the boundary.
  5.  *
  6.  * This solution uses Per's solution sketch approach to determine for any pair 
  7.  * of vertices (a, b), 
  8.  * (a) determine if a--b crosses the polygon boundary
  9.  * (b) extend the half lines outward until the polygon boundary is crossed.
  10.  *
  11.  * The results of certain boolean predicates must be computed exactly, 
  12.  * particularly the test whether the intersection point lies on the line
  13.  * or does not.
  14.  *
  15.  */
  16. #include <bits/stdc++.h>
  17. using namespace std;
  18.  
  19. //typedef long coord_t; // overflows
  20. //typedef __int128 coord_t;   // passes
  21. //typedef long double coord_t;   // passes
  22. typedef double coord_t; // passes on Kattis, too.
  23.  
  24. typedef std::tuple<coord_t, coord_t, coord_t> Hom;
  25.  
  26. Hom makePoint(coord_t x, coord_t y) {
  27.     return Hom { x, y, 1 };
  28. }
  29.  
  30. static coord_t 
  31. dot(const Hom &v, const Hom &w) {
  32.     coord_t v0, v1, v2, w0, w1, w2;
  33.     tie (v0, v1, v2) = v;
  34.     tie (w0, w1, w2) = w;
  35.     return v0*w0 + v1*w1 + v2*w2;
  36. }
  37.  
  38. static Hom 
  39. crossproduct(const Hom &u, const Hom &v) {
  40.     coord_t v0, v1, v2, u0, u1, u2;
  41.     tie (u0, u1, u2) = u;
  42.     tie (v0, v1, v2) = v;
  43.     return Hom { u1*v2 - u2*v1, u2*v0 - u0*v2, u0*v1 - u1*v0 };
  44. }
  45.  
  46. static Hom connect(const Hom &p0, const Hom &p1) { return crossproduct(p0, p1); }
  47.  
  48. static Hom intersect(const Hom &p0, const Hom &p1) { return crossproduct(p0, p1); }
  49.  
  50. static coord_t
  51. signum(coord_t x) {
  52.     return x > 0 ? 1 : x == 0 ? 0 : -1;
  53.     // return (x > 0) - (x < 0);
  54. }
  55.  
  56. static coord_t 
  57. isCCW(const Hom &p0, const Hom &p1, const Hom &p2) {
  58.     const Hom l0 = connect(p0, p1);
  59.     return dot(l0, p2);
  60. }
  61.  
  62. static pair<double, double> 
  63. getEuclidean(Hom &p) {
  64.     coord_t x, y, w;
  65.     tie (x, y, w) = p;
  66.     return make_pair(x/(double)w, y/(double)w);
  67. }
  68.  
  69. static coord_t
  70. isleft(const Hom &p0, const Hom &p1) {
  71.     coord_t x0, y0, w0, x1, y1, w1;
  72.     tie (x0, y0, w0) = p0;
  73.     tie (x1, y1, w1) = p1;
  74.     // see Ghali 16.6
  75.     return ((x0 * w1 - w0 * x1) * signum(w0 * w1));
  76. }
  77.  
  78. static coord_t
  79. isbelow(const Hom &p0, const Hom &p1) {
  80.     coord_t x0, y0, w0, x1, y1, w1;
  81.     tie (x0, y0, w0) = p0;
  82.     tie (x1, y1, w1) = p1;
  83.     return ((y0 * w1 - w0 * y1) * signum(w0 * w1));
  84. }
  85.  
  86. static coord_t
  87. isleftorbelow(const Hom &p0, const Hom &p1) {
  88.     coord_t c = isleft(p0, p1);
  89.     if (c == 0)
  90.         c = isbelow(p0, p1);
  91.     return signum(c);
  92. }
  93.  
  94. // is p in the rectangle spanned by q0-q1?
  95. static bool
  96. inBounds(const Hom &q0, const Hom& q1, const Hom &p)
  97. {
  98.     coord_t a = isleft(q0, p);
  99.     coord_t b = isleft(p, q1);
  100.     if (a * b < 0)
  101.         return false;
  102.  
  103.     a = isbelow(q0, p);
  104.     b = isbelow(p, q1);
  105.     return a * b >= 0;
  106. }
  107.  
  108. static bool
  109. pointOnSegment(const Hom& p0, const Hom& p1, const Hom& q) {
  110.     return inBounds(p0, p1, q) && isCCW(p0, p1, q) == 0;
  111. }
  112.  
  113. // standard winding number polygon contains
  114. static bool
  115. contains_wn(const vector<Hom> &points, const Hom &q) {
  116.     int wn = 0;
  117.  
  118.     const int n = points.size();
  119.     for (int i = 0; i < n; i++) {
  120.         const Hom &p = points[i];
  121.         const Hom &pn = points[(i+1)%n];
  122.         if (isbelow(p, q) <= 0) {
  123.             if (isbelow(pn, q) > 0)   
  124.                 if (isCCW(p, pn, q) > 0)
  125.                      wn++;
  126.         } else {
  127.             if (isbelow(pn, q) <= 0)
  128.                 if (isCCW(p, pn, q) < 0)
  129.                      wn--;
  130.         }
  131.     }
  132.  
  133.     return wn != 0;
  134. }
  135.  
  136. // contains: contains or on boundary.
  137. static bool
  138. contains(const vector<Hom> &points, Hom &q) {
  139.     if (contains_wn(points, q))
  140.         return true;
  141.  
  142.     const int n = points.size();
  143.     for (int i = 0; i < n; i++) {
  144.         const Hom &p = points[i];
  145.         const Hom &pn = points[(i+1)%n];
  146.         if (pointOnSegment(p, pn, q))
  147.             return true;
  148.     }
  149.     return false;
  150. }
  151.  
  152. enum IntersectType { INTERSECT_AT_VERTEX, INTERSECT_AT_POLYGON_BOUNDARY };
  153.  
  154. static Hom
  155. midpoint(const Hom &p0, const Hom &p1)
  156. {
  157.     coord_t x0, y0, w0, x1, y1, w1;
  158.     tie (x0, y0, w0) = p0;
  159.     tie (x1, y1, w1) = p1;
  160.     return Hom { x0*w1+x1*w0, y0*w1+y1*w0, 2*w0*w1 };
  161. }
  162.  
  163. struct sort_by_left_or_below
  164. {
  165.     inline bool operator() (const pair<IntersectType, Hom>& a, const pair<IntersectType, Hom>& b)
  166.     {
  167.         return isleftorbelow(get<1>(a), get<1>(b)) < 0;
  168.     }
  169. };
  170.  
  171. static Hom
  172. extend_half_line(vector<Hom> &p, vector<pair<IntersectType, Hom>> &relevantPoints, int i, int d, int bound) {
  173.     while (i + d != bound) {
  174.         IntersectType kind1, kind2;
  175.         Hom p1, p2;
  176.         tie(kind1, p1) = relevantPoints[i];
  177.         tie(kind2, p2) = relevantPoints[i+d];
  178.  
  179.         Hom mp = midpoint(p1, p2);
  180.         if (!contains(p, mp))
  181.             return p1;
  182.  
  183.         if (kind2 != INTERSECT_AT_VERTEX)
  184.             return p2;
  185.  
  186.         i += d;
  187.     }
  188.  
  189.     return get<1>(relevantPoints[i]);
  190. }
  191.  
  192. static double
  193. landingstrip(vector<Hom> &p, int a, int b)
  194. {
  195.     const int n = p.size();
  196.     vector<pair<IntersectType, Hom>> relevantPoints;
  197.  
  198.     for (int i = 0; i < n; i++) {
  199.         if (pointOnSegment(p[a], p[b], p[i])) {
  200.             relevantPoints.push_back(make_pair(INTERSECT_AT_VERTEX, p[i]));
  201.         } else
  202.         if (pointOnSegment(p[a], p[b], p[(i+1)%n])) {
  203.             relevantPoints.push_back(make_pair(INTERSECT_AT_VERTEX, p[(i+1)%n]));
  204.         } else {
  205.             Hom ls = connect(p[i], p[(i+1)%n]); // 42 bits
  206.             Hom ab = connect(p[a], p[b]);       // 42 bits
  207.             Hom is = intersect(ab, ls);         // 84 bits
  208.             coord_t x, y, w;
  209.             tie(x, y, w) = is;
  210.             if (w == 0) {   // parallel or coinciding
  211.                 if (x == 0 && y == 0) { // coinciding
  212.                     relevantPoints.push_back(make_pair(INTERSECT_AT_VERTEX, p[i]));
  213.                     relevantPoints.push_back(make_pair(INTERSECT_AT_VERTEX, p[(i+1)%n]));
  214.                 }   // ignore if parallel
  215.             } else {
  216.                 if (inBounds(p[i], p[(i+1)%n], is)) {
  217.                     relevantPoints.push_back(make_pair(INTERSECT_AT_POLYGON_BOUNDARY, is));
  218.                 }
  219.             }
  220.         }
  221.     }
  222.  
  223.     std::sort(relevantPoints.begin(), relevantPoints.end(), sort_by_left_or_below());
  224.  
  225.     int ai = std::distance(relevantPoints.begin(), std::find(relevantPoints.begin(), relevantPoints.end(), make_pair(INTERSECT_AT_VERTEX, p[a])));
  226.     int bi = std::distance(relevantPoints.begin(), std::find(relevantPoints.begin(), relevantPoints.end(), make_pair(INTERSECT_AT_VERTEX, p[b])));
  227.     int di = signum(bi-ai);
  228.  
  229.     // check that all points in [a, b] are inside polygon
  230.     for (int i = ai; i != bi; i += di) {
  231.         enum IntersectType kind;
  232.         Hom rp;
  233.         tie(kind, rp) = relevantPoints[i+di];
  234.  
  235.         if (kind != INTERSECT_AT_VERTEX)
  236.             return -1;
  237.  
  238.         Hom mp = midpoint(get<1>(relevantPoints[i]), rp);   // 42 bits since rp and relevantPoints[i] are vertices
  239.         if (!contains(p, mp))
  240.             return -1;
  241.     }
  242.  
  243.     Hom pmin = extend_half_line(p, relevantPoints, ai, -di, -di == -1 ? -1 : relevantPoints.size());
  244.     Hom pmax = extend_half_line(p, relevantPoints, bi,  di,  di == -1 ? -1 : relevantPoints.size());
  245.  
  246.     // now report distance as double
  247.     double xmin, ymin, xmax, ymax;
  248.     tie(xmin, ymin) = getEuclidean(pmin);
  249.     tie(xmax, ymax) = getEuclidean(pmax);
  250.  
  251.     return sqrt((double)(xmax-xmin)*(xmax-xmin) + (double)(ymax-ymin)*(ymax-ymin));
  252. }
  253.  
  254. int main(int ac, char *av[])
  255. {
  256.     int n;
  257.     cin >> n;
  258.     vector<Hom> p;
  259.     for (int i = 0; i < n; i++) {
  260.         long x, y;
  261.         cin >> x >> y;
  262.         p.push_back(makePoint(x, y));
  263.     }
  264.     double answer = -1.0;
  265.     for (int a = 0; a < n; a++)
  266.         for (int b = a+1; b < n; b++)
  267.             answer = max(answer, landingstrip(p, a, b));
  268.  
  269.     cout.precision(17);
  270.     cout << answer << endl;
  271. }
  272.  
Success #stdin #stdout 0s 15248KB
comments (?)
stdin
copy 
7
0 20
40 0
40 20
70 50
50 70
30 50
0 50
stdout
copy 
76.157731058639087
https://ideone.com/zy3uG9
language:
C++ (gcc 8.3)
created:
7 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!