const double pi = acos(-1.0);
const double eps = 1e-8;
// 判断 x 的大小，<0 返回 -1，>0 返回 1，==0 返回 0
int sgn(double x) {
if (fabs(x) < eps) return 0;
else return x < 0 ? -1 : 1;
}
// 比较两个浮点数
int dcmp(double x, double y) {
if (fabs(x - y) < eps) return 0;
else return x < y ? -1 : 1;
}
struct Point {
double x, y;
Point() {}
Point(double x, double y) : x(x), y(y) {}
Point operator + (const Point& B) const { return Point(x + B.x, y + B.y); }
Point operator - (const Point& B) const { return Point(x - B.x, y - B.y); }
Point operator * (double k) const { return Point(x * k, y * k); }
Point operator / (double k) const { return Point(x / k, y / k); }
bool operator == (const Point& B) const {
return sgn(x - B.x) == 0 && sgn(y - B.y) == 0;
}
bool operator < (const Point& B) const {
return sgn(x - B.x) < 0 || sgn(x - B.x) == 0 && sgb(y - B.y) < 0;
}
};
typedef Point Vector;
class Geometry {
public:
static double Distance(const Point& A, const Point& B) {
return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}
static double Dot(const Vector& A, const Vector& B) {
return A.x * B.x + A.y * B.y;
}
static int AngleJudge(const Vector& A, const Vector& B) {
return sgn(Dot(A, B));
}
static double Len(const Vector& A) {
return sqrt(Dot(A, A));
}
static double Len2(const Vector& A) {
return Dot(A, A);
}
static double Angle(const Vector& A, const Vector& B) {
return acos(Dot(A, B) / Len(A) / Len(B));
}
static double Cross(const Vector& A, const Vector& B) {
return A.x * B.y - A.y * B.x;
}
static double Area2(const Point& A, const Point& B, const Point& C) {
return Cross(B - A, C - A);
}
static double AreaTriangle(const Point& A, const Point& B, const Point& C) {
return Area2(A, B, C) / 2;
}
static Vector Rotate(const Vector& A, double rad) {
return Vector(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));
}
static Vector Normal(const Vector& A) {
double len = Len(A);
return Vector(-A.y / len, A.x / len);
}
static bool Parallel(const Vector& A, const Vector& B) {
return sgn(Cross(A, B)) == 0;
}
};
class Line {
public:
Point p1, p2;
Line() {}
Line(Point p1, Point p2) : p1(p1), p2(p2) {}
Line(Point p, double angle) {
p1 = p;
if (sgn(angle - pi / 2) == 0) {
p2 = p1 + Point(0, 1);
} else {
p2 = p1 + Point(1, tan(angle));
}
}
Line(double a, double b, double c) {
if (sgn(a) == 0) {
p1 = Point(0, -c / b);
p2 = Point(1, -c / b);
} else if (sgn(b) == 0) {
p1 = Point(-c / a, 0);
p2 = Point(-c / a, 1);
} else {
p1 = Point(0, -c / b);
p2 = Point(1, (-c - a) / b);
}
}
int PointLineRelation(const Point& p) const {
int c = sgn(Geometry::Cross(p - p1, p2 - p1));
if (c < 0) return 1; // p 在 v 的右侧
if (c > 0) return 2; // p 在 v 的左侧
return 0; // p 在 v 上
}
bool PointOnSegment(const Point& p) const {
return sgn(Geometry::Cross(p - p1, p2 - p1)) == 0 &&
sgn(Geometry::Dot(p - p1, p - p2)) <= 0;
}
double DistanceToPoint(const Point& p) const {
return fabs(Geometry::Cross(p - p1, p2 - p1)) / Geometry::Distance(p1, p2);
}
Point Projection(const Point& p) const {
double k = Geometry::Dot(p2 - p1, p - p1) / Geometry::Len2(p2 - p1);
return p1 + (p2 - p1) * k;
}
Point SymmetryPoint(const Point& p) const {
Point q = Projection(p);
return Point(2 * q.x - p.x, 2 * q.y - p.y);
}
static int LineRelation(const Line& v1, const Line& v2) {
if (sgn(Geometry::Cross(v1.p2 - v1.p1, v2.p2 - v2.p1)) == 0) {
if (v1.PointLineRelation(v2.p1) == 0) return 1; // 重合
else return 0; // 平行
}
return 2; // 相交
}
static Point CrossPoint(const Point& a, const Point& b, const Point& c, const Point& d) {
double s1 = Geometry::Cross(b - a, c - a);
double s2 = Geometry::Cross(b - a, d - a);
return Point(c.x * s2 - d.x * s1, c.y * s2 - d.y * s1) / (s2 - s1);
}
};
typedef Line Segment;
class SegmentOps {
public:
static bool CrossSegment(const Point& a, const Point& b, const Point& c, const Point& d) {
double c1 = Geometry::Cross(b - a, c - a);
double c2 = Geometry::Cross(b - a, d - a);
double d1 = Geometry::Cross(d - c, a - c);
double d2 = Geometry::Cross(d - c, b - c);
return sgn(c1) * sgn(c2) < 0 && sgn(d1) * sgn(d2) < 0;
}
static double DistanceToPoint(const Point& p, const Segment& v) {
if (sgn(Geometry::Dot(p - v.p1, v.p2 - v.p1)) < 0 || sgn(Geometry::Dot(p - v.p2, v.p1 - v.p2)) < 0)
return std::min(Geometry::Distance(p, v.p1), Geometry::Distance(p, v.p2));
return v.DistanceToPoint(p); // 点的投影在线段上
}
};
class PolygonOps {
public:
static int PointInPolygon(const Point& pt, const std::vector<Point>& p) {
int n = p.size();
for (int i = 0; i < n; i++) {
if (p[i] == pt) return 3; // 点在多边形顶点上
}
for (int i = 0; i < n; i++) {
Line v = Line(p[i], p[(i + 1) % n]);
if (v.PointOnSegment(pt)) return 2; // 点在边上
}
int num = 0;
for (int i = 0; i < n; i++) {
int j = (i + 1) % n;
int c = sgn(Geometry::Cross(pt - p[j], p[i] - p[j]));
int u = sgn(p[i].y - pt.y);
int v = sgn(p[j].y - pt.y);
if (c > 0 && u < 0 && v >= 0) num++;
if (c < 0 && u >= 0 && v < 0) num--;
}
return num != 0; // 1：点在内部；0：点在外部
}
static double Area(const std::vector<Point>& p) {
double area = 0;
int n = p.size();
for (int i = 0; i < n; i++)
area += Geometry::Cross(p[i], p[(i + 1) % n]);
return area / 2; // 面积有正负
}
static Point Center(const std::vector<Point>& p) {
Point ans(0, 0);
if (Area(p) == 0) return ans;
int n = p.size();
for (int i = 0; i < n; i++)
ans = ans + (p[i] + p[(i + 1) % n]) * Geometry::Cross(p[i], p[(i + 1) % n]);
return ans / Area(p) / 6;
}
};