diff options
| author | mzuenni <mzuenni@users.noreply.github.com> | 2024-07-28 22:54:40 +0200 |
|---|---|---|
| committer | GitHub <noreply@github.com> | 2024-07-28 22:54:40 +0200 |
| commit | 8d11c6c8213f46f0fa19826917c255edd5d43cb1 (patch) | |
| tree | 96d75baff33d5a04b5a60f1a41f514a26c716874 /content/geometry | |
| parent | 8c33b4e0d3030cfed17fc64b4fe41133339f6d87 (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.cpp | 12 | ||||
| -rw-r--r-- | content/geometry/circle.cpp | 33 | ||||
| -rw-r--r-- | content/geometry/closestPair.cpp | 27 | ||||
| -rw-r--r-- | content/geometry/convexHull.cpp | 18 | ||||
| -rw-r--r-- | content/geometry/delaunay.cpp | 124 | ||||
| -rw-r--r-- | content/geometry/formulas.cpp | 42 | ||||
| -rw-r--r-- | content/geometry/formulas3d.cpp | 53 | ||||
| -rw-r--r-- | content/geometry/geometry.tex | 62 | ||||
| -rw-r--r-- | content/geometry/hpi.cpp | 68 | ||||
| -rw-r--r-- | content/geometry/lines.cpp | 33 | ||||
| -rw-r--r-- | content/geometry/linesAndSegments.cpp | 89 | ||||
| -rw-r--r-- | content/geometry/polygon.cpp | 150 | ||||
| -rw-r--r-- | content/geometry/segmentIntersection.cpp | 63 | ||||
| -rw-r--r-- | content/geometry/sortAround.cpp | 11 | ||||
| -rw-r--r-- | content/geometry/spheres.cpp | 29 | ||||
| -rw-r--r-- | content/geometry/triangle.cpp | 43 | ||||
| -rw-r--r-- | content/geometry/triangle.tex | 41 |
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} |
