zoukankan      html  css  js  c++  java
  • LightOJ 1378 The Falling Circle (Basic Geometry)

    http://www.lightoj.com/volume_showproblem.php?problem=1378

      简单的计算几何,不过做这题要看清楚题目要求,对于自由下落的部分是求直线距离即可,可是对于在斜面上的运动,速度是多少就必须留意了。卡了一个晚上,出了好几组数据,结果就只是卡在没有看清楚题意上面。早上用手机提交,通过!

    带有几何模板以及部分debug的代码:

      1 #include <cstdio>
      2 #include <cstring>
      3 #include <cmath>
      4 #include <set>
      5 #include <vector>
      6 #include <iostream>
      7 #include <algorithm>
      8 
      9 using namespace std;
     10 
     11 // Point class
     12 struct Point {
     13     double x, y;
     14     Point() {}
     15     Point(double x, double y) : x(x), y(y) {}
     16 } ;
     17 template<class T> T sqr(T x) { return x * x;}
     18 inline double ptDis(Point a, Point b) { return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));}
     19 
     20 // basic calculations
     21 typedef Point Vec;
     22 Vec operator + (Vec a, Vec b) { return Vec(a.x + b.x, a.y + b.y);}
     23 Vec operator - (Vec a, Vec b) { return Vec(a.x - b.x, a.y - b.y);}
     24 Vec operator * (Vec a, double p) { return Vec(a.x * p, a.y * p);}
     25 Vec operator / (Vec a, double p) { return Vec(a.x / p, a.y / p);}
     26 
     27 const double EPS = 1e-8;
     28 const double PI = acos(-1.0);
     29 inline int sgn(double x) { return fabs(x) < EPS ? 0 : (x < 0 ? -1 : 1);}
     30 bool operator < (Point a, Point b) { return a.x < b.x || (a.x == b.x && a.y < b.y);}
     31 bool operator == (Point a, Point b) { return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0;}
     32 
     33 inline double dotDet(Vec a, Vec b) { return a.x * b.x + a.y * b.y;}
     34 inline double crossDet(Vec a, Vec b) { return a.x * b.y - a.y * b.x;}
     35 inline double crossDet(Point o, Point a, Point b) { return crossDet(a - o, b - o);}
     36 inline double vecLen(Vec x) { return sqrt(dotDet(x, x));}
     37 inline double toRad(double deg) { return deg / 180.0 * PI;}
     38 inline double angle(Vec v) { return atan2(v.y, v.x);}
     39 inline double angle(Vec a, Vec b) { return acos(dotDet(a, b) / vecLen(a) / vecLen(b));}
     40 inline double triArea(Point a, Point b, Point c) { return fabs(crossDet(b - a, c - a));}
     41 inline Vec vecUnit(Vec x) { return x / vecLen(x);}
     42 inline Vec rotate(Vec x, double rad) { return Vec(x.x * cos(rad) - x.y * sin(rad), x.x * sin(rad) + x.y * cos(rad));}
     43 Vec normal(Vec x) {
     44     double len = vecLen(x);
     45     return Vec(- x.y / len, x.x / len);
     46 }
     47 
     48 // Line class
     49 struct Line {
     50     Point s, t;
     51     Line() {}
     52     Line(Point s, Point t) : s(s), t(t) {}
     53     Point point(double x) {
     54         return s + (t - s) * x;
     55     }
     56     Line move(double x) { // while x > 0 move to (s->t)'s left
     57         Vec nor = normal(t - s);
     58         nor = nor * x;
     59         return Line(s + nor, t + nor);
     60     }
     61     Vec vec() { return t - s;}
     62 } ;
     63 typedef Line Seg;
     64 
     65 inline bool onLine(Point x, Point a, Point b) { return sgn(crossDet(a - x, b - x)) == 0;}
     66 inline bool onLine(Point x, Line l) { return onLine(x, l.s, l.t);}
     67 inline bool onSeg(Point x, Point a, Point b) { return sgn(crossDet(a - x, b - x)) == 0 && sgn(dotDet(a - x, b - x)) < 0;}
     68 inline bool onSeg(Point x, Seg s) { return onSeg(x, s.s, s.t);}
     69 
     70 // 0 : not intersect
     71 // 1 : proper intersect
     72 // 2 : improper intersect
     73 int segIntersect(Point a, Point c, Point b, Point d) {
     74     Vec v1 = b - a, v2 = c - b, v3 = d - c, v4 = a - d;
     75     int a_bc = sgn(crossDet(v1, v2));
     76     int b_cd = sgn(crossDet(v2, v3));
     77     int c_da = sgn(crossDet(v3, v4));
     78     int d_ab = sgn(crossDet(v4, v1));
     79     if (a_bc * c_da > 0 && b_cd * d_ab > 0) return 1;
     80     if (onSeg(b, a, c) && c_da) return 2;
     81     if (onSeg(c, b, d) && d_ab) return 2;
     82     if (onSeg(d, c, a) && a_bc) return 2;
     83     if (onSeg(a, d, b) && b_cd) return 2;
     84     return 0;
     85 }
     86 inline int segIntersect(Seg a, Seg b) { return segIntersect(a.s, a.t, b.s, b.t);}
     87 
     88 // point of the intersection of 2 lines
     89 Point lineIntersect(Point P, Vec v, Point Q, Vec w) {
     90     Vec u = P - Q;
     91     double t = crossDet(w, u) / crossDet(v, w);
     92     return P + v * t;
     93 }
     94 inline Point lineIntersect(Line a, Line b) { return lineIntersect(a.s, a.t - a.s, b.s, b.t - b.s);}
     95 
     96 // Warning: This is a DIRECTED Distance!!!
     97 double pt2Line(Point x, Point a, Point b) {
     98     Vec v1 = b - a, v2 = x - a;
     99     return crossDet(v1, v2) / vecLen(v1);
    100 }
    101 inline double pt2Line(Point x, Line L) { return pt2Line(x, L.s, L.t);}
    102 
    103 double pt2Seg(Point x, Point a, Point b) {
    104     if (a == b) return vecLen(x - a);
    105     Vec v1 = b - a, v2 = x - a, v3 = x - b;
    106     if (sgn(dotDet(v1, v2)) < 0) return vecLen(v2);
    107     if (sgn(dotDet(v1, v3)) > 0) return vecLen(v3);
    108     return fabs(crossDet(v1, v2)) / vecLen(v1);
    109 }
    110 inline double pt2Seg(Point x, Seg s) { return pt2Seg(x, s.s, s.t);}
    111 
    112 Point ptOnLine(Point p, Point a, Point b) {
    113     Vec v = b - a;
    114     return a + v * (dotDet(v, p - a) / dotDet(v, v));
    115 }
    116 inline Point ptOnLine(Point p, Line x) { return ptOnLine(p, x.s, x.t);}
    117 
    118 // Polygon class
    119 struct Poly {
    120     vector<Point> pt;
    121     Poly() { pt.clear();}
    122     ~Poly() {}
    123     Poly(vector<Point> pt) : pt(pt) {}
    124     Point operator [] (int x) const { return pt[x];}
    125     int size() { return pt.size();}
    126     double area() {
    127         double ret = 0.0;
    128         int sz = pt.size();
    129         for (int i = 1; i < sz; i++) {
    130             ret += crossDet(pt[i], pt[i - 1]);
    131         }
    132         return fabs(ret / 2.0);
    133     }
    134 } ;
    135 
    136 // Circle class
    137 struct Circle {
    138     Point c;
    139     double r;
    140     Circle() {}
    141     Circle(Point c, double r) : c(c), r(r) {}
    142     Point point(double a) {
    143         return Point(c.x + cos(a) * r, c.y + sin(a) * r);
    144     }
    145 } ;
    146 
    147 inline bool ptOnCircle(Point x, Circle c) { return sgn(ptDis(c.c, x) - c.r) == 0;}
    148 
    149 // Cirlce operations
    150 int lineCircleIntersect(Line L, Circle C, double &t1, double &t2, vector<Point> &sol) {
    151     double a = L.s.x, b = L.t.x - C.c.x, c = L.s.y, d = L.t.y - C.c.y;
    152     double e = sqr(a) + sqr(c), f = 2 * (a * b + c * d), g = sqr(b) + sqr(d) - sqr(C.r);
    153     double delta = sqr(f) - 4.0 * e * g;
    154     if (sgn(delta) < 0) return 0;
    155     if (sgn(delta) == 0) {
    156         t1 = t2 = -f / (2.0 * e);
    157         sol.push_back(L.point(t1));
    158         return 1;
    159     }
    160     t1 = (-f - sqrt(delta)) / (2.0 * e);
    161     sol.push_back(L.point(t1));
    162     t2 = (-f + sqrt(delta)) / (2.0 * e);
    163     sol.push_back(L.point(t2));
    164     return 2;
    165 }
    166 
    167 int lineCircleIntersect(Line L, Circle C, vector<Point> &sol) {
    168     Vec dir = L.t - L.s, nor = normal(dir);
    169     Point mid = lineIntersect(C.c, nor, L.s, dir);
    170     double len = sqr(C.r) - sqr(ptDis(C.c, mid));
    171     if (sgn(len) < 0) return 0;
    172     if (sgn(len) == 0) {
    173         sol.push_back(mid);
    174         return 1;
    175     }
    176     Vec dis = vecUnit(dir);
    177     len = sqrt(len);
    178     sol.push_back(mid + dis * len);
    179     sol.push_back(mid - dis * len);
    180     return 2;
    181 }
    182 
    183 // -1 : coincide
    184 int circleCircleIntersect(Circle C1, Circle C2, vector<Point> &sol) {
    185     double d = vecLen(C1.c - C2.c);
    186     if (sgn(d) == 0) {
    187         if (sgn(C1.r - C2.r) == 0) {
    188             return -1;
    189         }
    190         return 0;
    191     }
    192     if (sgn(C1.r + C2.r - d) < 0) return 0;
    193     if (sgn(fabs(C1.r - C2.r) - d) > 0) return 0;
    194     double a = angle(C2.c - C1.c);
    195     double da = acos((sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2.0 * C1.r * d));
    196     Point p1 = C1.point(a - da), p2 = C1.point(a + da);
    197     sol.push_back(p1);
    198     if (p1 == p2) return 1;
    199     sol.push_back(p2);
    200     return 2;
    201 }
    202 
    203 void circleCircleIntersect(Circle C1, Circle C2, vector<double> &sol) {
    204     double d = vecLen(C1.c - C2.c);
    205     if (sgn(d) == 0) return ;
    206     if (sgn(C1.r + C2.r - d) < 0) return ;
    207     if (sgn(fabs(C1.r - C2.r) - d) > 0) return ;
    208     double a = angle(C2.c - C1.c);
    209     double da = acos((sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2.0 * C1.r * d));
    210     sol.push_back(a - da);
    211     sol.push_back(a + da);
    212 }
    213 
    214 int tangent(Point p, Circle C, vector<Vec> &sol) {
    215     Vec u = C.c - p;
    216     double dist = vecLen(u);
    217     if (dist < C.r) return 0;
    218     if (sgn(dist - C.r) == 0) {
    219         sol.push_back(rotate(u, PI / 2.0));
    220         return 1;
    221     }
    222     double ang = asin(C.r / dist);
    223     sol.push_back(rotate(u, -ang));
    224     sol.push_back(rotate(u, ang));
    225     return 2;
    226 }
    227 
    228 // ptA : points of tangency on circle A
    229 // ptB : points of tangency on circle B
    230 int tangent(Circle A, Circle B, vector<Point> &ptA, vector<Point> &ptB) {
    231     if (A.r < B.r) {
    232         swap(A, B);
    233         swap(ptA, ptB);
    234     }
    235     double d2 = sqr(A.c.x - B.c.x) + sqr(A.c.y - B.c.y);
    236     double rdiff = A.r - B.r, rsum = A.r + B.r;
    237     if (d2 < sqr(rdiff)) return 0;
    238     double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
    239     if (d2 == 0 && A.r == B.r) return -1;
    240     if (d2 == sqr(rdiff)) {
    241         ptA.push_back(A.point(base));
    242         ptB.push_back(B.point(base));
    243         return 1;
    244     }
    245     double ang = acos((A.r - B.r) / sqrt(d2));
    246     ptA.push_back(A.point(base + ang));
    247     ptB.push_back(B.point(base + ang));
    248     ptA.push_back(A.point(base - ang));
    249     ptB.push_back(B.point(base - ang));
    250     if (d2 == sqr(rsum)) {
    251         ptA.push_back(A.point(base));
    252         ptB.push_back(B.point(PI + base));
    253     } else if (d2 > sqr(rsum)) {
    254         ang = acos((A.r + B.r) / sqrt(d2));
    255         ptA.push_back(A.point(base + ang));
    256         ptB.push_back(B.point(PI + base + ang));
    257         ptA.push_back(A.point(base - ang));
    258         ptB.push_back(B.point(PI + base - ang));
    259     }
    260     return (int) ptA.size();
    261 }
    262 
    263 // -1 : onside
    264 // 0 : outside
    265 // 1 : inside
    266 int ptInPoly(Point p, Poly &poly) {
    267     int wn = 0, sz = poly.size();
    268     for (int i = 0; i < sz; i++) {
    269         if (onSeg(p, poly[i], poly[(i + 1) % sz])) return -1;
    270         int k = sgn(crossDet(poly[(i + 1) % sz] - poly[i], p - poly[i]));
    271         int d1 = sgn(poly[i].y - p.y);
    272         int d2 = sgn(poly[(i + 1) % sz].y - p.y);
    273         if (k > 0 && d1 <= 0 && d2 > 0) wn++;
    274         if (k < 0 && d2 <= 0 && d1 > 0) wn--;
    275     }
    276     if (wn != 0) return 1;
    277     return 0;
    278 }
    279 
    280 // if DO NOT need a high precision
    281 /*
    282 int ptInPoly(Point p, Poly poly) {
    283     int sz = poly.size();
    284     double ang = 0.0, tmp;
    285     for (int i = 0; i < sz; i++) {
    286         if (onSeg(p, poly[i], poly[(i + 1) % sz])) return -1;
    287         tmp = angle(poly[i] - p) - angle(poly[(i + 1) % sz] - p) + PI;
    288         ang += tmp - floor(tmp / (2.0 * PI)) * 2.0 * PI - PI;
    289     }
    290     if (sgn(ang - PI) == 0) return -1;
    291     if (sgn(ang) == 0) return 0;
    292     return 1;
    293 }
    294 */
    295 
    296 
    297 // Convex Hull algorithms
    298 // return the number of points in convex hull
    299 
    300 // andwer's algorithm
    301 // if DO NOT want the points on the side of convex hull, change all "<" into "<="
    302 int andrew(Point *pt, int n, Point *ch) {
    303     sort(pt, pt + n);
    304     int m = 0;
    305     for (int i = 0; i < n; i++) {
    306         while (m > 1 && crossDet(ch[m - 1] - ch[m - 2], pt[i] - ch[m - 2]) <= 0) m--;
    307         ch[m++] = pt[i];
    308     }
    309     int k = m;
    310     for (int i = n - 2; i >= 0; i--) {
    311         while (m > k && crossDet(ch[m - 1] - ch[m - 2], pt[i] - ch[m - 2]) <= 0) m--;
    312         ch[m++] = pt[i];
    313     }
    314     if (n > 1) m--;
    315     return m;
    316 }
    317 
    318 // graham's algorithm
    319 // if DO NOT want the points on the side of convex hull, change all "<=" into "<"
    320 Point origin;
    321 inline bool cmpAng(Point p1, Point p2) { return crossDet(origin, p1, p2) > 0;}
    322 inline bool cmpDis(Point p1, Point p2) { return ptDis(p1, origin) > ptDis(p2, origin);}
    323 
    324 void removePt(Point *pt, int &n) {
    325     int idx = 1;
    326     for (int i = 2; i < n; i++) {
    327         if (sgn(crossDet(origin, pt[i], pt[idx]))) pt[++idx] = pt[i];
    328         else if (cmpDis(pt[i], pt[idx])) pt[idx] = pt[i];
    329     }
    330     n = idx + 1;
    331 }
    332 
    333 int graham(Point *pt, int n, Point *ch) {
    334     int top = -1;
    335     for (int i = 1; i < n; i++) {
    336         if (pt[i].y < pt[0].y || (pt[i].y == pt[0].y && pt[i].x < pt[0].x)) swap(pt[i], pt[0]);
    337     }
    338     origin = pt[0];
    339     sort(pt + 1, pt + n, cmpAng);
    340     removePt(pt, n);
    341     for (int i = 0; i < n; i++) {
    342         if (i >= 2) {
    343             while (!(crossDet(ch[top - 1], pt[i], ch[top]) <= 0)) top--;
    344         }
    345         ch[++top] = pt[i];
    346     }
    347     return top + 1;
    348 }
    349 
    350 
    351 // Half Plane
    352 // The intersected area is always a convex polygon.
    353 // Directed Line class
    354 struct DLine {
    355     Point p;
    356     Vec v;
    357     double ang;
    358     DLine() {}
    359     DLine(Point p, Vec v) : p(p), v(v) { ang = atan2(v.y, v.x);}
    360     bool operator < (const DLine &L) const { return ang < L.ang;}
    361     DLine move(double x) { // while x > 0 move to v's left
    362         Vec nor = normal(v);
    363         nor = nor * x;
    364         return DLine(p + nor, v);
    365     }
    366 
    367 } ;
    368 
    369 Poly cutPoly(Poly &poly, Point a, Point b) {
    370     Poly ret = Poly();
    371     int n = poly.size();
    372     for (int i = 0; i < n; i++) {
    373         Point c = poly[i], d = poly[(i + 1) % n];
    374         if (sgn(crossDet(b - a, c - a)) >= 0) ret.pt.push_back(c);
    375         if (sgn(crossDet(b - a, c - d)) != 0) {
    376             Point ip = lineIntersect(a, b - a, c, d - c);
    377             if (onSeg(ip, c, d)) ret.pt.push_back(ip);
    378         }
    379     }
    380     return ret;
    381 }
    382 inline Poly cutPoly(Poly &poly, DLine L) { return cutPoly(poly, L.p, L.p + L.v);}
    383 
    384 inline bool onLeft(DLine L, Point p) { return crossDet(L.v, p - L.p) > 0;}
    385 Point dLineIntersect(DLine a, DLine b) {
    386     Vec u = a.p - b.p;
    387     double t = crossDet(b.v, u) / crossDet(a.v, b.v);
    388     return a.p + a.v * t;
    389 }
    390 
    391 Poly halfPlane(DLine *L, int n) {
    392     Poly ret = Poly();
    393     sort(L, L + n);
    394     int fi, la;
    395     Point *p = new Point[n];
    396     DLine *q = new DLine[n];
    397     q[fi = la = 0] = L[0];
    398     for (int i = 1; i < n; i++) {
    399         while (fi < la && !onLeft(L[i], p[la - 1])) la--;
    400         while (fi < la && !onLeft(L[i], p[fi])) fi++;
    401         q[++la] = L[i];
    402         if (fabs(crossDet(q[la].v, q[la - 1].v)) < EPS) {
    403             la--;
    404             if (onLeft(q[la], L[i].p)) q[la] = L[i];
    405         }
    406         if (fi < la) p[la - 1] = dLineIntersect(q[la - 1], q[la]);
    407     }
    408     while (fi < la && !onLeft(q[fi], p[la - 1])) la--;
    409     if (la - fi <= 1) return ret;
    410     p[la] = dLineIntersect(q[la], q[fi]);
    411     for (int i = fi; i <= la; i++) ret.pt.push_back(p[i]);
    412     return ret;
    413 }
    414 
    415 
    416 // 3D Geometry
    417 void getCoor(double R, double lat, double lng, double &x, double &y, double &z) {
    418     lat = toRad(lat);
    419     lng = toRad(lng);
    420     x = R * cos(lat) * cos(lng);
    421     y = R * cos(lat) * sin(lng);
    422     z = R * sin(lat);
    423 }
    424 
    425 struct Point3 {
    426     double x, y, z;
    427     Point3() {}
    428     Point3(double x, double y, double z) : x(x), y(y), z(z) {}
    429 } ;
    430 typedef Point3 Vec3;
    431 
    432 Vec3 operator + (Vec3 a, Vec3 b) { return Vec3(a.x + b.x, a.y + b.y, a.z + b.z);}
    433 Vec3 operator - (Vec3 a, Vec3 b) { return Vec3(a.x - b.x, a.y - b.y, a.z - b.z);}
    434 Vec3 operator * (Vec3 a, double p) { return Vec3(a.x * p, a.y * p, a.z * p);}
    435 Vec3 operator / (Vec3 a, double p) { return Vec3(a.x / p, a.y / p, a.z / p);}
    436 
    437 bool operator < (Point3 a, Point3 b) {
    438     if (a.x != b.x) return a.x < b.x;
    439     if (a.y != b.y) return a.y < b.y;
    440     return a.z < b.z;
    441 }
    442 bool operator == (Point3 a, Point3 b) { return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0 && sgn(a.z - b.z) == 0;}
    443 
    444 inline double dotDet(Vec3 a, Vec3 b) { return a.x * b.x + a.y * b.y + a.z * b.z;}
    445 inline Vec3 crossDet(Vec3 a, Vec3 b) { return Vec3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);}
    446 inline double vecLen(Vec3 x) { return sqrt(dotDet(x, x));}
    447 inline double angle(Vec3 a, Vec3 b) { return acos(dotDet(a, b) / vecLen(a) / vecLen(b));}
    448 inline double triArea(Point3 a, Point3 b, Point3 c) { return vecLen(crossDet(b - a, c - a));}
    449 
    450 struct Plane {
    451     Point3 p;
    452     Vec3 n;
    453     Plane() {}
    454     Plane(Point3 p, Vec3 n) : p(p), n(n) {}
    455 } ;
    456 
    457 // Warning: This is a DIRECTED Distance!!!
    458 inline double pt2Plane(Point3 p, Point3 p0, Vec3 n) { return dotDet(p - p0, n) / vecLen(n);}
    459 inline double pt2Plane(Point3 p, Plane P) { return pt2Plane(p, P.p, P.n);}
    460 // get projection on plane
    461 inline Point3 ptOnPlane(Point3 p, Point3 p0, Vec3 n) { return p + n * pt2Plane(p, p0, n);}
    462 inline Point3 ptOnPlane(Point3 p, Plane P) { return ptOnPlane(p, P.p, P.n);}
    463 inline bool ptInPlane(Point3 p, Point3 p0, Vec3 n) { return sgn(dotDet(p - p0, n)) == 0;}
    464 inline bool ptInPlane(Point3 p, Plane P) { return ptInPlane(p, P.p, P.n);}
    465 
    466 struct Line3 {
    467     Point3 s, t;
    468     Line3() {}
    469     Line3(Point3 s, Point3 t) : s(s), t(t) {}
    470     Vec3 vec() { return t - s;}
    471 } ;
    472 typedef Line3 Seg3;
    473 
    474 double pt2Line(Point3 p, Point3 a, Point3 b) {
    475     Vec3 v1 = b - a, v2 = p - a;
    476     return vecLen(crossDet(v1, v2)) / vecLen(v1);
    477 }
    478 inline double pt2Line(Point3 p, Line3 l) { return pt2Line(p, l.s, l.t);}
    479 
    480 double pt2Seg(Point3 p, Point3 a, Point3 b) {
    481     if (a == b) return vecLen(p - a);
    482     Vec3 v1 = b - a, v2 = p - a, v3 = p - b;
    483     if (sgn(dotDet(v1, v2)) < 0) return vecLen(v2);
    484     else if (sgn(dotDet(v1, v3)) > 0) return vecLen(v3);
    485     else return vecLen(crossDet(v1, v2)) / vecLen(v1);
    486 }
    487 inline double pt2Seg(Point3 p, Seg3 s) { return pt2Seg(p, s.s, s.t);}
    488 
    489 struct Tri {
    490     Point3 a, b, c;
    491     Tri() {}
    492     Tri(Point3 a, Point3 b, Point3 c) : a(a), b(b), c(c) {}
    493 } ;
    494 
    495 bool ptInTri(Point3 p, Point3 a, Point3 b, Point3 c) {
    496     double area1 = triArea(p, a, b);
    497     double area2 = triArea(p, b, c);
    498     double area3 = triArea(p, c, a);
    499     double area = triArea(a, b, c);
    500     return sgn(area1 = area2 + area3 - area) == 0;
    501 }
    502 inline bool ptInTri(Point3 p, Tri t) { return ptInTri(p, t.a, t.b, t.c);}
    503 
    504 // 0 : not intersect
    505 // 1 : intersect at only 1 point
    506 // 2 : line in plane
    507 int linePlaneIntersect(Point3 s, Point3 t, Point3 p0, Vec3 n, Point3 &x) {
    508     double res1 = dotDet(n, p0 - s);
    509     double res2 = dotDet(n, t - s);
    510     if (sgn(res2)) {
    511         if (ptInPlane(s, p0, n)) return 2;
    512         return 0;
    513     }
    514     x = s + (t - s) * (res1 / res2);
    515     // use general form : a * x + b * y + c * z + d = 0
    516     // Vec3 v = s - t;
    517     // double k = (a * s.x + b * s.y + c * s.z + d) / (a * v.x + b * v.y + c * v.z);
    518     // x = s + (t - s) * k;
    519     return 1;
    520 }
    521 inline int linePlaneIntersect(Line3 l, Plane p, Point3 &x) { return linePlaneIntersect(l.s, l.t, p.p, p.n, x);}
    522 
    523 // ¡÷abc intersect with segment st
    524 bool triSegIntersect(Point3 a, Point3 b, Point3 c, Point3 s, Point3 t, Point3 &x) {
    525     Vec3 n = crossDet(b - a, c - a);
    526     if (sgn(dotDet(n, t - s)) == 0) return false;
    527     else {
    528         double k = dotDet(n, a - s) / dotDet(n, t - s);
    529         if (sgn(k) < 0 || sgn(k - 1) > 0) return false;
    530         x = s + (t - s) * k;
    531         return ptInTri(x, a, b, c);
    532     }
    533 }
    534 inline bool triSegIntersect(Tri t, Line3 l, Point3 &x) { return triSegIntersect(t.a, t.b, t.c, l.s, l.t, x);}
    535 
    536 // Warning: This is a DIRECTED Volume!!!
    537 inline double tetraVol(Point3 a, Point3 b, Point3 c, Point3 p) { return dotDet(p - a, crossDet(b - a, c - a)) / 6.0;}
    538 inline double tetraVol(Tri t, Point3 p) { return tetraVol(t.a, t.b, t.c, p);}
    539 
    540 /****************** template above *******************/
    541 
    542 int myTangent(Circle A, Circle B, vector<Point> &ptA, vector<Point> &ptB) {
    543     if (A.r < B.r) {
    544         swap(A, B);
    545         swap(ptA, ptB);
    546     }
    547     double d2 = sqr(A.c.x - B.c.x) + sqr(A.c.y - B.c.y);
    548     double rdiff = A.r - B.r, rsum = A.r + B.r;
    549     if (d2 < sqr(rdiff)) return 0;
    550     double base = atan2(B.c.y - A.c.y, B.c.x - A.c.x);
    551     if (d2 == 0 && A.r == B.r) return -1;
    552     if (d2 == sqr(rdiff)) {
    553         ptA.push_back(A.point(base));
    554         ptB.push_back(B.point(base));
    555         return 1;
    556     }
    557     double ang = acos((A.r - B.r) / sqrt(d2));
    558     ptA.push_back(A.point(base + ang));
    559     ptB.push_back(B.point(base + ang));
    560     ptA.push_back(A.point(base - ang));
    561     ptB.push_back(B.point(base - ang));
    562     return 2;
    563 }
    564 
    565 bool cmpY(Point a, Point b) {
    566     return a.y < b.y;
    567 }
    568 
    569 const double FINF = 1e8;
    570 
    571 // LOJ 1378
    572 int main() {
    573 //    freopen("in", "r", stdin);
    574 //    freopen("out", "w", stdout);
    575     Circle a, b, x;
    576     int T;
    577     cin >> T;
    578     for (int cas = 1; cas <= T; cas++) {
    579         vector<Point> A, B;
    580         A.clear(), B.clear();
    581         cin >> a.c.x >> a.c.y >> a.r;
    582         cin >> b.c.x >> b.c.y >> b.r;
    583         cin >> x.c.x >> x.c.y >> x.r;
    584         if (a.c.x > b.c.x) swap(a, b);
    585         int cnt = myTangent(a, b, A, B);
    586         if (cnt != 2) { puts("Shit!"); while (1) ;}
    587         if (!ptOnCircle(A[0], a)) swap(A, B);
    588 //        cout << cnt << endl;
    589         sort(A.begin(), A.end(), cmpY);
    590         sort(B.begin(), B.end(), cmpY);
    591 //        cout << A[0].x << ' ' << A[0].y << "= =" << B[0].x << ' ' << B[0].y << endl;
    592         cout << "Case " << cas << ": ";
    593         if (sgn(A[0].y - B[0].y) == 0) {
    594             double ans = x.c.y - x.r - A[0].y;
    595             printf("%.10f\n", ans);
    596         } else {
    597             Line L;
    598             if (A[0] < B[0]) {
    599                 L = Line(A[0], B[0]);
    600             } else {
    601                 L = Line(B[0], A[0]);
    602             }
    603 //            cout << "xr " << x.r << endl;
    604             L = L.move(x.r);
    605             Line ver = Line(Point(x.c.x, -FINF), Point(x.c.x, FINF));
    606 //            cout << ver.s.x << ' ' << ver.s.y << " ver " << ver.t.x << ' ' << ver.t.y << endl;
    607 //            cout << L.s.x << ' ' << L.s.y << " L " << L.t.x << ' ' << L.t.y << endl;
    608             Point ip = lineIntersect(ver, L);
    609 //            cout << "ip " << ip.x << ' ' << ip.y << endl;
    610             if (sgn(ip.x - x.c.x) || ip.y > x.c.y) { puts("Damn!"); while (1) ;}
    611             double ans1 = x.c.y - ip.y;
    612 //            cout << ans1 << endl;
    613             vector<Point> endPt;
    614             if (A[0].y < B[0].y) {
    615                 a.r += x.r;
    616                 lineCircleIntersect(L, a, endPt);
    617             } else {
    618                 b.r += x.r;
    619                 lineCircleIntersect(L, b, endPt);
    620             }
    621             double ans2 = FINF;
    622             for (int i = 0; i < endPt.size(); i++) {
    623                 ans2 = min(ans2, ptDis(endPt[i], ip));
    624             }
    625 //            cout << ans2 << endl;
    626             printf("%.10f\n", ans1 + ans2 / (2.0 * PI * x.r));
    627         }
    628     }
    629     return 0;
    630 }
    View Code

    ——written by Lyon

  • 相关阅读:
    流畅的python--函数
    流暢的python---函數闭包
    IDEA2019与Maven3.6以上版本不兼容
    java面试知识点汇总---第一章 重视基础
    《动手学深度学习》task10 文本分类;数据增强;模型微调 课后作业
    《动手学深度学习》task09 优化算法进阶;word2vec;词嵌入进阶 课后作业
    《动手学深度学习》task08 图像分类案例2;GAN;DCGAN 课后作业
    《动手学深度学习》task08GAN;DCGAN 笔记
    《动手学深度学习》task07 目标检测基础;图像风格迁移;图像分类案例1 课后作业
    《动手学深度学习》task07目标检测基础;图像风格迁移 笔记
  • 原文地址:https://www.cnblogs.com/LyonLys/p/LOJ_1378_Lyon.html
Copyright © 2011-2022 走看看