updated tcr

14 files changed, 483 insertions, 212 deletions

diff --git a/geometry/antipodalPoints.cpp b/geometry/antipodalPoints.cpp
new file mode 100644
index 0000000..c1921cc
--- /dev/null
+++ b/geometry/antipodalPoints.cpp
@@ -0,0 +1,13 @@
+vector<pair<int, int>> antipodalPoints(vector<pt>& h) {
+	vector<pair<int, int>> result;
+	int n = (int)h.size();
+	if (n < 2) return result;
+	for (int i = 0, j = 1; i < j; i++) {
+		while (true) {
+			result.push_back({i, j});
+			if (cross(h[(i + 1) % n] - h[i], 
+								h[(j + 1) % n] - h[j]) <= 0) break;
+			j = (j + 1) % n;
+	}}
+	return result;
+}
diff --git a/geometry/circle.cpp b/geometry/circle.cpp
new file mode 100644
index 0000000..f49ea27
--- /dev/null
+++ b/geometry/circle.cpp
@@ -0,0 +1,33 @@
+// berechnet die Schnittpunkte von zwei kreisen
+// (Kreise dürfen nicht gleich sein!)
+vector<pt> circleIntersection(pt c1, double r1, 
+															pt c2, double r2) {
+	double d = abs(c1 - c2);
+	if (d < abs(r1 - r2) || d > abs(r1 + r2)) return {};
+	double a = (r1 * r1 - r2 * r2 + d * d) / (2 * d);
+	pt p = (c2 - c1) * a / d + c1;
+	if (d == abs(r1 - r2) || d == abs(r1 + r2)) return {p};
+	double h = sqrt(r1 * r1 - a * a);
+	return {p + pt{0, 1} * (c2 - c1) * h / d,
+					p - pt{0, 1} * (c2 - c1) * h / d};
+}
+
+// berechnet die Schnittpunkte zwischen
+// einem Kreis(Kugel) und einer Grade 2d und 3d
+vector<pt> circleRayIntersection(pt center, double r, 
+																 pt orig, pt dir) {
+	vector<pt> result;
+	double a = dot(dir, dir);
+	double b = 2 * dot(dir, orig - center);
+	double c = dot(orig - center, orig - center) - r * r;
+	double discr = b * b - 4 * a * c;
+	if (discr >= 0) {
+		//t in [0, 1] => schnitt mit segment [orig, orig + dir]
+		double t1 = -(b + sqrt(discr)) / (2 * a);
+		double t2 = -(b - sqrt(discr)) / (2 * a);
+		if (t1 >= 0) result.push_back(t1 * dir + orig);
+		if (t2 >= 0 && abs(t1 - t2) > EPS) {
+			result.push_back(t2 * dir + orig);
+	}}
+	return result;
+}
diff --git a/geometry/closestPair.cpp b/geometry/closestPair.cpp
index cac544b..a7662b6 100644
--- a/geometry/closestPair.cpp
+++ b/geometry/closestPair.cpp
@@ -1,35 +1,40 @@
 double squaredDist(pt a, pt b) {
-	return (a.fst-b.fst) * (a.fst-b.fst) + (a.snd-b.snd) * (a.snd-b.snd);
+	return real(conj(a-b) * (a-b));
 }
 
 bool compY(pt a, pt b) {
-	if (a.snd == b.snd) return a.fst < b.fst;
-	return a.snd < b.snd;
+	return (imag(a) == imag(b)) ? real(a) < real(b) 
+															: imag(a) < imag(b);
 }
 
