merged

6 files changed, 229 insertions, 21 deletions

diff --git a/geometry/delaunay.cpp b/geometry/delaunay.cpp
new file mode 100644
index 0000000..1008b39
--- /dev/null
+++ b/geometry/delaunay.cpp
@@ -0,0 +1,124 @@
+using lll = __int128;
+using pt = complex<ll>;
+
+constexpr pt INF_PT = pt(1e18, 1e18);
+
+bool circ(pt p, pt a, pt b, pt c) {// p in circle(A,B,C)
+	return imag((c-b)*conj(p-c)*(a-p)*conj(b-a)) < 0;
+}
+
+struct QuadEdge {
+	QuadEdge* rot = nullptr;
+	QuadEdge* onext = nullptr;
+	pt orig = INF_PT;
+	bool used = false;
+	QuadEdge* rev() const {return rot->rot;}
+	QuadEdge* lnext() const {return rot->rev()->onext->rot;}
+	QuadEdge* oprev() const {return rot->onext->rot;}
+	pt dest() const {return rev()->orig;}
+};
+
+deque<QuadEdge> edgeData;
+
+QuadEdge* makeEdge(pt from, pt to) {
+	for (int j : {0,1,2,3}) edgeData.push_back({});
+	auto e = edgeData.end() - 4;
+	for (int j : {0,1,2,3}) e[j].onext = e[j^3].rot = &e[j^(j>>1)];
+	e[0].orig = from;
+	e[1].orig = to;
+	return &e[0];
+}
+
+void splice(QuadEdge* a, QuadEdge* b) {
+	swap(a->onext->rot->onext, b->onext->rot->onext);
+	swap(a->onext, b->onext);
+}
+
+QuadEdge* connect(QuadEdge* a, QuadEdge* b) {
+	QuadEdge* e = makeEdge(a->dest(), b->orig);
+	splice(e, a->lnext());
+	splice(e->rev(), b);
+	return e;
+}
+
+bool valid(QuadEdge* e, QuadEdge* base) {
+	return cross(e->dest(), base->orig, base->dest()) < 0;
+}
+
+template<bool ccw>
+QuadEdge* deleteAll(QuadEdge* e, QuadEdge* base) {
+	if (valid(e, base)) {
+		while (circ(base->dest(), base->orig, e->dest(), (ccw ? e->onext : e->oprev())->dest())) {
+			QuadEdge* t = ccw ? e->onext : e->oprev();
+			splice(e, e->oprev());
+			splice(e->rev(), e->rev()->oprev());
+			e = t;
+	}}
+	return e;
+}
+
+template<typename IT>
+pair<QuadEdge*, QuadEdge*> rec(IT l, IT r) {
+	int n = distance(l, r);
+	if (n <= 3) {
+		QuadEdge* a = makeEdge(l[0], l[1]);
+		if (n == 2) return {a, a->rev()};
+		QuadEdge* b = makeEdge(l[1], l[2]);
+		splice(a->rev(), b);
+		int side = cross(l[0], l[1], l[2]);
+		QuadEdge* c = nullptr;
+		if (side != 0) c = connect(b, a);
+		if (side >= 0) return {a, b->rev()};
+		else return {c->rev(), c};
+	}
+	auto m = l + (n / 2);
+	auto [ldo, ldi] = rec(l, m);
+	auto [rdi, rdo] = rec(m, r);
+	while (true) {
+		if (cross(rdi->orig, ldi->orig, ldi->dest()) > 0) {
+			ldi = ldi->lnext();
+		} else if (cross(ldi->orig, rdi->orig, rdi->dest()) < 0) {
+			rdi = rdi->rev()->onext;
+		} else break;
+	}
+	QuadEdge* base = connect(rdi->rev(), ldi);
+	if (ldi->orig == ldo->orig)	ldo = base->rev();
+	if (rdi->orig == rdo->orig)	rdo = base;
+	while (true) {
+		QuadEdge* lcand = deleteAll<true>(base->rev()->onext, base);
+		QuadEdge* rcand = deleteAll<false>(base->oprev(), base);
+		if (!valid(lcand, base) && !valid(rcand, base))	break;
+		if (!valid(lcand, base) || (valid(rcand, base) &&
+		    circ(lcand->dest(), lcand->orig, rcand->orig, rcand->dest()))) {
+			base = connect(rcand, base->rev());
+		} else {
+			base = connect(base->rev(), lcand->rev());
+	}}
+	return {ldo, rdo};
+}
+
+vector<pt> delaunay(vector<pt> pts) {
+	if (sz(pts) <= 2) return {};
+	sort(all(pts), [](const pt& a, const pt& b) {
+		if (real(a) != real(b)) return real(a) < real(b);
+		return imag(a) < imag(b);
+	});
+	QuadEdge* r = rec(all(pts)).first;
+	vector<QuadEdge*> edges = {r};
+	while (cross(r->onext->dest(), r->dest(), r->orig) < 0) r = r->onext;
+	auto add = [&](QuadEdge* e){
+		QuadEdge* cur = e;
+		do {
+			cur->used = true;
+			pts.push_back(cur->orig);
+			edges.push_back(cur->rev());
+			cur = cur->lnext();
+		} while (cur != e);
+	};
+	add(r);
+	pts.clear();
+	for (int i = 0; i < sz(edges); i++) {
+		if (!edges[i]->used) add(edges[i]);
+	}
+	return pts;
+}
diff --git a/geometry/geometry.tex b/geometry/geometry.tex
index a027de4..d92bad1 100644
--- a/geometry/geometry.tex
+++ b/geometry/geometry.tex
@@ -7,14 +7,6 @@
 	\sourcecode{geometry/closestPair.cpp}
 \end{algorithm}
 
