首页 > ACM题库 > HDU-杭电 > HDU 1140 War on Weather-计算几何-[解题报告] C++
2013
11-27

HDU 1140 War on Weather-计算几何-[解题报告] C++

War on Weather

问题描述 :

After an unprovoked hurricane attack on the south shore, Glorious Warrior has declared war on weather. The first salvo in this campaign will be a coordinated pre-emptive attack on as many tropical depressions as possible. GW reckons that the attack will neutralize the tropical depressions before they become storms, and dissuade others from forming.
GW has at his disposal k space-to-earth killer satellites at various locations in space. m tropical depressions are known to exist at various locations on the earth’s surface. Each satellite can attack any number of targets on the earth provided there is line of sight between the satellite and each target. How many different targets can be hit?

输入:

The input consists of several test cases. Each case begins with a line containing integers 0 < k, m &le 100 as defined above. k lines follow, each giving x,y,z – the location in space of a satellite at the scheduled time of attack. m lines then follow, each giving x,y,z – the location of a target tropical depression. Assume the earth is a sphere centred at (0,0,0) with circumference 40,000 km. All targets will be on the surface of the earth (within 10-9 km) and all satellites will be at least 50 km above the surface. A line containing 0 0 follows the last test case.

输出:

For each test case, output a line giving the total number of targets that can be hit. If a particular target falls within 10-8 km of the boundary between being within line-of-sight and not, it may be counted either way. (That is, you need not consider rounding error so long as it does not exceed 10-8 km.)

样例输入:

3 2
-10.82404031 -1594.10929753 -6239.77925152
692.58497298 -5291.64700245 4116.92402298
3006.49210582 2844.61925179 5274.03201053
2151.03635167 2255.29684503 5551.13972186
-1000.08700886 -4770.25497971 4095.48127333
3 4
0 0 6466.197723676
0 6466.197723676 0
6466.197723676 0 0
6366.197723676 0 0
6365.197723676 112.833485488 0
0 0 6366.197723676
0 -6366.197723676 0
0 0

样例输出:

2
3

Problem – 1140

  简单的三维几何题,对于每一个点,搜索是否存在卫星直接可见,也就是直线与球体是否交于较近的点。其中,套用的主要是直线与平面交点的函数。

模板以及代码:

