/*
博主效率很低，可能要很久才能完成一篇，所以-_-#
*/
向量与点的表示 向量的基本运算
我们之所以采用向量代数而非解析几何的方法来描述二维平面，是因为解析几何中过多除法运算会带来精度问题，以及斜率是否存在等需要分类讨论的问题会使程序累赘。
在二维平面上，向量可由一个二维坐标表示，而点可以看作是从原点出发的向量，所以表示方法和向量一样。
struct Vector{ double x, y; Vector(double x_ = 0, double y_ = 0):x(x_),y(y_){} }; typedef Vector Point;
对双精度实数，判断大小的时候要考虑浮点误差。
const double EPS = 1e-8; int sgn(double x) {if (-EPS < x && x < EPS) return 0; return x > EPS ? 1 : -1;} bool operator==(Vector p1, Vector p2){return !sgn(p1.x - p2.x) && !sgn(p1.y - p2.y);}
向量的运算主要包括线性运算、点乘运算和叉乘运算。
线性运算包括加法与数乘，对运算符重载即可实现。
点乘运算（Dot）可用于求向量模长（len）、求点与点的距离（dis）、求夹角（angle）、求投影长度（project）、求镜面对称（sym）等。
叉乘运算（Cross）非常重要，主要用于求左右和求面积。向量a叉乘向量b的几何意义为以a、b为邻边的平行四边形的有向面积，当a以小于180度的角逆时针（左）转到b时结果为正，反之（右）为负。叉乘为0当且仅当两个向量平行。
double Dot(Vector p1, Vector p2){return p1.x * p2.x + p1.y * p2.y;} double Cross(Vector p1, Vector p2){return p1.x * p2.y - p1.y * p2.x;}
/* Linear operations */ Vector operator+(Vector p1, Vector p2){return Vector(p1.x + p2.x, p1.y + p2.y);} Vector operator-(Vector p1, Vector p2){return Vector(p1.x - p2.x, p1.y - p2.y);} Vector operator*(Vector p, double k){return Vector(k * p.x , k * p.y);} Vector operator*(double k, Vector p){return p * k;} Vector operator/(Vector p, double k){return Vector(p.x / k, p.y / k);} bool operator==(Vector p1, Vector p2){return !sgn(p1.x - p2.x) && !sgn(p1.y - p2.y);} /* Dot product operations */ double Dot(Vector p1, Vector p2){return p1.x * p2.x + p1.y * p2.y;} double len(Vector p){return sqrt(Dot(p, p));} double dis(Point p1, Point p2){return len(p1 - p2);} double angle(Vector p1, Vector p2){return acos(Dot(p1, p2) / len(p1) / len(p2));} double project(Vector a, Vector b){return Dot(a, b) / len(b);} Vector sym(Vector b, Vector a){return 2 * Dot(b, a) / Dot(a, a) * a - b;} /* Cross product operations */ double Cross(Vector p1, Vector p2){return p1.x * p2.y - p1.y * p2.x;} double area(Point p1, Point p2, Point p3){return fabs(Cross(p2 - p1, p3 - p1)) / 2;}
利用叉乘求任意多边形的面积：设多边形顶点按照逆时针顺序为$P_0, P_1, ..., P_{n-1}$，则所求面积
$S=\frac12 \sum_{i=0}^{n-1} \mathrm{Cross}\ (P_i, P_{i+1}), \mathrm{assuming}\ P_n=P_0.$
如果选取一个顶点（如$P_0$）作为所有向量的起点，面积也可以写成
$S=\frac12 \sum_{i=1}^{n-2} \mathrm{Cross}\ (P_i-P_0, P_{i+1}-P_0),$
即用$P_0$对多边形进行三角剖分，然后求出所有三角形的有向面积之和。如果顶点顺序未知，取绝对值即可。
对于格点多边形的面积，有Pick定理：S = in + on/2 - 1，其中in和on分别表示多边形内部和边界上的格点数。由于多边形每条边（dx, dy）上的格点数目为 gcd(|dx|, |dy|)+1，于是on的数量可求，所以可以通过先用叉乘求出面积，再求出on的数量，反过来求出in的数量。
线的表示 点与线的位置关系
线通常用两个点A、B表示，表示线的方向为A->B。
关于线的操作细节比较多，因为有时需要区分直线（line）和线段（segment），例如求点到线段的距离就比求点到直线的距离复杂一些。
点与线的位置关系涉及点的旋转、求垂足、求垂线、判断线线相交情况、求交点等操作。
struct Segment{ Point a, b;
Segment(){} Segment(Point a_, Point b_):a(a_),b(b_){} }; typedef Segment Line;
/* Distance between Points and Lines */ double len(Segment l){return len(l.a - l.b);} double disToLine(Point p, Line l){return fabs(Cross(p - l.a, p - l.b) / len(l));} double disToSeg(Point p, Segment l){ Vector v1 = p - l.a, v2 = p - l.b, v3 = l.b - l.a; if (sgn(Dot(v1, v3)) < 0) return len(v1); else if (sgn(Dot(v2, v3)) > 0) return len(v2); else return fabs(Cross(v1, v2) / len(v3)); } /* Angles and Verticals */ Point rotate(Point b, Point a, double alpha){ Vector p = b - a; return Point(a.x + (p.x * cos(alpha) - p.y * sin(alpha)), a.y + (p.x * sin(alpha) + p.y * cos(alpha))); } Point foot(Point p, Line l){return l.a + project(p - l.a, l.b - l.a) * (l.b - l.a) / len(l);} Line Vertical(Point p, Line l){return Line(p, foot(p, l));}
判断线线相交情况分为三类问题：
I. 判断直线与直线的相交情况
II. 判断直线与线段的相交情况
III. 判断线段与线段的相交情况
均可利用叉乘解决。对问题I检查叉积是否为0，对问题II先检查叉积是否为0，再检查线段的端点是否在直线的异侧，对问题III先检查叉积是否为0，再检查l1的端点是否在l2异侧且l2的端点是否在l1异侧。注：这里的相交包括交点为线段的一个端点的情形，但不包括共线且有公共点的情形。如果要求必须规范相交，则检查叉积符号从 <=0 改为 <0 即可。如果只要求有公共点而未必是相交，则后两个问题中不需要检查叉积是否为0。
/* Positional Relations */ bool Intersect1(Line l1, Line l2) { Point p1 = l1.a, p2 = l1.b, p3 = l2.a, p4 = l2.b; return sgn(Cross(p2 - p1, p4 - p3)) != 0; } bool Intersect2(Line l1, Segment l2) { Point p1 = l1.a, p2 = l1.b, p3 = l2.a, p4 = l2.b; return sgn(Cross(p2 - p1, p4 - p3)) != 0 && sgn(Cross(p2 - p1, p3 - p1) * Cross(p2 - p1, p4 - p1)) <= 0; } bool Intersect3(Segment l1, Segment l2) { Point p1 = l1.a, p2 = l1.b, p3 = l2.a, p4 = l2.b; return sgn(Cross(p2 - p1, p4 - p3)) != 0 && sgn(Cross(p2 - p1, p3 - p1) * Cross(p2 - p1, p4 - p1)) <= 0 && sgn(Cross(p4 - p3, p1 - p3) * Cross(p4 - p3, p2 - p3)) <= 0; }
还有两个比较常用的操作：求两条直线的交点（Intersection）以及判断点是否在线段上（PointInSeg）。这里求交点的算法对于延长后相交的线段也适用，思想是利用一条对角线截另一条对角线所得线段长度之比等于相应三角形面积之比。
Point Intersection(Line l1, Line l2) { Point p1 = l1.a, p2 = l1.b, p3 = l2.a, p4 = l2.b; double a1 = Cross(p3 - p1, p4 - p1) / 2; double a2 = Cross(p4 - p2, p3 - p2) / 2; return p1 + (p2 - p1) * (a1 / (a1 + a2)); } bool PointInSeg(Point p, Segment l) { if (sgn(Cross(p - l.a, p - l.b))) return 0; return sgn((p.x - l.a.x) * (p.x - l.b.x)) <= 0 && sgn((p.y - l.a.y) * (p.y - l.b.y)) <= 0; }
如果要对问题III的情形进行彻底分类，则需要考虑以下所有共10种情形：
| A. l1与l2所在直线相交 |
| a. 交点在线段上 |
| 0. 规范相交 |
| 1. l1端点在l2内部 |
| 2. l2端点在l1内部 |
| 3. l1与l2一个端点重合 |
| b. 交点不在线段上 |
| 4. l1与l2的延长线相交 |
| 5. l2与l1的延长线相交 |
| 6. l1与l2均延长后相交 |
| B. l1与l2所在直线平行或共线 |
| 7. 平行 |
| 8. 共线且无公共点 |
| 9. 共线且有公共点 |
下面的Intersect3plus函数实现了对两条线段所有可能位置关系的彻底分类，并返回一个pair<int, Point>对象，int值取值0~9按照上表分别对应十种情况，Point值（若存在）则为两条线段的公共点。
pair<int, Point> Intersect3plus(Segment l1, Segment l2) { Point p1 = l1.a, p2 = l1.b, p3 = l2.a, p4 = l2.b; if (sgn(Cross(p2 - p1, p4 - p3))) // 所在直线相交 { Point p = Intersection(l1, l2); if (sgn(Cross(p2 - p1, p3 - p1) * Cross(p2 - p1, p4 - p1)) < 0 && sgn(Cross(p4 - p3, p1 - p3) * Cross(p4 - p3, p2 - p3)) < 0) return make_pair(0, p); // 规范相交 if (p1 == p3 || p1 == p4 || p2 == p3 || p2 == p4) return make_pair(3, p); if (PointInSeg(p1, l2) || PointInSeg(p2, l2)) return make_pair(1, p); if (PointInSeg(p3, l1) || PointInSeg(p4, l1)) return make_pair(2, p); if (PointInSeg(p, l1)) return make_pair(4, p); if (PointInSeg(p, l2)) return make_pair(5, p); return make_pair(6, p); } else // 所在直线平行或共线 { if (sgn(disToLine(p1, l2))) return make_pair(7, Point(0, 0)); if (PointInSeg(p1, l2)) return make_pair(9, p1); if (PointInSeg(p2, l2)) return make_pair(9, p2); if (PointInSeg(p3, l1)) return make_pair(9, p3); if (PointInSeg(p4, l1)) return make_pair(9, p4); return make_pair(8, Point(0, 0)); } }
有了以上这些模板，就可以着手开始解决计算几何的问题了。
Problems
题意：将点按照逆时针排序，排序结果保证能和原点一起构成凸多边形。
由于P1<P2等价于P1绕原点逆时针旋转到P2，等价于Cross(P1, P2)>0，所以利用叉乘排序即可。
1 // Problem: poj2007 - Scrambled Polygon 2 // Category: Computational geometry 3 // Author: Niwatori 4 // Date: 2016/07/27 5 6 #include <stdio.h> 7 #include <algorithm> 8 #define MAXN 60 9 10 struct Vector{ 11 int x, y; 12 Vector(int x_ = 0, int y_ = 0):x(x_),y(y_){} 13 }; typedef Vector Point; 14 15 Vector operator-(Vector p1, Vector p2){return Vector(p1.x - p2.x, p1.y - p2.y);} 16 int Cross(Vector p1, Vector p2){return p1.x * p2.y - p1.y * p2.x;} 17 bool operator<(const Point & p1, const Point & p2){return Cross(p1, p2) > 0;} 18 19 int main() 20 { 21 Point p[MAXN]; 22 int cnt = 0; 23 while (scanf("%d%d", &p[cnt].x, &p[cnt].y) == 2) ++cnt; 24 25 std::sort(p + 1, p + cnt); 26 for (int i = 0; i < cnt; ++i) 27 printf("(%d,%d)\n", p[i].x, p[i].y); 28 return 0; 29 }
题意：给定一个凸多边形，求其内部格点数、边界上的格点数以及多边形的面积。
用叉乘求面积，用gcd求边界上格点数，用Pick定理求内部格点数（然而不知何故G++WA，C++AC）。
提醒：浮点数的绝对值函数fabs在<math.h>头文件里，int和long long的绝对值函数abs、llabs在头文件<stdlib.h>里。
1 // Problem: poj1265 - Area 2 // Category: Computational geometry 3 // Author: Niwatori 4 // Date: 2016/07/28 5 6 #include <stdio.h> 7 #include <stdlib.h> 8 9 struct Vector{ 10 int x, y; 11 Vector(int x_ = 0, int y_ = 0):x(x_),y(y_){} 12 }; typedef Vector Point; 13 14 Vector operator+(Vector p1, Vector p2){return Vector(p1.x + p2.x, p1.y + p2.y);} 15 int Cross(Vector p1, Vector p2){return p1.x * p2.y - p1.y * p2.x;} 16 int gcd(int a, int b){return b ? gcd(b, a % b) : a;} 17 18 int main() 19 { 20 int t; scanf("%d\n", &t); 21 for (int i = 1; i <= t; ++i) 22 { 23 int n, dx, dy, ans = 0, on = 0; 24 scanf("%d", &n); 25 Point p1, p2; 26 while (n--) 27 { 28 scanf("%d%d", &dx, &dy); 29 p1 = p2; 30 p2 = p1 + Vector(dx, dy); 31 ans += Cross(p1, p2); 32 on += gcd(abs(dx), abs(dy)); 33 } 34 printf("Scenario #%d:\n", i); 35 printf("%d %d %.1f\n\n", (ans - on) / 2 + 1, on, ans / 2.0); 36 } 37 return 0; 38 }
题意：给定n条线段（n<=100），判断是否存在这样一条直线，使得所有线段在其上的投影有至少一个公共点。
假设存在这样的一条直线l，使得所有线段在其上的投影有公共点P，那么过点P作l的垂线m，必有m与所有线段有公共点，因此只需判断是否存在一条直线与所有线段有公共点。而如果这样的直线存在，则将其适当平移及绕点旋转后，必可让它碰到所有线段中的某两个端点，因而只需枚举任两个（不重合的）端点构造直线，然后检查这条直线是否与所有线段有公共点即可。注意：本题中共线且有公共点的情形是合法的。
1 // Problem: poj3304 - Segments 2 // Category: Computational geometry 3 // Author: Niwatori 4 // Date: 2016/07/30 5 6 #include <cstdio> 7 const double EPS = 1e-8; 8 inline int sgn(double x) {if (-EPS < x && x < EPS) return 0; return x > EPS ? 1 : -1;} 9 10 struct Vector{ 11 double x, y; 12 Vector(double x_ = 0, double y_ = 0):x(x_),y(y_){} 13 }; typedef Vector Point; 14 15 struct Segment{ 16 Point a, b; 17 Segment(){} 18 Segment(Point a_, Point b_):a(a_),b(b_){} 19 }; typedef Segment Line; 20 21 Vector operator-(Vector p1, Vector p2){return Vector(p1.x - p2.x, p1.y - p2.y);} 22 bool operator!=(Vector p1, Vector p2){return sgn(p1.x - p2.x) || sgn(p1.y - p2.y);} 23 double Cross(Vector p1, Vector p2){return p1.x * p2.y - p1.y * p2.x;} 24 25 bool Intersect(Line l1, Segment l2) 26 { 27 Point p1 = l1.a, p2 = l1.b, p3 = l2.a, p4 = l2.b; 28 return sgn(Cross(p2 - p1, p3 - p1) * Cross(p2 - p1, p4 - p1)) <= 0; 29 } 30 31 int main() 32 { 33 int t; scanf("%d", &t); 34 while (t--) 35 { 36 Segment seg[105]; Point p[210]; 37 int n, cnt = 0, ok = 0; 38 scanf("%d", &n); 39 for (int i = 0; i < n; ++i) 40 scanf("%lf%lf%lf%lf", &seg[i].a.x, &seg[i].a.y, &seg[i].b.x, &seg[i].b.y), 41 p[cnt++] = seg[i].a, p[cnt++] = seg[i].b; 42 for (int i = 0; i < cnt; ++i) 43 for (int j = i + 1; j < cnt; ++j) 44 if (p[i] != p[j]) 45 { 46 int tmp = 1; 47 for (int k = 0; k < n; ++k) 48 tmp = tmp && Intersect(Line(p[i], p[j]), seg[k]); 49 ok = ok || tmp; 50 } 51 printf("%s\n", ok ? "Yes!" : "No!"); 52 } 53 return 0; 54 }
题意：给定一根折线管道（上下边界均由n条线段组成，n<=20），上边界总比下边界高1个单位，求从入口处射入的光线最远能射到哪儿，或者能穿过整根管道。
能射到最远的光线必然会擦到一个上端点和一个下端点，因而只需枚举一个上端点和一个下端点构成直线l，然后一方面检查光线前半部分是否都在管道内，另一方面计算光线后半部分最远能打到哪儿。要检查光线是否在管道内，只需检查直线是否与管道每一个拐处的竖直线段（(xi, yi-1), (xi, yi)）相交即可。
1 // Problem: poj1039 - Pipe 2 // Category: Computational geometry 3 // Author: Niwatori 4 // Date: 2016/08/02 5 6 #include <stdio.h> 7 #define MAX(a,b) ((a)>(b)?(a):(b)) 8 #define INF 1e30 9 10 const double EPS = 1e-8; 11 int sgn(double x) {if (-EPS < x && x < EPS) return 0; return x > EPS ? 1 : -1;} 12 13 struct Vector{ 14 double x, y; 15 Vector(double x_ = 0, double y_ = 0):x(x_),y(y_){} 16 }; typedef Vector Point; 17 18 struct Segment{ 19 Point a, b; 20 Segment(){} 21 Segment(Point a_, Point b_):a(a_),b(b_){} 22 }; typedef Segment Line; 23 24 Vector operator+(Vector p1, Vector p2){return Vector(p1.x + p2.x, p1.y + p2.y);} 25 Vector operator-(Vector p1, Vector p2){return Vector(p1.x - p2.x, p1.y - p2.y);} 26 Vector operator*(Vector p, double k){return Vector(k * p.x , k * p.y);} 27 double Cross(Vector p1, Vector p2){return p1.x * p2.y - p1.y * p2.x;} 28 29 bool Intersect(Line l1, Segment l2) 30 { 31 Point p1 = l1.a, p2 = l1.b, p3 = l2.a, p4 = l2.b; 32 return sgn(Cross(p2 - p1, p4 - p3)) != 0 33 && sgn(Cross(p2 - p1, p3 - p1) * Cross(p2 - p1, p4 - p1)) <= 0; 34 } 35 36 Point Intersection(Line l1, Line l2) 37 { 38 Point p1 = l1.a, p2 = l1.b, p3 = l2.a, p4 = l2.b; 39 double a1 = Cross(p3 - p1, p4 - p1) / 2; 40 double a2 = Cross(p4 - p2, p3 - p2) / 2; 41 return p1 + (p2 - p1) * (a1 / (a1 + a2)); 42 } 43 44 int main() 45 { 46 int n; 47 while (scanf("%d", &n) == 1) 48 { 49 if (!n) break; 50 Segment seg[30]; 51 for (int i = 0; i < n; ++i) 52 { 53 scanf("%lf%lf", &seg[i].b.x, &seg[i].b.y); 54 seg[i].a.x = seg[i].b.x; 55 seg[i].a.y = seg[i].b.y - 1; 56 } 57 58 double ans = -INF; 59 int ok = 0; 60 for (int i = 0; i < n; ++i) 61 for (int j = 0; j < n; ++j) 62 { 63 int check = 1; 64 Line l = Line(seg[i].a, seg[j].b); 65 for (int k = 0; k <= MAX(i, j); ++k) 66 check = check && Intersect(l, seg[k]); 67 if (!check) continue; 68 69 for (int k = MAX(i, j); k < n; ++k) 70 if (!Intersect(l, seg[k])) 71 { 72 check = 0; 73 Segment up = Segment(seg[k - 1].b, seg[k].b); 74 Segment down = Segment(seg[k - 1].a, seg[k].a); 75 if (Intersect(l, up)) 76 ans = MAX(ans, Intersection(l, up).x); 77 if (Intersect(l, down)) 78 ans = MAX(ans, Intersection(l, down).x); 79 break; 80 } 81 ok = ok || check; 82 } 83 84 if (ok) printf("Through all the pipe.\n"); 85 else printf("%.2f\n", ans); 86 } 87 return 0; 88 }
题意：正方形区域内有n堵贯穿区域的墙（0<=n<=30，没有三线共点），问从区域外走到区域内的指定点至少需要穿越多少堵墙（穿墙时只能从墙的中点处穿过）。
一方面，容易知道只需从正方形边上出发求解即可，进一步也容易知道只需从所有墙与区域边的交点出发求解即可（从这些点出发可达区域边上任意一点）。另一方面，需要穿越的最少墙数与路径无关，因而只需枚举墙与区域边的交点，将其与指定点相连，检查它与多少堵墙相交即可。注意这题中要求的是线段规范相交。由于穿墙只从墙的中点处穿过，所以不需要考虑从墙的交点处穿过的情形。对n=0的情形注意特判。
1 // Problem: poj1066 - Treasure Hunt 2 // Category: Computational geometry 3 // Author: Niwatori 4 // Date: 2016/08/06 5 6 #include <stdio.h> 7 #define INF 100 8 #define MIN(a,b) ((a)<(b)?(a):(b)) 9 10 const double EPS = 1e-8; 11 int sgn(double x) {if (-EPS < x && x < EPS) return 0; return x > EPS ? 1 : -1;} 12 13 struct Vector{ 14 double x, y; 15 Vector(double x_ = 0, double y_ = 0):x(x_),y(y_){} 16 }; typedef Vector Point; 17 18 struct Segment{ 19 Point a, b; 20 Segment(){} 21 Segment(Point a_, Point b_):a(a_),b(b_){} 22 }; typedef Segment Line; 23 24 Vector operator-(Vector p1, Vector p2){return Vector(p1.x - p2.x, p1.y - p2.y);} 25 double Cross(Vector p1, Vector p2){return p1.x * p2.y - p1.y * p2.x;} 26 27 bool Intersect(Segment l1, Segment l2) 28 { 29 Point p1 = l1.a, p2 = l1.b, p3 = l2.a, p4 = l2.b; 30 return sgn(Cross(p2 - p1, p3 - p1) * Cross(p2 - p1, p4 - p1)) < 0 31 && sgn(Cross(p4 - p3, p1 - p3) * Cross(p4 - p3, p2 - p3)) < 0; 32 } 33 34 int main() 35 { 36 int n; scanf("%d", &n); 37 Point p[80], dest; 38 for (int i = 0; i < n; ++i) 39 scanf("%lf%lf%lf%lf", &p[2*i].x, &p[2*i].y, &p[2*i+1].x, &p[2*i+1].y); 40 scanf("%lf%lf", &dest.x, &dest.y); 41 42 int ans = INF; 43 if (!n) ans = 0; 44 for (int i = 0; i < 2 * n; ++i) 45 { 46 Segment l = Segment(p[i], dest); 47 int cnt = 0; 48 for (int j = 0; j < 2 * n; j += 2) 49 cnt += Intersect(l, Segment(p[j], p[j + 1])); 50 ans = MIN(ans, cnt); 51 } 52 printf("Number of doors = %d\n", ans + 1); 53 return 0; 54 }
251
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