-\begin{algorithm}{Rotating calipers}
-	\begin{methods}
-		\method{antipodalPoints}{berechnet antipodale Punkte}{n}
-	\end{methods}
-	\textbf{WICHTIG:} Punkte müssen gegen den Uhrzeigersinn Sortiert sein und konvexes Polygon bilden!
-	\sourcecode{geometry/antipodalPoints.cpp}
-\end{algorithm}
-
 \begin{algorithm}{Konvexe Hülle}
 	\begin{methods}
 		\method{convexHull}{berechnet Konvexehülle}{n\*\log(n)}
@@ -27,23 +19,43 @@
 	\sourcecode{geometry/convexHull.cpp}
 \end{algorithm}
 
+\begin{algorithm}{Rotating calipers}
+	\begin{methods}
+		\method{antipodalPoints}{berechnet antipodale Punkte}{n}
+	\end{methods}
+	\textbf{WICHTIG:} Punkte müssen gegen den Uhrzeigersinn Sortiert sein und konvexes Polygon bilden!
+	\sourcecode{geometry/antipodalPoints.cpp}
+\end{algorithm}
+
 \subsection{Formeln~~--~\texttt{std::complex}}
 \sourcecode{geometry/formulars.cpp}
 \sourcecode{geometry/linesAndSegments.cpp}
+\sourcecode{geometry/sortAround.cpp}
 \input{geometry/triangle}
 \sourcecode{geometry/triangle.cpp}
 \sourcecode{geometry/polygon.cpp}
-\sourcecode{geometry/sortAround.cpp}
 \sourcecode{geometry/circle.cpp}
 
 \subsection{Formeln - 3D}
 \sourcecode{geometry/formulars3d.cpp}
 
 \optional{
-\subsection{3D-Kugeln}
-\sourcecode{geometry/spheres.cpp}
+	\subsection{3D-Kugeln}
+	\sourcecode{geometry/spheres.cpp}
 }
 