#include <cstdio>
 #include <cstring>
 #include <cmath>
 #include <set>
 #include <vector>
 #include <iostream>
 #include <algorithm>
 
 using namespace std;
 
 // Point class
 struct Point {
     double x, y;
     Point() {}
     Point(double x, double y) : x(x), y(y) {}
 } ;
 template<class T> T sqr(T x) { return x * x;}
 inline double ptDis(Point a, Point b) { return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));}
 
 // basic calculations
 typedef Point Vec;
 Vec operator + (Vec a, Vec b) { return Vec(a.x + b.x, a.y + b.y);}
 Vec operator - (Vec a, Vec b) { return Vec(a.x - b.x, a.y - b.y);}
 Vec operator * (Vec a, double p) { return Vec(a.x * p, a.y * p);}
 Vec operator / (Vec a, double p) { return Vec(a.x / p, a.y / p);}
 
 const double EPS = 1e-8;
 const double PI = acos(-1.0);
 inline int sgn(double x) { return fabs(x) < EPS ? 0 : (x < 0 ? -1 : 1);}
 bool operator < (Point a, Point b) { return a.x < b.x || (a.x == b.x && a.y < b.y);}
 bool operator == (Point a, Point b) { return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0;}
 
 inline double dotDet(Vec a, Vec b) { return a.x * b.x + a.y * b.y;}
 inline double crossDet(Vec a, Vec b) { return a.x * b.y - a.y * b.x;}
 inline double crossDet(Point o, Point a, Point b) { return crossDet(a - o, b - o);}
 inline double vecLen(Vec x) { return sqrt(dotDet(x, x));}
 inline double toRad(double deg) { return deg / 180.0 * PI;}
 inline double angle(Vec v) { return atan2(v.y, v.x);}
 inline double angle(Vec a, Vec b) { return acos(dotDet(a, b) / vecLen(a) / vecLen(b));}
 inline double triArea(Point a, Point b, Point c) { return fabs(crossDet(b - a, c - a));}
 inline Vec vecUnit(Vec x) { return x / vecLen(x);}
 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));}
 Vec normal(Vec x) {
     double len = vecLen(x);
     return Vec(- x.y / len, x.x / len);
 }
 
 // Line class
 struct Line {
     Point s, t;
     Line() {}
     Line(Point s, Point t) : s(s), t(t) {}
     Point point(double x) {
         return s + (t - s) * x;
     }
     Line move(double x) { // while x > 0 move to (s->t)'s left
         Vec nor = normal(t - s);
         nor = nor * x;
         return Line(s + nor, t + nor);
     }
     Vec vec() { return t - s;}
 } ;
 typedef Line Seg;
 
 inline bool onLine(Point x, Point a, Point b) { return sgn(crossDet(a - x, b - x)) == 0;}
 inline bool onLine(Point x, Line l) { return onLine(x, l.s, l.t);}
 inline bool onSeg(Point x, Point a, Point b) { return sgn(crossDet(a - x, b - x)) == 0 && sgn(dotDet(a - x, b - x)) < 0;}
 inline bool onSeg(Point x, Seg s) { return onSeg(x, s.s, s.t);}
 
 // 0 : not intersect
 // 1 : proper intersect
 // 2 : improper intersect
 int segIntersect(Point a, Point c, Point b, Point d) {
     Vec v1 = b - a, v2 = c - b, v3 = d - c, v4 = a - d;
     int a_bc = sgn(crossDet(v1, v2));
     int b_cd = sgn(crossDet(v2, v3));
     int c_da = sgn(crossDet(v3, v4));
     int d_ab = sgn(crossDet(v4, v1));
     if (a_bc * c_da > 0 && b_cd * d_ab > 0) return 1;
     if (onSeg(b, a, c) && c_da) return 2;
     if (onSeg(c, b, d) && d_ab) return 2;
     if (onSeg(d, c, a) && a_bc) return 2;
     if (onSeg(a, d, b) && b_cd) return 2;
     return 0;
 }
 inline int segIntersect(Seg a, Seg b) { return segIntersect(a.s, a.t, b.s, b.t);}
 
 // point of the intersection of 2 lines
 Point lineIntersect(Point P, Vec v, Point Q, Vec w) {
     Vec u = P - Q;
     double t = crossDet(w, u) / crossDet(v, w);
     return P + v * t;
 }
 inline Point lineIntersect(Line a, Line b) { return lineIntersect(a.s, a.t - a.s, b.s, b.t - b.s);}
 
 // Warning: This is a DIRECTED Distance!!!
 double pt2Line(Point x, Point a, Point b) {
     Vec v1 = b - a, v2 = x - a;
     return crossDet(v1, v2) / vecLen(v1);
 }
 inline double pt2Line(Point x, Line L) { return pt2Line(x, L.s, L.t);}
 
 double pt2Seg(Point x, Point a, Point b) {
     if (a == b) return vecLen(x - a);
     Vec v1 = b - a, v2 = x - a, v3 = x - b;
     if (sgn(dotDet(v1, v2)) < 0) return vecLen(v2);
     if (sgn(dotDet(v1, v3)) > 0) return vecLen(v3);
     return fabs(crossDet(v1, v2)) / vecLen(v1);
 }
 inline double pt2Seg(Point x, Seg s) { return pt2Seg(x, s.s, s.t);}
 
 Point ptOnLine(Point p, Point a, Point b) {
     Vec v = b - a;
     return a + v * (dotDet(v, p - a) / dotDet(v, v));
 }
 inline Point ptOnLine(Point p, Line x) { return ptOnLine(p, x.s, x.t);}
 
 // Polygon class
 struct Poly {
     vector<Point> pt;
     Poly() { pt.clear();}
     ~Poly() {}
     Poly(vector<Point> pt) : pt(pt) {}
     Point operator [] (int x) const { return pt[x];}
     int size() { return pt.size();}
     double area() {
         double ret = 0.0;
         int sz = pt.size();
         for (int i = 1; i < sz; i++) {
             ret += crossDet(pt[i], pt[i - 1]);
         }
         return fabs(ret / 2.0);
     }
 } ;
 
 // Circle class
 struct Circle {
     Point c;
     double r;
     Circle() {}
     Circle(Point c, double r) : c(c), r(r) {}
     Point point(double a) {
         return Point(c.x + cos(a) * r, c.y + sin(a) * r);
     }
 } ;
 
 inline bool ptOnCircle(Point x, Circle c) { return sgn(ptDis(c.c, x) - c.r) == 0;}
 
 // Cirlce operations
 int lineCircleIntersect(Line L, Circle C, double &t1, double &t2, vector<Point> &sol) {
     double a = L.s.x, b = L.t.x - C.c.x, c = L.s.y, d = L.t.y - C.c.y;
     double e = sqr(a) + sqr(c), f = 2 * (a * b + c * d), g = sqr(b) + sqr(d) - sqr(C.r);
     double delta = sqr(f) - 4.0 * e * g;
     if (sgn(delta) < 0) return 0;
     if (sgn(delta) == 0) {
         t1 = t2 = -f / (2.0 * e);
         sol.push_back(L.point(t1));
         return 1;
     }
     t1 = (-f - sqrt(delta)) / (2.0 * e);
     sol.push_back(L.point(t1));
     t2 = (-f + sqrt(delta)) / (2.0 * e);
     sol.push_back(L.point(t2));
     return 2;
 }
 
 int lineCircleIntersect(Line L, Circle C, vector<Point> &sol) {
     Vec dir = L.t - L.s, nor = normal(dir);
     Point mid = lineIntersect(C.c, nor, L.s, dir);
     double len = sqr(C.r) - sqr(ptDis(C.c, mid));
     if (sgn(len) < 0) return 0;
     if (sgn(len) == 0) {
         sol.push_back(mid);
         return 1;
     }
     Vec dis = vecUnit(dir);
     len = sqrt(len);
     sol.push_back(mid + dis * len);
     sol.push_back(mid - dis * len);
     return 2;
 }
 
 // -1 : coincide
 int circleCircleIntersect(Circle C1, Circle C2, vector<Point> &sol) {
     double d = vecLen(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((sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2.0 * 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;
 }
 
 void circleCircleIntersect(Circle C1, Circle C2, vector<double> &sol) {
     double d = vecLen(C1.c - C2.c);
     if (sgn(d) == 0) return ;
     if (sgn(C1.r + C2.r - d) < 0) return ;
     if (sgn(fabs(C1.r - C2.r) - d) > 0) return ;
     double a = angle(C2.c - C1.c);
     double da = acos((sqr(C1.r) + sqr(d) - sqr(C2.r)) / (2.0 * C1.r * d));
     sol.push_back(a - da);
     sol.push_back(a + da);
 }
 
 int tangent(Point p, Circle C, vector<Vec> &sol) {
     Vec u = C.c - p;
     double dist = vecLen(u);
     if (dist < C.r) return 0;
     if (sgn(dist - C.r) == 0) {
         sol.push_back(rotate(u, PI / 2.0));
         return 1;
     }
     double ang = asin(C.r / dist);
     sol.push_back(rotate(u, -ang));
     sol.push_back(rotate(u, ang));
     return 2;
 }
 
 // ptA : points of tangency on circle A
 // ptB : points of tangency on circle B
 int tangent(Circle A, Circle B, vector<Point> &ptA, vector<Point> &ptB) {
     if (A.r < B.r) {
         swap(A, B);
         swap(ptA, ptB);
     }
     double d2 = sqr(A.c.x - B.c.x) + sqr(A.c.y - B.c.y);
     double rdiff = A.r - B.r, rsum = A.r + B.r;
     if (d2 < sqr(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 == sqr(rdiff)) {
         ptA.push_back(A.point(base));
         ptB.push_back(B.point(base));
         return 1;
     }
     double ang = acos((A.r - B.r) / sqrt(d2));
     ptA.push_back(A.point(base + ang));
     ptB.push_back(B.point(base + ang));
     ptA.push_back(A.point(base - ang));
     ptB.push_back(B.point(base - ang));
     if (d2 == sqr(rsum)) {
         ptA.push_back(A.point(base));
         ptB.push_back(B.point(PI + base));
         return 3;
     } else if (d2 > sqr(rsum)) {
         ang = acos((A.r + B.r) / sqrt(d2));
         ptA.push_back(A.point(base + ang));
         ptB.push_back(B.point(PI + base + ang));
         ptA.push_back(A.point(base - ang));
         ptB.push_back(B.point(PI + base - ang));
         return 4;
     }
     return 2;
 }
 
 // -1 : onside
 // 0 : outside
 // 1 : inside
 int ptInPoly(Point p, Poly &poly) {
     int wn = 0, sz = poly.size();
     for (int i = 0; i < sz; i++) {
         if (onSeg(p, poly[i], poly[(i + 1) % sz])) return -1;
         int k = sgn(crossDet(poly[(i + 1) % sz] - poly[i], p - poly[i]));
         int d1 = sgn(poly[i].y - p.y);
         int d2 = sgn(poly[(i + 1) % sz].y - p.y);
         if (k > 0 && d1 <= 0 && d2 > 0) wn++;
         if (k < 0 && d2 <= 0 && d1 > 0) wn--;
     }
     if (wn != 0) return 1;
     return 0;
 }
 
 // if DO NOT need a high precision
 /*
 int ptInPoly(Point p, Poly poly) {
     int sz = poly.size();
     double ang = 0.0, tmp;
     for (int i = 0; i < sz; i++) {
         if (onSeg(p, poly[i], poly[(i + 1) % sz])) return -1;
         tmp = angle(poly[i] - p) - angle(poly[(i + 1) % sz] - p) + PI;
         ang += tmp - floor(tmp / (2.0 * PI)) * 2.0 * PI - PI;
     }
     if (sgn(ang - PI) == 0) return -1;
     if (sgn(ang) == 0) return 0;
     return 1;
 }
 */
 
 
 // Convex Hull algorithms
 // return the number of points in convex hull
 
 // andwer's algorithm
 // if DO NOT want the points on the side of convex hull, change all "<" into "<="
 int andrew(Point *pt, int n, Point *ch) {
     sort(pt, pt + n);
     int m = 0;
     for (int i = 0; i < n; i++) {
         while (m > 1 && crossDet(ch[m - 1] - ch[m - 2], pt[i] - ch[m - 2]) <= 0) m--;
         ch[m++] = pt[i];
     }
     int k = m;
     for (int i = n - 2; i >= 0; i--) {
         while (m > k && crossDet(ch[m - 1] - ch[m - 2], pt[i] - ch[m - 2]) <= 0) m--;
         ch[m++] = pt[i];
     }
     if (n > 1) m--;
     return m;
 }
 
 // graham's algorithm
 // if DO NOT want the points on the side of convex hull, change all "<=" into "<"
 Point origin;
 inline bool cmpAng(Point p1, Point p2) { return crossDet(origin, p1, p2) > 0;}
 inline bool cmpDis(Point p1, Point p2) { return ptDis(p1, origin) > ptDis(p2, origin);}
 
 void removePt(Point *pt, int &n) {
     int idx = 1;
     for (int i = 2; i < n; i++) {
         if (sgn(crossDet(origin, pt[i], pt[idx]))) pt[++idx] = pt[i];
         else if (cmpDis(pt[i], pt[idx])) pt[idx] = pt[i];
     }
     n = idx + 1;
 }
 
 int graham(Point *pt, int n, Point *ch) {
     int top = -1;
     for (int i = 1; i < n; i++) {
         if (pt[i].y < pt[0].y || (pt[i].y == pt[0].y && pt[i].x < pt[0].x)) swap(pt[i], pt[0]);
     }
     origin = pt[0];
     sort(pt + 1, pt + n, cmpAng);
     removePt(pt, n);
     for (int i = 0; i < n; i++) {
         if (i >= 2) {
             while (!(crossDet(ch[top - 1], pt[i], ch[top]) <= 0)) top--;
         }
         ch[++top] = pt[i];
     }
     return top + 1;
 }
 
 
 // Half Plane
 // The intersected area is always a convex polygon.
 // Directed Line class
 struct DLine {
     Point p;
     Vec v;
     double ang;
     DLine() {}
     DLine(Point p, Vec v) : p(p), v(v) { ang = atan2(v.y, v.x);}
     bool operator < (const DLine &L) const { return ang < L.ang;}
     DLine move(double x) { // while x > 0 move to v's left
         Vec nor = normal(v);
         nor = nor * x;
         return DLine(p + nor, v);
     }
 
 } ;
 
 Poly cutPoly(Poly &poly, Point a, Point b) {
     Poly ret = Poly();
     int n = poly.size();
     for (int i = 0; i < n; i++) {
         Point c = poly[i], d = poly[(i + 1) % n];
         if (sgn(crossDet(b - a, c - a)) >= 0) ret.pt.push_back(c);
         if (sgn(crossDet(b - a, c - d)) != 0) {
             Point ip = lineIntersect(a, b - a, c, d - c);
             if (onSeg(ip, c, d)) ret.pt.push_back(ip);
         }
     }
     return ret;
 }
 inline Poly cutPoly(Poly &poly, DLine L) { return cutPoly(poly, L.p, L.p + L.v);}
 
 inline bool onLeft(DLine L, Point p) { return crossDet(L.v, p - L.p) > 0;}
 Point dLineIntersect(DLine a, DLine b) {
     Vec u = a.p - b.p;
     double t = crossDet(b.v, u) / crossDet(a.v, b.v);
     return a.p + a.v * t;
 }
 
 Poly halfPlane(DLine *L, int n) {
     Poly ret = Poly();
     sort(L, L + n);
     int fi, la;
     Point *p = new Point[n];
     DLine *q = new DLine[n];
     q[fi = la = 0] = L[0];
     for (int i = 1; i < n; i++) {
         while (fi < la && !onLeft(L[i], p[la - 1])) la--;
         while (fi < la && !onLeft(L[i], p[fi])) fi++;
         q[++la] = L[i];
         if (fabs(crossDet(q[la].v, q[la - 1].v)) < EPS) {
             la--;
             if (onLeft(q[la], L[i].p)) q[la] = L[i];
         }
         if (fi < la) p[la - 1] = dLineIntersect(q[la - 1], q[la]);
     }
     while (fi < la && !onLeft(q[fi], p[la - 1])) la--;
     if (la - fi <= 1) return ret;
     p[la] = dLineIntersect(q[la], q[fi]);
     for (int i = fi; i <= la; i++) ret.pt.push_back(p[i]);
     return ret;
 }
 
 
 // 3D Geometry
 void getCoor(double R, double lat, double lng, double &x, double &y, double &z) {
     lat = toRad(lat);
     lng = toRad(lng);
     x = R * cos(lat) * cos(lng);
     y = R * cos(lat) * sin(lng);
     z = R * sin(lat);
 }
 
 struct Point3 {
     double x, y, z;
     Point3() {}
     Point3(double x, double y, double z) : x(x), y(y), z(z) {}
 } ;
 typedef Point3 Vec3;
 
 Vec3 operator + (Vec3 a, Vec3 b) { return Vec3(a.x + b.x, a.y + b.y, a.z + b.z);}
 Vec3 operator - (Vec3 a, Vec3 b) { return Vec3(a.x - b.x, a.y - b.y, a.z - b.z);}
 Vec3 operator * (Vec3 a, double p) { return Vec3(a.x * p, a.y * p, a.z * p);}
 Vec3 operator / (Vec3 a, double p) { return Vec3(a.x / p, a.y / p, a.z / p);}
 
 bool operator < (Point3 a, Point3 b) {
     if (a.x != b.x) return a.x < b.x;
     if (a.y != b.y) return a.y < b.y;
     return a.z < b.z;
 }
 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;}
 
 inline double dotDet(Vec3 a, Vec3 b) { return a.x * b.x + a.y * b.y + a.z * b.z;}
 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);}
 inline double vecLen(Vec3 x) { return sqrt(dotDet(x, x));}
 inline double angle(Vec3 a, Vec3 b) { return acos(dotDet(a, b) / vecLen(a) / vecLen(b));}
 inline double ptDis(Point3 a, Point3 b) { return sqrt(vecLen(a - b));}
 inline double triArea(Point3 a, Point3 b, Point3 c) { return vecLen(crossDet(b - a, c - a));}
 Vec3 vecUnit(Vec3 x) { return x / vecLen(x);}
 
 struct Plane {
     Point3 p;
     Vec3 n;
     Plane() {}
     Plane(Point3 p, Vec3 n) : p(p), n(n) {}
 } ;
 
 // Warning: This is a DIRECTED Distance!!!
 inline double pt2Plane(Point3 p, Point3 p0, Vec3 n) { return dotDet(p - p0, n) / vecLen(n);}
 inline double pt2Plane(Point3 p, Plane P) { return pt2Plane(p, P.p, P.n);}
 // get projection on plane
 inline Point3 ptOnPlane(Point3 p, Point3 p0, Vec3 n) { return p + n * pt2Plane(p, p0, n);}
 inline Point3 ptOnPlane(Point3 p, Plane P) { return ptOnPlane(p, P.p, P.n);}
 inline bool ptInPlane(Point3 p, Point3 p0, Vec3 n) { return sgn(dotDet(p - p0, n)) == 0;}
 inline bool ptInPlane(Point3 p, Plane P) { return ptInPlane(p, P.p, P.n);}
 
 struct Line3 {
     Point3 s, t;
     Line3() {}
     Line3(Point3 s, Point3 t) : s(s), t(t) {}
     Vec3 vec() { return t - s;}
 } ;
 typedef Line3 Seg3;
 
 double pt2Line(Point3 p, Point3 a, Point3 b) {
     Vec3 v1 = b - a, v2 = p - a;
     return vecLen(crossDet(v1, v2)) / vecLen(v1);
 }
 inline double pt2Line(Point3 p, Line3 l) { return pt2Line(p, l.s, l.t);}
 
 double pt2Seg(Point3 p, Point3 a, Point3 b) {
     if (a == b) return vecLen(p - a);
     Vec3 v1 = b - a, v2 = p - a, v3 = p - b;
     if (sgn(dotDet(v1, v2)) < 0) return vecLen(v2);
     else if (sgn(dotDet(v1, v3)) > 0) return vecLen(v3);
     else return vecLen(crossDet(v1, v2)) / vecLen(v1);
 }
 inline double pt2Seg(Point3 p, Seg3 s) { return pt2Seg(p, s.s, s.t);}
 
 struct Tri {
     Point3 a, b, c;
     Tri() {}
     Tri(Point3 a, Point3 b, Point3 c) : a(a), b(b), c(c) {}
 } ;
 
 bool ptInTri(Point3 p, Point3 a, Point3 b, Point3 c) {
     double area1 = triArea(p, a, b);
     double area2 = triArea(p, b, c);
     double area3 = triArea(p, c, a);
     double area = triArea(a, b, c);
     return sgn(area1 = area2 + area3 - area) == 0;
 }
 inline bool ptInTri(Point3 p, Tri t) { return ptInTri(p, t.a, t.b, t.c);}
 
 // 0 : not intersect
 // 1 : intersect at only 1 point
 // 2 : line in plane
 int linePlaneIntersect(Point3 s, Point3 t, Point3 p0, Vec3 n, Point3 &x) {
     double res1 = dotDet(n, p0 - s);
     double res2 = dotDet(n, t - s);
     if (sgn(res2) == 0) {
         if (ptInPlane(s, p0, n)) return 2;
         return 0;
     }
     x = s + (t - s) * (res1 / res2);
     // use general form : a * x + b * y + c * z + d = 0
     // Vec3 v = s - t;
     // double k = (a * s.x + b * s.y + c * s.z + d) / (a * v.x + b * v.y + c * v.z);
     // x = s + (t - s) * k;
     return 1;
 }
 inline int linePlaneIntersect(Line3 l, Plane p, Point3 &x) { return linePlaneIntersect(l.s, l.t, p.p, p.n, x);}
 
 Point3 ptOnLine(Point3 p, Line3 L) {
     Plane P = Plane(p, L.t - L.s);
     Point3 ret;
     if (linePlaneIntersect(L, P, ret) != 1) puts("Shit!");
     return ret;
 }
 
 // ¡÷abc intersect with segment st
 bool triSegIntersect(Point3 a, Point3 b, Point3 c, Point3 s, Point3 t, Point3 &x) {
     Vec3 n = crossDet(b - a, c - a);
     if (sgn(dotDet(n, t - s)) == 0) return false;
     else {
         double k = dotDet(n, a - s) / dotDet(n, t - s);
         if (sgn(k) < 0 || sgn(k - 1) > 0) return false;
         x = s + (t - s) * k;
         return ptInTri(x, a, b, c);
     }
 }
 inline bool triSegIntersect(Tri t, Line3 l, Point3 &x) { return triSegIntersect(t.a, t.b, t.c, l.s, l.t, x);}
 
 // Warning: This is a DIRECTED Volume!!!
 inline double tetraVol(Point3 a, Point3 b, Point3 c, Point3 p) { return dotDet(p - a, crossDet(b - a, c - a)) / 6.0;}
 inline double tetraVol(Tri t, Point3 p) { return tetraVol(t.a, t.b, t.c, p);}
 
 struct Ball {
     Point3 c;
     double r;
     Ball() {}
     Ball(Point3 c, double r) : c(c), r(r) {}
 } ;
 
 int lineBallIntersect(Ball B, Line3 L, vector<Point3> &sol) {
     double dis = pt2Line(B.c, L);
     if (dis > B.r) return 0;
     Point3 ip = ptOnLine(B.c, L);
     if (sgn(dis - B.r)) {
         double ds = sqrt(sqr(B.r) - sqr(dis));
         sol.push_back(ip + vecUnit(L.t - L.s) * ds);
         sol.push_back(ip - vecUnit(L.t - L.s) * ds);
         return 2;
     } else {
         sol.push_back(ip);
         return 1;
     }
 }
 
 /****************** template above *******************/
 
 const int N = 111;
 const double FINF = 1e100;
 Point3 sate[N];
 
 bool test(Point3 x, Point3 st) {
     vector<Point3> ips;
     ips.clear();
     int c = lineBallIntersect(Ball(Point3(0.0, 0.0, 0.0), vecLen(x)), Line3(st, x), ips);
     if (c == 0) return false;
     if (c == 1) return true;
     double mini = FINF;
     for (int i = 0; i < 2; i++) {
         mini = min(mini, ptDis(st, ips[i]));
     }
     return sgn(ptDis(st, x) - mini) == 0;
 }
 
 bool check(Point3 x, int n) {
     for (int i = 0; i < n; i++) {
         if (test(x, sate[i])) return true;
     }
     return false;
 }
 
 int main() {
 //    freopen("in", "r", stdin);
     int n, m;
     while (cin >> n >> m && (n || m)) {
         double x, y, z;
         for (int i = 0; i < n; i++) {
             cin >> x >> y >> z;
             sate[i] = Point3(x, y, z);
         }
         int cnt = 0;
         for (int i = 0; i < m; i++) {
             cin >> x >> y >> z;
             if (check(Point3(x, y, z), n)) cnt++;
         }
         cout << cnt << endl;
     }
     return 0;
 }

  至此,hdu分类中的几道基础几何已经切完。

P.S.: 模板有更新,之前的模板在点到直线上的投影一函数中有错。

——written by Lyon


  1. #include <cstdio>

    int main() {
    //answer must be odd
    int n, u, d;
    while(scanf("%d%d%d",&n,&u,&d)==3 && n>0) {
    if(n<=u) { puts("1"); continue; }
    n-=u; u-=d; n+=u-1; n/=u;
    n<<=1, ++n;
    printf("%dn",n);
    }
    return 0;
    }

  2. #include <cstdio>
    #include <cstring>

    const int MAXSIZE=256;
    //char store[MAXSIZE];
    char str1[MAXSIZE];
    /*
    void init(char *store) {
    int i;
    store['A']=’V', store['B']=’W',store['C']=’X',store['D']=’Y',store['E']=’Z';
    for(i=’F';i<=’Z';++i) store =i-5;
    }
    */
    int main() {
    //freopen("input.txt","r",stdin);
    //init(store);
    char *p;
    while(fgets(str1,MAXSIZE,stdin) && strcmp(str1,"STARTn")==0) {
    if(p=fgets(str1,MAXSIZE,stdin)) {
    for(;*p;++p) {
    //*p=store[*p]
    if(*p<’A’ || *p>’Z') continue;
    if(*p>’E') *p=*p-5;
    else *p=*p+21;
    }
    printf("%s",str1);
    }
    fgets(str1,MAXSIZE,stdin);
    }
    return 0;
    }

  3. 在方法1里面:

    //遍历所有的边,计算入度
    for(int i=0; i<V; i++)
    {
    degree = 0;
    for (j = adj .begin(); j != adj .end(); ++j)
    {
    degree[*j]++;
    }
    }

    为什么每遍历一条链表,要首先将每个链表头的顶点的入度置为0呢?
    比如顶点5,若在顶点1、2、3、4的链表中出现过顶点5,那么要增加顶点5的入度,但是在遍历顶点5的链表时,又将顶点5的入度置为0了,那之前的从顶点1234到顶点5的边不是都没了吗?