-// points.size() > 1 und alle Punkte müssen verschieden sein!
-double shortestDist(vector<pt> &points) {
+bool compX(pt a, pt b) {
+	return (real(a) == real(b)) ? imag(a) < imag(b)
+															: real(a) < real(b);
+}
+
+double shortestDist(vector<pt>& points) { // points.size() > 1
 	set<pt, bool(*)(pt, pt)> status(compY);
-	sort(points.begin(), points.end());
-	double opt = 1e30, sqrtOpt = 1e15;
+	sort(points.begin(), points.end(), compX);
+	double opt = 1.0/0.0, sqrtOpt = 1.0/0.0;
 	auto left = points.begin(), right = points.begin();
 	status.insert(*right); right++;
-	
+
 	while (right != points.end()) {
-		if (fabs(left->fst - right->fst) >= sqrtOpt) {
-			status.erase(*(left++));
+		if (left != right &&
+				abs(real(*left - *right)) >= sqrtOpt) {
+			status.erase(*left);
+			left++;
 		} else {
-			auto lower = status.lower_bound(pt(-1e20, right->snd - sqrtOpt));
-			auto upper = status.upper_bound(pt(-1e20, right->snd + sqrtOpt));
-			while (lower != upper) {
+			auto lower = status.lower_bound({-1.0/0.0, imag(*right) - sqrtOpt});
+			auto upper = status.upper_bound({-1.0/0.0, imag(*right) + sqrtOpt});
+			for (;lower != upper; lower++) {
 				double cand = squaredDist(*right, *lower);
 				if (cand < opt) {
 					opt = cand;
 					sqrtOpt = sqrt(opt);
-				}
-				++lower;
-			}
-			status.insert(*(right++));
+			}}
+			status.insert(*right);
+			right++;
 	}}
 	return sqrtOpt;
 }
diff --git a/geometry/convexHull.cpp b/geometry/convexHull.cpp
index 9c14c45..e988977 100644
--- a/geometry/convexHull.cpp
+++ b/geometry/convexHull.cpp
@@ -1,24 +1,19 @@
-// Laufzeit: O(n*log(n))
-ll cross(const pt p, const pt a, const pt b) {
-  return (a.x - p.x) * (b.y - p.y) - (a.y - p.y) * (b.x - p.x);
-}
-
-// Punkte auf der konvexen Hülle, gegen den Uhrzeigersinn sortiert.
-// Kollineare Punkte nicht enthalten, entferne dafür "=" im CCW-Test.
-// Achtung: Der erste und letzte Punkt im Ergebnis sind gleich.
-// Achtung: Alle Punkte müssen verschieden sein.
 vector<pt> convexHull(vector<pt> p){
-  int n = p.size(), k = 0;
-  vector<pt> h(2 * n);
-  sort(p.begin(), p.end());
-  for (int i = 0; i < n; i++) { // Untere Hülle.
-    while (k >= 2 && cross(h[k - 2], h[k - 1], p[i]) <= 0.0) k--;
-    h[k++] = p[i];
-  }
-  for (int i = n - 2, t = k + 1; i >= 0; i--) { // Obere Hülle.
-    while (k >= t && cross(h[k - 2], h[k - 1], p[i]) <= 0.0) k--;
-    h[k++] = p[i];
-  }
-  h.resize(k);
-  return h;
+	sort(p.begin(), p.end(), [](const pt& a, const pt& b){
+		return real(a) == real(b) ? imag(a) < imag(b)
+															: real(a) < real(b);
+	});
+	p.erase(unique(p.begin(), p.end()), p.end());
+	int k = 0;
+	vector<pt> h(2 * sz(p));
+	for (int i = 0; i < sz(p); i++) {// Untere Hülle.
+		while (k > 1 && cross(h[k-2], h[k-1], p[i]) <= 0) k--;
+		h[k++] = p[i];
+	}
+	for (int i = sz(p)-2, t = k; i >= 0; i--) {// Obere Hülle.
+		while (k > t && cross(h[k-2], h[k-1], p[i]) <= 0) k--;
+		h[k++] = p[i];
+	}
+	h.resize(k);
+	return h;
 }
diff --git a/geometry/formulars.cpp b/geometry/formulars.cpp
index ed44b8b..43affbb 100644
--- a/geometry/formulars.cpp
+++ b/geometry/formulars.cpp
@@ -1,165 +1,37 @@
-// Komplexe Zahlen als Darstellung für Punkte. Wenn immer möglich
-// complex<int> verwenden. Funktionen wie abs() geben dann int zurück.
-typedef complex<double> pt;
+// Komplexe Zahlen als Darstellung für Punkte. 
+// Wenn immer möglich complex<int> verwenden.
+// Funktionen wie abs() geben dann int zurück.
+using pt = complex<double>;
 
-// Winkel zwischen Punkt und x-Achse in [0, 2 * PI).
-double angle = arg(a);
+// PIL < PI < PIU
+constexpr double PIU = acos(-1.0l);
+constexpr double PIL = PIU-2e-19l;
 
-// Punkt rotiert um Winkel theta.
-pt a_rotated = a * exp(pt(0, theta));
+// Winkel zwischen Punkt und x-Achse in [-PI, PI].
+double angle(pt a) {return arg(a);}
 
-// Mittelpunkt des Dreiecks abc.
-pt centroid = (a + b + c) / 3.0;
+// rotiert Punkt im Uhrzeigersinn um den Ursprung.
+pt rotate(pt a, double theta) {return a * exp(pt(0.0, theta));}
 
 // Skalarprodukt.
-double dot(pt a, pt b) { return real(conj(a) * b); }
+double dot(pt a, pt b) {return real(conj(a) * b);}
 
-// Kreuzprodukt, 0, falls kollinear.
-double cross(pt a, pt b) { return imag(conj(a) * b); }
-
-// Flächeninhalt eines Dreicks bei bekannten Eckpunkten.
-double areaOfTriangle(pt a, pt b, pt c) {
-	return abs(cross(b - a, c - a)) / 2.0;
-}
+// abs()^2.(pre c++20)
+double norm(pt a) {return dot(a, a);}
 
-// Flächeninhalt eines Dreiecks bei bekannten Seitenlängen.
-double areaOfTriangle(double a, double b, double c) {
-	double s = (a + b + c) / 2;
-	return sqrt(s * (s-a) * (s-b) * (s-c));
-}
-
-// Sind die Dreiecke a1, b1, c1, and a2, b2, c2 ähnlich?
-// Erste Zeile testet Ähnlichkeit mit gleicher Orientierung,
-// zweite Zeile testet Ähnlichkeit mit unterschiedlicher 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)
-	);
-}
+// Kreuzprodukt, 0, falls kollinear.
+double cross(pt a, pt b) {return imag(conj(a) * b);}
+double cross(pt p, pt a, pt b) {return cross(a - p, b - p);}
 
