Geometry Template (By Arpa) - YessineJallouli/Competitive-Programming GitHub Wiki

// warning : Always use fabs (not abs) to compare
double INF = 1e100;
double EPS = 1e-6;

struct pt { 
	double x, y; 
	pt() {}
	pt(double x, double y) : x(x), y(y) {}
	pt(const pt &p) : x(p.x), y(p.y)    {}
	pt operator + (const pt &p)  const { return pt(x+p.x, y+p.y); }
	pt operator - (const pt &p)  const { return pt(x-p.x, y-p.y); }
	pt operator * (double c)     const { return pt(x*c,   y*c  ); }
	pt operator / (double c)     const { return pt(x/c,   y/c  ); }
};

double dot(pt p, pt q)     { return p.x*q.x+p.y*q.y; }
double dist2(pt p, pt q)   { return dot(p-q,p-q); }
double cross(pt p, pt q)   { return p.x*q.y-p.y*q.x; }
ostream &operator<<(ostream &os, const pt &p) {
	return os << "(" << p.x << "," << p.y << ")"; 
}

// rotate a point CCW or CW around the origin
pt RotateCCW90(pt p)   { return pt(-p.y,p.x); }
pt RotateCW90(pt p)    { return pt(p.y,-p.x); }
pt RotateCCW(pt p, double t) { 
	return pt(p.x*cos(t)-p.y*sin(t), p.x*sin(t)+p.y*cos(t)); 
}

// sort around origin
// use this : sort(pts.begin(), pts.end(), AngleCmp);
bool AngleCmp(const pt &a, const pt &b) {
    double A = atan2(a.y, a.x);
    double B = atan2(b.y, b.x);
    if (fabs(A - B) > EPS) return A < B;       // smaller angle first
    return dist2(a, pt(0,0)) < dist2(b, pt(0,0)); // tie-breaker
}

// usage
// PIVOT = pts[0];        // or any pivot
// sort(pts.begin(), pts.end(), AngleCmpPivot);
pt PIVOT;
bool AngleCmpPivot(const pt &a, const pt &b) {
    double A = atan2(a.y - PIVOT.y, a.x - PIVOT.x);
    double B = atan2(b.y - PIVOT.y, b.x - PIVOT.x);
    if (fabs(A - B) > EPS) return A < B;
    return dist2(a, PIVOT) < dist2(b, PIVOT); // closer to pivot first
}

// project point c onto line through a and b
// assuming a != b
pt ProjectPointLine(pt a, pt b, pt c) {
	return a + (b-a)*dot(c-a, b-a)/dot(b-a, b-a);
}

// project point c onto line segment through a and b
pt ProjectPointSegment(pt a, pt b, pt c) {
	double r = dot(b-a,b-a);
	if (fabs(r) < EPS) return a;
	r = dot(c-a, b-a)/r;
	if (r < 0) return a;
	if (r > 1) return b;
	return a + (b-a)*r;
}

// compute distance from c to segment between a and b
double DistancePointSegment(pt a, pt b, pt c) {
	return sqrt(dist2(c, ProjectPointSegment(a, b, c)));
}

// compute distance between point (x,y,z) and plane ax+by+cz=d
double DistancePointPlane(double x, double y, double z,
                          double a, double b, double c, double d)
{
	return fabs(a*x+b*y+c*z-d)/sqrt(a*a+b*b+c*c);
}

// reflect p on the line (l1,l2)
pt reflect(pt p, pt l1, pt l2) {
    pt z = p - l1, w = l2 - l1;
    double w2 = dot(w, w);
    double a = (z.x * w.x + z.y * w.y) / w2;
    double b = (z.y * w.x - z.x * w.y) / w2;
    return l1 + pt(a*w.x + b*w.y, a*w.y - b*w.x);
}

// returns a point D on the bisector of angle BAC
pt Bisector(pt a, pt b, pt c) {
    pt u = b - a;
    pt v = c - a;
    double lu = sqrt(dot(u, u));
    double lv = sqrt(dot(v, v));
    if (lu < EPS || lv < EPS) return a;  // degenerate
    u = u / lu;   // make them unit vectors
    v = v / lv;
    // direction of angle bisector = u + v
    return a + u + v;   // this point D lies on the bisector
}

// determine if lines from a to b and c to d are parallel or collinear
bool LinesParallel(pt a, pt b, pt c, pt d) { 
	return fabs(cross(b-a, c-d)) < EPS; 
}

