From 9bd7638dff1ca1c93953596ab7b44b25dea47748 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebasti=C3=A1n=20Torrealba?= Date: Sun, 1 Mar 2026 17:26:56 -0300 Subject: [PATCH 1/3] Add fast prime sieve (min_25) to WF and normal tracker --- content/maths/number-theory/fast-sieve.cpp | 94 ++++++++++++++++++++++ tracker-wf.yaml | 1 + tracker.yaml | 1 + 3 files changed, 96 insertions(+) create mode 100644 content/maths/number-theory/fast-sieve.cpp diff --git a/content/maths/number-theory/fast-sieve.cpp b/content/maths/number-theory/fast-sieve.cpp new file mode 100644 index 0000000..b5b51c3 --- /dev/null +++ b/content/maths/number-theory/fast-sieve.cpp @@ -0,0 +1,94 @@ +/* + *Author:* min\_25 + *Description:* Fast prime sieve up to $N$. + Returns a sorted vector of all primes $<= N$. + Uses a segmented sieve with wheel factorization (mod 30). + *Time:* ~0.5s for $N = 10^9$. + *Status:* Tested on Yosupo (enumerate\_primes) +*/ +vector sieve(int N, int Q = 17, int L = 1 << 15) { + static const int rs[] = {1, 7, 11, 13, 17, 19, 23, 29}; + struct P { + P(int p) : p(p) {} + int p, pos[8]; + }; + + auto prime_count = [](int N) -> int { + return N > 60184 ? N / (log(N) - 1.1) : max(1., N / (log(N) - 1.11)) + 1; + }; + + int v = sqrt(N), vv = sqrt(v); + vector isp(v + 1, true); + for (int i = 2; i <= vv; i++) + if (isp[i]) + for (int j = i * i; j <= v; j += i) + isp[j] = false; + + int rsize = prime_count(N + 30); + vector primes = {2, 3, 5}; + int psize = 3; + primes.resize(rsize); + + vector

sprimes; + size_t pbeg = 0; + int prod = 1; + + for (int p = 7; p <= v; p++) { + if (!isp[p]) + continue; + if (p <= Q) + prod *= p, pbeg++, primes[psize++] = p; + auto pp = P(p); + for (int t = 0; t < 8; t++) { + int j = (p <= Q) ? p : p * p; + while (j % 30 != rs[t]) + j += p << 1; + pp.pos[t] = j / 30; + } + sprimes.push_back(pp); + } + + vector pre(prod, 0xFF); + for (size_t pi = 0; pi < pbeg; pi++) { + auto pp = sprimes[pi]; + int p = pp.p; + for (int t = 0; t < 8; t++) { + unsigned char m = ~(1 << t); + for (int i = pp.pos[t]; i < prod; i += p) + pre[i] &= m; + } + } + + int block_size = (L + prod - 1) / prod * prod; + vector block(block_size); + unsigned char *pblock = block.data(); + int M = (N + 29) / 30; + + for (int beg = 0; beg < M; beg += block_size, pblock -= block_size) { + int end = min(M, beg + block_size); + for (int i = beg; i < end; i += prod) + copy(pre.begin(), pre.end(), pblock + i); + if (beg == 0) + pblock[0] &= 0xFE; + for (size_t pi = pbeg; pi < sprimes.size(); pi++) { + auto &pp = sprimes[pi]; + int p = pp.p; + for (int t = 0; t < 8; t++) { + int i = pp.pos[t]; + unsigned char m = ~(1 << t); + for (; i < end; i += p) + pblock[i] &= m; + pp.pos[t] = i; + } + } + for (int i = beg; i < end; i++) + for (int m = pblock[i]; m > 0; m &= m - 1) + primes[psize++] = i * 30 + rs[__builtin_ctz(m)]; + } + + assert(psize <= rsize); + while (psize > 0 && primes[psize - 1] > N) + psize--; + primes.resize(psize); + return primes; +} diff --git a/tracker-wf.yaml b/tracker-wf.yaml index 4b3ae4a..f7c354d 100644 --- a/tracker-wf.yaml +++ b/tracker-wf.yaml @@ -99,6 +99,7 @@ content: - counting-divisors.cpp - discrete-log.cpp - erathostenes-sieve.cpp + - fast-sieve.cpp - euler-phi.cpp - euler-phi-sieve.cpp - mobius-sieve.cpp diff --git a/tracker.yaml b/tracker.yaml index 0a10bcb..513385a 100644 --- a/tracker.yaml +++ b/tracker.yaml @@ -76,6 +76,7 @@ content: - divisors.cpp - discrete-log.cpp - erathostenes-sieve.cpp + - fast-sieve.cpp - euler-phi.cpp - euler-phi-sieve.cpp - mobius-sieve.cpp From 3aea1d57844bedc4592ba74d43c40041272f04ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebasti=C3=A1n=20Torrealba?= Date: Mon, 2 Mar 2026 18:06:56 -0300 Subject: [PATCH 2/3] Add Delaunay triangulation and Delaunay graph templates - Add delaunay.cpp: Delaunay triangulation via paraboloid lifting - Add delaunay-graph.cpp: planar graph with O(sqrt(n)) greedy walk - Update convex-hull-3d.cpp: exact integer arithmetic, degenerate flag, build_hull() method - Update point-3d.cpp: add nonzero() method - Update both trackers --- .../geometry/algorithms/convex-hull-3d.cpp | 97 ++++++++----------- .../geometry/algorithms/delaunay-graph.cpp | 39 ++++++++ content/geometry/algorithms/delaunay.cpp | 37 +++++++ content/geometry/primitives/point-3d.cpp | 1 + tracker-wf.yaml | 2 + tracker.yaml | 2 + 6 files changed, 123 insertions(+), 55 deletions(-) create mode 100644 content/geometry/algorithms/delaunay-graph.cpp create mode 100644 content/geometry/algorithms/delaunay.cpp diff --git a/content/geometry/algorithms/convex-hull-3d.cpp b/content/geometry/algorithms/convex-hull-3d.cpp index 61017d7..f021fb1 100644 --- a/content/geometry/algorithms/convex-hull-3d.cpp +++ b/content/geometry/algorithms/convex-hull-3d.cpp @@ -1,57 +1,52 @@ /* - *Description:* 3D Convex Hull using randomized incremental algorithm + *Description:* 3D Convex Hull using randomized incremental algorithm. Exact integer arithmetic. Sets `degenerate=true` if less than 4 non-coplanar points. *Time:* $O(n log n)$ expected *Status:* Tested */ template struct convex_hull_3d { using pt = point3d; - const ld EPS = 1e-9; struct edge; struct face { int a, b, c; - pt pa, pb, pc, q; + pt q; edge *e1, *e2, *e3; vector points; int dead = 1e9; face(int a, int b, int c, pt pa, pt pb, pt pc) - : a(a), b(b), c(c), pa(pa), pb(pb), pc(pc) { - q = ((pb - pa)^(pc - pa)); + : a(a), b(b), c(c) { + q = (pb - pa) ^ (pc - pa); e1 = e2 = e3 = nullptr; } }; - - struct edge { - edge *rev; face *f; - }; + struct edge { edge *rev; face *f; }; vector p; vector f; + bool degenerate = false; static void glue(face *F1, face *F2, edge *&e1, edge *&e2) { - e1 = new edge, e2 = new edge; - e1->rev = e2, e2->rev = e1; - e1->f = F2, e2->f = F1; + e1 = new edge; e2 = new edge; + e1->rev = e2; e2->rev = e1; + e1->f = F2; e2->f = F1; } void prepare() { int n = (int)p.size(); - mt19937 rng(chrono::steady_clock::now().time_since_epoch().count()); + mt19937 rng(42); shuffle(p.begin(), p.end(), rng); vector ve = {0}; - for (int i = 1; i < n; i++) { + for (int i = 1; i < n && (int)ve.size() < 4; i++) { if (ve.size() == 1) { - if ((p[ve[0]] - p[i]).abs() > EPS) ve.push_back(i); + if ((p[ve[0]] - p[i]).nonzero()) ve.push_back(i); } else if (ve.size() == 2) { - if ((p[ve[1]] - p[ve[0]]).cross(p[i] - p[ve[0]]).abs() > EPS) + if (((p[ve[1]] - p[ve[0]]) ^ (p[i] - p[ve[0]])).nonzero()) + ve.push_back(i); + } else { + if ((p[i] - p[ve[0]]).dot((p[ve[1]] - p[ve[0]]) ^ (p[ve[2]] - p[ve[0]])) != 0) ve.push_back(i); - } else if (abs((p[i] - p[ve[0]]) - .dot((p[ve[1]] - p[ve[0]]) - .cross(p[ve[2]] - p[ve[0]]))) > EPS) { - ve.push_back(i); - break; } } - assert((int)ve.size() == 4); + if ((int)ve.size() < 4) { degenerate = true; return; } vector ve2; for (int i : ve) ve2.push_back(p[i]); reverse(ve.begin(), ve.end()); @@ -59,43 +54,40 @@ template struct convex_hull_3d { p.insert(p.begin(), ve2.begin(), ve2.end()); } - convex_hull_3d(vector points) : p(move(points)) { - int n = p.size(); + convex_hull_3d() {} + convex_hull_3d(vector points) : p(move(points)) { build_hull(); } + + void build_hull() { + int n = (int)p.size(); prepare(); + if (degenerate) return; vector new_face(n, nullptr); vector> conflict(n); - auto add_face = [&](int a, int b, int c) { face *F = new face(a, b, c, p[a], p[b], p[c]); - f.push_back(F); - return F; + f.push_back(F); return F; }; - face *F1 = add_face(0, 1, 2); face *F2 = add_face(0, 2, 1); glue(F1, F2, F1->e1, F2->e3); glue(F1, F2, F1->e2, F2->e2); glue(F1, F2, F1->e3, F2->e1); - for (int i = 3; i < n; i++) { for (face *F : {F1, F2}) { T Q = (p[i] - p[F->a]).dot(F->q); - if (Q > EPS) conflict[i].push_back(F); - if (Q >= -EPS) F->points.push_back(i); + if (Q > 0) conflict[i].push_back(F); + if (Q >= 0) F->points.push_back(i); } } - for (int i = 3; i < n; i++) { - for (face *F : conflict[i]) - F->dead = min(F->dead, i); + for (face *F : conflict[i]) F->dead = min(F->dead, i); int v = -1; for (face *F : conflict[i]) { - if (F->dead != i) - continue; + if (F->dead != i) continue; int parr[3] = {F->a, F->b, F->c}; edge *earr[3] = {F->e1, F->e2, F->e3}; for (int j = 0; j < 3; j++) { - int j2 = (j + 1)%3; + int j2 = (j + 1) % 3; if (earr[j]->f->dead > i) { face *Fn = new_face[parr[j]] = add_face(parr[j], parr[j2], i); set_union(F->points.begin(), F->points.end(), @@ -103,10 +95,7 @@ template struct convex_hull_3d { back_inserter(Fn->points)); Fn->points.erase( stable_partition(Fn->points.begin(), Fn->points.end(), - [&](int k) { - return k > i and - (p[k] - p[Fn->a]).dot(Fn->q) > EPS; - }), + [&](int k) { return k > i && (p[k] - p[Fn->a]).dot(Fn->q) > 0; }), Fn->points.end()); for (int k : Fn->points) conflict[k].push_back(Fn); earr[j]->rev->f = Fn; @@ -122,26 +111,24 @@ template struct convex_hull_3d { v = u; } } - f.erase(remove_if(f.begin(), f.end(), [&](face *F) { return F->dead < n; }), f.end()); + f.erase(remove_if(f.begin(), f.end(), + [&](face *F) { return F->dead < n; }), f.end()); } - int num_faces() const { return f.size(); } - T area() { - T res = 0; - for (face *F : f) - res += F->q.abs() / 2.0; + int num_faces() const { return (int)f.size(); } + double area() const { + double res = 0; + for (face *F : f) res += F->q.abs() / 2.0; return res; } - T volume() { - T vol = 0; - for (face *F : f) - vol += F->pa.dot(F->q); - return abs(vol) / 6.0; + double volume() const { + double vol = 0; + for (face *F : f) vol += p[F->a].dot(F->q); + return std::abs(vol) / 6.0; } - vector> faces() { + vector> faces() const { vector> res; - for (face *F : f) - res.push_back({F->a, F->b, F->c}); + for (face *F : f) res.push_back({F->a, F->b, F->c}); return res; } }; \ No newline at end of file diff --git a/content/geometry/algorithms/delaunay-graph.cpp b/content/geometry/algorithms/delaunay-graph.cpp new file mode 100644 index 0000000..8efa298 --- /dev/null +++ b/content/geometry/algorithms/delaunay-graph.cpp @@ -0,0 +1,39 @@ +/* + *Description:* Planar graph of a Delaunay Triangulation with greedy walk. Walk finds nearest site by the empty-circle property. + *Time:* Build $O(n log n)$ expected | Walk $O(sqrt(n))$ expected + *Status:* Tested + *Dependencies:* point3d, convex_hull_3d, delaunay + */ +template struct delaunay_graph : delaunay { + vector> adj; + + delaunay_graph() {} + delaunay_graph(vector &px, vector &py) { build(px, py); } + + void build(vector &px, vector &py) { + delaunay::build(px, py); + adj.assign(this->n, {}); + for (auto &t : this->tri) + for (int i = 0; i < 3; i++) for (int j = i+1; j < 3; j++) { + adj[t[i]].push_back(t[j]); + adj[t[j]].push_back(t[i]); + } + for (auto &v : adj) { + sort(v.begin(), v.end()); + v.erase(unique(v.begin(), v.end()), v.end()); + } + } + + int walk(vector &px, vector &py, T qx, T qy) const { + auto d2 = [&](int i) { return (qx-px[i])*(qx-px[i]) + (qy-py[i])*(qy-py[i]); }; + int cur = 0; T best = d2(0); + for (bool ok = true; ok; ) { + ok = false; + for (int nb : adj[cur]) { + T d = d2(nb); + if (d < best) { best = d; cur = nb; ok = true; break; } + } + } + return cur; + } +}; \ No newline at end of file diff --git a/content/geometry/algorithms/delaunay.cpp b/content/geometry/algorithms/delaunay.cpp new file mode 100644 index 0000000..383250e --- /dev/null +++ b/content/geometry/algorithms/delaunay.cpp @@ -0,0 +1,37 @@ +/* + *Description:* Delaunay Triangulation via lifting to paraboloid $z = x^2 + y^2$. Exact integer arithmetic. `degenerate=true` if all points are collinear. + *Time:* $O(n log n)$ expected + *Status:* Tested + *Dependencies:* point3d, convex_hull_3d + */ +template struct delaunay : convex_hull_3d { + using pt = point3d; + int n; + vector> tri; + + delaunay() : n(0) {} + delaunay(vector &px, vector &py) : n(0) { build(px, py); } + + void build(vector &px, vector &py) { + n = (int)px.size(); + tri.clear(); + if (n < 3) return; + vector lifted(n); + for (int i = 0; i < n; i++) + lifted[i] = {px[i], py[i], px[i] * px[i] + py[i] * py[i]}; + this->p = move(lifted); + this->build_hull(); + if (this->degenerate) return; + map, int> idx; + for (int i = 0; i < n; i++) idx[{px[i], py[i]}] = i; + int hn = (int)this->p.size(); + vector re(hn); + for (int i = 0; i < hn; i++) + re[i] = idx[{this->p[i].x, this->p[i].y}]; + for (auto &f : this->faces()) { + pt nm = (this->p[f[1]] - this->p[f[0]]) ^ (this->p[f[2]] - this->p[f[0]]); + if (nm.z < 0) + tri.push_back({re[f[0]], re[f[2]], re[f[1]]}); + } + } +}; \ No newline at end of file diff --git a/content/geometry/primitives/point-3d.cpp b/content/geometry/primitives/point-3d.cpp index 9eac166..1e81e11 100644 --- a/content/geometry/primitives/point-3d.cpp +++ b/content/geometry/primitives/point-3d.cpp @@ -51,6 +51,7 @@ template struct point3d { return dot(a ^ b); } + bool nonzero() const { return x != 0 || y != 0 || z != 0; } T norm() const { return dot(*this); } T sq_dist(const point3d &b) const { return (*this - b).norm(); } double abs() const { return sqrt((double)norm()); } diff --git a/tracker-wf.yaml b/tracker-wf.yaml index f7c354d..faec849 100644 --- a/tracker-wf.yaml +++ b/tracker-wf.yaml @@ -157,6 +157,8 @@ content: - polygon-center.cpp - polygon-cut.cpp - hull-diameter.cpp + - delaunay.cpp + - delaunay-graph.cpp - minkowski.cpp - minimum-enclosing-sphere.cpp - strings: diff --git a/tracker.yaml b/tracker.yaml index 513385a..ddda096 100644 --- a/tracker.yaml +++ b/tracker.yaml @@ -192,6 +192,8 @@ content: - polygon-center.cpp - polygon-cut.cpp - hull-diameter.cpp + - delaunay.cpp + - delaunay-graph.cpp - dynamic-programming: - classic: - lis.cpp From f730c9bbba899480777fef4e96a689bef0f93a76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebasti=C3=A1n=20Torrealba?= Date: Tue, 3 Mar 2026 18:09:38 -0300 Subject: [PATCH 3/3] Replace Delaunay lifting (chull3d) with divide-and-conquer quad-edge implementation, add locate to delaunay_graph --- .../geometry/algorithms/delaunay-graph.cpp | 38 +++--- content/geometry/algorithms/delaunay.cpp | 126 ++++++++++++++---- 2 files changed, 117 insertions(+), 47 deletions(-) diff --git a/content/geometry/algorithms/delaunay-graph.cpp b/content/geometry/algorithms/delaunay-graph.cpp index 8efa298..c03e4d8 100644 --- a/content/geometry/algorithms/delaunay-graph.cpp +++ b/content/geometry/algorithms/delaunay-graph.cpp @@ -1,39 +1,41 @@ /* - *Description:* Planar graph of a Delaunay Triangulation with greedy walk. Walk finds nearest site by the empty-circle property. - *Time:* Build $O(n log n)$ expected | Walk $O(sqrt(n))$ expected + *Description:* Planar graph of a Delaunay Triangulation with greedy walk and locate. `walk` finds nearest site. `locate` finds all equidistant nearest sites via BFS: 1 site $=$ Voronoi region, 2 $=$ edge, 3+ $=$ vertex (cocircular). + *Time:* Build $O(n log n)$ | Walk $O(sqrt(n))$ expected | Locate $O(sqrt(n))$ expected *Status:* Tested - *Dependencies:* point3d, convex_hull_3d, delaunay + *Dependencies:* point2d, delaunay */ template struct delaunay_graph : delaunay { vector> adj; - delaunay_graph() {} delaunay_graph(vector &px, vector &py) { build(px, py); } - void build(vector &px, vector &py) { delaunay::build(px, py); adj.assign(this->n, {}); for (auto &t : this->tri) for (int i = 0; i < 3; i++) for (int j = i+1; j < 3; j++) { - adj[t[i]].push_back(t[j]); - adj[t[j]].push_back(t[i]); + adj[t[i]].push_back(t[j]); adj[t[j]].push_back(t[i]); } - for (auto &v : adj) { - sort(v.begin(), v.end()); - v.erase(unique(v.begin(), v.end()), v.end()); - } + for (auto &v : adj) { sort(v.begin(), v.end()); v.erase(unique(v.begin(), v.end()), v.end()); } } - int walk(vector &px, vector &py, T qx, T qy) const { auto d2 = [&](int i) { return (qx-px[i])*(qx-px[i]) + (qy-py[i])*(qy-py[i]); }; int cur = 0; T best = d2(0); - for (bool ok = true; ok; ) { - ok = false; - for (int nb : adj[cur]) { - T d = d2(nb); - if (d < best) { best = d; cur = nb; ok = true; break; } - } + for (bool ok = true; ok; ) { ok = false; + for (int nb : adj[cur]) { T d = d2(nb); if (d < best) { best=d; cur=nb; ok=true; break; } } } return cur; } + vector locate(vector &px, vector &py, T qx, T qy) const { + auto d2 = [&](int i) { return (qx-px[i])*(qx-px[i]) + (qy-py[i])*(qy-py[i]); }; + int v = walk(px, py, qx, qy); + T best = d2(v); + vector res = {v}; + vector vis(this->n, false); + vis[v] = true; + for (int qi = 0; qi < (int)res.size(); qi++) + for (int nb : adj[res[qi]]) + if (!vis[nb]) { vis[nb] = true; if (d2(nb) == best) res.push_back(nb); } + sort(res.begin(), res.end()); + return res; + } }; \ No newline at end of file diff --git a/content/geometry/algorithms/delaunay.cpp b/content/geometry/algorithms/delaunay.cpp index 383250e..866ff50 100644 --- a/content/geometry/algorithms/delaunay.cpp +++ b/content/geometry/algorithms/delaunay.cpp @@ -1,37 +1,105 @@ /* - *Description:* Delaunay Triangulation via lifting to paraboloid $z = x^2 + y^2$. Exact integer arithmetic. `degenerate=true` if all points are collinear. - *Time:* $O(n log n)$ expected + *Description:* Delaunay Triangulation using divide and conquer with quad-edges. Each circumcircle contains none of the input points. If all points are collinear, `tri` is empty. Uses `__int128` for circumcircle test (coords up to $10^9$). + *Time:* $O(n log n)$ *Status:* Tested - *Dependencies:* point3d, convex_hull_3d + *Dependencies:* point2d */ -template struct delaunay : convex_hull_3d { - using pt = point3d; +template struct delaunay { + using pt = point2d; + using lll = __int128_t; + struct Quad { + Quad *rot, *o; + pt p; int id; bool mark; + pt& F() { return rot->rot->p; } + int& Fi() { return rot->rot->id; } + Quad*& r() { return rot->rot; } + Quad* prev() { return rot->o->rot; } + Quad* next() { return r()->prev(); } + }; + Quad* pool; + bool circ(pt p, pt a, pt b, pt c) { + lll p2=p.norm(), A=a.norm()-p2, B=b.norm()-p2, C=c.norm()-p2; + return p.cross(a,b)*C + p.cross(b,c)*A + p.cross(c,a)*B > 0; + } + Quad* makeEdge(pt orig, int oi, pt dest, int di) { + auto nq = [&]() { return new Quad{nullptr,nullptr,{}, -1, false}; }; + Quad* r; + if (pool) { r = pool; pool = r->o; } + else { r = nq(); r->rot = nq(); r->rot->rot = nq(); r->rot->rot->rot = nq(); } + r->r()->r() = r; r->rot->rot->rot->rot = r; + pt arb{(T)1e18, (T)1e18}; + Quad* c = r; + for (int i = 0; i < 4; i++) { + c = c->rot; c->p = arb; c->id = -1; + c->o = (i & 1) ? c : c->r(); + } + r->p = orig; r->id = oi; r->F() = dest; r->Fi() = di; + return r; + } + void splice(Quad* a, Quad* b) { + swap(a->o->rot->o, b->o->rot->o); swap(a->o, b->o); + } + Quad* connect(Quad* a, Quad* b) { + Quad* q = makeEdge(a->F(), a->Fi(), b->p, b->id); + splice(q, a->next()); splice(q->r(), b); return q; + } + pair rec(vector& s, vector& ids, int lo, int hi) { + int sz = hi - lo; + if (sz <= 3) { + Quad* a = makeEdge(s[lo],ids[lo], s[lo+1],ids[lo+1]); + Quad* b = makeEdge(s[lo+1],ids[lo+1], s[lo+sz-1],ids[lo+sz-1]); + if (sz == 2) return {a, a->r()}; + splice(a->r(), b); + auto side = s[lo].cross(s[lo+1], s[lo+2]); + Quad* c = side ? connect(b, a) : nullptr; + return {side<0 ? c->r() : a, side<0 ? c : b->r()}; + } + int mid = lo + sz/2; + auto [ra, A] = rec(s, ids, lo, mid); + auto [B, rb] = rec(s, ids, mid, hi); + while ((B->p.cross(A->F(), A->p)<0 && (A=A->next())) || + (A->p.cross(B->F(), B->p)>0 && (B=B->r()->o))); + Quad* base = connect(B->r(), A); + if (A->p.x==ra->p.x && A->p.y==ra->p.y) ra = base->r(); + if (B->p.x==rb->p.x && B->p.y==rb->p.y) rb = base; + auto valid = [&](Quad* e) { return e->F().cross(base->F(), base->p) > 0; }; + for (;;) { + Quad* LC = base->r()->o; + if (valid(LC)) while (circ(LC->o->F(), base->F(), base->p, LC->F())) + { Quad* t=LC->o; splice(LC,LC->prev()); splice(LC->r(),LC->r()->prev()); LC->o=pool; pool=LC; LC=t; } + Quad* RC = base->prev(); + if (valid(RC)) while (circ(RC->prev()->F(), base->F(), base->p, RC->F())) + { Quad* t=RC->prev(); splice(RC,RC->prev()); splice(RC->r(),RC->r()->prev()); RC->o=pool; pool=RC; RC=t; } + if (!valid(LC) && !valid(RC)) break; + if (!valid(LC) || (valid(RC) && circ(RC->F(), base->F(), base->p, LC->F()))) + base = connect(RC, base->r()); + else base = connect(base->r(), LC->r()); + } + return {ra, rb}; + } int n; vector> tri; - - delaunay() : n(0) {} - delaunay(vector &px, vector &py) : n(0) { build(px, py); } - + delaunay() : pool(nullptr), n(0) {} + delaunay(vector &px, vector &py) : pool(nullptr), n(0) { build(px, py); } void build(vector &px, vector &py) { - n = (int)px.size(); - tri.clear(); - if (n < 3) return; - vector lifted(n); - for (int i = 0; i < n; i++) - lifted[i] = {px[i], py[i], px[i] * px[i] + py[i] * py[i]}; - this->p = move(lifted); - this->build_hull(); - if (this->degenerate) return; - map, int> idx; - for (int i = 0; i < n; i++) idx[{px[i], py[i]}] = i; - int hn = (int)this->p.size(); - vector re(hn); - for (int i = 0; i < hn; i++) - re[i] = idx[{this->p[i].x, this->p[i].y}]; - for (auto &f : this->faces()) { - pt nm = (this->p[f[1]] - this->p[f[0]]) ^ (this->p[f[2]] - this->p[f[0]]); - if (nm.z < 0) - tri.push_back({re[f[0]], re[f[2]], re[f[1]]}); - } + n = (int)px.size(); tri.clear(); pool = nullptr; + if (n < 2) return; + vector ord(n); iota(ord.begin(), ord.end(), 0); + sort(ord.begin(), ord.end(), [&](int a, int b) { + return px[a] sp(n); vector si(n); + for (int i = 0; i < n; i++) { sp[i] = {px[ord[i]], py[ord[i]]}; si[i] = ord[i]; } + Quad* e = rec(sp, si, 0, n).first; + vector q = {e}; + int qi = 0; + while (e->o->F().cross(e->F(), e->p) < 0) e = e->o; + auto add = [&](Quad* e) { + Quad* c = e; vector face; + do { c->mark=1; face.push_back(c->id); q.push_back(c->r()); c=c->next(); } while (c!=e); + if ((int)face.size()==3) tri.push_back({face[0],face[1],face[2]}); + }; + add(e); + while (qi < (int)q.size()) if (!(e=q[qi++])->mark) add(e); } }; \ No newline at end of file