-// -1 => gegen den Uhrzeigersinn, 0 => kolliniear, 1 => im Uhrzeigersinn.
-// Einschränken der Rückgabe auf [-1,1] ist sicherer gegen Overflows.
-double orientation(pt a, pt b, pt c) {
+// -1 => gegen den Uhrzeigersinn
+//  0 => kolliniear
+//  1 => im Uhrzeigersinn.
+int orientation(pt a, pt b, pt c) {
 	double orien = cross(b - a, c - a);
-	if (abs(orien) < EPSILON) return 0; // Braucht großes EPSILON: ~1e-6
-	return orien < 0 ? -1 : 1;
-}
-
-// Test auf Streckenschnitt zwischen a-b und c-d.
-bool lineSegmentIntersection(pt a, pt b, pt c, pt d) {
-	if (orientation(a, b, c) == 0 && orientation(a, b, d) == 0) {
-		double dist = abs(a - b);
-		return (abs(a - c) <= dist && abs(b - c) <= dist) ||
-					 (abs(a - d) <= dist && abs(b - d) <= dist);
-	}
-	return orientation(a, b, c) * orientation(a, b, d) <= 0 &&
-				 orientation(c, d, a) * orientation(c, d, b) <= 0;
-}
-
-// Berechnet die Schnittpunkte der Strecken a-b und c-d. Enthält entweder
-// keinen Punkt, den einzigen Schnittpunkt oder die Endpunkte der
-// Schnittstrecke. operator<, min, max müssen noch geschrieben werden!
-vector<pt> lineSegmentIntersection(pt a, pt b, pt c, pt d) {
-	vector<pt> result;
-	if (orientation(a, b, c) == 0 && orientation(a, b, d) == 0 &&
-			orientation(c, d, a) == 0 && orientation(c, d, b) == 0) {
-		pt minAB = min(a, b), maxAB = max(a, b);
-		pt minCD = min(c, d), maxCD = max(c, d);
-		if (minAB < minCD && maxAB < minCD) return result;
-		if (minCD < minAB && maxCD < minAB) return result;
-		pt start = max(minAB, minCD), end = min(maxAB, maxCD);
-		result.push_back(start);
-		if (start != end) result.push_back(end); 
-		return result;
-	}
-	double x1 = real(b) - real(a), y1 = imag(b) - imag(a);
-	double x2 = real(d) - real(c), y2 = imag(d) - imag(c);
-	double u1 = (-y1 * (real(a) - real(c)) + x1 * (imag(a) - imag(c))) /
-			(-x2 * y1 + x1 * y2);
-	double u2 = (x2 * (imag(a) - imag(c)) - y2 * (real(a) - real(c))) /
-			(-x2 * y1 + x1 * y2);
-	if (u1 >= 0 && u1 <= 1 && u2 >= 0 && u2 <= 1) {
-		double x = real(a) + u2 * x1, y = imag(a) + u2 * y1;
-		result.push_back(pt(x, y));
-	}
-	return result;
-}
-
-// Entfernung von Punkt p zur Gearden durch a-b.
-double distToLine(pt a, pt b, pt p) {
-	return abs(cross(p - a, b - a)) / abs(b - a);
-}
-
-// Liegt p auf der Geraden a-b?
-bool pointOnLine(pt a, pt b, pt p) {
-	return orientation(a, b, p) == 0;
-}
-
-// Liegt p auf der Strecke a-b?
-bool pointOnLineSegment(pt a, pt b, pt p) {
-	if (orientation(a, b, p) != 0) return false;
-	return real(p) >= min(real(a), real(b)) &&
-				 real(p) <= max(real(a), real(b)) &&
-				 imag(p) >= min(imag(a), imag(b)) &&
-				 imag(p) <= max(imag(a), imag(b));
-}
-
-// Entfernung von Punkt p zur Strecke a-b.
-double distToSegment(pt a, pt b, pt p) {
-  if (a == b) return abs(p - a);
-  double segLength = abs(a - b);
-	double u = ((real(p) - real(a)) * (real(b) - real(a)) +
-			(imag(p) - imag(a)) * (imag(b) - imag(a))) /
-			(segLength * segLength);
-  pt projection(real(a) + u * (real(b) - real(a)),
-  		imag(a) + u * (imag(b) - imag(a)));
-  double projectionDist = abs(p - projection);
-  if (!pointOnLineSegment(a, b, projection)) projectionDist = 1e30;
-  return min(projectionDist, min(abs(p - a), abs(p - b)));
-}
-
-// Kürzeste Entfernung zwischen den Strecken a-b und c-d.
-double distBetweenSegments(pt a, pt b, pt c, pt d) {
-	if (lineSegmentIntersection(a, b, c, d)) return 0.0;
-	double result = distToSegment(a, b, c);
-	result = min(result, distToSegment(a, b, d));
-	result = min(result, distToSegment(c, d, a));
-	return min(result, distToSegment(c, d, b));
+	return (orien > EPS) - (orien < -EPS);
 }
 
 // Liegt d in der gleichen Ebene wie a, b, und c?
 bool isCoplanar(pt a, pt b, pt c, pt d) {
-	return abs((b - a) * (c - a) * (d - a)) < EPSILON;
-}
-
-// Berechnet den Flächeninhalt eines Polygons (nicht selbstschneidend).
-// Punkte gegen den Uhrzeigersinn: positiv, sonst negativ.
-double areaOfPolygon(vector<pt> &polygon) { // Jeder Eckpunkt nur einmal.
-	double res = 0; int n = polygon.size();
-	for (int i = 0; i < n; i++)
-		res += real(polygon[i]) * imag(polygon[(i + 1) % n]) -
-					 real(polygon[(i + 1) % n]) * imag(polygon[i]);
-	return 0.5 * res;
-}
-
-// Schneiden sich (p1, p2) und (p3, p4) (gegenüberliegende Ecken).
-bool rectIntersection(pt p1, pt p2, pt p3, pt p4) {
-	double minx12=min(real(p1), real(p2)), maxx12=max(real(p1), real(p2));
-	double minx34=min(real(p3), real(p4)), maxx34=max(real(p3), real(p4));
-	double miny12=min(imag(p1), imag(p2)), maxy12=max(imag(p1), imag(p2));
-	double miny34=min(imag(p3), imag(p4)), maxy34=max(imag(p3), imag(p4));
-	return (maxx12 >= minx34) && (maxx34 >= minx12) &&
-				 (maxy12 >= miny34) && (maxy34 >= miny12);
-}
-
-// Testet, ob ein Punkt im Polygon liegt (beliebige Polygone).
-bool pointInPolygon(pt p, vector<pt> &polygon) { // Punkte nur einmal.
-	pt rayEnd = p + pt(1, 1000000);
-	int counter = 0, n = polygon.size();
-	for (int i = 0; i < n; i++) {
-		pt start = polygon[i], end = polygon[(i + 1) % n];
-		if (lineSegmentIntersection(p, rayEnd, start, end)) counter++;
-	}
-	return counter & 1;
+	return abs((b - a) * (c - a) * (d - a)) < EPS;
 }
diff --git a/geometry/formulars3d.cpp b/geometry/formulars3d.cpp
new file mode 100644
index 0000000..f527d57
--- /dev/null
+++ b/geometry/formulars3d.cpp
@@ -0,0 +1,53 @@
+// Skalarprodukt
+double operator|(pt3 a, pt3 b) {
+	return a.x * b.x + a.y*b.y + a.z*b.z;
+}
+double dot(pt3 a, pt3 b) {return a|b;}
+
+// Kreuzprodukt
+pt3 operator*(pt3 a, pt3 b) {return {a.y*b.z - a.z*b.y,
+																		 a.z*b.x - a.x*b.z,
+																		 a.x*b.y - a.y*b.x};}
+pt3 cross(pt3 a, pt3 b) {return a*b;}
+
+// Länge von a
+double abs(pt3 a) {return sqrt(dot(a, a));}
+double abs(pt3 a, pt3 b) {return abs(b - a);}
+
+// Mixedprodukt
+double mixed(pt3 a, pt3 b, pt3 c) {return a*b|c;};
+
+// orientierung von p zu der Ebene durch a, b, c
+// -1 => gegen den Uhrzeigersinn,
+//  0 => kolliniear,
+//  1 => im Uhrzeigersinn.
+int orientation(pt3 a, pt3 b, pt3 c, pt3 p) {
+	double orien = mixed(b - a, c - a, p - a);
+	return (orien > EPS) - (orien < -EPS);
+}
+
+// Entfernung von Punkt p zur Ebene a,b,c.
+double distToPlane(pt3 a, pt3 b, pt3 c, pt3 p) {
+	pt3 n = cross(b-a, c-a);
+	return (abs(dot(n, p)) - dot(n, a)) / abs(n);
+}
+
+// Liegt p in der Ebene a,b,c?
+bool pointOnPlane(pt3 a, pt3 b, pt3 c, pt3 p) {
+	return orientation(a, b, c, p) == 0;
+}
+
+// Schnittpunkt von der Grade a-b und der Ebene c,d,e
+// die Grade darf nicht parallel zu der Ebene sein!
+pt3 linePlaneIntersection(pt3 a, pt3 b, pt3 c, pt3 d, pt3 e) {
+	pt3 n = cross(d-c, e-c);
+	pt3 d = b - a;
+	return a - d * (dot(n, a) - dot(n, c)) / dot(n, d);
+}
+
+// Abstand zwischen der Grade a-b und c-d
+double lineLineDist(pt3 a, pt3 b, pt3 c, pt3 d) {
+	pt3 n = cross(b - a, d - c);
+	if (abs(n) < EPS) return distToLine(a, b, c);
+	return abs(dot(a - c, n)) / abs(n);
+}
\ No newline at end of file
diff --git a/geometry/geometry.tex b/geometry/geometry.tex
index de7b2f6..1201917 100644
--- a/geometry/geometry.tex
+++ b/geometry/geometry.tex
@@ -1,13 +1,47 @@
 \section{Geometrie}
 
-\subsection{Closest Pair}
-\lstinputlisting{geometry/closestPair.cpp}
+\begin{algorithm}{Closest Pair}
+	\begin{methods}
+		\method{shortestDist}{kürzester Abstand zwischen Punkten}{n\*\log(n)}
+	\end{methods}
+	\sourcecode{geometry/closestPair.cpp}
+\end{algorithm}
 
-\subsection{Geraden}
-\lstinputlisting{geometry/lines.cpp}
+\begin{algorithm}{Konvexe Hülle}
+	\begin{methods}
+		\method{convexHull}{berechnet Konvexehülle}{n\*\log(n)}
+	\end{methods}
+	\begin{itemize}
+		\item Konvexehülle gegen den Uhrzeigersinn Sortiert
+		\item nur Eckpunkte enthalten(für alle Punkte = im CCW Test entfernen)
+		\item Erster und Letzter Punkt sind identisch
+	\end{itemize}
+	\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/triangle.cpp}
+\sourcecode{geometry/polygon.cpp}
+\sourcecode{geometry/circle.cpp}
+\sourcecode{geometry/sortAround.cpp}
 