bool LinesCollinear(pt a, pt b, pt c, pt d) { 
	return LinesParallel(a, b, c, d)
		&& fabs(cross(a-b, a-c)) < EPS
					   && fabs(cross(c-d, c-a)) < EPS; 
}

// determine if line segment from a to b intersects with 
// line segment from c to d
bool SegmentsIntersect(pt a, pt b, pt c, pt d) {
	if (LinesCollinear(a, b, c, d)) {
		if (dist2(a, c) < EPS || dist2(a, d) < EPS ||
		    dist2(b, c) < EPS || dist2(b, d) < EPS) return true;
		if (dot(c-a, c-b) > 0 && dot(d-a, d-b) > 0 && dot(c-b, d-b) > 0)
			return false;
		return true;
	}
	if (cross(d-a, b-a) * cross(c-a, b-a) > 0) return false;
	if (cross(a-c, d-c) * cross(b-c, d-c) > 0) return false;
	return true;
}

// compute intersection of line passing through a and b
// with line passing through c and d, assuming that unique
// intersection exists; for segment intersection, check if
// segments intersect first
pt ComputeLineIntersection(pt a, pt b, pt c, pt d) {
	b=b-a; d=c-d; c=c-a;
	assert(dot(b, b) > EPS && dot(d, d) > EPS);
	return a + b*cross(c, d)/cross(b, d);
}

// compute center of circle given three points
pt ComputeCircleCenter(pt a, pt b, pt c) {
	b=(a+b)/2;
	c=(a+c)/2;
	return ComputeLineIntersection(b, b+RotateCW90(a-b), c, c+RotateCW90(a-c));
}

// determine if point is in a possibly non-convex polygon (by William
// Randolph Franklin); returns 1 for strictly interior points, 0 for
// strictly exterior points, and 0 or 1 for the remaining points.
// Note that it is possible to convert this into an *exact* test using
// integer arithmetic by taking care of the division appropriately
// (making sure to deal with signs properly) and then by writing exact
// tests for checking point on polygon boundary
bool PointInPolygon(const vector<pt> &p, pt q) {
	bool c = 0;
	for (int i = 0; i < p.size(); i++){
		int j = (i+1)%p.size();
		if ((p[i].y <= q.y && q.y < p[j].y || 
		     p[j].y <= q.y && q.y < p[i].y) &&
		    q.x < p[i].x + (p[j].x - p[i].x) * (q.y - p[i].y) / (p[j].y - p[i].y))
			c = !c;
	}
	return c;
}

// determine if point is on the boundary of a polygon
bool PointOnPolygon(const vector<pt> &p, pt q) {
	for (int i = 0; i < p.size(); i++)
		if (dist2(ProjectPointSegment(p[i], p[(i+1)%p.size()], q), q) < EPS)
			return true;
	return false;
}

// compute intersection of line through points a and b with
// circle centered at c with radius r > 0
// going from a to b, t[1] is the first intersection and t[0] is the second
vector<pt> CircleLineIntersection(pt a, pt b, pt c, double r) {
	vector<pt> ret;
	b = b-a;
	a = a-c;
	double A = dot(b, b);
	double B = dot(a, b);
	double C = dot(a, a) - r*r;
	double D = B*B - A*C;
	if (D < -EPS) return ret;
	ret.push_back(c+a+b*(-B+sqrt(D+EPS))/A);
	if (D > EPS)
		ret.push_back(c+a+b*(-B-sqrt(D))/A);
	return ret;
}

// compute intersection of circle centered at a with radius r
// with circle centered at b with radius R
// order is counter clock wise
vector<pt> CircleCircleIntersection(pt a, pt b, double r, double R) {
	vector<pt> ret;
	double d = sqrt(dist2(a, b));
	if (d > r+R || d+min(r, R) < max(r, R)) return ret;
	double x = (d*d-R*R+r*r)/(2*d);
	double y = sqrt(r*r-x*x);
	pt v = (b-a)/d;
	ret.push_back(a+v*x + RotateCCW90(v)*y);
	if (y > 0)
		ret.push_back(a+v*x - RotateCCW90(v)*y);
	return ret;
}

// This code computes the area or centroid of a (possibly nonconvex)
// polygon, assuming that the coordinates are listed in a clockwise or
// counterclockwise fashion.  Note that the centroid is often known as
// the "center of gravity" or "center of mass".

// use ll for this 
ll ComputeSignedArea(const vector<pt> &p) {
	ll area = 0;
	for(int i = 0; i < p.size(); i++) {
		int j = (i+1) % p.size();
		area += p[i].x*p[j].y - p[j].x*p[i].y;
	}
	return area;
}

