Update Geometry section

c++/Geometry/closest_pair.cpp

@@ -0,0 +1,51 @@
+/*	
+	Recibe una lista ordenada de puntos, orden lexicografico.
+	Retorna el cuadrado de la distancia más corta entre dos puntos (d*d)
+
+	Existe otro algoritmo pero es más complejo de entender.
+*/
+
+long long closest_pair(vector<pt> pts){
+	set<pair<ll, ll>> cur;
+	ll dist = 1e18;
+	int j = 0;
+	for(int i=0; i<(int) pts.size(); i++){
+		ll d = ceil(sqrt(dist));
+		while(j <= i && pts[i].x - pts[j].x >= d){
+			cur.erase({pts[j].y, pts[j].x});
+			j++;
+		}
+		auto it1 = cur.lower_bound({pts[i].y - d, pts[i].x});
+		auto it2 = cur.upper_bound({pts[i].y + d, pts[i].x});
+
+		for(auto it = it1; it != it2; it++){
+			ll dx = pts[i].x - (*it).second;
+			ll dy = pts[i].y - (*it).first;
+			dist = min(dist, dx*dx + dy*dy);
+		}
+		cur.insert({pts[i].y, pts[i].x});
+	}
+	return dist;
+}
+
+T closest_pair_double_ver(vector<pt> &pts){
+	set<pair<T, T>> cur;
+	T dist = 1e18;
+	int j = 0;
+	for(int i=0; i<(int) pts.size(); i++){
+		while(j <= i && (pts[i].x - pts[j].x) > dist-EPS){
+			cur.erase({pts[j].y, pts[j].x});
+			j++;
+		}
+		auto it1 = cur.lower_bound({pts[i].y - dist, pts[i].x});
+		auto it2 = cur.upper_bound({pts[i].y + dist, pts[i].x});
+
+		for(auto it = it1; it != it2; it++){
+			T dx = pts[i].x - (*it).second;
+			T dy = pts[i].y - (*it).first;
+			dist = min(dist, sqrt(dx*dx + dy*dy));
+		}
+		cur.insert({pts[i].y, pts[i].x});
+	}
+	return dist;
+}

c++/Geometry/diameter_polygon.cpp

@@ -0,0 +1,33 @@
+double area(pt &a, pt &b, pt &c) {
+    return abs((a.x - c.x)*(b.y - a.y) - (a.x - b.x)*(c.y - a.y));
+}
+
+//Compute the diameter of a polygon using rotating-caliper method to generate
+//all antipodal pairs. Return the max distance between an antipodal pair extremes.
+
+double diameter(vector<pt> &h) {
+    int m = h.size();
+
+    //Return trivial results for polygons with less than 3 points.
+    if (m <= 1)
+        return 0;
+    if (m == 2)
+        return abs(h[0] - h[1]);
+
+    //Search the point k which is antipodal of point 0
+    int k = 1;
+    while (area(h[m - 1], h[0], h[(k + 1) % m]) > area(h[m - 1], h[0], h[k]))
+        ++k;
+
+    double res = 0;
+    //Generate all antipodal pairs in i,j
+    for (int i = 0, j = k; i <= k && j < m; ++i) {
+        res = max(res, abs(h[i] - h[j]));
+        //Advance point j while j is antipodal pair of i
+        while (j < m && area(h[i], h[(i + 1) % m], h[(j + 1) % m]) > area(h[i], h[(i + 1) % m], h[j])) {
+            res = max(res, abs(h[i] - h[(j + 1) % m]));
+            ++j;
+        }
+    }
+    return res;
+}

c++/Geometry/geo_dlc.cpp

@@ -0,0 +1,94 @@
+struct line {
+	pt v; T c;
+	// From direction vector v and offset c
+	line(pt v, T c) : v(v), c(c) {} 
+	// From equation ax+by=c
+	line(T a, T b, T c) : v({b,-a}), c(c) {}
+	// From points P and Q
+	line(pt p, pt q) : v(q-p), c(cross(v,p)) {}
+	// Will be defined later:
+	// - these work with T = int
+	T side(pt p) {return cross(v,p)-c;} //Algo parecido a orient
+	double dist(pt p) {return abs(side(p)) / abs(v);} //Distancia de P a la linea
+	line perpThrough(pt p) {return {p, p + perp(v)};} //Linea perpendicular a l
+
+	bool cmpProj(pt p, pt q) { //Comparador para ordenar puntos sobre la linea
+		return dot(v,p) < dot(v,q);
+	}
+
+	line translate(pt t) {return {v, c + cross(v,t)};} //Traslada la linea con base en un vector
+	line shiftLeft(double dist) {return {v, c + dist*abs(v)};} //Mueve la linea perpendicularmente una distancia
+	// - these require T = double
+	pt proj(pt p) {return p - perp(v)*side(p)/sq(v);} //Retorna el punto mas cercano de la linea a p
+	pt refl(pt p) {return p - perp(v)*2*side(p)/sq(v);} //Retorna el punto reflejo de p con respecto a la linea
+
+};
+
+bool inter(line l1, line l2, pt &out) { 
+	T d = cross(l1.v, l2.v);
+	if (d == 0) return false;
+	out = (l2.v*l1.c - l1.v*l2.c) / d; // requires floating-point coordinates
+	return true;
+	//Retorna true si se pudo hallar un punto de cruce
+	//Determina en que punto se cruzan dos lineas
+}
+
+
+line bisector(line l1, line l2, bool interior) {
+	assert(cross(l1.v, l2.v) != 0); // l1 and l2 cannot be parallel!
+	double sign = interior ? 1 : -1;
+	return {l2.v/abs(l2.v) + l1.v/abs(l1.v) * sign,
+	l2.c/abs(l2.v) + l1.c/abs(l1.v) * sign};
+	//Retorna una linea que forma angulos iguales con las otras dos lines
+	//La interna es la que esta entre el vector l1 y l2
+	//La externa es la otra
+}
+
+double segPoint(pt a, pt b, pt p, pt &ans) {
+	if (a != b) {
+		line l(a,b);
+		if (l.cmpProj(a,p) && l.cmpProj(p,b)){// if closest to projection
+			ans = l.proj(p);
+			return l.dist(p); // output distance to line
+		} 		
+	}
+	if(abs(p-a) < abs(p-b)) ans = a;
+	else ans = b;
+	return min(abs(p-a), abs(p-b)); // otherwise distance to A or B
+}
+
+double segSeg(pt a, pt b, pt c, pt d) {
+	pt ax;
+	return min({segPoint(a,b,c,ax), segPoint(a,b,d,ax),
+	segPoint(c,d,a,ax), segPoint(c,d,b,ax)});
+}
+
+bool inDisk(pt a, pt b, pt p) {
+	return dot(a-p, b-p) <= 0;
+}
+
+bool onSegment(pt a, pt b, pt p) {
+	return orient(a,b,p) == 0 && inDisk(a,b,p);
+}
+
+bool properInter(pt a, pt b, pt c, pt d, pt &out) {
+	double oa = orient(c,d,a), // basado en la orientacion sabemos si se cruzan
+	ob = orient(c,d,b),
+	oc = orient(a,b,c),
+	od = orient(a,b,d);
+	// Proper intersection exists iff opposite signs
+	if (oa*ob < 0 && oc*od < 0) {
+		out = (a*ob - b*oa) / (ob-oa);
+		return true;
+	}
+	return false;
+}
+
+//Circulo dado tres puntos
+//Es la interseccion del segmento perpendicular
+//al segmento AB y del segmento perpendicular a AC
+pt circumCenter(pt a, pt b, pt c) {
+	b = b-a, c = c-a; // consider coordinates relative to A
+	assert(cross(b,c) != 0); // no circumcircle if A,B,C aligned
+	return a + perp(b*sq(c) - c*sq(b))/cross(b,c)/2;
+}