+\begin{algorithm}{Half-plane intersection}
+	\sourcecode{geometry/hpi.cpp}
+\end{algorithm}
+
+\begin{algorithm}[optional]{Delaunay Triangulierung}
+	\begin{methods}
+		\method{delaunay}{berechnet Triangulierung}{n\*\log(n)}
+	\end{methods}
+	\textbf{WICHTIG:} Wenn alle Punkte kollinear sind gibt es keine Traingulierung! Wenn 4 Punkte auf einem Kreis liegen ist die Triangulierung nicht eindeutig.
+	\sourcecode{geometry/delaunay.cpp}
+\end{algorithm}
+
 \optional{
 \subsection{Geraden}
 \sourcecode{geometry/lines.cpp}
diff --git a/geometry/hpi.cpp b/geometry/hpi.cpp
new file mode 100644
index 0000000..3509e0e
--- /dev/null
+++ b/geometry/hpi.cpp
@@ -0,0 +1,68 @@
+constexpr ll inf = 0x1FFF'FFFF'FFFF'FFFF;//THIS CODE IS WIP
+
+bool left(pt p) {return real(p) < 0 || 
+					   (real(p) == 0 && imag(p) < 0);}
+struct hp {
+	pt from, to;
+
+	hp(pt a, pt b) : from(a), to(b) {}
+	hp(pt dummy) : hp(dummy, dummy) {}
+
+	bool dummy() const {return from == to;}
+	pt dir() const {return dummy() ? to : to - from;}
+	bool operator<(const hp& o) const {
+		if (left(dir()) != left(o.dir())) 
+			return left(dir()) > left(o.dir());
+		return cross(dir(), o.dir()) > 0;
+	}
+
+	using lll = __int128;
+	using ptl = complex<lll>;
+	ptl mul(lll m, ptl p) const {return m*p;}//ensure 128bit
+
+	bool check(const hp& a, const hp& b) const {
+		if (dummy() || b.dummy()) return false;
+		if (a.dummy()) {
+			ll ort = sgn(cross(b.dir(), dir()));
+			if (ort == 0) return cross(from, to, a.from) < 0;
+			return cross(b.dir(), a.dir()) * ort > 0;
+		}
+		ll y = cross(a.dir(), b.dir());
+		ll z = cross(b.from - a.from, b.dir());
+		ptl i = mul(y, a.from) + mul(z, a.dir()); //intersect a and b
+		// check if i is outside/right of x
+		return imag(conj(mul(sgn(y),dir()))*(i-mul(y,from))) < 0;
+	}
+};
+
+constexpr ll lim = 2e9+7;
+
+deque<hp> intersect(vector<hp> hps) {
+	hps.push_back(hp(pt{lim+1,-1}));
+	hps.push_back(hp(pt{lim+1,1}));
+	sort(all(hps));
+
+	deque<hp> dq = {hp(pt{-lim, 1})};
+	for (auto x : hps) {
+		while (sz(dq) > 1 && x.check(dq.end()[-1], dq.end()[-2]))
+			dq.pop_back();
+		while (sz(dq) > 1 && x.check(dq[0], dq[1]))
+			dq.pop_front();
+
+		if (cross(x.dir(), dq.back().dir()) == 0) {
+			if (dot(x.dir(), dq.back().dir()) < 0) return {};
+			if (cross(x.from, x.to, dq.back().from) < 0)
+				dq.pop_back();
+			else continue;
+		}
+		dq.push_back(x);
+	}
+
+	while (sz(dq) > 2 && dq[0].check(dq.end()[-1], dq.end()[-2]))
+		dq.pop_back();
+	while (sz(dq) > 2 && dq.end()[-1].check(dq[0], dq[1]))
+		dq.pop_front();
+
+	if (sz(dq) < 3) return {};
+	return dq;
+}
diff --git a/geometry/linesAndSegments.cpp b/geometry/linesAndSegments.cpp
index ba6c468..c86b331 100644
--- a/geometry/linesAndSegments.cpp
+++ b/geometry/linesAndSegments.cpp
@@ -17,8 +17,7 @@ vector<pt> lineSegmentIntersection(pt p0, pt p1, pt p2, pt p3) {
 	double b = cross(p2 - p0, p3 - p2);
 	double c = cross(p1 - p0, p0 - p2);
 	if (a < 0) {a = -a; b = -b; c = -c;}
-	if (b < -EPS || a-b < -EPS ||
-	    c < -EPS || a-c < -EPS) return {};
+	if (b < -EPS || b-a > EPS || c < -EPS || c-a > EPS) return {};
 	if (a > EPS) return {p0 + b/a*(p1 - p0)};
 	vector<pt> result;
 	auto insertUnique = [&](pt p) {
diff --git a/geometry/polygon.cpp b/geometry/polygon.cpp
index f562b7a..409f765 100644
--- a/geometry/polygon.cpp
+++ b/geometry/polygon.cpp
@@ -10,13 +10,13 @@ double area(const vector<pt>& poly) { //poly[0] == poly.back()
 // Anzahl drehungen einer Polyline um einen Punkt
 // p nicht auf rand und poly[0] == poly.back()
 // res != 0 or (res & 1) != 0 um inside zu prüfen bei
-// selbstschneidenden polygonen (definitions sache)
+// selbstschneidenden Polygonen (definitions Sache)
 ll windingNumber(pt p, const vector<pt>& poly) {
 	ll res = 0;
 	for (int i = 0; i + 1 < sz(poly); i++) {
 		pt a = poly[i], b = poly[i + 1];
-		if (real(a) > real(b)) swap(a,b);
-		if (real(a) <= real(p) &&real(p) < real(b) &&
+		if (real(a) > real(b)) swap(a, b);
+		if (real(a) <= real(p) && real(p) < real(b) &&
 		    cross(p, a, b) < 0) {
 			res += orientation(p, poly[i], poly[i + 1]);
 	}}
diff --git a/geometry/triangle.cpp b/geometry/triangle.cpp
index a00eb56..ebbbc88 100644
--- a/geometry/triangle.cpp
+++ b/geometry/triangle.cpp
@@ -28,12 +28,17 @@ pt outCenter(pt a, pt b, pt c) {
 	        c*conj(c)*conj(a-b)) / d;
 }
 
+//  1 => p außerhalb Kreis durch a,b,c
+//  0 => p auf Kreis durch a,b,c
+// -1 => p im Kreis durch a,b,c
+int insideOutCenter(pt a, pt b, pt c, pt p) {// braucht lll
+	return sgn(imag((c-b)*conj(p-c)*(a-p)*conj(b-a)));
+}
+
 // Sind die Dreiecke a1, b1, c1, and a2, b2, c2 ähnlich?
 // Erste Zeile testet Ähnlichkeit mit gleicher Orientierung,
 // zweite Zeile testet Ähnlichkeit mit verschiedener Orientierung
-bool similar (pt a1, pt b1, pt c1, pt a2, pt b2, pt c2) {
-	return ((b2-a2) * (c1-a1) == (b1-a1) * (c2-a2) ||
-	        (b2-a2) * (conj(c1)-conj(a1)) == (conj(b1)-conj(a1))
-			* (c2-a2)
-	);
+bool similar(pt a1, pt b1, pt c1, pt a2, pt b2, pt c2) {
+	return (b2-a2) * (c1-a1) == (b1-a1) * (c2-a2) ||
+	       (b2-a2) * conj(c1-a1) == conj(b1-a1) * (c2-a2);
 }