ll ComputeArea(const vector<pt> &p) {
	return abs(ComputeSignedArea(p));
}

pt ComputeCentroid(const vector<pt> &p) {
	pt c(0,0);
	double scale = 6.0 * ComputeSignedArea(p);
	for (int i = 0; i < p.size(); i++){
		int j = (i+1) % p.size();
		c = c + (p[i]+p[j])*(p[i].x*p[j].y - p[j].x*p[i].y);
	}
	return c / scale;
}

// tests whether or not a given polygon (in CW or CCW order) is simple
bool IsSimple(const vector<pt> &p) {
	for (int i = 0; i < p.size(); i++) {
		for (int k = i+1; k < p.size(); k++) {
			int j = (i+1) % p.size();
			int l = (k+1) % p.size();
			if (i == l || j == k) continue;
			if (SegmentsIntersect(p[i], p[j], p[k], p[l])) 
				return false;
		}
	}
	return true;
}

// finds the position of b on the segment from a to c.
double coefOnLine(pt a, pt b, pt c){
	if(abs(a.x - c.x) < EPS)
		return (b.y - a.y) / (c.y - a.y);
	return (b.x - a.x) / (c.x - a.x);
}
void Union(vector<pair<double, double >> &segs){
	sort(segs.begin(), segs.end());
	int sz = 0;
	for(auto [l, r] : segs)
		if(l <= r)
			if(!sz || l > segs[sz - 1].second + EPS)
				segs[sz++] = {l, r};
			else
				segs[sz - 1].second = max(segs[sz - 1].second, r);
	segs.resize(sz);
}
vector<pair<double, double > > PolygonSegmentIntersection(vector<pt> &pol, pt a, pt b){
	vector<pair<double, double> > segs;
	vector<pt> impos({a, b});
	for(int k = 0; k < pol.size(); k++)
		if(SegmentsIntersect(a, b, pol[k], pol[(k + 1) % pol.size()]))
			impos.push_back(ComputeLineIntersection(a, b, pol[k], pol[(k + 1) % pol.size()]));
	sort(impos.begin(), impos.end(), [&](pt x, pt y){
			return coefOnLine(a, x, b) < coefOnLine(a, y, b);
		});
	for(int k = 0; k < impos.size() - 1; k++) {
		pt mid = (impos[k] + impos[k + 1]) / 2;
		if(PointInPolygon(pol, mid))
			segs.emplace_back(coefOnLine(a, impos[k], b), coefOnLine(a, impos[k + 1], b));
	}
	return segs;
}
pair<double, double> CircleSegmentIntersection(pt a, pt b, pt c, double r) {
	vector<pt> ret = CircleLineIntersection(a, b, c, r);
	if(ret.size() < 2)
		return {0, 0};
	return {max<double>(0, min(coefOnLine(a, ret[0], b), coefOnLine(a, ret[1], b))),
			min<double>(1, max(coefOnLine(a, ret[0], b), coefOnLine(a, ret[1], b)))};
}

bool cmpx(const pt& a, const pt& b) {
    if (a.x == b.x) return a.y < b.y;
    return a.x < b.x;
}

bool cmpy(const pt& a, const pt& b) {
    if (a.y == b.y) return a.x < b.x;
    return a.y < b.y;
}

// change everything to long long
ll closest_pair_rec(int l, int r, vector<pt>& a, vector<pt>& tmp) {
    if (r - l <= 3) {
        ll d = INF;
        for (int i = l; i < r; ++i)
            for (int j = i + 1; j < r; ++j)
                d = min<ll>(d, dist2(a[i], a[j]));
        sort(a.begin() + l, a.begin() + r, cmpy); // sort by y in this segment
        return d;
    }
    int m = (l + r) >> 1;
    ll midx = a[m].x;
    ll d = min(
        closest_pair_rec(l, m, a, tmp),
        closest_pair_rec(m, r, a, tmp)
    );
    // merge by y-coordinate
    merge(a.begin() + l, a.begin() + m,
          a.begin() + m, a.begin() + r,
          tmp.begin(), cmpy);
    copy(tmp.begin(), tmp.begin() + (r - l), a.begin() + l);
    // build vertical strip |x - midx|^2 < d
    vector<pt> strip;
    strip.reserve(r - l);
    for (int i = l; i < r; ++i) {
        ll dx = a[i].x - midx;
        if (dx*dx < d) strip.push_back(a[i]);
    }
    // only check next points with small y-diff
    for (int i = 0; i < (int)strip.size(); ++i) {
        for (int j = i + 1; j < (int)strip.size(); ++j) {
            ll dy = strip[j].y - strip[i].y;
            if (dy*dy >= d) break;
            d = min<ll>(d, dist2(strip[i], strip[j]));
        }
    }
    return d;
}

