diff --git a/geo/euclidean_mst/checker.cpp b/geo/euclidean_mst/checker.cpp
new file mode 100644
index 00000000..02199846
--- /dev/null
+++ b/geo/euclidean_mst/checker.cpp
@@ -0,0 +1,91 @@
+#include <algorithm>
+#include <vector>
+#include <utility>
+#include "testlib.h"
+
+using namespace std;
+using ll = long long;
+
+
+struct DsuFast{
+private:
+    std::vector<int> w;
+public:
+    DsuFast(int n = 0) : w(n, -1) {}
+    int leader(int u){
+        if(w[u] < 0) return u;
+        return w[u] = leader(w[u]);
+    }
+    int operator[](int u){ return leader(u); }
+    int merge(int u, int v){
+        u = leader(u);
+        v = leader(v);
+        if(u == v) return u;
+        if(-w[u] < -w[v]) std::swap(u, v);
+        w[u] += w[v];
+        w[v] = u;
+        return u;
+    }
+    int size(int u){ return -w[leader(u)]; }
+    bool same(int u, int v){ return leader(u) == leader(v); }
+};
+
+int main(int argc, char* argv[]) {
+    registerTestlibCmd(argc, argv);
+
+
+    int N = inf.readInt();
+
+    vector<ll> X(N), Y(N);
+    for(int i=0; i<N; i++){
+        X[i] = inf.readInt();
+        Y[i] = inf.readInt();
+    }
+
+    auto weight = [&](int a, int b){
+        ll dx = X[a] - X[b];
+        ll dy = Y[a] - Y[b];
+        return dx * dx + dy * dy;
+    };
+
+    DsuFast dsu_ans(N);
+    vector<ll> weights_ans;
+    for(int i=0; i<N-1; i++){
+        int a = ans.readInt(0, N - 1);
+        int b = ans.readInt(0, N - 1);
+        dsu_ans.merge(a, b);
+        weights_ans.push_back(weight(a, b));
+    }
+
+    if(dsu_ans.size(0) != N){
+        quit(_fail, "the output of writer\'s solution is not connected");
+    }
+
+    sort(weights_ans.begin(), weights_ans.end());
+
+    DsuFast dsu_out(N);
+    vector<ll> weights_out;
+    for(int i=0; i<N-1; i++){
+        int a = ouf.readInt(0, N - 1);
+        int b = ouf.readInt(0, N - 1);
+        dsu_out.merge(a, b);
+        weights_out.push_back(weight(a, b));
+    }
+
+    if(dsu_out.size(0) != N){
+        quit(_wa, "not connected");
+    }
+
+    sort(weights_out.begin(), weights_out.end());
+
+    for(int i=0; i<N-1; i++){
+        if(weights_ans[i] > weights_out[i]){
+            quit(_fail, "your output showed the writer solution doesn\'t output the minumum spanning tree");
+        }
+        if(weights_ans[i] < weights_out[i]){
+            quit(_wa, "not the minumum cost");
+        }
+    }
+    
+    quitf(_ok, "OK");
+}
diff --git a/geo/euclidean_mst/gen/all_same.cpp b/geo/euclidean_mst/gen/all_same.cpp
new file mode 100644
index 00000000..35d84943
--- /dev/null
+++ b/geo/euclidean_mst/gen/all_same.cpp
@@ -0,0 +1,33 @@
+#include <cstdio>
+#include <cstdlib>
+#include <utility>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using P = pair<int, int>;
+
+P random_point(Random& gen, int LIM) {
+    int x = gen.uniform(-LIM, LIM);
+    int y = gen.uniform(-LIM, LIM);
+    return {x, y};
+}
+
+void out(vector<P> S) {
+    int n = S.size();
+    printf("%d\n", n);
+    for (auto& [x, y]: S) { printf("%d %d\n", x, y); }
+}
+
+int main(int, char* argv[]) {
+    long long seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    int n = N_MAX;
+    auto [x, y] = random_point(gen, X_AND_Y_ABS_MAX);
+    vector<P> S(n, { x,y });
+
+    out(S);
+    return 0;
+}
diff --git a/geo/euclidean_mst/gen/example_00.in b/geo/euclidean_mst/gen/example_00.in
new file mode 100644
index 00000000..bc09ec33
--- /dev/null
+++ b/geo/euclidean_mst/gen/example_00.in
@@ -0,0 +1,6 @@
+5
+-1 -1
+-6 4
+-9 -7
+2 5
+-7 6
diff --git a/geo/euclidean_mst/gen/max_colinear.cpp b/geo/euclidean_mst/gen/max_colinear.cpp
new file mode 100644
index 00000000..3436a13e
--- /dev/null
+++ b/geo/euclidean_mst/gen/max_colinear.cpp
@@ -0,0 +1,35 @@
+#include <cstdio>
+#include <cstdlib>
+#include <utility>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using P = pair<int, int>;
+
+void out(vector<P> S) {
+    int n = S.size();
+    printf("%d\n", n);
+    for (auto& [x, y]: S) { printf("%d %d\n", x, y); }
+}
+
+int main(int, char* argv[]) {
+    long long seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    int x0 = gen.uniform<int>(1, 10);
+    int y0 = gen.uniform<int>(1, 10);
+    int lim = X_AND_Y_ABS_MAX / max(x0, y0);
+
+    int n = N_MAX;
+    vector<P> S;
+    for (int i = 0; i < n; ++i) {
+        int t = gen.uniform<int>(-lim, lim);
+        int x = x0 * t, y = y0 * t;
+        S.emplace_back(x, y);
+    }
+
+    out(S);
+    return 0;
+}
diff --git a/geo/euclidean_mst/gen/max_random.cpp b/geo/euclidean_mst/gen/max_random.cpp
new file mode 100644
index 00000000..8965b752
--- /dev/null
+++ b/geo/euclidean_mst/gen/max_random.cpp
@@ -0,0 +1,37 @@
+#include <cstdio>
+#include <cstdlib>
+#include <utility>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using P = pair<int, int>;
+
+P random_point(Random& gen, int LIM) {
+    int x = gen.uniform(-LIM, LIM);
+    int y = gen.uniform(-LIM, LIM);
+    return {x, y};
+}
+
+void out(vector<P> S) {
+    int n = S.size();
+    printf("%d\n", n);
+    for (auto& [x, y]: S) { printf("%d %d\n", x, y); }
+}
+
+int main(int, char* argv[]) {
+    long long seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    int lims[] = { 10, 1000, X_AND_Y_ABS_MAX, X_AND_Y_ABS_MAX };
+
+    int LIM = lims[seed % 4];
+
+    int n = N_MAX;
+    vector<P> S;
+    for (int i = 0; i < n; ++i) S.emplace_back(random_point(gen, LIM));
+
+    out(S);
+    return 0;
+}
diff --git a/geo/euclidean_mst/gen/max_slender.cpp b/geo/euclidean_mst/gen/max_slender.cpp
new file mode 100644
index 00000000..85c342c0
--- /dev/null
+++ b/geo/euclidean_mst/gen/max_slender.cpp
@@ -0,0 +1,52 @@
+#include <cstdio>
+#include <cstdlib>
+#include <utility>
+#include <algorithm>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using P = pair<int, int>;
+
+void out(vector<P> S) {
+    int n = S.size();
+    printf("%d\n", n);
+    for (auto& [x, y]: S) { printf("%d %d\n", x, y); }
+}
+
+int main(int, char* argv[]) {
+    long long seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    int n = N_MAX;
+    int h = 100;
+    int w = X_AND_Y_ABS_MAX - h;
+
+    vector<P> S;
+    for(int i=0; i<n; i++){
+        int x = gen.uniform<int>(-w, w);
+        int y = gen.uniform<int>(-h, h);
+        S.push_back(make_pair(x, y));
+    }
+
+    if(seed == 1){
+        for(auto& a : S){
+            swap(a.first, a.second);
+        }
+    }
+    else if(seed == 2){
+        for(auto& a : S){
+            a.second += a.first;
+        }
+    }
+    else if(seed == 3){
+        for(auto& a : S){
+            a.second -= a.first;
+        }
+    }
+
+
+    out(S);
+    return 0;
+}
diff --git a/geo/euclidean_mst/gen/mid_random.cpp b/geo/euclidean_mst/gen/mid_random.cpp
new file mode 100644
index 00000000..04bc1573
--- /dev/null
+++ b/geo/euclidean_mst/gen/mid_random.cpp
@@ -0,0 +1,37 @@
+#include <cstdio>
+#include <cstdlib>
+#include <utility>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using P = pair<int, int>;
+
+P random_point(Random& gen, int LIM) {
+    int x = gen.uniform(-LIM, LIM);
+    int y = gen.uniform(-LIM, LIM);
+    return {x, y};
+}
+
+void out(vector<P> S) {
+    int n = S.size();
+    printf("%d\n", n);
+    for (auto& [x, y]: S) { printf("%d %d\n", x, y); }
+}
+
+int main(int, char* argv[]) {
+    long long seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    int lims[] = { 10, 1000, X_AND_Y_ABS_MAX, X_AND_Y_ABS_MAX };
+
+    int LIM = lims[seed % 4];
+
+    int n = gen.uniform<int>(100, 1000);
+    vector<P> S;
+    for (int i = 0; i < n; ++i) S.emplace_back(random_point(gen, LIM));
+
+    out(S);
+    return 0;
+}
diff --git a/geo/euclidean_mst/gen/min_00.in b/geo/euclidean_mst/gen/min_00.in
new file mode 100644
index 00000000..c3db04f7
--- /dev/null
+++ b/geo/euclidean_mst/gen/min_00.in
@@ -0,0 +1,2 @@
+1
+0 0
diff --git a/geo/euclidean_mst/gen/min_01.in b/geo/euclidean_mst/gen/min_01.in
new file mode 100644
index 00000000..87b4f446
--- /dev/null
+++ b/geo/euclidean_mst/gen/min_01.in
@@ -0,0 +1,3 @@
+2
+0 0
+0 1
diff --git a/geo/euclidean_mst/gen/min_02.in b/geo/euclidean_mst/gen/min_02.in
new file mode 100644
index 00000000..256e46ca
--- /dev/null
+++ b/geo/euclidean_mst/gen/min_02.in
@@ -0,0 +1,3 @@
+2
+0 0
+0 0
diff --git a/geo/euclidean_mst/gen/near_circle.cpp b/geo/euclidean_mst/gen/near_circle.cpp
new file mode 100644
index 00000000..984e09b7
--- /dev/null
+++ b/geo/euclidean_mst/gen/near_circle.cpp
@@ -0,0 +1,41 @@
+#include <cstdio>
+#include <cstdlib>
+#include <utility>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using P = pair<int, int>;
+
+P random_point(Random& gen) {
+    long long R = X_AND_Y_ABS_MAX;
+    long long x = gen.uniform(-R, R);
+    long long ok = 0, ng = R + 1;
+    while (ok + 1 < ng) {
+        long long mi = (ok + ng) / 2;
+        (x * x + mi * mi <= R * R ? ok : ng) = mi;
+    }
+    long long y = ok;
+    if (gen.uniform<int>(0, 1)) y = -y;
+    if (gen.uniform<int>(0, 1)) swap(x, y);
+    return {x, y};
+}
+
+void out(vector<P> S) {
+    int n = S.size();
+    printf("%d\n", n);
+    for (auto& [x, y]: S) { printf("%d %d\n", x, y); }
+}
+
+int main(int, char* argv[]) {
+    long long seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    int n = N_MAX;
+    vector<P> S;
+    for (int i = 0; i < n; ++i) S.emplace_back(random_point(gen));
+
+    out(S);
+    return 0;
+}
diff --git a/geo/euclidean_mst/gen/near_grid.cpp b/geo/euclidean_mst/gen/near_grid.cpp
new file mode 100644
index 00000000..d79efe76
--- /dev/null
+++ b/geo/euclidean_mst/gen/near_grid.cpp
@@ -0,0 +1,49 @@
+#include <cstdio>
+#include <cstdlib>
+#include <utility>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using P = pair<int, int>;
+
+void out(vector<P> S) {
+    int n = S.size();
+    printf("%d\n", n);
+    for (auto& [x, y]: S) { printf("%d %d\n", x, y); }
+}
+
+int main(int, char* argv[]) {
+    long long seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    vector<vector<P>> ALL;
+    int n = 1;
+    while ((n + 1) * (n + 1) <= N_MAX) ++n;
+
+    int eps = 5;
+    /*
+    x[0],x[1],...,x[N-1]
+    */
+    int LIM = X_AND_Y_ABS_MAX;
+    int d = LIM / (n - 1);
+
+    auto get = [&](int i) -> int {
+        int x = i * d - (LIM / 2);
+        return x + gen.uniform<int>(-eps, eps);
+    };
+
+    vector<P> S;
+    for (int i = 0; i < n; ++i) {
+        for (int j = 0; j < n; ++j) {
+        int x = get(i), y = get(j);
+        S.emplace_back(x, y);
+        }
+    }
+
+    gen.shuffle(S.begin(), S.end());
+
+    out(S);
+    return 0;
+}
diff --git a/geo/euclidean_mst/gen/small_chunks.cpp b/geo/euclidean_mst/gen/small_chunks.cpp
new file mode 100644
index 00000000..7f5ad6e1
--- /dev/null
+++ b/geo/euclidean_mst/gen/small_chunks.cpp
@@ -0,0 +1,55 @@
+#include <cstdio>
+#include <cstdlib>
+#include <utility>
+#include "random.h"
+#include "../params.h"
+
+using namespace std;
+
+using P = pair<int, int>;
+
+void out(vector<P> S) {
+    int n = S.size();
+    printf("%d\n", n);
+    for (auto& [x, y]: S) { printf("%d %d\n", x, y); }
+}
+
+int sqrt_int_ceil(int nn){
+    int n = 0;
+    while((n+1) * (n+1) <= nn) n++;
+    return n;
+}
+
+int main(int, char* argv[]) {
+    long long seed = atoll(argv[1]);
+    auto gen = Random(seed);
+
+    int chunk_width = 10;
+    int n = N_MAX;
+    int chunk_count = n / 20;
+
+    vector<int> chunk_n(chunk_count);
+    for(int i=0; i<n; i++) chunk_n[gen.uniform<int>(0, chunk_count-1)]++;
+
+    const int global_coord_min = int(-X_AND_Y_ABS_MAX);
+    const int global_coord_width = int(X_AND_Y_ABS_MAX) * 2 + 1;
+    int chunk_count_side = sqrt_int_ceil(chunk_count);
+    int chunk_interval = (global_coord_width - chunk_width) / (chunk_count_side - 1);
+
+    vector<P> S;
+    
+    for(int chunk_i = 0; chunk_i < chunk_count; chunk_i++){
+        int offset_x = global_coord_min + chunk_i / chunk_count_side * chunk_interval;
+        int offset_y = global_coord_min + chunk_i % chunk_count_side * chunk_interval;
+        for(int i = 0; i < chunk_n[i]; i++){
+            int dx = gen.uniform<int>(0, chunk_width - 1);
+            int dy = gen.uniform<int>(0, chunk_width - 1);
+            S.push_back(make_pair(offset_x + dx, offset_y + dy));
+        }
+    }
+
+    gen.shuffle(S.begin(), S.end());
+
+    out(S);
+    return 0;
+}
diff --git a/geo/euclidean_mst/hash.json b/geo/euclidean_mst/hash.json
new file mode 100644
index 00000000..b44e95ff
--- /dev/null
+++ b/geo/euclidean_mst/hash.json
@@ -0,0 +1,60 @@
+{
+  "all_same_00.in": "15819da0571c986fb677dcd0a3b6eb03c84327250434defe2f100f05956766ca",
+  "all_same_00.out": "a9d04e964c2573e049281626f8b9c7abb7360dfb49f258348eb01b26f61711e1",
+  "example_00.in": "c2a911d5e776235fc515736900a5206e4ebd9329735f32127e4b68b305a356c8",
+  "example_00.out": "1bbfda1800d0a8236bab40348f624252ae0cfbf4166bdb2fb8b773554728bc28",
+  "max_colinear_00.in": "eb16e81d0788a2d332f8f710435102553f5e4edb0299cf91e38053e458980367",
+  "max_colinear_00.out": "33f0a34b994284f7ef487895201d058fb6f83c776af6ff0acd19e126838f8a38",
+  "max_colinear_01.in": "d26d19ca74a216567a80d7067467d44dedeacf82da0a4827458fbdb889ecaacf",
+  "max_colinear_01.out": "38a02c3a278b7a9ceafbf37db829b7271ef80d149499f751360dd8c74557d73a",
+  "max_random_00.in": "0888114ee54313b2eea183c0cd91596d5779d3d279cb3fca2459e8d34c21251f",
+  "max_random_00.out": "df7e9bc3c8d2f708a9fc8518e3496abb1648d70518ac833c0ea63e8997d855b6",
+  "max_random_01.in": "ae5f9bf50c5d5bb33f5cfe0cc8802ec63012cad9671f0eaca61a6951ff4c42cc",
+  "max_random_01.out": "1f81bb598e9816f210bf5241d9b5429fbde673f946ba734036975d085dacab9d",
+  "max_random_02.in": "0d70b4b5fd328292a700d32ea7cbe45be5d97f1b0092bec7ae6f91d419ea62e1",
+  "max_random_02.out": "68db4f26ddf31f28d00c47ee85ea50ca4054f824388611cc3a2c6af924c18aa5",
+  "max_random_03.in": "1caac0b4abf5cfab6b4e9f1abfbef1bff051a45962dc8c3233b87f8d7465083f",
+  "max_random_03.out": "c7e444586b85aee27e549216c71820d2cda02a1cf5265dfe30db42ee8161e751",
+  "max_random_04.in": "62832801fdb905fe2147f2f34c30c00a83cc8b37e23ff38a80508de747ee93d0",
+  "max_random_04.out": "d28588687324375aa2a8ce2ce43d25f4e4b805accdf0d1d0c587213752f10dbb",
+  "max_random_05.in": "e6edc485c8e59a9af54ad4cbf8e1b6b42379d165d25b50f12991d06c64ffe1dd",
+  "max_random_05.out": "7fd8b2440931b495060cb244b153030499ef90240af4ccedd36e0aa9c383ba94",
+  "max_random_06.in": "c070996f6f925d29a10a64af4c17e9be585d1512299ca4581e0810d8a3d9fb2d",
+  "max_random_06.out": "3187a4725eb03962e84336e298c3e1a9b5a9f0f5544c767c8716aee75d2de0d3",
+  "max_random_07.in": "69fed19041b7b6b1dafdbb695543603f9fea3685efbf1b64eee5103e6d61bd2d",
+  "max_random_07.out": "05171527b47d839789b2e1ac0768995c94e20337abaca101a5042575880513c4",
+  "max_slender_00.in": "e6cfabd781ed24d3c82f08ae6e5e522f77cfc8b9d01394a0e7cda7f0ef514297",
+  "max_slender_00.out": "7d0cf4e83b9b4f73167991312f69d9c45a4818acf9e27d454c2e58dd0eba1bd7",
+  "max_slender_01.in": "eedcae1a2763c588c4e5159f3447d0fe2d37630871d511c72b325e377f5f63dc",
+  "max_slender_01.out": "60738f2976d60d316f83195fa5ba8317f2c13ee2ad9bc96546e1ac351be9981f",
+  "max_slender_02.in": "2d5b7a68c27ed517a400537a18af0e6156c38f0bdc6e80a400af82f97d0bc3cc",
+  "max_slender_02.out": "fa7d86fcb017304425de079ffdae31e6fc56a40a0e49a345c5feb2ab4bb0702e",
+  "max_slender_03.in": "6614626ef89fa96d2b95d354ec0625b3db0f2fb7ad7dc5599cdceca0740d98be",
+  "max_slender_03.out": "b801f7af82ac164557130cbe69ebc69aec73e6487e7aede5e629ae89b3ca6720",
+  "mid_random_00.in": "396c8b427a7a2ba551c4d1762f6ff44d55f34d66cc8ff572bcad6d678600c106",
+  "mid_random_00.out": "a588f51abf00bc4caf0fd8ed033ba47cdf96c8720b7c39c9418ef4c2c985a500",
+  "mid_random_01.in": "322befea098b561c5ca8fbf19f6f61d82263acbb9b288a64d2e3e9cf3aaa9f5e",
+  "mid_random_01.out": "b773808fb1b445e2a1801ead038fccb408ff63d0777e024b016dd106633d445a",
+  "mid_random_02.in": "961c99e231a293a574a1be2504ea7d3696f5f31b10cd836939df063cf88534bd",
+  "mid_random_02.out": "4750e5179f4178d158da637107d3cf6496bb3bb70bcf1cd56a94aa5be7c54b08",
+  "mid_random_03.in": "9fa4a55f12ddac77310ad8b3fe8583f556daba516d347f30bad1efbbdbb4b290",
+  "mid_random_03.out": "9ef502b043dbc01f1efd7b4d1ab76374f4bfb7f69e85dacd7611793ccdde214c",
+  "min_00.in": "cfabffaab1ba46d5f6db92ce5fa5c78fc6446f39cb4d6c627d630715ac48f92c",
+  "min_00.out": "e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855",
+  "min_01.in": "55bc6ea841574879ce514f4e6c9761c42c1158c536a10c6cbb625d85d2de1afc",
+  "min_01.out": "a79122992d53d358e6bbbbb98883d64fa0c15df3bcb08ff7b65a0580870af424",
+  "min_02.in": "a5a9ddc641534c3d4bd22d72fe09c3a0a84f2a6e2fe0b464a94526011a6648c3",
+  "min_02.out": "a79122992d53d358e6bbbbb98883d64fa0c15df3bcb08ff7b65a0580870af424",
+  "near_circle_00.in": "c685c137c4f30dfabe14db630457d7aeabf92bb0e01a0cefe6d7039c3a6bf498",
+  "near_circle_00.out": "6b0513cbfeba19ea3c687f31134c3b983a3785bc25df3d5be7b244c3fab47c15",
+  "near_circle_01.in": "c2e31ed77955da1f7265c41aa8c6c33c993abcf7ecf670a7bbf71892ad583bc1",
+  "near_circle_01.out": "6aa2e523cbde3c265615458f8d85f6a18f52e74aff19f900bf1960f9b547613e",
+  "near_grid_00.in": "17d7049955cb3aef4618656cf37ec5d5601e61f43c68f18dcf720c1722747517",
+  "near_grid_00.out": "57841d6dcad27be76e5a8bb2389c9e6c1912a4e1e3f83412c06fb30d05599290",
+  "near_grid_01.in": "0f916038d08619bd159d197aba42cb698ed02dd2b14d2ab48665675e71b5974a",
+  "near_grid_01.out": "ecceafe41ddcb140b132ce83a4a398e03aee189a1046e188aaaceeaf18534263",
+  "small_chunks_00.in": "23b5789c8ed8f433cd8114c831a7e4a3ac1c9b1fabb0dfaf4de2b1a3a89cd052",
+  "small_chunks_00.out": "eb01548baa5abc1c99597298dbb3d055f7b71ffa9fc462bcd1381f5624093a35",
+  "small_chunks_01.in": "aafd58498ce03e6d31d89a41c7b93ad33af581c4e55d002c877d8854c7a90916",
+  "small_chunks_01.out": "a1792895c9238c621d19b4e28794eacdee36fbca9e319d3bbddab610c35ca716"
+}
\ No newline at end of file
diff --git a/geo/euclidean_mst/info.toml b/geo/euclidean_mst/info.toml
new file mode 100644
index 00000000..90264dbc
--- /dev/null
+++ b/geo/euclidean_mst/info.toml
@@ -0,0 +1,51 @@
+title = 'Euclidean MST'
+timelimit = 5.0
+forum = "https://github.com/yosupo06/library-checker-problems/issues/916"
+
+[[tests]]
+    name = "example.in"
+    number = 1
+    
+[[tests]]
+    name = "all_same.cpp"
+    number = 1
+
+[[tests]]
+    name = "mid_random.cpp"
+    number = 4
+
+[[tests]]
+    name = "min.in"
+    number = 3
+
+[[tests]]
+    name = "max_random.cpp"
+    number = 8
+
+[[tests]]
+    name = "max_colinear.cpp"
+    number = 2
+
+[[tests]]
+    name = "max_slender.cpp"
+    number = 4
+
+[[tests]]
+    name = "near_circle.cpp"
+    number = 2
+
+[[tests]]
+    name = "near_grid.cpp"
+    number = 2
+
+[[tests]]
+    name = "small_chunks.cpp"
+    number = 2
+
+[[solutions]]
+    name = "naive_prim.cpp"
+    expect = "TLE"
+
+[params]
+    N_MAX = 200_000
+    X_AND_Y_ABS_MAX = 10000
diff --git a/geo/euclidean_mst/sol/correct.cpp b/geo/euclidean_mst/sol/correct.cpp
new file mode 100644
index 00000000..3eedd914
--- /dev/null
+++ b/geo/euclidean_mst/sol/correct.cpp
@@ -0,0 +1,452 @@
+#include <utility>
+#include <vector>
+#include <tuple>
+#include <algorithm>
+
+template<class Int = long long, class Int2 = long long>
+struct VecI2 {
+    Int x, y;
+    VecI2() : x(0), y(0) {}
+    VecI2(std::pair<Int, Int> _p) : x(std::move(_p.first)), y(std::move(_p.second)) {}
+    VecI2(Int _x, Int _y) : x(std::move(_x)), y(std::move(_y)) {}
+    VecI2& operator+=(VecI2 r){ x+=r.x; y+=r.y; return *this; }
+    VecI2& operator-=(VecI2 r){ x-=r.x; y-=r.y; return *this; }
+    VecI2& operator*=(Int r){ x*=r; y*=r; return *this; }
+    VecI2 operator+(VecI2 r) const { return VecI2(x+r.x, y+r.y); }
+    VecI2 operator-(VecI2 r) const { return VecI2(x-r.x, y-r.y); }
+    VecI2 operator*(Int r) const { return VecI2(x*r, y*r); }
+    VecI2 operator-() const { return VecI2(-x, -y); }
+    Int2 operator*(VecI2 r) const { return Int2(x) * Int2(r.x) + Int2(y) * Int2(r.y); }
+    Int2 operator^(VecI2 r) const { return Int2(x) * Int2(r.y) - Int2(y) * Int2(r.x); }
+    bool operator<(VecI2 r) const { return x < r.x || (!(r.x < x) && y < r.y); }
+    Int2 norm() const { return Int2(x) * Int2(x) + Int2(y) * Int2(y); }
+    static bool compareYX(VecI2 a, VecI2 b){ return a.y < b.y || (!(b.y < a.y) && a.x < b.x); }
+    static bool compareXY(VecI2 a, VecI2 b){ return a.x < b.x || (!(b.x < a.x) && a.y < b.y); }
+    bool operator==(VecI2 r) const { return x == r.x && y == r.y; }
+    bool operator!=(VecI2 r) const { return x != r.x || y != r.y; }
+};
+
+template<class Elem>
+class CsrArray{
+public:
+    struct ListRange{
+        using iterator = typename std::vector<Elem>::iterator;
+        iterator begi, endi;
+        iterator begin() const { return begi; }
+        iterator end() const { return endi; }
+        int size() const { return (int)std::distance(begi, endi); }
+        Elem& operator[](int i) const { return begi[i]; }
+    };
+    struct ConstListRange{
+        using iterator = typename std::vector<Elem>::const_iterator;
+        iterator begi, endi;
+        iterator begin() const { return begi; }
+        iterator end() const { return endi; }
+        int size() const { return (int)std::distance(begi, endi); }
+        const Elem& operator[](int i) const { return begi[i]; }
+    };
+private:
+    int m_n;
+    std::vector<Elem> m_list;
+    std::vector<int> m_pos;
+public:
+    CsrArray() : m_n(0), m_list(), m_pos() {}
+    static CsrArray Construct(int n, std::vector<std::pair<int, Elem>> items){
+        CsrArray res;
+        res.m_n = n;
+        std::vector<int> buf(n+1, 0);
+        for(auto& [u,v] : items){ ++buf[u]; }
+        for(int i=1; i<=n; i++) buf[i] += buf[i-1];
+        res.m_list.resize(buf[n]);
+        for(int i=(int)items.size()-1; i>=0; i--){
+            res.m_list[--buf[items[i].first]] = std::move(items[i].second);
+        }
+        res.m_pos = std::move(buf);
+        return res;
+    }
+    static CsrArray FromRaw(std::vector<Elem> list, std::vector<int> pos){
+        CsrArray res;
+        res.m_n = pos.size() - 1;
+        res.m_list = std::move(list);
+        res.m_pos = std::move(pos);
+        return res;
+    }
+    ListRange operator[](int u) { return ListRange{ m_list.begin() + m_pos[u], m_list.begin() + m_pos[u+1] }; }
+    ConstListRange operator[](int u) const { return ConstListRange{ m_list.begin() + m_pos[u], m_list.begin() + m_pos[u+1] }; }
+    int size() const { return m_n; }
+    int fullSize() const { return (int)m_list.size(); }
+};
+
+// Int3 must be able to handle the value range :
+//  |x| <= | (any input - any input) ** 4 * 12 |
+
+template<class Int = long long, class Int2 = long long, class Int3 = Int2>
+class DelaunayTriangulation {
+public:
+
+    using GPos2 = VecI2<Int, Int2>;
+
+    struct Edge {
+        int to;
+        int ccw;
+        int cw;
+        int rev;
+        bool enabled = false;
+    };
+    
+private:
+
+    static int isDinOABC(GPos2 a, GPos2 b, GPos2 c, GPos2 d){
+        a = a - d;
+        b = b - d;
+        c = c - d;
+        auto val = Int3(b^c) * Int3(a.norm()) + Int3(c^a) * Int3(b.norm()) + Int3(a^b) * Int3(c.norm());
+        return val > Int3(0) ? 1 : 0;
+    }
+
+    int getOpenAddress(){
+        if(openAddress.empty()){
+            edges.push_back({});
+            return (int)edges.size() - 1;
+        }
+        int res = openAddress.back();
+        openAddress.pop_back();
+        return res;
+    }
+
+    std::pair<int, int> newEdge(int u, int v){
+        int euv = getOpenAddress();
+        int evu = getOpenAddress();
+        edges[euv].ccw = edges[euv].cw = euv;
+        edges[evu].ccw = edges[evu].cw = evu;
+        edges[euv].to = v;
+        edges[evu].to = u;
+        edges[euv].rev = evu;
+        edges[evu].rev = euv;
+        edges[euv].enabled = true;
+        edges[evu].enabled = true;
+        return { euv, evu };
+    }
+
+    void eraseSingleEdge(int e){
+        int eccw = edges[e].ccw;
+        int ecw = edges[e].cw;
+        edges[eccw].cw = ecw;
+        edges[ecw].ccw = eccw;
+        edges[e].enabled = false;
+    }
+
+    void eraseEdgeBidirectional(int e){
+        int ex = edges[e].rev;
+        eraseSingleEdge(e);
+        eraseSingleEdge(ex);
+        openAddress.push_back(e);
+        openAddress.push_back(ex);
+    }
+
+    void insertCcwAfter(int e, int x){
+        int xccw = edges[x].ccw;
+        edges[e].ccw = xccw;
+        edges[xccw].cw = e;
+        edges[e].cw = x;
+        edges[x].ccw = e;
+    }
+
+    void insertCwAfter(int e, int x){
+        int xcw = edges[x].cw;
+        edges[e].cw = xcw;
+        edges[xcw].ccw = e;
+        edges[e].ccw = x;
+        edges[x].cw = e;
+    }
+
+    // move from ab to ac ... is this ccw?
+    int isCcw(int a, int b, int c) const {
+        auto ab = pos[b] - pos[a];
+        auto ac = pos[c] - pos[a];
+        auto cp = ab ^ ac;
+        if(0 < cp) return 1;
+        if(cp < 0) return -1;
+        return 0;
+    }
+
+    std::pair<int, int> goNext(int , int ea){
+        int ap = edges[ea].to;
+        int eap = edges[edges[ea].rev].ccw;
+        return { ap, eap };
+    }
+
+    std::pair<int, int> goPrev(int , int ea){
+        int ap = edges[edges[ea].cw].to;
+        int eap = edges[edges[ea].cw].rev;
+        return { ap, eap };
+    }
+
+    std::tuple<int, int, int, int> goBottom(int a, int ea, int b, int eb){
+        while(true){
+            auto [ap, eap] = goPrev(a, ea);
+            if(isCcw(b, a, ap) > 0){
+                std::tie(a, ea) = { ap, eap };
+                continue;
+            }
+            auto [bp, ebp] = goNext(b, eb);
+            if(isCcw(a, b, bp) < 0){
+                std::tie(b, eb) = { bp, ebp };
+                continue;
+            }
+            break;
+        }
+        return { a, ea, b, eb };
+    }
+
+    std::pair<int, int> getMaximum(int a, int ea, bool toMin){
+        std::pair<int, int> ans = { a, ea };
+        int p = a, ep = ea;
+        do {
+            std::tie(p, ep) = goNext(p, ep);
+            if(toMin) ans = std::min(ans, std::make_pair(p, ep));
+            else ans = std::max(ans, std::make_pair(p, ep));
+        } while(ep != ea);
+        return ans;
+    }
+
+    bool isDinOABC(int a, int b, int c, int d){
+        return isDinOABC(pos[a], pos[b], pos[c], pos[d]);
+    }
+
+    std::pair<int, int> dfs(int a, int ea, int b, int eb){
+        std::tie(a, ea) = getMaximum(a, ea, false);
+        std::tie(b, eb) = getMaximum(b, eb, true);
+        auto [al, eal, bl, ebl] = goBottom(a, ea, b, eb);
+        auto [bu, ebu, au, eau] = goBottom(b, eb, a, ea);
+        ebl = edges[ebl].cw;
+        ebu = edges[ebu].cw;
+
+        auto [abl, bal] = newEdge(al, bl);
+        insertCwAfter(abl, eal);
+        insertCcwAfter(bal, ebl);
+        if(al == au) eau = abl;
+        if(bl == bu) ebu = bal;
+
+        int ap = al, eap = eal;
+        int bp = bl, ebp = ebl;
+        while(ap != au || bp != bu){
+            int a2 = edges[eap].to;
+            int b2 = edges[ebp].to;
+            int nxeap = edges[eap].ccw;
+            int nxebp = edges[ebp].cw;
+
+            if(eap != eau && nxeap != abl){
+                int a1 = edges[nxeap].to;
+                if(isDinOABC(ap, bp, a2, a1)){
+                    eraseEdgeBidirectional(eap);
+                    eap = nxeap;
+                    continue;
+                }
+            }
+
+            if(ebp != ebu && nxebp != bal){
+                int b1 = edges[nxebp].to;
+                if(isDinOABC(b2, ap, bp, b1)){
+                    eraseEdgeBidirectional(ebp);
+                    ebp = nxebp;
+                    continue;
+                }
+            }
+
+            bool chooseA = ebp == ebu;
+            if(eap != eau && ebp != ebu){
+				if(isCcw(ap, bp, b2) < 0) chooseA = true;
+				else if(isCcw(a2, ap, bp) < 0) chooseA = false;
+				else chooseA = isDinOABC(ap, bp, b2, a2);
+            }
+
+            if(chooseA){
+                nxeap = edges[edges[eap].rev].ccw;
+                auto [hab, hba] = newEdge(a2, bp);
+                insertCwAfter(hab, nxeap);
+                insertCcwAfter(hba, ebp);
+                eap = nxeap; ap = a2;
+            }
+            else {
+                nxebp = edges[edges[ebp].rev].cw;
+                auto [hba, hab] = newEdge(b2, ap);
+                insertCcwAfter(hba, nxebp);
+                insertCwAfter(hab, eap);
+                ebp = nxebp; bp = b2;
+            }
+        }
+
+        return { al, abl };
+    }
+
+    std::pair<int, int> solveRange(int l, int r){
+        if(r - l == 2){
+            int u = l;
+            int v = l + 1;
+            auto [uv, vu] = newEdge(u, v);
+            return { u, uv };
+        }
+        if(r - l == 3){
+            int u = l;
+            int v = l + 1;
+            int w = l + 2;
+            auto [uv, vu] = newEdge(u, v);
+            auto [vw, wv] = newEdge(v, w);
+            int ccw = isCcw(u, v, w);
+            if(ccw == 0){
+                insertCcwAfter(vu, vw);
+            }
+            if(ccw > 0){
+                auto [uw, wu] = newEdge(u, w);
+                insertCwAfter(uv, uw);
+                insertCwAfter(vw, vu);
+                insertCwAfter(wu, wv);
+                return { u, uv };
+            }
+            if(ccw < 0){
+                auto [uw, wu] = newEdge(u, w);
+                insertCcwAfter(uv, uw);
+                insertCcwAfter(vw, vu);
+                insertCcwAfter(wu, wv);
+                return { v, vu };
+            }
+            return { u, uv };
+        }
+        int m = (l + r) / 2;
+
+        auto [a, ea] = solveRange(l, m);
+        auto [b, eb] = solveRange(m, r);
+
+        return dfs(a, ea, b, eb);
+    }
+
+    void solve(){
+        int sz = (int)pos.size();
+        if(sz <= 1) return;
+
+        std::vector<int> pi(pos.size());
+        for(int i=0; i<(int)pi.size(); i++) pi[i] = i;
+        std::stable_sort(
+            pi.begin(), pi.end(),
+            [&](int l, int r){
+                return pos[l].x != pos[r].x ?
+                    pos[l].x < pos[r].x : pos[l].y < pos[r].y;
+            }
+        );
+        auto posbuf = pos;
+        int posptr = 0;
+        mappings.assign(sz, 0);
+        for(int i=0; i<sz; i++){
+            int v = pi[i];
+            if(i == 0 || !(posbuf[pi[posptr-1]] == posbuf[v])){
+                pi[posptr] = v;
+                pos[posptr++] = posbuf[v];
+                mappings[v] = v;
+            } else {
+                mappings[v] = pi[posptr-1];
+            }
+        }
+        
+        if(posptr >= 2) outerOneEdge = solveRange(0, posptr).second;
+        std::swap(pos, posbuf);
+        for(auto& e : edges) e.to = pi[e.to];
+    }
+    
+    std::vector<int> openAddress;
+    std::vector<GPos2> pos;
+    std::vector<Edge> edges;
+    std::vector<int> mappings;
+    int outerOneEdge = -1;
+
+public:
+
+    DelaunayTriangulation()
+        : pos()
+    { solve(); }
+
+    DelaunayTriangulation(std::vector<GPos2> x_points)
+        : pos(std::move(x_points))
+    {
+        solve();
+    }
+
+    std::vector<std::pair<int, int>> getEdges() const {
+        std::vector<std::pair<int, int>> res;
+        for(int e=0; e<(int)edges.size(); e++) if(edges[e].enabled){
+            int re = edges[e].rev;
+            if(e < re) continue;
+            res.push_back({ edges[e].to, edges[re].to });
+        }
+        for(int v=0; v<int(mappings.size()); v++){
+            if(mappings[v] != v) res.push_back({ v, mappings[v] });
+        }
+        return res;
+    }
+
+};
+
+struct DsuFast{
+private:
+    std::vector<int> w;
+public:
+    DsuFast(int n = 0) : w(n, -1) {}
+    int leader(int u){
+        if(w[u] < 0) return u;
+        return w[u] = leader(w[u]);
+    }
+    int operator[](int u){ return leader(u); }
+    int merge(int u, int v){
+        u = leader(u);
+        v = leader(v);
+        if(u == v) return u;
+        if(-w[u] < -w[v]) std::swap(u, v);
+        w[u] += w[v];
+        w[v] = u;
+        return u;
+    }
+    int size(int u){ return -w[leader(u)]; }
+    bool same(int u, int v){ return leader(u) == leader(v); }
+};
+
+#include <cstdio>
+
+int main(){
+    using namespace std;
+
+    using i64 = long long;
+    using Point = VecI2<i64, i64>;
+
+    int N; scanf("%d", &N);
+    vector<Point> A(N);
+    for(int i=0; i<N; i++){
+        int x, y; scanf("%d%d", &x, &y);
+        A[i].x = x;
+        A[i].y = y;
+    }
+    auto weight = [&](pair<int,int> p){ return (A[p.first] - A[p.second]).norm(); };
+
+    // if coord value bounding is like 0 <= x, y <= 10^9
+    //    DelaunayTriangulation<i64, i64, __int128_t>
+    auto tri = DelaunayTriangulation(A).getEdges();
+    stable_sort(tri.begin(), tri.end(),
+        [&](pair<int,int> l, pair<int,int> r) -> bool {
+            return weight(l) < weight(r); });
+
+    auto dsu = DsuFast(N);
+    vector<pair<int,int>> ans;
+    for(auto a : tri){
+        if(dsu.same(a.first, a.second)) continue;
+        dsu.merge(a.first, a.second);
+        if(a.first > a.second) swap(a.first, a.second);
+        ans.push_back(a);
+    }
+
+    stable_sort(ans.begin(), ans.end());
+
+    for(auto a : ans){
+        printf("%d %d\n", a.first, a.second);
+    }
+
+    return 0;
+}
diff --git a/geo/euclidean_mst/sol/naive_prim.cpp b/geo/euclidean_mst/sol/naive_prim.cpp
new file mode 100644
index 00000000..4a6f7c4f
--- /dev/null
+++ b/geo/euclidean_mst/sol/naive_prim.cpp
@@ -0,0 +1,52 @@
+#include <vector>
+#include <cstdio>
+
+using namespace std;
+
+
+int main() {
+    using ll = long long;
+
+    int N;
+    scanf("%d", &N);
+
+    vector<ll> X(N), Y(N);
+    for(int i=0; i<N; i++){
+        int x, y; scanf("%d%d", &x, &y);
+        X[i] = x;
+        Y[i] = y;
+    }
+
+    auto weight = [&](int a, int b){
+        ll dx = X[a] - X[b];
+        ll dy = Y[a] - Y[b];
+        return dx * dx + dy * dy;
+    };
+
+    vector<pair<ll, int>> nx(N, { ll(1) << 60, -1 });
+    vector<int> vis(N);
+    vector<pair<int,int>> ans;
+
+    vis[0] = 1;
+    for(int i=1; i<N; i++) nx[i] = { weight(0,i), 0 };
+
+    for(int t=0; t<N-1; t++){
+        int v = -1;
+        for(int w=0; w<N; w++) if(vis[w] == 0){
+            if(v == -1 || nx[w] < nx[v]) v = w;
+        }
+        vis[v] = 1;
+        ans.push_back({ v, nx[v].second });
+
+        for(int w=0; w<N; w++){
+            auto d = make_pair(weight(v,w), v);
+            if(d < nx[w]) nx[w] = d;
+        }
+    }
+
+    for(auto a : ans){
+        printf("%d %d\n", a.first, a.second);
+    }
+
+    return 0;
+}
diff --git a/geo/euclidean_mst/task.md b/geo/euclidean_mst/task.md
new file mode 100644
index 00000000..e315d4cb
--- /dev/null
+++ b/geo/euclidean_mst/task.md
@@ -0,0 +1,53 @@
+## @{keyword.statement}
+
+@{lang.en}
+Given $N$ 2D-points. $i$-th is $(x _ i, y _ i)$.
+
+For each $i$ , $j$, we add edge between them and whose weight is $\sqrt{(x _ i - x _ j)^2 + (y _ i - y _ j)^2}$.
+
+Calculate a minimum spanning tree of this graph.
+
+@{lang.ja}
+2次元平面上に $N$ 個の点が与えられる。$i$ 個目の頂点の座標は $(x _ i, y _ i)$ である。
+
+2点間の距離をユークリッド距離、つまり $\sqrt{(x _ i - x _ j)^2 + (y _ i - y _ j)^2}$ で定義するときの、最小全域木を求めよ。
+@{lang.end}
+
+## @{keyword.constraints}
+
+- $1 \leq N \leq @{param.N_MAX}$
+- $|x _ i|, |y _ i| \leq @{param.X_AND_Y_ABS_MAX}$
+
+## @{keyword.input}
+
+~~~
+$N$
+$x _ 0$ $y _ 0$
+$x _ 1$ $y _ 1$
+:
+$x _ {N - 1}$ $y _ {N - 1}$
+~~~
+
+## @{keyword.output}
+
+~~~
+$u _ 0$ $v _ 0$
+$u _ 1$ $v _ 1$
+:
+$u _ {N - 2}$ $v _ {N - 2}$
+~~~
+
+@{lang.en}
+Output all those $N-1$ pairs of $u _ i$ , $v _ i$ as the indices of the points that are the terminals of the edges in a minimum spanning tree.
+
+If there are multiple solutions, print any of them.
+@{lang.ja}
+最小全域木の辺が結ぶ点の番号を $u _ i$ , $v _ i$ として、それら $N-1$ 組をすべて出力せよ。
+
+解が複数存在する場合、どれを返しても構わない。
+@{lang.end}
+
+
+## @{keyword.sample}
+
+@{example.example_00}
diff --git a/geo/euclidean_mst/verifier.cpp b/geo/euclidean_mst/verifier.cpp
new file mode 100644
index 00000000..378272b8
--- /dev/null
+++ b/geo/euclidean_mst/verifier.cpp
@@ -0,0 +1,19 @@
+#include "testlib.h"
+#include "params.h"
+
+int main() {
+    registerValidation();
+
+    int n = inf.readInt(1, N_MAX);
+    inf.readChar('\n');
+    for (int i = 0; i < n; i++) {
+        inf.readInt(-X_AND_Y_ABS_MAX, X_AND_Y_ABS_MAX);
+        inf.readSpace();
+        inf.readInt(-X_AND_Y_ABS_MAX, X_AND_Y_ABS_MAX);
+        inf.readChar('\n');
+    }
+    
+    inf.readEof();
+
+    return 0;
+}