summaryrefslogtreecommitdiff
path: root/content/geometry
diff options
context:
space:
mode:
authormzuenni <mzuenni@users.noreply.github.com>2024-07-28 22:54:40 +0200
committerGitHub <noreply@github.com>2024-07-28 22:54:40 +0200
commit8d11c6c8213f46f0fa19826917c255edd5d43cb1 (patch)
tree96d75baff33d5a04b5a60f1a41f514a26c716874 /content/geometry
parent8c33b4e0d3030cfed17fc64b4fe41133339f6d87 (diff)
Test (#4)
* update * moved content in subdir * rename file * add test setup * add test setup * add github action * automaticly test all cpp files * timeout after 10s * setulimit and dont zero memory * test build pdf * install latexmk * update * update * ngerman * fonts * removed old code * add first test * added tests * test in sorted order * more tests * simplified test * more tests * fix suffix tree * fixes and improvements * done ust lst directly * fix swap * add links to pdf * fix constants * add primorial * add comment * various improvements * more tests * added missing stuf * more tests * fix tests * more tests * more tests * more tests * fix recursion? * test trie * more tests * only use python temporarily for listings * only use python temporarily for listings * more tests * fix longestCommonSubstring * more tests * more tests * made code more similiar * fix? * more tests * more tests * more tests * add ahoCorasick test + limit 4GB stack size * more tests * fix test * add additional test * more tests * more tests * fix? * better fix * fix virtual tree * more tests * more tests * recursive closest pair * more tests * decrease limit * new tests * more tests * fix name * more tests * add test * new test * more tests * more tests * more tests * more tests * new test and content * new code * new code * larger tests * fix and test * new test * new test * update pdf * remove comments * new test * more tests * more testcases * more tests * increased limit * more tests * more tests * more tests * new tests * more tests * shortened code * new test * add basic tests for bigint * more tests * removed old files * new test * ignore some files * more auto more ccw * fix test * more tests * fix * new tests * more tests * more tests * stronger test * actually verify delaunay... * more tests * fix header * more tests * run tests parallel? * test parralel? * add --missing * separate workflows * test * is the pdf checked? * separate workflows * fix workflow * more workflows --------- Co-authored-by: Yidi <noob999noob999@gmail.com>
Diffstat (limited to 'content/geometry')
-rw-r--r--content/geometry/antipodalPoints.cpp12
-rw-r--r--content/geometry/circle.cpp33
-rw-r--r--content/geometry/closestPair.cpp27
-rw-r--r--content/geometry/convexHull.cpp18
-rw-r--r--content/geometry/delaunay.cpp124
-rw-r--r--content/geometry/formulas.cpp42
-rw-r--r--content/geometry/formulas3d.cpp53
-rw-r--r--content/geometry/geometry.tex62
-rw-r--r--content/geometry/hpi.cpp68
-rw-r--r--content/geometry/lines.cpp33
-rw-r--r--content/geometry/linesAndSegments.cpp89
-rw-r--r--content/geometry/polygon.cpp150
-rw-r--r--content/geometry/segmentIntersection.cpp63
-rw-r--r--content/geometry/sortAround.cpp11
-rw-r--r--content/geometry/spheres.cpp29
-rw-r--r--content/geometry/triangle.cpp43
-rw-r--r--content/geometry/triangle.tex41
17 files changed, 898 insertions, 0 deletions
diff --git a/content/geometry/antipodalPoints.cpp b/content/geometry/antipodalPoints.cpp
new file mode 100644
index 0000000..110cc74
--- /dev/null
+++ b/content/geometry/antipodalPoints.cpp
@@ -0,0 +1,12 @@
+vector<pair<int, int>> antipodalPoints(vector<pt>& h) {
+ if (sz(h) < 2) return {};
+ vector<pair<int, int>> result;
+ for (int i = 0, j = 1; i < j; i++) {
+ while (true) {
+ result.push_back({i, j});
+ if (cross(h[(i + 1) % sz(h)] - h[i],
+ h[(j + 1) % sz(h)] - h[j]) <= 0) break;
+ j = (j + 1) % sz(h);
+ }}
+ return result;
+}
diff --git a/content/geometry/circle.cpp b/content/geometry/circle.cpp
new file mode 100644
index 0000000..6789c52
--- /dev/null
+++ b/content/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 einem Strahl (2D und 3D)
+vector<pt> circleRayIntersection(pt center, double r,
+ pt orig, pt dir) {
+ vector<pt> result;
+ double a = norm(dir);
+ double b = 2 * dot(dir, orig - center);
+ double c = norm(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/content/geometry/closestPair.cpp b/content/geometry/closestPair.cpp
new file mode 100644
index 0000000..9b115f3
--- /dev/null
+++ b/content/geometry/closestPair.cpp
@@ -0,0 +1,27 @@
+ll rec(vector<pt>::iterator a, int l, int r) {
+ if (r - l < 2) return INF;
+ int m = (l + r) / 2;
+ ll midx = a[m].real();
+ ll ans = min(rec(a, l, m), rec(a, m, r));
+
+ inplace_merge(a+l, a+m, a+r, [](const pt& x, const pt& y) {
+ return x.imag() < y.imag();
+ });
+
+ pt tmp[8];
+ fill(all(tmp), a[l]);
+ for (int i = l + 1, next = 0; i < r; i++) {
+ if (ll x = a[i].real() - midx; x * x < ans) {
+ for (pt& p : tmp) ans = min(ans, norm(p - a[i]));
+ tmp[next++ & 7] = a[i];
+ }
+ }
+ return ans;
+}
+
+ll shortestDist(vector<pt> a) { // sz(pts) > 1
+ sort(all(a), [](const pt& x, const pt& y) {
+ return x.real() < y.real();
+ });
+ return rec(a.begin(), 0, sz(a));
+}
diff --git a/content/geometry/convexHull.cpp b/content/geometry/convexHull.cpp
new file mode 100644
index 0000000..6d89e05
--- /dev/null
+++ b/content/geometry/convexHull.cpp
@@ -0,0 +1,18 @@
+vector<pt> convexHull(vector<pt> pts){
+ sort(all(pts), [](const pt& a, const pt& b){
+ return real(a) == real(b) ? imag(a) < imag(b)
+ : real(a) < real(b);
+ });
+ pts.erase(unique(all(pts)), pts.end());
+ int k = 0;
+ vector<pt> h(2 * sz(pts));
+ auto half = [&](auto begin, auto end, int t) {
+ for (auto it = begin; it != end; it++) {
+ while (k > t && cross(h[k-2], h[k-1], *it) <= 0) k--;
+ h[k++] = *it;
+ }};
+ half(all(pts), 1);// Untere Hülle.
+ half(next(pts.rbegin()), pts.rend(), k);// Obere Hülle.
+ h.resize(k);
+ return h;
+}
diff --git a/content/geometry/delaunay.cpp b/content/geometry/delaunay.cpp
new file mode 100644
index 0000000..c813892
--- /dev/null
+++ b/content/geometry/delaunay.cpp
@@ -0,0 +1,124 @@
+using lll = __int128;
+using pt = complex<lll>;
+
+constexpr pt INF_PT = pt(2e18, 2e18);
+
+bool circ(pt p, pt a, pt b, pt c) {// p in circle(A,B,C), ABC must be ccw
+ 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 _ : {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);
+ auto 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/content/geometry/formulas.cpp b/content/geometry/formulas.cpp
new file mode 100644
index 0000000..5d4e10d
--- /dev/null
+++ b/content/geometry/formulas.cpp
@@ -0,0 +1,42 @@
+// Komplexe Zahlen als Punkte. Wenn immer möglich complex<ll>
+// verwenden. Funktionen wie abs() geben dann aber ll zurück.
+using pt = complex<double>;
+
+constexpr double PIU = acos(-1.0l); // PIL < PI < PIU
+constexpr double PIL = PIU-2e-19l;
+
+// Winkel zwischen Punkt und x-Achse in [-PI, PI].
+double angle(pt a) {return arg(a);}
+
+// rotiert Punkt im Uhrzeigersinn um den Ursprung.
+pt rotate(pt a, double theta) {return a * polar(1.0, theta);}
+
+// Skalarprodukt.
+auto dot(pt a, pt b) {return real(conj(a) * b);}
+
+// abs()^2.(pre c++20)
+auto norm(pt a) {return dot(a, a);}
+
+// Kreuzprodukt, 0, falls kollinear.
+auto cross(pt a, pt b) {return imag(conj(a) * b);}
+auto cross(pt p, pt a, pt b) {return cross(a - p, b - p);}
+
+// 1 => c links von a->b
+// 0 => a, b und c kolliniear
+// -1 => c rechts von a->b
+int ccw(pt a, pt b, pt c) {
+ auto orien = cross(b - a, c - a);
+ 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)) < EPS;
+}
+
+// charakterisiert winkel zwischen Vektoren u und v
+pt uniqueAngle(pt u, pt v) {
+ pt tmp = v * conj(u);
+ ll g = abs(gcd(real(tmp), imag(tmp)));
+ return tmp / g;
+}
diff --git a/content/geometry/formulas3d.cpp b/content/geometry/formulas3d.cpp
new file mode 100644
index 0000000..dee3ce8
--- /dev/null
+++ b/content/geometry/formulas3d.cpp
@@ -0,0 +1,53 @@
+// Skalarprodukt
+auto operator|(pt3 a, pt3 b) {
+ return a.x * b.x + a.y*b.y + a.z*b.z;
+}
+auto 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
+auto 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 ccw(pt3 a, pt3 b, pt3 c, pt3 p) {
+ auto 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 ccw(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);
+}
diff --git a/content/geometry/geometry.tex b/content/geometry/geometry.tex
new file mode 100644
index 0000000..92285c4
--- /dev/null
+++ b/content/geometry/geometry.tex
@@ -0,0 +1,62 @@
+\section{Geometrie}
+
+\begin{algorithm}{Closest Pair}
+ \begin{methods}
+ \method{shortestDist}{kürzester Abstand zwischen Punkten}{n\*\log(n)}
+ \end{methods}
+ \sourcecode{geometry/closestPair.cpp}
+\end{algorithm}
+
+\begin{algorithm}{Konvexehülle}
+ \begin{methods}
+ \method{convexHull}{berechnet konvexe Hülle}{n\*\log(n)}
+ \end{methods}
+ \begin{itemize}
+ \item konvexe Hü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/formulas.cpp}
+\sourcecode{geometry/linesAndSegments.cpp}
+\sourcecode{geometry/sortAround.cpp}
+\input{geometry/triangle}
+\sourcecode{geometry/triangle.cpp}
+\sourcecode{geometry/polygon.cpp}
+\sourcecode{geometry/circle.cpp}
+
+\subsection{Formeln -- 3D}
+\sourcecode{geometry/formulas3d.cpp}
+
+\optional{
+ \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/content/geometry/hpi.cpp b/content/geometry/hpi.cpp
new file mode 100644
index 0000000..3509e0e
--- /dev/null
+++ b/content/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/content/geometry/lines.cpp b/content/geometry/lines.cpp
new file mode 100644
index 0000000..95536a4
--- /dev/null
+++ b/content/geometry/lines.cpp
@@ -0,0 +1,33 @@
+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 (abs(real(p1 - p2)) < EPS) {
+ l.a = 1; l.b = 0.0; l.c = -real(p1);
+ } else {
+ l.a = -imag(p1 - p2) / real(p1 - p2);
+ l.b = 1.0;
+ l.c = -(l.a * real(p1)) - imag(p1);
+ }
+ return l;
+}
+
+bool parallel(line l1, line l2) {
+ return (abs(l1.a - l2.a) < EPS) && (abs(l1.b - l2.b) < EPS);
+}
+
+bool same(line l1, line l2) {
+ return parallel(l1, l2) && (abs(l1.c - l2.c) < EPS);
+}
+
+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/content/geometry/linesAndSegments.cpp b/content/geometry/linesAndSegments.cpp
new file mode 100644
index 0000000..1e21cba
--- /dev/null
+++ b/content/geometry/linesAndSegments.cpp
@@ -0,0 +1,89 @@
+// Liegt p auf der Geraden a-b? 2d und 3d
+bool pointOnLine(pt a, pt b, pt p) {
+ return ccw(a, b, p) == 0;
+}
+
+// Test auf Linienschnitt zwischen a-b und c-d. (nicht identisch)
+bool lineIntersection(pt a, pt b, pt c, pt d) {
+ return abs(cross(a - b, c - d)) < EPS;
+}
+
+// Berechnet den Schnittpunkt der Graden a-b und c-d.
+// die Graden dürfen nicht parallel sein!
+pt lineIntersection2(pt a, pt b, pt c, pt d) {
+ double x = cross(b - a, d - c);
+ double y = cross(c - a, d - c);
+ return a + y/x*(b - a);
+}
+
+// Entfernung von Punkt p zur Geraden durch a-b. 2d und 3d
+double distToLine(pt a, pt b, pt p) {
+ return abs(cross(p - a, b - a)) / abs(b - a);
+}
+
+// Projiziert p auf die Gerade a-b
+pt projectToLine(pt a, pt b, pt p) {
+ return a + (b - a) * dot(p - a, b - a) / norm(b - a);
+}
+
+// sortiert alle Punkte pts auf einer Linie entsprechend dir
+void sortLine(pt dir, vector<pt>& pts) { // (2d und 3d)
+ sort(all(pts), [&](pt a, pt b){
+ return dot(dir, a) < dot(dir, b);
+ });
+}
+
+// Liegt p auf der Strecke a-b? (nutze < für inberhalb)
+bool pointOnSegment(pt a, pt b, pt p) {
+ if (ccw(a, b, p) != 0) return false;
+ auto 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);
+ if (dot(p - a, b - a) <= 0) return abs(p - a);
+ if (dot(p - b, b - a) >= 0) return abs(p - b);
+ return distToLine(a, b, p);
+}
+
+// Test auf Streckenschnitt zwischen a-b und c-d.
+bool segmentIntersection(pt a, pt b, pt c, pt d) {
+ if (ccw(a, b, c) == 0 && ccw(a, b, d) == 0)
+ return pointOnSegment(a,b,c) ||
+ pointOnSegment(a,b,d) ||
+ pointOnSegment(c,d,a) ||
+ pointOnSegment(c,d,b);
+ return ccw(a, b, c) * ccw(a, b, d) <= 0 &&
+ ccw(c, d, a) * ccw(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.
+vector<pt> segmentIntersection2(pt a, pt b, pt c, pt d) {
+ double x = cross(b - a, d - c);
+ double y = cross(c - a, d - c);
+ double z = cross(b - a, a - c);
+ if (x < 0) {x = -x; y = -y; z = -z;}
+ if (y < -EPS || y-x > EPS || z < -EPS || z-x > EPS) return {};
+ if (x > EPS) return {a + y/x*(b - a)};
+ vector<pt> result;
+ auto insertUnique = [&](pt p) {
+ for (auto q : result) if (abs(p - q) < EPS) return;
+ result.push_back(p);
+ };
+ if (dot(c-a, d-a) < EPS) insertUnique(a);
+ if (dot(c-b, d-b) < EPS) insertUnique(b);
+ if (dot(a-c, b-c) < EPS) insertUnique(c);
+ if (dot(a-d, b-d) < EPS) insertUnique(d);
+ return result;
+}
+
+// Kürzeste Entfernung zwischen den Strecken a-b und c-d.
+double distBetweenSegments(pt a, pt b, pt c, pt d) {
+ if (segmentIntersection(a, b, c, d)) return 0.0;
+ return min({distToSegment(a, b, c), distToSegment(a, b, d),
+ distToSegment(c, d, a), distToSegment(c, d, b)});
+}
diff --git a/content/geometry/polygon.cpp b/content/geometry/polygon.cpp
new file mode 100644
index 0000000..3178290
--- /dev/null
+++ b/content/geometry/polygon.cpp
@@ -0,0 +1,150 @@
+// Flächeninhalt eines Polygons (nicht selbstschneidend).
+// Punkte gegen den Uhrzeigersinn: positiv, sonst negativ.
+double area(const vector<pt>& poly) { //poly[0] == poly.back()
+ ll res = 0;
+ for (int i = 0; i + 1 < sz(poly); i++)
+ res += cross(poly[i], poly[i + 1]);
+ return 0.5 * res;
+}
+
+// Anzahl ccw 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)
+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) &&
+ cross(p, a, b) < 0) {
+ res += ccw(p, poly[i], poly[i + 1]);
+ }}
+ return res;
+}
+
+// Testet, ob ein Punkt im Polygon liegt (beliebige Polygone).
+// Ändere Zeile 32 falls rand zählt, poly[0] == poly.back()
+bool inside(pt p, const vector<pt>& poly) {
+ bool in = false;
+ for (int i = 0; i + 1 < sz(poly); i++) {
+ pt a = poly[i], b = poly[i + 1];
+ if (pointOnLineSegment(a, b, p)) return false;
+ if (real(a) > real(b)) swap(a,b);
+ if (real(a) <= real(p) && real(p) < real(b) &&
+ cross(p, a, b) < 0) {
+ in ^= 1;
+ }}
+ return in;
+}
+
+// convex hull without duplicates, h[0] != h.back()
+// apply comments if border counts as inside
+bool insideConvex(pt p, const vector<pt>& hull) {
+ int l = 0, r = sz(hull) - 1;
+ if (cross(hull[0], hull[r], p) >= 0) return false; // > 0
+ while (l + 1 < r) {
+ int m = (l + r) / 2;
+ if (cross(hull[0], hull[m], p) > 0) l = m; // >= 0
+ else r = m;
+ }
+ return cross(hull[l], hull[r], p) > 0; // >= 0
+}
+
+void rotateMin(vector<pt>& hull) {
+ auto mi = min_element(all(hull), [](const pt& a, const pt& b){
+ return real(a) == real(b) ? imag(a) < imag(b)
+ : real(a) < real(b);
+ });
+ rotate(hull.begin(), mi, hull.end());
+}
+
+// convex hulls without duplicates, h[0] != h.back()
+vector<pt> minkowski(vector<pt> ps, vector<pt> qs) {
+ rotateMin(ps);
+ rotateMin(qs);
+ ps.push_back(ps[0]);
+ qs.push_back(qs[0]);
+ ps.push_back(ps[1]);
+ qs.push_back(qs[1]);
+ vector<pt> res;
+ for (ll i = 0, j = 0; i + 2 < sz(ps) || j + 2 < sz(qs);) {
+ res.push_back(ps[i] + qs[j]);
+ auto c = cross(ps[i + 1] - ps[i], qs[j + 1] - qs[j]);
+ if(c >= 0) i++;
+ if(c <= 0) j++;
+ }
+ return res;
+}
+
+// convex hulls without duplicates, h[0] != h.back()
+double dist(const vector<pt>& ps, vector<pt> qs) {
+ for (pt& q : qs) q *= -1;
+ auto p = minkowski(ps, qs);
+ p.push_back(p[0]);
+ double res = INF;
+ bool intersect = true;
+ for (ll i = 0; i + 1 < sz(p); i++) {
+ intersect &= cross(p[i], p[i+1]) >= 0;
+ res = min(res, distToSegment(p[i], p[i+1], 0));
+ }
+ return intersect ? 0 : res;
+}
+
+bool left(pt of, pt p) {return cross(p, of) < 0 ||
+ (cross(p, of) == 0 && dot(p, of) > 0);}
+
+// convex hulls without duplicates, hull[0] == hull.back() and
+// hull[0] must be a convex point (with angle < pi)
+// returns index of corner where dot(dir, corner) is maximized
+int extremal(const vector<pt>& hull, pt dir) {
+ dir *= pt(0, 1);
+ int l = 0, r = sz(hull) - 1;
+ while (l + 1 < r) {
+ int m = (l + r) / 2;
+ pt dm = hull[m+1]-hull[m];
+ pt dl = hull[l+1]-hull[l];
+ if (left(dl, dir) != left(dl, dm)) {
+ if (left(dl, dm)) l = m;
+ else r = m;
+ } else {
+ if (cross(dir, dm) < 0) l = m;
+ else r = m;
+ }}
+ return r % (sz(hull) - 1);
+}
+
+// convex hulls without duplicates, hull[0] == hull.back() and
+// hull[0] must be a convex point (with angle < pi)
+// {} if no intersection
+// {x} if corner is only intersection
+// {i, j} segments (i,i+1) and (j,j+1) intersected (if only the
+// border is intersected corners i and j are the start and end)
+vector<int> intersectLine(const vector<pt>& hull, pt a, pt b) {
+ int endA = extremal(hull, (a-b) * pt(0, 1));
+ int endB = extremal(hull, (b-a) * pt(0, 1));
+ // cross == 0 => line only intersects border
+ if (cross(hull[endA], a, b) > 0 ||
+ cross(hull[endB], a, b) < 0) return {};
+
+ int n = sz(hull) - 1;
+ vector<int> res;
+ for (auto _ : {0, 1}) {
+ int l = endA, r = endB;
+ if (r < l) r += n;
+ while (l + 1 < r) {
+ int m = (l + r) / 2;
+ if (cross(hull[m % n], a, b) <= 0 &&
+ cross(hull[m % n], a, b) != cross(hull[endB], a, b))
+ l = m;
+ else r = m;
+ }
+ if (cross(hull[r % n], a, b) == 0) l++;
+ res.push_back(l % n);
+ swap(endA, endB);
+ swap(a, b);
+ }
+ if (res[0] == res[1]) res.pop_back();
+ return res;
+}
+
diff --git a/content/geometry/segmentIntersection.cpp b/content/geometry/segmentIntersection.cpp
new file mode 100644
index 0000000..4262ddc
--- /dev/null
+++ b/content/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 = ccw(a, b, o.a);
+ return (s > 0 || (s == 0 && imag(a) < imag(o.a)));
+ } else if (real(a) > real(o.a)) {
+ int s = ccw(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};
+}
diff --git a/content/geometry/sortAround.cpp b/content/geometry/sortAround.cpp
new file mode 100644
index 0000000..98d17a8
--- /dev/null
+++ b/content/geometry/sortAround.cpp
@@ -0,0 +1,11 @@
+bool left(pt p) {return real(p) < 0 ||
+ (real(p) == 0 && imag(p) < 0);}
+
+// counter clockwise, starting with "11:59"
+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;
+ });
+}
diff --git a/content/geometry/spheres.cpp b/content/geometry/spheres.cpp
new file mode 100644
index 0000000..ec22262
--- /dev/null
+++ b/content/geometry/spheres.cpp
@@ -0,0 +1,29 @@
+// Great Circle 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 Circle 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/content/geometry/triangle.cpp b/content/geometry/triangle.cpp
new file mode 100644
index 0000000..534bb10
--- /dev/null
+++ b/content/geometry/triangle.cpp
@@ -0,0 +1,43 @@
+// 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(a, b, c)) / 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; //unpräzise
+ return sqrt(s * (s-a) * (s-b) * (s-c));
+}
+
+// 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) / (x+y+z);
+}
+
+// Zentrum des Kreises durch alle Eckpunkte
+// a, b und c nicht kollinear
+pt circumCenter(pt a, pt b, pt c) {
+ b -= a, c -= a;
+ pt d = b * norm(c) - c * norm(b);
+ d = {-d.imag(), d.real()};
+ return a + d / cross(b, c) / 2.0;
+}
+
+// -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 ccw(a,b,c) * 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-a1) == conj(b1-a1) * (c2-a2);
+}
diff --git a/content/geometry/triangle.tex b/content/geometry/triangle.tex
new file mode 100644
index 0000000..3decd54
--- /dev/null
+++ b/content/geometry/triangle.tex
@@ -0,0 +1,41 @@
+
+\begin{minipage}[T]{0.27\linewidth}
+ Generell:
+ \begin{itemize}
+ \item $\cos(\gamma)=\frac{a^2+b^2-c^2}{2ab}$
+ \item $b=\frac{a}{\sin(\alpha)}\sin(\beta)$
+ %\item $b=\frac{a}{\sin(\pi-\beta-\gamma)}\sin(\beta)$
+ %\item $\sin(\beta)=\frac{b\sin(\alpha)}{a}$ %asin is not uniquely invertible
+ \item $\Delta=\frac{bc}{2}\sin(\alpha)$
+ \end{itemize}
+\end{minipage}
+\hfill
+\begin{minipage}[B]{0.5\linewidth}
+ \centering
+ \begin{tikzpicture}[line cap=round,minimum size=0,x=.7cm,y=0.7cm]
+ \node[circle,inner sep=0] (AA) at (0,0) {$A$};
+ \node[circle,inner sep=0] (BB) at (3,-1) {$B$};
+ \node[circle,inner sep=0] (CC) at (3.666667,1) {$C$};
+
+ \coordinate (A) at (AA.0);
+ \coordinate (B) at (BB.100);
+ \coordinate (C) at (CC.210);
+
+ \pic[draw,angle radius=15,pic text=$\gamma$]{angle = A--C--B};
+ \pic[draw,angle radius=15,pic text=$\beta$]{angle = C--B--A};
+ \pic[draw,angle radius=20,pic text=$\alpha$]{angle = B--A--C};
+
+ \draw (A) to[edge label={$b$},inner sep=1] (C);
+ \draw (A) to[edge label'={$c$},inner sep=1.3] (B);
+ \draw (B) to[edge label'={$a$},inner sep=0.6] (C);
+ \end{tikzpicture}
+\end{minipage}
+\hfill
+\begin{minipage}[T]{0.16\linewidth}
+ $\beta=90^\circ$:
+ \begin{itemize}
+ \item $\sin(\alpha)=\frac{a}{b}$
+ \item $\cos(\alpha)=\frac{c}{b}$
+ \item $\tan(\alpha)=\frac{a}{c}$
+ \end{itemize}
+\end{minipage}