From 19436370fb48bc0a5e356d1ba713dc39b1a35d3a Mon Sep 17 00:00:00 2001 From: mzuenni Date: Tue, 6 Aug 2024 16:54:13 +0200 Subject: upper case INF --- content/geometry/hpi.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'content/geometry') diff --git a/content/geometry/hpi.cpp b/content/geometry/hpi.cpp index 3509e0e..bd3ab1e 100644 --- a/content/geometry/hpi.cpp +++ b/content/geometry/hpi.cpp @@ -1,4 +1,4 @@ -constexpr ll inf = 0x1FFF'FFFF'FFFF'FFFF;//THIS CODE IS WIP +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);} -- cgit v1.2.3 From fdb4d0c1b54c0987367069d57a0ee56d30243431 Mon Sep 17 00:00:00 2001 From: mzuenni Date: Sat, 10 Aug 2024 21:29:15 +0200 Subject: more tests --- content/geometry/formulas3d.cpp | 16 +-- test/geometry/formulas3d.cpp | 236 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 244 insertions(+), 8 deletions(-) create mode 100644 test/geometry/formulas3d.cpp (limited to 'content/geometry') diff --git a/content/geometry/formulas3d.cpp b/content/geometry/formulas3d.cpp index dee3ce8..63de2ce 100644 --- a/content/geometry/formulas3d.cpp +++ b/content/geometry/formulas3d.cpp @@ -26,23 +26,23 @@ int ccw(pt3 a, pt3 b, pt3 c, pt3 p) { return (orien > EPS) - (orien < -EPS); } -// Entfernung von Punkt p zur Ebene a,b,c. +// 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); + pt3 n = cross(b - a, c - a); + return abs(dot(n, a - p)) / abs(n); } -// Liegt p in der Ebene a,b,c? +// 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 +// 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); + pt3 n = cross(d - c, e - c); + pt3 dir = b - a; + return a + dir * dot(n, c - a) / dot(n, dir); } // Abstand zwischen der Grade a-b und c-d diff --git a/test/geometry/formulas3d.cpp b/test/geometry/formulas3d.cpp new file mode 100644 index 0000000..f6c04b7 --- /dev/null +++ b/test/geometry/formulas3d.cpp @@ -0,0 +1,236 @@ +#include "../util.h" +constexpr double EPS = 1e-6; +struct pt3 { + double x, y, z; + pt3 operator-(const pt3 o) const { + return {x-o.x, y-o.y, z-o.z}; + } + pt3 operator+(const pt3 o) const { + return {x+o.x, y+o.y, z+o.z}; + } + + pt3 operator*(double m) const { + return {x*m, y*m, z*m}; + } + + pt3 operator/(double d) const { + return {x/d, y/d, z/d}; + } + + friend ostream& operator<<(ostream& os, pt3 p) { + return os << "(" << p.x << ", " << p.y << ", " << p.z << ")"; + } +}; + +pt3 cross(pt3 a, pt3 b); +double abs(pt3 p); +double distToLine(pt3 a, pt3 b, pt3 p) { + return abs(cross(p - a, b - a)) / abs(b - a); +} + +#include + +namespace Random { + pt3 point3d(ll l, ll r) { + return { + (double)integer(l, r), + (double)integer(l, r), + (double)integer(l, r), + }; + } + + template + std::array distinct(ll l, ll r) { + std::array res = {}; + for (size_t i = 0; i < C; i++) { + bool unqiue; + do { + res[i] = point3d(l, r); + unqiue = true; + for (size_t j = 0; j < i; j++) { + unqiue &= res[j].x != res[i].x || + res[j].y != res[i].y || + res[j].z != res[i].z; + } + } while (!unqiue); + } + return res; + } +} + +void test_dot(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto p = Random::point3d(-range, range); + auto q = Random::point3d(-range, range); + + auto expected = p.x*q.x + p.y*q.y + p.z*q.z; + auto got = dot(p, q); + + if (got != expected) cerr << "got: " << got << ", expected: " << expected << FAIL; + queries++; + } + cerr << "tested dot: " << queries << endl; +} + +void test_cross(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto p = Random::point3d(-range, range); + auto q = Random::point3d(-range, range); + + auto expected = pt3{p.y*q.z - p.z*q.y, + p.z*q.x - p.x*q.z, + p.x*q.y - p.y*q.x}; + auto got = cross(p, q); + + if (got.x != expected.x) cerr << "error: x" << FAIL; + if (got.y != expected.y) cerr << "error: y" << FAIL; + if (got.z != expected.z) cerr << "error: z" << FAIL; + queries++; + } + cerr << "tested cross: " << queries << endl; +} + +void test_abs(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto p = Random::point3d(-range, range); + + auto expected = sqrt(p.x*p.x + p.y*p.y + p.z*p.z); + auto got = abs(p); + + if (abs(got - expected) > 1e-6) cerr << "got: " << got << ", expected: " << expected << FAIL; + queries++; + } + cerr << "tested abs: " << queries << endl; +} + +void test_mixed(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto a = Random::point3d(-range, range); + auto b = Random::point3d(-range, range); + auto c = Random::point3d(-range, range); + + auto expected = dot(cross(a, b), c); + auto got = mixed(a, b, c); + + if (got != expected) cerr << "got: " << got << ", expected: " << expected << FAIL; + queries++; + } + cerr << "tested mixed: " << queries << endl; +} + +void test_ccw(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto a = Random::point3d(-range, range); + auto b = Random::point3d(-range, range); + auto c = Random::point3d(-range, range); + auto d = Random::point3d(-range, range); + + ll expected = mixed(b - a, c - a, d - a); + if (expected < 0) expected = -1; + if (expected > 0) expected = 1; + auto got = ccw(a, b, c, d); + + if (got != expected) cerr << "got: " << got << ", expected: " << expected << FAIL; + queries++; + } + cerr << "tested ccw: " << queries << endl; +} + +void test_distToPlane(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto [p, q] = Random::distinct<2>(-range, range); + auto [a, b, c] = Random::distinct<3>(-range, range); + + auto norm = cross(a - c, b - c); + norm = norm / abs(norm); + + auto gotA = distToPlane(a, b, c, p); + auto gotB = distToPlane(a, b, c, q); + auto got = ccw(a, b, c, p) * ccw(a, b, c, q) < 0 ? gotA + gotB : abs(gotA - gotB); + + auto expected = abs(dot(norm, p - q)); + + if (abs(got - expected) > 1e-6) cerr << "got: " << got << ", expected: " << expected << FAIL; + queries++; + } + cerr << "tested distToPlane: " << queries << endl; +} + +void stress_pointOnPlane(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto p = Random::point3d(-range, range); + auto [a, b, c] = Random::distinct<3>(-range, range); + + bool expected = ccw(a, b, c, p) == 0; + bool got = pointOnPlane(a, b, c, p); + + if (got != expected) cerr << "error" << FAIL; + queries++; + } + cerr << "tested pointOnPlane: " << queries << endl; +} + +void test_linePlaneIntersection(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto [p, q] = Random::distinct<2>(-range, range); + auto [a, b, c] = Random::distinct<3>(-range, range); + + if (abs(dot(p - q, cross(a-c, b-c))) < 1e-6) continue; + + auto got = linePlaneIntersection(p, q, a, b, c); + + if (distToLine(p, q, got) > 1e-6) cerr << "error: 1" << FAIL; + if (distToPlane(a, b, c, got) > 1e-6) cerr << "error: 2" << FAIL; + queries++; + } + cerr << "tested linePlaneIntersection: " << queries << endl; +} + +void test_lineLineDist(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto [p, q] = Random::distinct<2>(-range, range); + auto [a, b] = Random::distinct<2>(-range, range); + + double expected = 0; + if (ccw(a, b, p, q) != 0) { + pt3 c = a + p - q; + expected = distToPlane(a, b, c, p); + } + auto got = lineLineDist(p, q, a, b); + + if (abs(got - expected) > 1e-6) cerr << "got: " << got << ", expected: " << expected << FAIL; + queries++; + } + cerr << "tested lineLineDist: " << queries << endl; +} + +int main() { + test_dot(100); + test_dot(1'000'000); + test_cross(100); + test_cross(1'000'000); + test_abs(100); + test_abs(1'000'000); + test_mixed(100); + test_mixed(1'000'000); + test_ccw(100); + test_ccw(1'000'000); + + test_distToPlane(100); + test_distToPlane(1'000'000); + stress_pointOnPlane(100); + stress_pointOnPlane(1'000'000); + test_linePlaneIntersection(100); + test_linePlaneIntersection(1'000);//can be far away + test_lineLineDist(100); + test_lineLineDist(1'000'000); +} -- cgit v1.2.3 From 6285c377d819cbcf1883d3496e7f7d461b29e171 Mon Sep 17 00:00:00 2001 From: mzuenni Date: Sat, 10 Aug 2024 22:41:19 +0200 Subject: more tests --- content/geometry/lines.cpp | 16 ++--- content/geometry/linesAndSegments.cpp | 2 +- tcr.pdf | Bin 691160 -> 691136 bytes test/geometry/lines.cpp | 107 ++++++++++++++++++++++++++++++++++ test/geometry/linesAndSegments.cpp | 2 +- 5 files changed, 117 insertions(+), 10 deletions(-) create mode 100644 test/geometry/lines.cpp (limited to 'content/geometry') diff --git a/content/geometry/lines.cpp b/content/geometry/lines.cpp index 95536a4..36de1db 100644 --- a/content/geometry/lines.cpp +++ b/content/geometry/lines.cpp @@ -1,10 +1,10 @@ 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)) {} + double a, b, c; // ax + by + c = 0; vertikale Gerade: b = 0 + line(pt p, pt q) : a(imag(p-q)), b(real(q-p)), c(cross({-b, a}, p)) {} }; -line pointsToLine(pt p1, pt p2) { - line l; +line pointsToLine(pt p1, pt p2) { // vertikale Gerade: b = 1 oder b = 0 + line l(0, 0); if (abs(real(p1 - p2)) < EPS) { l.a = 1; l.b = 0.0; l.c = -real(p1); } else { @@ -15,19 +15,19 @@ line pointsToLine(pt p1, pt p2) { return l; } -bool parallel(line l1, line l2) { +bool parallel(line l1, line l2) { // braucht b in {0, 1} return (abs(l1.a - l2.a) < EPS) && (abs(l1.b - l2.b) < EPS); } -bool same(line l1, line l2) { +bool same(line l1, line l2) { // braucht b in {0, 1} return parallel(l1, l2) && (abs(l1.c - l2.c) < EPS); } -bool intersect(line l1, line l2, pt& p) { +bool intersect(line l1, line l2, pt& res) { // braucht b in {0, 1} 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}; + res = {x, y}; return true; } diff --git a/content/geometry/linesAndSegments.cpp b/content/geometry/linesAndSegments.cpp index 1e21cba..ddab554 100644 --- a/content/geometry/linesAndSegments.cpp +++ b/content/geometry/linesAndSegments.cpp @@ -5,7 +5,7 @@ bool pointOnLine(pt a, pt b, pt p) { // 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; + return abs(cross(a - b, c - d)) > EPS; } // Berechnet den Schnittpunkt der Graden a-b und c-d. diff --git a/tcr.pdf b/tcr.pdf index 649aa24..5ebfcb9 100644 Binary files a/tcr.pdf and b/tcr.pdf differ diff --git a/test/geometry/lines.cpp b/test/geometry/lines.cpp new file mode 100644 index 0000000..7b1b99a --- /dev/null +++ b/test/geometry/lines.cpp @@ -0,0 +1,107 @@ +#include "../util.h" +constexpr double EPS = 1e-6; +#define ll double +double gcd(double x, double /**/) {return x;} //hacky +#include +#undef ll +#include +#include + +#include "../geometry.h" + +void stress_pointsToLine(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto [a, b] = Random::line(range); + + line gotA(a, b); + if (real(a) != real(b)) { + gotA.a /= gotA.b; + gotA.c /= gotA.b; + gotA.b /= gotA.b; + } else { + gotA.c /= gotA.a; + gotA.a /= gotA.a; + } + line gotB = pointsToLine(a, b); + + if (!same(gotA, gotB)) cerr << "error" << FAIL; + queries++; + } + cerr << "tested pointsToLine: " << queries << endl; +} + +void stress_same(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto [a, b] = Random::line(range); + auto [c, d] = Random::line(range); + + line lAB = pointsToLine(a, b); + line lCD = pointsToLine(c, d); + + auto got = same(lAB, lCD); + auto expected = pointOnLine(a, b, c) && pointOnLine(a, b, d); + + if (got != expected) cerr << "error" << FAIL; + queries++; + } + cerr << "tested same: " << queries << endl; +} + +void stress_parallel(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto [a, b] = Random::line(range); + auto [c, d] = Random::line(range); + + line lAB = pointsToLine(a, b); + line lCD = pointsToLine(c, d); + + auto got = parallel(lAB, lCD); + auto expected = cross(b-a, d-c) == 0; + + if (got != expected) cerr << "error" << FAIL; + queries++; + } + cerr << "tested parallel: " << queries << endl; +} + +void stress_intersect(ll range) { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto [a, b] = Random::line(range); + auto [c, d] = Random::line(range); + + line lAB = pointsToLine(a, b); + line lCD = pointsToLine(c, d); + + if (same(lAB, lCD)) continue; + + pt gotPT; + auto got = intersect(lAB, lCD, gotPT); + auto expected = lineIntersection(a, b, c, d); + + if (got != expected) cerr << "error: 1" << FAIL; + if (got) { + pt expectedPt = lineIntersection2(a, b, c, d); + if (float_error(real(gotPT), real(expectedPt)) > 1e-6) cerr << "error: 2" << FAIL; + if (float_error(imag(gotPT), imag(expectedPt)) > 1e-6) cerr << "error: 2" << FAIL; + } + queries++; + } + cerr << "tested intersect: " << queries << endl; +} + +int main() { + stress_pointsToLine(100); + stress_pointsToLine(100'000); + stress_same(10); + stress_same(100); + stress_same(1'000'000'000);//no precision issue since this will alwas be false... + stress_parallel(10); + stress_parallel(100); + stress_parallel(1'000'000'000); + stress_intersect(100); + stress_intersect(1'000'000'000); +} diff --git a/test/geometry/linesAndSegments.cpp b/test/geometry/linesAndSegments.cpp index 233546b..a2da3ba 100644 --- a/test/geometry/linesAndSegments.cpp +++ b/test/geometry/linesAndSegments.cpp @@ -30,7 +30,7 @@ void stress_lineIntersection(ll range) { auto [c, d] = Random::line(range); if (ccw(a, b, c) == 0 && ccw(a, b, d) == 0) continue; - bool expected = ccw(0, a-b, c-d) == 0; + bool expected = ccw(0, a-b, c-d) != 0; bool got = lineIntersection(a, b, c, d); if (got != expected) cerr << "error" << FAIL; -- cgit v1.2.3 From 497dc5e137b908e694c55bdd7a18842484939e7b Mon Sep 17 00:00:00 2001 From: mzuenni Date: Sat, 10 Aug 2024 23:38:14 +0200 Subject: more tests --- content/geometry/spheres.cpp | 30 +++++++++++++++--------------- test/geometry/spheres.cpp | 28 ++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+), 15 deletions(-) create mode 100644 test/geometry/spheres.cpp (limited to 'content/geometry') diff --git a/content/geometry/spheres.cpp b/content/geometry/spheres.cpp index ec22262..d34bca9 100644 --- a/content/geometry/spheres.cpp +++ b/content/geometry/spheres.cpp @@ -1,3 +1,16 @@ +// 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); + } +}; + // Great Circle Distance mit Längen- und Breitengrad. double gcDist(double pLat, double pLon, double qLat, double qLon, double radius) { @@ -11,19 +24,6 @@ double gcDist(double pLat, double pLon, } // Great Circle Distance mit kartesischen Koordinaten. -double gcDist(point p, point q) { +double gcDist(point p, point q) { // radius = 1 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); - } -}; +} \ No newline at end of file diff --git a/test/geometry/spheres.cpp b/test/geometry/spheres.cpp new file mode 100644 index 0000000..16e0ebd --- /dev/null +++ b/test/geometry/spheres.cpp @@ -0,0 +1,28 @@ +#include "../util.h" +constexpr double PI = acos(-1.0); +#pragma GCC diagnostic ignored "-Wshadow" +#include + +void test_consistent() { + ll queries = 0; + for (int tries = 0; tries < 100'000; tries++) { + auto pLat = Random::real(-180, 180); + auto pLon = Random::real(0, 360); + auto qLat = Random::real(-180, 180); + auto qLon = Random::real(0, 360); + + point p(pLat, pLon); + point q(qLat, qLon); + + auto gotA = gcDist(pLat, pLon, qLat, qLon, 1); + auto gotB = gcDist(p, q); + + if (abs(gotA - gotB) > 1e-6) cerr << "gotA: " << gotA << ", gotB: " << gotB << FAIL; + queries++; + } + cerr << "tested random: " << queries << endl; +} + +int main() { + test_consistent(); +} -- cgit v1.2.3 From f32a00178f0d3b2152a6fc1dc492c987aaede85f Mon Sep 17 00:00:00 2001 From: MZuenni Date: Tue, 13 Aug 2024 17:14:17 +0200 Subject: small improvements --- content/datastructures/dynamicConvexHull.cpp | 16 ++++++++-------- content/datastructures/pbds.cpp | 4 ++-- content/geometry/convexHull.cpp | 4 ++-- content/geometry/hpi.cpp | 2 +- content/graph/bitonicTSP.cpp | 2 +- content/graph/bitonicTSPsimple.cpp | 2 +- content/graph/pushRelabel.cpp | 2 +- tcr.pdf | Bin 691284 -> 691101 bytes test/datastructures/fenwickTree2.cpp | 2 +- test/geometry/delaunay.cpp | 6 +++--- test/math/gauss.cpp | 2 +- test/math/inversionsMerge.cpp | 2 +- 12 files changed, 22 insertions(+), 22 deletions(-) (limited to 'content/geometry') diff --git a/content/datastructures/dynamicConvexHull.cpp b/content/datastructures/dynamicConvexHull.cpp index 27ec898..abae1d7 100644 --- a/content/datastructures/dynamicConvexHull.cpp +++ b/content/datastructures/dynamicConvexHull.cpp @@ -1,22 +1,22 @@ struct Line { - mutable ll m, b, p; + mutable ll m, c, p; bool operator<(const Line& o) const {return m < o.m;} bool operator<(ll x) const {return p < x;} }; struct HullDynamic : multiset> { // max über Geraden - // (for doubles, use INF = 1/.0, div(a,b) = a/b) - ll div(ll a, ll b) {return a / b - ((a ^ b) < 0 && a % b);} + // (for doubles, use INF = 1/.0, div(a,c) = a/c) + ll div(ll a, ll c) {return a / b - ((a ^ c) < 0 && a % c);} bool isect(iterator x, iterator y) { if (y == end()) {x->p = INF; return false;} - if (x->m == y->m) x->p = x->b > y->b ? INF : -INF; - else x->p = div(y->b - x->b, x->m - y->m); + if (x->m == y->m) x->p = x->c > y->c ? INF : -INF; + else x->p = div(y->c - x->c, x->m - y->m); return x->p >= y->p; } - void add(ll m, ll b) { - auto x = insert({m, b, 0}); + void add(ll m, ll c) { + auto x = insert({m, c, 0}); while (isect(x, next(x))) erase(next(x)); if (x != begin()) { x--; @@ -31,6 +31,6 @@ struct HullDynamic : multiset> { // max über Geraden ll query(ll x) { auto l = *lower_bound(x); - return l.m * x + l.b; + return l.m * x + l.c; } }; diff --git a/content/datastructures/pbds.cpp b/content/datastructures/pbds.cpp index f0889a2..de0ace6 100644 --- a/content/datastructures/pbds.cpp +++ b/content/datastructures/pbds.cpp @@ -6,11 +6,11 @@ using Tree = tree, rb_tree_tag, // T.order_of_key(x): number of elements strictly less than x // *T.find_by_order(k): k-th element +constexpr uint64_t RNG = ll(2e18 * acos(-1)) | 199; // random odd template struct chash { - static const uint64_t C = ll(2e18 * acos(-1)) | 199; // random odd size_t operator()(T o) const { - return __builtin_bswap64(hash()(o) * C); + return __builtin_bswap64(hash()(o) * RNG); }}; template using hashMap = gp_hash_table>; diff --git a/content/geometry/convexHull.cpp b/content/geometry/convexHull.cpp index 6d89e05..1173924 100644 --- a/content/geometry/convexHull.cpp +++ b/content/geometry/convexHull.cpp @@ -11,8 +11,8 @@ vector convexHull(vector pts){ 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. + 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/hpi.cpp b/content/geometry/hpi.cpp index bd3ab1e..c58a6e7 100644 --- a/content/geometry/hpi.cpp +++ b/content/geometry/hpi.cpp @@ -1,4 +1,4 @@ -constexpr ll INF = 0x1FFF'FFFF'FFFF'FFFF;//THIS CODE IS WIP +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);} diff --git a/content/graph/bitonicTSP.cpp b/content/graph/bitonicTSP.cpp index 6470232..eee5082 100644 --- a/content/graph/bitonicTSP.cpp +++ b/content/graph/bitonicTSP.cpp @@ -27,5 +27,5 @@ auto bitonicTSP() { (lt.back() == 1 ? lt : ut).push_back(0); reverse(all(lt)); lt.insert(lt.end(), all(ut)); - return lt;// Enthält Knoten 0 zweimal. An erster und letzter Position. + return lt; // Enthält Knoten 0 zweimal. An erster und letzter Position. } diff --git a/content/graph/bitonicTSPsimple.cpp b/content/graph/bitonicTSPsimple.cpp index 8b6e6c5..cacfb9c 100644 --- a/content/graph/bitonicTSPsimple.cpp +++ b/content/graph/bitonicTSPsimple.cpp @@ -23,5 +23,5 @@ auto bitonicTSP() { rl.push_back(v); p2 = v; }} lr.insert(lr.end(), rl.rbegin(), rl.rend()); - return lr;// Enthält Knoten 0 zweimal. An erster und letzter Position. + return lr; // Enthält Knoten 0 zweimal. An erster und letzter Position. } diff --git a/content/graph/pushRelabel.cpp b/content/graph/pushRelabel.cpp index 73a9eae..ec36026 100644 --- a/content/graph/pushRelabel.cpp +++ b/content/graph/pushRelabel.cpp @@ -29,7 +29,7 @@ ll maxFlow(int s, int t) { cur.assign(n, 0); H.assign(n, 0); H[s] = n; - ec[t] = 1;//never set t to active... + ec[t] = 1; //never set t to active... vector co(2*n); co[0] = n - 1; for (Edge& e : adj[s]) addFlow(e, e.c); diff --git a/tcr.pdf b/tcr.pdf index b096438..57691c5 100644 Binary files a/tcr.pdf and b/tcr.pdf differ diff --git a/test/datastructures/fenwickTree2.cpp b/test/datastructures/fenwickTree2.cpp index aa99576..89d5b0f 100644 --- a/test/datastructures/fenwickTree2.cpp +++ b/test/datastructures/fenwickTree2.cpp @@ -9,7 +9,7 @@ void stress_test() { ll queries = 0; for (int tries = 0; tries < 100; tries++) { int n = Random::integer(10, 100); - vector naive(n);// = Random::integers(n, -1000, 1000); + vector naive = Random::integers(n, -1000, 1000); init(naive); for (int operations = 0; operations < 1000; operations++) { { diff --git a/test/geometry/delaunay.cpp b/test/geometry/delaunay.cpp index 7f8ec30..5740b95 100644 --- a/test/geometry/delaunay.cpp +++ b/test/geometry/delaunay.cpp @@ -16,11 +16,11 @@ vector convexHull(vector pts){ vector 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--;//allow collinear points! + while (k > t && cross(h[k-2], h[k-1], *it) < 0) k--; //allow collinear points! h[k++] = *it; }}; - half(all(pts), 1);// Untere Hülle. - half(next(pts.rbegin()), pts.rend(), k);// Obere Hülle. + 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/test/math/gauss.cpp b/test/math/gauss.cpp index 37bacce..6e4759d 100644 --- a/test/math/gauss.cpp +++ b/test/math/gauss.cpp @@ -14,7 +14,7 @@ vector> inverseMat(const vector>& m) { mat[i].resize(2*n); mat[i][n+i] = 1; } - gauss(n);//the unique cetc. checks are not usefull since we dont solve an lgs... + gauss(n); //the unique cetc. checks are not usefull since we dont solve an lgs... vector> res(m); for (int i = 0; i < n; i++) { res[i] = vector(mat[i].begin() + n, mat[i].end()); diff --git a/test/math/inversionsMerge.cpp b/test/math/inversionsMerge.cpp index 85ab0d2..acdc555 100644 --- a/test/math/inversionsMerge.cpp +++ b/test/math/inversionsMerge.cpp @@ -16,7 +16,7 @@ void stress_test() { for (ll i = 0; i < 100'000; i++) { int n = Random::integer(1, 100); vector v(n); - for (ll j = 0; j < n; j++) v[j] = (j-10) * 100000 + Random::integer(0, 10000);//values must be unique ): + for (ll j = 0; j < n; j++) v[j] = (j-10) * 100000 + Random::integer(0, 10000); //values must be unique ): shuffle(all(v), Random::rng); ll expected = naive(v); ll got = mergeSort(v); -- cgit v1.2.3 From e1803aee4e150088b0865954de476edf182a9ffe Mon Sep 17 00:00:00 2001 From: mzuenni Date: Thu, 22 Aug 2024 22:20:14 +0200 Subject: improve segmentIntersection --- content/geometry/segmentIntersection.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'content/geometry') diff --git a/content/geometry/segmentIntersection.cpp b/content/geometry/segmentIntersection.cpp index 4262ddc..acce7cc 100644 --- a/content/geometry/segmentIntersection.cpp +++ b/content/geometry/segmentIntersection.cpp @@ -18,8 +18,8 @@ struct event { 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); + if (imag(p) != imag(o.p)) return imag(p) < imag(o.p); + return type > o.type; } }; -- cgit v1.2.3 From 961706d93af4daf58d5991deea88981a01f962b0 Mon Sep 17 00:00:00 2001 From: mzuenni Date: Fri, 23 Aug 2024 21:36:51 +0200 Subject: more stuff --- content/geometry/segmentIntersection.cpp | 2 +- content/graph/dinic.cpp | 55 +++++++++++++++++++++++++++ content/string/duval.cpp | 8 ++-- tcr.pdf | Bin 691111 -> 691029 bytes test/graph/dinic.cpp | 62 +++++++++++++++++++++++++++++++ 5 files changed, 121 insertions(+), 6 deletions(-) create mode 100644 content/graph/dinic.cpp create mode 100644 test/graph/dinic.cpp (limited to 'content/geometry') diff --git a/content/geometry/segmentIntersection.cpp b/content/geometry/segmentIntersection.cpp index acce7cc..afc01b2 100644 --- a/content/geometry/segmentIntersection.cpp +++ b/content/geometry/segmentIntersection.cpp @@ -29,7 +29,7 @@ bool lessPT(const pt& a, const pt& b) { } bool intersect(const seg& a, const seg& b) { - return lineSegmentIntersection(a.a, a.b, b.a, b.b); + return segmentIntersection(a.a, a.b, b.a, b.b); //@\sourceref{geometry/linesAndSegments.cpp}@ } pair intersect(vector& segs) { diff --git a/content/graph/dinic.cpp b/content/graph/dinic.cpp new file mode 100644 index 0000000..2e58a2d --- /dev/null +++ b/content/graph/dinic.cpp @@ -0,0 +1,55 @@ +struct Edge { + int to, rev; + ll f, c; +}; + +vector> adj; +int s, t; +vector pt, dist; + +void addEdge(int u, int v, ll c) { + adj[u].push_back({v, (int)sz(adj[v]), 0, c}); + adj[v].push_back({u, (int)sz(adj[u]) - 1, 0, 0}); +} + +bool bfs() { + dist.assign(sz(adj), -1); + dist[s] = 0; + queue q({s}); + while (!q.empty() && dist[t] < 0) { + int v = q.front(); q.pop(); + for (Edge& e : adj[v]) { + if (dist[e.to] < 0 && e.c - e.f > 0) { + dist[e.to] = dist[v] + 1; + q.push(e.to); + }}} + return dist[t] >= 0; +} + +ll dfs(int v, ll flow = INF) { + if (v == t || flow == 0) return flow; + for (; pt[v] < sz(adj[v]); pt[v]++) { + Edge& e = adj[v][pt[v]]; + if (dist[e.to] != dist[v] + 1) continue; + ll cur = dfs(e.to, min(e.c - e.f, flow)); + if (cur > 0) { + e.f += cur; + adj[e.to][e.rev].f -= cur; + return cur; + }} + return 0; +} + +ll maxFlow(int source, int target) { + s = source, t = target; + ll flow = 0; + while (bfs()) { + pt.assign(sz(adj), 0); + ll cur; + do { + cur = dfs(s); + flow += cur; + } while (cur > 0); + } + return flow; +} diff --git a/content/string/duval.cpp b/content/string/duval.cpp index bf36cce..253bae1 100644 --- a/content/string/duval.cpp +++ b/content/string/duval.cpp @@ -6,9 +6,8 @@ vector> duval(const string& s) { if (s[k] < s[j]) k = i; else k++; } - while (i <= k) { + for (; i <= k; i += j - k) { res.push_back({i, i + j - k}); - i += j - k; }} return res; } @@ -16,6 +15,5 @@ vector> duval(const string& s) { int minrotation(const string& s) { auto parts = duval(s+s); for (auto [l, r] : parts) { - if (l < sz(s) && r >= sz(s)) { - return l; -}}} + if (r >= sz(s)) return l; +}} diff --git a/tcr.pdf b/tcr.pdf index 199c3de..2f24510 100644 Binary files a/tcr.pdf and b/tcr.pdf differ diff --git a/test/graph/dinic.cpp b/test/graph/dinic.cpp new file mode 100644 index 0000000..5af7c6f --- /dev/null +++ b/test/graph/dinic.cpp @@ -0,0 +1,62 @@ +#include "../util.h" +constexpr ll INF = LL::INF; +namespace dinic { +#include +} + +namespace pushRelabel { +#include +} + +void stress_test() { + ll queries = 0; + for (int tries = 0; tries < 20'000; tries++) { + int n = Random::integer(2, 30); + int m = Random::integer(n-1, max(n, min(500, n*(n-1) / 2 + 1))); + + dinic::adj.assign(n, {}); + pushRelabel::adj.assign(n, {}); + + Graph g(n); + g.erdosRenyi(m); + g.forEdges([](int a, int b){ + ll w = Random::integer(1, 1'000'000'000'000ll); + dinic::addEdge(a, b, w); + pushRelabel::addEdge(a, b, w); + }); + + ll got = dinic::maxFlow(0, n - 1); + ll expected = pushRelabel::maxFlow(0, n - 1); + + if (got != expected) cerr << "got: " << got << ", expected: " << expected << FAIL; + queries += n; + } + cerr << "tested random queries: " << queries << endl; +} + +constexpr int N = 50000; +constexpr int M = 200000; +void performance_test() { + using namespace dinic; + timer t; + Graph g(N); + g.erdosRenyi(M); + adj.assign(N, {}); + g.forEdges([](int a, int b){ + ll w1 = Random::integer(1, 1'000'000'000'000ll); + ll w2 = Random::integer(1, 1'000'000'000'000ll); + addEdge(a, b, w1); + addEdge(b, a, w2); + }); + + t.start(); + hash_t hash = maxFlow(0, N - 1); + t.stop(); + if (t.time > 2000) cerr << "too slow: " << t.time << FAIL; + cerr << "tested performance: " << t.time << "ms (hash: " << hash << ")" << endl; +} + +int main() { + stress_test(); + performance_test(); +} -- cgit v1.2.3 From 0257f0b3c61f203f64c3817dfe19a08f6191517c Mon Sep 17 00:00:00 2001 From: mzuenni Date: Wed, 11 Sep 2024 00:29:27 +0200 Subject: moved stuff --- content/geometry/hpi.cpp | 16 +++++++------- content/math/linearRecurrence.cpp | 14 ++++++------ content/math/math.tex | 7 +++++- content/math/recover.cpp | 13 +++++++++++ content/other/other.tex | 18 ++++++---------- content/other/recover.cpp | 13 ----------- tcr.pdf | Bin 698834 -> 695781 bytes test/geometry.h | 6 +++--- test/math/recover.cpp | 44 ++++++++++++++++++++++++++++++++++++++ test/other/recover.cpp | 44 -------------------------------------- 10 files changed, 87 insertions(+), 88 deletions(-) create mode 100644 content/math/recover.cpp delete mode 100644 content/other/recover.cpp create mode 100644 test/math/recover.cpp delete mode 100644 test/other/recover.cpp (limited to 'content/geometry') diff --git a/content/geometry/hpi.cpp b/content/geometry/hpi.cpp index c58a6e7..f3dc08d 100644 --- a/content/geometry/hpi.cpp +++ b/content/geometry/hpi.cpp @@ -27,22 +27,22 @@ struct hp { 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; + ll x = cross(a.dir(), b.dir()); + ll y = cross(b.from - a.from, b.dir()); + ptl i = mul(x, a.from) + mul(y, a.dir()); //intersect a and b + // check if i is outside/right of this + return imag(conj(mul(sgn(x),dir()))*(i-mul(x,from))) < 0; } }; constexpr ll lim = 2e9+7; deque intersect(vector hps) { - hps.push_back(hp(pt{lim+1,-1})); - hps.push_back(hp(pt{lim+1,1})); + hps.push_back(hp(pt{lim + 1, -1})); + hps.push_back(hp(pt{lim + 1, 1})); sort(all(hps)); - deque dq = {hp(pt{-lim, 1})}; + deque dq = {hp(pt{-lim - 1, 1})}; for (auto x : hps) { while (sz(dq) > 1 && x.check(dq.end()[-1], dq.end()[-2])) dq.pop_back(); diff --git a/content/math/linearRecurrence.cpp b/content/math/linearRecurrence.cpp index c15c25c..ab86f71 100644 --- a/content/math/linearRecurrence.cpp +++ b/content/math/linearRecurrence.cpp @@ -10,21 +10,21 @@ // return c; // } -ll kthTerm(const vector& f, const vector& c, ll k){ +ll kthTerm(const vector& f, const vector& c, ll k) { int n = sz(c); - vector q(n+1, 1); - for(int i = 1; i <= n; i++) q[i] = (mod-c[i-1])%mod; + vector q(n + 1, 1); + for (int i = 0; i < n; i++) q[i + 1] = (mod - c[i])%mod; vector p = mul(f, q); p.resize(n); p.push_back(0); - do{ + do { vector q2 = q; - for(int i = 1; i <= n; i += 2) q2[i] = (mod - q2[i]) % mod; + for (int i = 1; i <= n; i += 2) q2[i] = (mod - q2[i]) % mod; vector x = mul(p, q2), y = mul(q, q2); - for(int i = 0; i <= n; i++){ + for (int i = 0; i <= n; i++){ p[i] = i == n ? 0 : x[2*i + (k&1)]; q[i] = y[2*i]; } - }while(k /= 2); + } while (k /= 2); return p[0]; } \ No newline at end of file diff --git a/content/math/math.tex b/content/math/math.tex index fb66110..4ac6c9e 100644 --- a/content/math/math.tex +++ b/content/math/math.tex @@ -544,6 +544,11 @@ Wenn man $k$ Spiele in den Zuständen $X_1, \ldots, X_k$ hat, dann ist die \text \subsection{Wichtige Zahlen} \input{math/tables/composite} +\subsection{Recover $\boldsymbol{x}$ and $\boldsymbol{y}$ from $\boldsymbol{y}$ from $\boldsymbol{x\*y^{-1}}$ } +\method{recover}{findet $x$ und $y$ für $x=x\*y^{-1}\bmod m$}{\log(m)} +\textbf{WICHTIG:} $x$ und $y$ müssen kleiner als $\sqrt{\nicefrac{m}{2}}$ sein! +\sourcecode{math/recover.cpp} + \optional{ \subsection{Primzahlzählfunktion $\boldsymbol{\pi}$} \begin{methods} @@ -552,10 +557,10 @@ Wenn man $k$ Spiele in den Zuständen $X_1, \ldots, X_k$ hat, dann ist die \text \method{pi}{zählt Primzahlen $\leq n$ ($n < N^2$)}{n^{2/3}} \end{methods} \sourcecode{math/piLehmer.cpp} -} \subsection{Primzahlzählfunktion $\boldsymbol{\pi}$} \sourcecode{math/piLegendre.cpp} +} \begin{algorithm}[optional]{Big Integers} \sourcecode{math/bigint.cpp} diff --git a/content/math/recover.cpp b/content/math/recover.cpp new file mode 100644 index 0000000..1a593f0 --- /dev/null +++ b/content/math/recover.cpp @@ -0,0 +1,13 @@ +ll sq(ll x) {return x*x;} + +array recover(ll c, ll m) { + array u = {m, 0}, v = {c, 1}; + while (m <= 2 * sq(v[0])) { + ll q = u[0] / v[0]; + u[0] -= q * v[0]; + u[1] -= q * v[1]; + swap(u, v); + } + if (v[1] <= 0 || 2 * sq(v[1]) >= m) return {-1, -1}; + return v; +} diff --git a/content/other/other.tex b/content/other/other.tex index 368d0b3..191a6da 100644 --- a/content/other/other.tex +++ b/content/other/other.tex @@ -72,6 +72,12 @@ \sourcecode{other/sos.cpp} \end{algorithm} +\begin{algorithm}{Fast Subset Sum} + \method{fastSubsetSum}{findet maximale subset $\mathit{sum}\leq t$}{n \cdot A} + Die Laufzeit hängt vom maximalen Wert $A$ in der Menge ab. + \sourcecode{other/fastSubsetSum.cpp} +\end{algorithm} + \begin{algorithm}{Parallel Binary Search} \sourcecode{other/pbs.cpp} \end{algorithm} @@ -95,18 +101,6 @@ \textbf{Beachte bei der Ausgabe, dass die Personen im ersten Fall von $\boldsymbol{1, \ldots, n}$ nummeriert sind, im zweiten Fall von $\boldsymbol{0, \ldots, n-1}$!} \end{algorithm} -\vfill\null\columnbreak - -\subsection{Recover $\boldsymbol{x}$ and $\boldsymbol{y}$ from $\boldsymbol{y}$ from $\boldsymbol{x\*y^{-1}}$ } -\method{recover}{findet $x$ und $y$ für $x=x\*y^{-1}\bmod m$}{\log(m)} -\textbf{WICHTIG:} $x$ und $y$ müssen kleiner als $\sqrt{\nicefrac{m}{2}}$ sein! -\sourcecode{other/recover.cpp} - -\subsection{Fast Subset Sum} -\method{fastSubsetSum}{findet maximale subset sum $\leq t$}{n \cdot A} -Die Laufzeit hängt vom maximalen Wert $A$ in der Menge ab. -\sourcecode{other/fastSubsetSum.cpp} - \begin{algorithm}[optional]{Zeileneingabe} \sourcecode{other/split.cpp} \end{algorithm} diff --git a/content/other/recover.cpp b/content/other/recover.cpp deleted file mode 100644 index 1a593f0..0000000 --- a/content/other/recover.cpp +++ /dev/null @@ -1,13 +0,0 @@ -ll sq(ll x) {return x*x;} - -array recover(ll c, ll m) { - array u = {m, 0}, v = {c, 1}; - while (m <= 2 * sq(v[0])) { - ll q = u[0] / v[0]; - u[0] -= q * v[0]; - u[1] -= q * v[1]; - swap(u, v); - } - if (v[1] <= 0 || 2 * sq(v[1]) >= m) return {-1, -1}; - return v; -} diff --git a/tcr.pdf b/tcr.pdf index d383dd9..6f156e4 100644 Binary files a/tcr.pdf and b/tcr.pdf differ diff --git a/test/geometry.h b/test/geometry.h index 7886fe2..0167d5c 100644 --- a/test/geometry.h +++ b/test/geometry.h @@ -1,6 +1,6 @@ -#include - namespace details { + #include + // Liegt p auf der Strecke a-b? bool pointInLineSegment(pt a, pt b, pt p) { if (cross(a, b, p) != 0) return false; @@ -59,7 +59,7 @@ namespace Random { for (size_t i = 0; i < dirs.size(); i++) { dirs[i] = pt(x[i], y[i]); } - sortAround(0, dirs); + details::sortAround(0, dirs); vector res = {{0, 0}}; ll maxX = 0; diff --git a/test/math/recover.cpp b/test/math/recover.cpp new file mode 100644 index 0000000..72853e5 --- /dev/null +++ b/test/math/recover.cpp @@ -0,0 +1,44 @@ +#include "../util.h" +#include +#include + +void stress_test() { + ll queries = 0; + timer t; + for (int i = 0; i < 500; i++) { + ll p = Random::prime(10000); + for (ll j = 0; 2*j*j < p; j++) { + for (ll b = 1; 2*b*b < p; b++) { + if (gcd(j, b) != 1) continue; + for (ll a : {-j, j}) { + ll c = a * multInv(b, p); + + t.start(); + auto [x, y] = recover(c, p); + t.stop(); + + if (a != x || b != y) cerr << "got: " << x << "/" << y << ", expected: " << a << "/" << b << FAIL; + queries++; + } + } + } + for (ll c = 0; c < p; c++) { + t.start(); + auto [x, y] = recover(c, p); + t.stop(); + + if (y < 0) continue; + if (y == 0) cerr << "error: y=0" << FAIL; + ll got = (((x * multInv(y, p)) % p) + p) % p; + if (got != c) cerr << "got: " << got << ", expected: " << c << FAIL; + queries++; + } + } + cerr << "tested random queries: " << queries << endl; + if (t.time > 500) cerr << "too slow: " << t.time << FAIL; + cerr << "tested performance: " << t.time << "ms" << endl; +} + +int main() { + stress_test(); +} diff --git a/test/other/recover.cpp b/test/other/recover.cpp deleted file mode 100644 index 72853e5..0000000 --- a/test/other/recover.cpp +++ /dev/null @@ -1,44 +0,0 @@ -#include "../util.h" -#include -#include - -void stress_test() { - ll queries = 0; - timer t; - for (int i = 0; i < 500; i++) { - ll p = Random::prime(10000); - for (ll j = 0; 2*j*j < p; j++) { - for (ll b = 1; 2*b*b < p; b++) { - if (gcd(j, b) != 1) continue; - for (ll a : {-j, j}) { - ll c = a * multInv(b, p); - - t.start(); - auto [x, y] = recover(c, p); - t.stop(); - - if (a != x || b != y) cerr << "got: " << x << "/" << y << ", expected: " << a << "/" << b << FAIL; - queries++; - } - } - } - for (ll c = 0; c < p; c++) { - t.start(); - auto [x, y] = recover(c, p); - t.stop(); - - if (y < 0) continue; - if (y == 0) cerr << "error: y=0" << FAIL; - ll got = (((x * multInv(y, p)) % p) + p) % p; - if (got != c) cerr << "got: " << got << ", expected: " << c << FAIL; - queries++; - } - } - cerr << "tested random queries: " << queries << endl; - if (t.time > 500) cerr << "too slow: " << t.time << FAIL; - cerr << "tested performance: " << t.time << "ms" << endl; -} - -int main() { - stress_test(); -} -- cgit v1.2.3 From a6f428e91457dae36efda91fbdc1eb4f0617132d Mon Sep 17 00:00:00 2001 From: mzuenni Date: Thu, 12 Sep 2024 00:25:38 +0200 Subject: improve halfplane intersection tests --- content/geometry/hpi.cpp | 2 +- tcr.pdf | Bin 695518 -> 695511 bytes test/geometry/hpi.cpp | 240 +++++++++++++++++++++++++++++++++++++++++------ 3 files changed, 212 insertions(+), 30 deletions(-) (limited to 'content/geometry') diff --git a/content/geometry/hpi.cpp b/content/geometry/hpi.cpp index f3dc08d..02c71e3 100644 --- a/content/geometry/hpi.cpp +++ b/content/geometry/hpi.cpp @@ -60,7 +60,7 @@ deque intersect(vector hps) { 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])) + while (sz(dq) > 2 && dq.back().check(dq[0], dq[1])) dq.pop_front(); if (sz(dq) < 3) return {}; diff --git a/tcr.pdf b/tcr.pdf index 47a96a2..7147158 100644 Binary files a/tcr.pdf and b/tcr.pdf differ diff --git a/test/geometry/hpi.cpp b/test/geometry/hpi.cpp index e61178a..a2326bc 100644 --- a/test/geometry/hpi.cpp +++ b/test/geometry/hpi.cpp @@ -11,30 +11,168 @@ ll sgn(ll x) { } #include -using ptd = complex; -ptd convert(pt x) {return ptd(real(x), imag(x));} -auto cross(ptd a, ptd b) {return imag(conj(a) * b);} -auto cross(ptd p, ptd a, ptd b) {return cross(a - p, b - p);} -ptd lineIntersection2(ptd a, ptd b, ptd c, ptd d) { - ld x = cross(b - a, d - c); - ld y = cross(c - a, d - c); - return a + y/x*(b - a); +//https://cp-algorithms.com/geometry/halfplane-intersection.html +namespace cpalgo { + // Redefine epsilon and infinity as necessary. Be mindful of precision errors. + const long double eps = 1e-9, inf = 1e9; + + // Basic point/vector struct. + struct Point { + + long double x, y; + explicit Point(long double x_ = 0, long double y_ = 0) : x(x_), y(y_) {} + Point(pt p) : x(real(p)), y(imag(p)) {} + + // Addition, substraction, multiply by constant, dot product, cross product. + + friend Point operator + (const Point& p, const Point& q) { + return Point(p.x + q.x, p.y + q.y); + } + + friend Point operator - (const Point& p, const Point& q) { + return Point(p.x - q.x, p.y - q.y); + } + + friend Point operator * (const Point& p, const long double& k) { + return Point(p.x * k, p.y * k); + } + + friend long double dot(const Point& p, const Point& q) { + return p.x * q.x + p.y * q.y; + } + + friend long double cross(const Point& p, const Point& q) { + return p.x * q.y - p.y * q.x; + } + + friend std::ostream& operator<<(std::ostream& os, const Point& p) { + return os << "(" << p.x << ", " << p.y << ")"; + } + + + }; + + // Basic half-plane struct. + struct Halfplane { + + // 'p' is a passing point of the line and 'pq' is the direction vector of the line. + Point p, pq; + long double angle; + + Halfplane() {} + Halfplane(const Point& a, const Point& b) : p(a), pq(b - a) { + angle = atan2l(pq.x, -pq.y);//patched to get same sort + } + Halfplane(array ps) : Halfplane(ps[0], ps[1]) {} + Halfplane(hp h) : Halfplane(h.from, h.to) {} + + // Check if point 'r' is outside this half-plane. + // Every half-plane allows the region to the LEFT of its line. + bool out(const Point& r) { + return cross(pq, r - p) < -eps; + } + + // Comparator for sorting. + bool operator < (const Halfplane& e) const { + return angle < e.angle; + } + + // Intersection point of the lines of two half-planes. It is assumed they're never parallel. + friend Point inter(const Halfplane& s, const Halfplane& t) { + long double alpha = cross((t.p - s.p), t.pq) / cross(s.pq, t.pq); + return s.p + (s.pq * alpha); + } + + friend std::ostream& operator<<(std::ostream& os, const Halfplane& hp) { + return os << hp.p << " " << hp.p+hp.pq; + } + }; + + // Actual algorithm + vector hp_intersect(vector& H) { + + /*Point box[4] = { // Bounding box in CCW order + Point(inf, inf), + Point(-inf, inf), + Point(-inf, -inf), + Point(inf, -inf) + }; + + for(int i = 0; i<4; i++) { // Add bounding box half-planes. + Halfplane aux(box[i], box[(i+1) % 4]); + H.push_back(aux); + }*/ + // Add bounding box half-planes, fixed: + H.push_back({pt( 1e9, 1e9), pt( 1e9-1, 1e9 )}); + H.push_back({pt(-1e9, 1e9), pt(-1e9 , 1e9-1)}); + H.push_back({pt(-1e9, -1e9), pt(-1e9+1, -1e9 )}); + H.push_back({pt( 1e9, -1e9), pt( 1e9 , -1e9+1)}); + + // Sort by angle and start algorithm + sort(H.begin(), H.end()); + deque dq; + int len = 0; + for(int i = 0; i < int(H.size()); i++) { + + // Remove from the back of the deque while last half-plane is redundant + while (len > 1 && H[i].out(inter(dq[len-1], dq[len-2]))) { + dq.pop_back(); + --len; + } + + // Remove from the front of the deque while first half-plane is redundant + while (len > 1 && H[i].out(inter(dq[0], dq[1]))) { + dq.pop_front(); + --len; + } + + // Special case check: Parallel half-planes + if (len > 0 && fabsl(cross(H[i].pq, dq[len-1].pq)) < eps) { + // Opposite parallel half-planes that ended up checked against each other. + if (dot(H[i].pq, dq[len-1].pq) < 0.0) + return vector(); + + // Same direction half-plane: keep only the leftmost half-plane. + if (H[i].out(dq[len-1].p)) { + dq.pop_back(); + --len; + } + else continue; + } + + // Add new half-plane + dq.push_back(H[i]); + ++len; + } + + // Final cleanup: Check half-planes at the front against the back and vice-versa + while (len > 2 && dq[0].out(inter(dq[len-1], dq[len-2]))) { + dq.pop_back(); + --len; + } + + while (len > 2 && dq[len-1].out(inter(dq[0], dq[1]))) { + dq.pop_front(); + --len; + } + + // Report empty intersection if necessary + if (len < 3) return vector(); + + // Reconstruct the convex polygon from the remaining half-planes. + vector ret(len); + for(int i = 0; i+1 < len; i++) { + ret[i] = inter(dq[i], dq[i+1]); + } + ret.back() = inter(dq[len-1], dq[0]); + return ret; + } } //check if a is dominated by b and c -bool naiveCheck(array a, array b, array c) { - //https://cp-algorithms.com/geometry/halfplane-intersection.html - //intersect b and c - //check cross of a and intersection - ptd intersection = lineIntersection2( - convert(b[0]), - convert(b[1]), - convert(c[0]), - convert(c[1]) - ); - return cross(convert(a[0]), convert(a[1]), intersection) <= -1e-12; +bool naiveCheck(cpalgo::Halfplane a, cpalgo::Halfplane b, cpalgo::Halfplane c) { + return a.out(inter(b, c)); } -//real(a - p)*imag(b - p)-imag(a - p)*real(b - p) void test_check(ll range) { ll queries = 0; @@ -42,10 +180,7 @@ void test_check(ll range) { auto a = Random::line(range); auto b = Random::line(range); auto c = b; - while ( - cross(b[0] - b[1], c[0] - c[1]) == 0 - //||cross(a[0] - a[1], b[0] - b[1], c[0] - c[1]) >= 0 - ) c = Random::line(range); + while (cross(b[0] - b[1], c[0] - c[1]) == 0) c = Random::line(range); bool got = hp(a[0], a[1]).check(hp(b[0], b[1]), hp(c[0], c[1])); bool expected = naiveCheck(a, b, c); @@ -63,14 +198,61 @@ void test_check(ll range) { cerr << "tested random queries: " << queries << endl; } -/*void stress_test(ll range) { +void stress_test(ll range) { ll queries = 0; for (int tries = 0; tries < 100'000; tries++) { - auto centre = Random::point(n, -range / 2, range / 2); + auto centre = Random::point(-range / 2, range / 2); + int n = Random::integer(3, 30); + + vector hps1 = {//+cpalgo bounding box + {pt( 1e9, 1e9), pt(-1e9, 1e9)}, + {pt(-1e9, 1e9), pt(-1e9, -1e9)}, + {pt(-1e9, -1e9), pt( 1e9, -1e9)}, + {pt( 1e9, -1e9), pt( 1e9, 1e9)}, + }; + + vector hps2; + for (int i = 0; i < n; i++) { + auto [a, b] = Random::line(range); + if (cross(a, b, centre) < 0) swap(a, b); + hps1.emplace_back(a, b); + hps2.emplace_back(a, b); + } + + auto expected = cpalgo::hp_intersect(hps2); + auto gotHp = intersect(hps1); + if (sz(gotHp) == 1) cerr << "WHAT?" << FAIL; + for (hp h : gotHp) { + if (h.dummy()) cerr << "DUMMY!" << FAIL;//we added a bounding box... + } + vector got; + if (!gotHp.empty()) { + gotHp.push_back(gotHp[0]); + for (int i = 0; i + 1 < sz(gotHp); i++) { + got.push_back(inter(cpalgo::Halfplane(gotHp[i]), cpalgo::Halfplane(gotHp[i + 1]))); + } + } + + bool eq = sz(got) == sz(expected); + for (ll i = 0; eq && i < sz(got); i++) { + eq &= float_error(got[i].x, expected[i].x) < 1e-6; + eq &= float_error(got[i].y, expected[i].y) < 1e-6; + } + if (!eq) { + cerr << tries << endl; + cerr << setprecision(20) << endl; + for (auto tmp : hps2) cerr << tmp << endl; + cerr << endl; + for (auto tmp : expected) cerr << tmp << endl; + cerr << endl; + for (auto tmp : got) cerr << tmp << endl; + cerr << FAIL; + } + queries += n; } cerr << "tested random queries: " << queries << endl; -}*/ +} /*constexpr int N = 1'000'000; void performance_test() { @@ -103,7 +285,7 @@ void performance_test() { int main() { test_check(10); test_check(100); - test_check(1000); - test_check(10000); + stress_test(10); + stress_test(100); //performance_test(); } -- cgit v1.2.3 From 35d485bcf6a9ed0a9542628ce4aa94a3326d0884 Mon Sep 17 00:00:00 2001 From: Yidi Date: Wed, 25 Sep 2024 20:18:42 +0200 Subject: fix function name --- content/geometry/polygon.cpp | 2 +- test/geometry/polygon.cpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) (limited to 'content/geometry') diff --git a/content/geometry/polygon.cpp b/content/geometry/polygon.cpp index 3178290..064d81f 100644 --- a/content/geometry/polygon.cpp +++ b/content/geometry/polygon.cpp @@ -29,7 +29,7 @@ bool inside(pt p, const vector& 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 (pointOnSegment(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) { diff --git a/test/geometry/polygon.cpp b/test/geometry/polygon.cpp index 1dd46ca..643ea70 100644 --- a/test/geometry/polygon.cpp +++ b/test/geometry/polygon.cpp @@ -10,7 +10,7 @@ double abs(pt p) { return hypot(real(p), imag(p)); } // Liegt p auf der Strecke a-b? -bool pointOnLineSegment(pt a, pt b, pt p) { +bool pointOnSegment(pt a, pt b, pt p) { if (cross(a, b, p) != 0) return false; auto dist = norm(a - b); return norm(a - p) <= dist && norm(b - p) <= dist; @@ -65,7 +65,7 @@ void test_windingNumber(ll range) { int cur = details::lineSegmentIntersection(p, p + pt(1, 2'000'000'007), ps[j], ps[j+1]); if (ptLess(ps[j], ps[j+1])) expected -= cur; else expected += cur; - onBorder |= pointOnLineSegment(ps[j], ps[j+1], p); + onBorder |= pointOnSegment(ps[j], ps[j+1], p); } if (onBorder) continue; if (area(ps) < 0) expected = -expected; @@ -93,7 +93,7 @@ void test_inside(ll range) { bool onBorder = false; for (int j = 0; j < n; j++) { count += details::lineSegmentIntersection(p, p + pt(1, 2'000'000'007), ps[j], ps[j+1]); - onBorder |= pointOnLineSegment(ps[j], ps[j+1], p); + onBorder |= pointOnSegment(ps[j], ps[j+1], p); } bool expected = (count % 2) && !onBorder; bool got = inside(p, ps); -- cgit v1.2.3