1 #include <cstdio>  2 #include <algorithm>  3 #include <cmath>  4 using namespace std;  5 //lrj计算几何模板  6 struct Point  7 {  8     double x, y;  9     Point(double x=0, double y=0) :x(x),y(y) {} 10 }; 11 typedef Point Vector; 12 const double EPS = 1e-10; 13  14 //向量+向量=向量 点+向量=点 15 Vector operator + (Vector A, Vector B)    { return Vector(A.x + B.x, A.y + B.y); } 16  17 //向量-向量=向量 点-点=向量 18 Vector operator - (Vector A, Vector B)    { return Vector(A.x - B.x, A.y - B.y); } 19  20 //向量*数=向量 21 Vector operator * (Vector A, double p)    { return Vector(A.x*p, A.y*p); } 22  23 //向量/数=向量 24 Vector operator / (Vector A, double p)    { return Vector(A.x/p, A.y/p); } 25  26 bool operator < (const Point& a, const Point& b) 27 { return a.x < b.x || (a.x == b.x && a.y < b.y); } 28  29 int dcmp(double x) 30 { if(fabs(x) < EPS) return 0; else x < 0 ? -1 : 1; } 31  32 bool operator == (const Point& a, const Point& b) 33 { return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0; } 34  35 /**********************基本运算**********************/ 36  37 //点积 38 double Dot(Vector A, Vector B) 39 { return A.x*B.x + A.y*B.y; } 40 //向量的模 41 double Length(Vector A)    { return sqrt(Dot(A, A)); } 42  43 //向量的夹角，返回值为弧度 44 double Angle(Vector A, Vector B) 45 { return acos(Dot(A, B) / Length(A) / Length(B)); } 46  47 //叉积 48 double Cross(Vector A, Vector B) 49 { return A.x*B.y - A.y*B.x; } 50  51 //向量AB叉乘AC的有向面积 52 double Area2(Point A, Point B, Point C) 53 { return Cross(B-A, C-A); } 54  55 //向量A旋转rad弧度 56 Vector VRotate(Vector A, double rad) 57 { 58     return Vector(A.x*cos(rad) - A.y*sin(rad), A.x*sin(rad) + A.y*cos(rad)); 59 } 60  61 //将B点绕A点旋转rad弧度 62 Point PRotate(Point A, Point B, double rad) 63 { 64     return A + VRotate(B-A, rad); 65 } 66  67 //求向量A向左旋转90°的单位法向量,调用前确保A不是零向量 68 Vector Normal(Vector A) 69 { 70     double l = Length(A); 71     return Vector(-A.y/l, A.x/l); 72 } 73  74 /**********************点和直线**********************/ 75  76 //求直线P + tv 和 Q + tw的交点，调用前要确保两条直线有唯一交点 77 Point GetLineIntersection(Point P, Vector v, Point Q, Vector w) 78 { 79     Vector u = P - Q; 80     double t = Cross(w, u) / Cross(v, w); 81     return P + v*t; 82 }//在精度要求极高的情况下，可以自定义分数类 83  84 //P点到直线AB的距离 85 double DistanceToLine(Point P, Point A, Point B) 86 { 87     Vector v1 = B - A, v2 = P - A; 88     return fabs(Cross(v1, v2)) / Length(v1);    //不加绝对值是有向距离 89 } 90  91 //点到线段的距离 92 double DistanceToSegment(Point P, Point A, Point B) 93 { 94     if(A == B)    return Length(P - A); 95     Vector v1 = B - A, v2 = P - A, v3 = P - B; 96     if(dcmp(Dot(v1, v2)) < 0)    return Length(v2); 97     else if(dcmp(Dot(v1, v3)) > 0)    return Length(v3); 98     else return fabs(Cross(v1, v2)) / Length(v1); 99 }100 101 //点在直线上的射影102 Point GetLineProjection(Point P, Point A, Point B)103 {104     Vector v = B - A;105     return A + v * (Dot(v, P - A) / Dot(v, v));106 }107 108 //线段“规范”相交判定109 bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point b2)110 {111     double c1 = Cross(a2-a1, b1-a1), c2 = Cross(a2-a1, b2-a1);112     double c3 = Cross(b2-b1, a1-b1), c4 = Cross(b2-b1, a2-b1);113     return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0;114 }115 116 //判断点是否在线段上117 bool OnSegment(Point P, Point a1, Point a2)118 {119     Vector v1 = a1 - P, v2 = a2 - P;120     return dcmp(Cross(v1, v2)) == 0 && dcmp(Dot(v1, v2)) < 0;121 }122 123 //求多边形面积124 double PolygonArea(Point* P, int n)125 {126     double ans = 0.0;127     for(int i = 1; i < n - 1; ++i)128         ans += Cross(P[i]-P[0], P[i+1]-P[0]);129     return ans/2;130 }