c++/Geometry/inPolygon.cpp

@@ -0,0 +1,18 @@
+bool above(pt a, pt p) {
+	return p.y >= a.y;
+}
+// check if [PQ] crosses ray from A
+bool crossesRay(pt a, pt p, pt q) {
+	return (above(a,q) - above(a,p)) * orient(a,p,q) > 0;
+}
+bool inPolygon(vector<pt> p, pt a, bool strict = true) {
+	int numCrossings = 0;
+	for (int i = 0, n = p.size(); i < n; i++) {
+		if (onSegment(p[i], p[(i+1)%n], a))
+			return !strict; // strict es para determinar si esta en el borde o no
+								 // (si para el problema estar en el borde es estar dentro
+								 // return strict, sino return !strict)
+		numCrossings += crossesRay(a, p[i], p[(i+1)%n]);
+	}
+	return numCrossings & 1;
+}

c++/Geometry/minkowski_sum.cpp

@@ -0,0 +1,28 @@
+void ordMinY(vector<pt> &pts){
+	int pos = 0;
+	for(int i=0; i<(int)pts.size(); i++){
+		if(pts[i].y < pts[pos].y || (pts[i].y == pts[pos].y && pts[i].x < pts[pos].x)){
+			pos = i;
+		}
+	}
+	rotate(pts.begin(), pts.begin() + pos, pts.end());
+}
+
+vector<pt> minkowski(vector<pt> a, vector<pt> b){
+	ordMinY(a); ordMinY(b);
+	vector<pt> sum;
+	a.push_back(a[0]); a.push_back(a[1]);
+	b.push_back(b[0]); b.push_back(b[1]);
+	int i = 0, j = 0;
+	while(i < (int)a.size() - 2 || j < (int)b.size() - 2){
+		sum.push_back(a[i] + b[j]);
+		ll crp = cross(a[i+1] - a[i], b[j+1] - b[j]);
+		if(crp >= 0 && i < (int)a.size() - 2){
+			i++;
+		} 
+		if(crp <= 0 && j < (int)b.size() - 2){
+			j++;
+		} 
+	}
+	return sum;
+}

c++/Geometry/pointInConvexPolygon.cpp

@@ -0,0 +1,31 @@
+bool pointInTriangle(pt a, pt b, pt c, pt p){
+	ll s1 = abs(cross(b - a, c - a));
+	ll s2 = abs(cross(b - a, p - a)) + abs(cross(p - a, c - a)) + abs(cross(b - p, c - p));
+	return s1 == s2;
+}
+
+bool pointInConvexPolygon(vector<pt> &pts, pt p){
+	//Podrias mover todos los puntos con respecto a 
+	//pts[0] antes de hacer esto, eso reduce el codigo
+	int n = pts.size();
+	if(cross(pts[1] - pts[0], p - pts[0]) >= 0 && cross(p - pts[0], pts[n-1] - pts[0]) >= 0){
+		if(cross(pts[1] - pts[0], p - pts[0]) == 0){
+			return sq(p - pts[0]) <= sq(pts[1] - pts[0]);
+		}
+		if(cross(pts[n-1] - pts[0], p - pts[0]) == 0){
+			return sq(p - pts[0]) <= sq(pts[n-1] - pts[0]);
+		}	
+		int l = 1, r = n-1;
+		while(l <= r){
+			int mid = (l+r)/2;
+			if(cross(pts[mid] - pts[0], p - pts[0]) >= 0){
+				l = mid + 1;
+			}else{
+				r = mid - 1;
+			}
+		}
+		if(r == 0) return false;
+		return pointInTriangle(pts[0], pts[r], pts[r+1], p); 
+	}
+	return false;
+}