-\subsection{Konvexe Hülle}
-\lstinputlisting{geometry/convexHull.cpp}
+\subsection{Formeln - 3D}
+\sourcecode{geometry/formulars3d.cpp}
 
-\subsection{Formeln - \lstinline{std::complex}}
-\lstinputlisting{geometry/formulars.cpp}
+\subsection{3D-Kugeln}
+\sourcecode{geometry/spheres.cpp}
+
+\optional{
+\subsection{Geraden}
+\sourcecode{geometry/lines.cpp}
+}
diff --git a/geometry/lines.cpp b/geometry/lines.cpp
index 2594ab4..95536a4 100644
--- a/geometry/lines.cpp
+++ b/geometry/lines.cpp
@@ -1,32 +1,33 @@
-// Nicht complex<double> benutzen. Eigene struct schreiben.
 struct line {
 	double a, b, c; // ax + by + c = 0; vertikale Line: b = 0, sonst: b = 1
+	line(pt p, pt q) : a(-imag(q-p)), b(real(q-p)), c(cross({b, -a},p)) {}
 };
 
 line pointsToLine(pt p1, pt p2) {
 	line l;
-	if (fabs(p1.x - p2.x) < EPSILON) {
-		l.a = 1; l.b = 0.0; l.c = -p1.x;
+	if (abs(real(p1 - p2)) < EPS) {
+		l.a = 1; l.b = 0.0; l.c = -real(p1);
 	} else {
-		l.a = -(double)(p1.y - p2.y) / (p1.x - p2.x);
+		l.a = -imag(p1 - p2) / real(p1 - p2);
 		l.b = 1.0;
-		l.c = -(double)(l.a * p1.x) - p1.y;
+		l.c = -(l.a * real(p1)) - imag(p1);
 	}
 	return l;
 }
 
-bool areParallel(line l1, line l2) {
-	return (fabs(l1.a - l2.a) < EPSILON) && (fabs(l1.b - l2.b) < EPSILON);
+bool parallel(line l1, line l2) {
+	return (abs(l1.a - l2.a) < EPS) && (abs(l1.b - l2.b) < EPS);
 }
 
-bool areSame(line l1, line l2) {
-	return areParallel(l1, l2) && (fabs(l1.c - l2.c) < EPSILON);
+bool same(line l1, line l2) {
+	return parallel(l1, l2) && (abs(l1.c - l2.c) < EPS);
 }
 
-bool areIntersect(line l1, line l2, pt &p) {
-	if (areParallel(l1, l2)) return false;
-	p.x = (l2.b * l1.c - l1.b * l2.c) / (l2.a * l1.b - l1.a * l2.b);
-	if (fabs(l1.b) > EPSILON) p.y = -(l1.a * p.x + l1.c);
-	else p.y = -(l2.a * p.x + l2.c);
+bool intersect(line l1, line l2, pt& p) {
+	if (parallel(l1, l2)) return false;
+	double y, x = (l2.b * l1.c - l1.b * l2.c) / (l2.a * l1.b - l1.a * l2.b);
+	if (abs(l1.b) > EPS) y = -(l1.a * x + l1.c);
+	else y = -(l2.a * x + l2.c);
+	p = {x, y};
 	return true;
 }
diff --git a/geometry/linesAndSegments.cpp b/geometry/linesAndSegments.cpp
new file mode 100644
index 0000000..6c53cee
--- /dev/null
+++ b/geometry/linesAndSegments.cpp
@@ -0,0 +1,90 @@
+// Test auf Streckenschnitt zwischen a-b und c-d.
+bool lineSegmentIntersection(pt a, pt b, pt c, pt d) {
+	if (orientation(a, b, c) == 0 && orientation(a, b, d) == 0)
+			return pointOnLineSegment(a,b,c) ||
+						 pointOnLineSegment(a,b,d) ||
+						 pointOnLineSegment(c,d,a) ||
+						 pointOnLineSegment(c,d,b);
+	return orientation(a, b, c) * orientation(a, b, d) <= 0 &&
+				 orientation(c, d, a) * orientation(c, d, b) <= 0;
+}
+
+// Berechnet die Schnittpunkte der Strecken p0-p1 und p2-p3.
+// Enthält entweder keinen Punkt, den einzigen Schnittpunkt 
+// oder die Endpunkte der Schnittstrecke.
+vector<pt> lineSegmentIntersection(pt p0, pt p1, pt p2, pt p3) {
+	double a = cross(p1 - p0, p3 - p2);
+	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 (a > EPS) return {p0 + b/a*(p1 - p0)};
+	vector<pt> result;
+	auto insertUnique = [&](pt p) {
+		for (auto q: result) if (abs(p - q) < EPS) return;
+		result.push_back(p);
+	};
+	if (dot(p2-p0, p3-p0) < EPS) insertUnique(p0);
+	if (dot(p2-p1, p3-p1) < EPS) insertUnique(p1);
+	if (dot(p0-p2, p1-p2) < EPS) insertUnique(p2);
+	if (dot(p0-p3, p1-p3) < EPS) insertUnique(p3);
+	return result;
+}
+
+// Entfernung von Punkt p zur Gearden durch a-b. 2d und 3d
+double distToLine(pt a, pt b, pt p) {
+	return abs(cross(p - a, b - a)) / abs(b - a);
+}
+
+// Liegt p auf der Geraden a-b? 2d und 3d
+bool pointOnLine(pt a, pt b, pt p) {
+	return orientation(a, b, p) == 0;
+}
+
+// Test auf Linienschnitt zwischen a-b und c-d.
+bool lineIntersection(pt a, pt b, pt c, pt d) {
+	return abs(cross(a - b, c - d)) < EPS;
+}
+
+// Berechnet den Schnittpunkt der Graden p0-p1 und p2-p3.
+// die Graden dürfen nicht parallel sein! 
+pt lineIntersection(pt p0, pt p1, pt p2, pt p3) {
+	double a = cross(p1 - p0, p3 - p2);
+	double b = cross(p2 - p0, p3 - p2);
+	return {p0 + b/a*(p1 - p0)};
+}
+
+// Liegt p auf der Strecke a-b?
+bool pointOnLineSegment(pt a, pt b, pt p) {
+	if (orientation(a, b, p) != 0) return false;
+	ld dist = norm(a - b);
+	return norm(a - p) <= dist && norm(b - p) <= dist;
+}
+
+// Entfernung von Punkt p zur Strecke a-b.
+double distToSegment(pt a, pt b, pt p) {
+	if (a == b) return abs(p - a);
+	pt dir = b - a;
+	if (dot(dir, a) <= dot(dir, p) && dot(dir, p) <= dot(dir, b)) {
+		return distToLine(a, b, p);
+	} else {
+		return min(abs(p - a), abs(p - b));
+}}
+
+// Kürzeste Entfernung zwischen den Strecken a-b und c-d.
+double distBetweenSegments(pt a, pt b, pt c, pt d) {
+	if (lineSegmentIntersection(a, b, c, d)) return 0.0;
+	double result = distToSegment(a, b, c);
+	result = min(result, distToSegment(a, b, d));
+	result = min(result, distToSegment(c, d, a));
+	return min(result, distToSegment(c, d, b));
+}
+
+// sortiert alle Punkte pts auf einer Linie 
+// entsprechend der richtung dir 2d und 3d
+void sortLine(pt dir, vector<pt>& pts) {
+	sort(pts.begin(), pts.end(), [&](pt a, pt b){
+		return dot(dir, a) < dot(dir, b);
+	});
+}
\ No newline at end of file
diff --git a/geometry/polygon.cpp b/geometry/polygon.cpp
new file mode 100644
index 0000000..6a8080d
--- /dev/null
+++ b/geometry/polygon.cpp
@@ -0,0 +1,33 @@
+// Berechnet den Flächeninhalt eines Polygons
+// (nicht selbstschneidend).
+// Punkte gegen den Uhrzeigersinn: positiv, sonst negativ.
+double area(const vector<pt>& polygon) {//points unique
+	double res = 0; int n = sz(polygon);
+	for (int i = 0; i < n; i++)
+		res += cross(polygon[i], polygon[(i + 1) % n]);
+	return 0.5 * res;
+}
+
+// Testet, ob ein Punkt im Polygon liegt (beliebige Polygone).
+bool inside(pt p, const vector<pt>& polygon) {//points unique
+	pt rayEnd = p + pt(1, 1000000);
+	int counter = 0, n = sz(polygon);
+	for (int i = 0; i < n; i++) {
+		pt start = polygon[i], end = polygon[(i + 1) % n];
+		if (lineSegmentIntersection(p, rayEnd, start, end)) {
+			counter++;
+	}}
+	return counter & 1;
+}
+
+//convex hull without duplicates, h[0] == h.back()
+bool inside(pt p, const vector<pt>& hull) {
+	if (orientation(hull[0], hull.back(), p) > 0) return false;
+	ll l = 0, r = sz(hull) - 1;
+	while (l + 1 < r) {
+		ll m = (l + r) / 2;
+		if (orientation(hull[0], hull[m], p) >= 0) l = m;
+		else r = m;
+	}
+	return orientation(hull[l], hull[r], p) >= 0;
+}
diff --git a/geometry/segmentIntersection.cpp b/geometry/segmentIntersection.cpp
new file mode 100644
index 0000000..c2961da
--- /dev/null
+++ b/geometry/segmentIntersection.cpp
@@ -0,0 +1,63 @@
+struct seg {
+	pt a, b;
+	int id;
+	bool operator<(const seg& o) const {
+		if (real(a) < real(o.a)) {
+			int s = orientation(a, b, o.a);
+			return (s > 0 || (s == 0 && imag(a) < imag(o.a)));//
+		} else if (real(a) > real(o.a)) {
+			int s = orientation(o.a, o.b, a);
+			return (s < 0 || (s == 0 && imag(a) < imag(o.a)));
+		}
+		return imag(a) < imag(o.a);
+	}
+};
+ 
+struct event {
+	pt p;
+	int id, type;
+	bool operator<(const event& o) const {
+		if (real(p) != real(o.p)) return real(p) < real(o.p);
+		if (type != o.type) return type > o.type;
+		return imag(p) < imag(o.p);
+	}
+};
+ 
+bool lessPT(const pt& a, const pt& b) {
+	return real(a) != real(b) ? real(a) < real(b)
+														: imag(a) < imag(b);
+}
+ 
+bool intersect(const seg& a, const seg& b) {
+	return lineSegmentIntersection(a.a, a.b, b.a, b.b);
+}
+ 
+pair<int, int> intersect(vector<seg>& segs) {
+	vector<event> events;
+	for (seg& s : segs) {
+		if (lessPT(s.b, s.a)) swap(s.b, s.a);
+		events.push_back({s.a, s.id, 1});
+		events.push_back({s.b, s.id, -1});
+	}
+	sort(all(events));
+ 
+	set<seg> q;
+	vector<set<seg>::iterator> where(sz(segs));
+	for (auto e : events) {
+		int id = e.id;
+		if (e.type > 0) {
+			auto it = q.lower_bound(segs[id]);
+			if (it != q.end() && intersect(*it, segs[id]))
+				return {it->id, segs[id].id};
+			if (it != q.begin() && intersect(*prev(it), segs[id]))
+				return {prev(it)->id, segs[id].id};
+			where[id] = q.insert(it, segs[id]);
+		} else {
+			auto it = where[id];
+			if (it != q.begin() && next(it) != q.end() && intersect(*next(it), *prev(it)))
+				return {next(it)->id, prev(it)->id};
+			q.erase(it);
+		}
+	}
+	return {-1, -1};
+}
\ No newline at end of file
diff --git a/geometry/sortAround.cpp b/geometry/sortAround.cpp
new file mode 100644
index 0000000..a95d224
--- /dev/null
+++ b/geometry/sortAround.cpp
@@ -0,0 +1,10 @@
+bool left(pt p) {return real(p) < 0 || 

+											 (real(p) == 0 && imag(p) < 0);}

+

+void sortAround(pt p, vector<pt>& ps) {

+	sort(all(ps), [&](const pt& a, const pt& b){

+		if (left(a - p) != left(b - p)) 

+			return left(a - p) > left(b - p);

+		return cross(p, a, b) > 0;

+	});

+}
\ No newline at end of file
diff --git a/geometry/spheres.cpp b/geometry/spheres.cpp
new file mode 100644
index 0000000..0657ab8
--- /dev/null
+++ b/geometry/spheres.cpp
@@ -0,0 +1,29 @@
+// Great Cirlce Distance mit Längen- und Breitengrad.
+double gcDist(double pLat, double pLon, 
+							double qLat, double qLon, double radius) {
+	pLat *= PI / 180; pLon *= PI / 180;
+	qLat *= PI / 180; qLon *= PI / 180;
+	return radius * acos(cos(pLat) * cos(pLon) *
+											 cos(qLat) * cos(qLon) +
+											 cos(pLat) * sin(pLon) *
+											 cos(qLat) * sin(qLon) +
+											 sin(pLat) * sin(qLat));
+}
+
+// Great Cirlce Distance mit kartesischen Koordinaten.
+double gcDist(point p, point q) {
+	return acos(p.x * q.x + p.y * q.y + p.z * q.z);
+}
+
+// 3D Punkt in kartesischen Koordinaten.
+struct point{
+	double x, y, z;
+	point() {}
+	point(double x, double y, double z) : x(x), y(y), z(z) {}
+	point(double lat, double lon) {
+		lat *= PI / 180.0; lon *= PI / 180.0;
+		x = cos(lat) * sin(lon);
+		y = cos(lat) * cos(lon);
+		z = sin(lat);
+	}
+};
diff --git a/geometry/triangle.cpp b/geometry/triangle.cpp
new file mode 100644
index 0000000..3a39302
--- /dev/null
+++ b/geometry/triangle.cpp
@@ -0,0 +1,40 @@
+// Mittelpunkt des Dreiecks abc.
+pt centroid(pt a, pt b, pt c) {return (a + b + c) / 3.0;}
+
+// Flächeninhalt eines Dreicks bei bekannten Eckpunkten.
+double area(pt a, pt b, pt c) {
+	return abs(cross(b - a, c - a)) / 2.0;
+}
+
+// Flächeninhalt eines Dreiecks bei bekannten Seitenlängen.
+double area(double a, double b, double c) {
+	double s = (a + b + c) / 2.0;
+	return sqrt(s * (s-a) * (s-b) * (s-c));
+}
+
+// Zentrum des Kreises durch alle Eckpunkte
+pt outCenter(pt a, pt b, pt c) {
+	double d = 2.0 * (real(a) * imag(b-c) +
+										real(b) * imag(c-a) +
+										real(c) * imag(a-b));
+	return (a*conj(a)*conj(b-c) +
+					b*conj(b)*conj(c-a) +
+					c*conj(c)*conj(a-b)) / d;
+}
+
+// Zentrum des größten Kreises im Dreiecke
+pt inCenter(pt a, pt b, pt c) {
+	double x = abs(a-b), y = abs(b-c), z = abs(a-c);
+	return (y*a + z*b + x*c) / (a+b+c);
+}
+
+
+// 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)
+	);
+}