ll closest_pair(vector<pt>& a) {
    int n = (int)a.size();
    if (n < 2) return 0; // or INF, depending on problem statement
    sort(a.begin(), a.end(), cmpx);   // sort by x
    vector<pt> tmp(n);
    return closest_pair_rec(0, n, a, tmp);
}

double orient(const pt &a, const pt &b, const pt &c) {
    // >0 : left turn (CCW), <0 : right turn (CW), 0 : collinear
    return cross(b - a, c - a);
}
 
// Monotone chain convex hull
// - returns hull in CCW order, without repeating the first point
// - if include_collinear = true, keeps collinear boundary points
vector<pt> ConvexHull(vector<pt> pts, bool include_collinear = false) {
    int n = (int)pts.size();
    if (n <= 1) return pts;
 
    // sort by x, then y
    sort(pts.begin(), pts.end(), [](const pt &a, const pt &b) {
        if (fabs(a.x - b.x) > EPS) return a.x < b.x;
        return a.y < b.y;
    });
 
    // remove duplicates
    pts.erase(unique(pts.begin(), pts.end(), [](const pt &a, const pt &b) {
        return fabs(a.x - b.x) < EPS && fabs(a.y - b.y) < EPS;
    }), pts.end());
 
    if ((int)pts.size() <= 1) return pts;
 
    vector<pt> lower, upper;
    // build lower hull
    for (const pt &p : pts) {
        while ((int)lower.size() >= 2) {
            pt a = lower[(int)lower.size() - 2];
            pt b = lower[(int)lower.size() - 1];
            double cr = orient(a, b, p);
            if (include_collinear) {
                if (cr < -EPS) lower.pop_back();
                else break;
            } else {
                if (cr <= EPS) lower.pop_back();
                else break;
            }
        }
        lower.push_back(p);
    }
    for (int i = (int)pts.size() - 1; i >= 0; --i) {
        const pt &p = pts[i];
        while ((int)upper.size() >= 2) {
            pt a = upper[(int)upper.size() - 2];
            pt b = upper[(int)upper.size() - 1];
            double cr = orient(a, b, p);
            if (include_collinear) {
                if (cr < -EPS) upper.pop_back();
                else break;
            } else {
                if (cr <= EPS) upper.pop_back();
                else break;
            }
        }
        upper.push_back(p);
    }
    lower.pop_back();
    upper.pop_back();
    vector<pt> hull = lower;
    hull.insert(hull.end(), upper.begin(), upper.end());
 
    return hull;
}

struct Circle {
    pt c;
    double r;
    Circle() : c(0, 0), r(0) {}
    Circle(const pt &c, double r) : c(c), r(r) {}
};

bool InsideCircle(const Circle &C, const pt &p) {
    return dist2(C.c, p) <= (C.r + EPS) * (C.r + EPS);
}

Circle CircleFromTwoPoints(const pt &a, const pt &b) {
    pt c = (a + b) / 2.0;
    double r = sqrt(dist2(a, b)) / 2.0;
    return Circle(c, r);
}

Circle CircleFromThreePoints(const pt &a, const pt &b, const pt &c) {
    pt center = ComputeCircleCenter(a, b, c);
    double r = sqrt(dist2(center, a));
    return Circle(center, r);
}

// O(n)
Circle MinimumEnclosingCircle(vector<pt> pts) {
    int n = (int)pts.size();
    if (n == 0) return Circle(pt(0, 0), 0.0);
    if (n == 1) return Circle(pts[0], 0.0);

    static mt19937_64 mt(chrono::steady_clock::now().time_since_epoch().count());
    shuffle(pts.begin(), pts.end(), mt);

    Circle C(pts[0], 0.0);
    for (int i = 1; i < n; ++i) {
        if (InsideCircle(C, pts[i])) continue;

        C = Circle(pts[i], 0.0);
        for (int j = 0; j < i; ++j) {
            if (InsideCircle(C, pts[j])) continue;

            C = CircleFromTwoPoints(pts[i], pts[j]);
            for (int k = 0; k < j; ++k) {
                if (InsideCircle(C, pts[k])) continue;

                C = CircleFromThreePoints(pts[i], pts[j], pts[k]);
            }
        }
    }
    return C;
}