diff --git a/cpp/graph_gcm.hpp b/cpp/graph_gcm.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b53a8acb6625f4f23ac6b1c5fc4ad33ee8c464a3
--- /dev/null
+++ b/cpp/graph_gcm.hpp
@@ -0,0 +1,117 @@
+#include "graph.hpp"
+#include <algorithm>
+#include <iostream>
+#include <random>
+#include <vector>
+
+//
+// Assumes degree sequence does NOT contain any zeroes!
+//
+// method2 = true  -> take highest degree and finish its pairing completely
+// method2 = false -> take new highest degree after every pairing
+template <typename RNG>
+bool greedyConfigurationModel(DegreeSequence& ds, Graph& g, RNG& rng, bool method2) {
+    // Similar to Havel-Hakimi but instead of pairing up with the highest ones
+    // that remain, simply pair up with random ones
+    unsigned int n = ds.size();
+
+    // degree, vertex index
+    std::vector<std::pair<unsigned int, unsigned int>> degrees(n);
+    for (unsigned int i = 0; i < n; ++i) {
+        degrees[i].first = ds[i];
+        degrees[i].second = i;
+    }
+
+    std::vector<decltype(degrees.begin())> available;
+    available.reserve(n);
+
+    // Clear the graph
+    g.reset(n);
+
+    while (!degrees.empty()) {
+        std::shuffle(degrees.begin(), degrees.end(), rng);
+        // Get the highest degree:
+        // If there are multiple highest ones, we pick a random one,
+        // ensured by the shuffle.
+        // The shuffle is needed anyway for the remaining part
+        unsigned int dmax = 0;
+        auto uIter = degrees.begin();
+        for (auto iter = degrees.begin(); iter != degrees.end(); ++iter) {
+            if (iter->first >= dmax) {
+                dmax = iter->first;
+                uIter = iter;
+            }
+        }
+
+        if (dmax > degrees.size() - 1)
+            return false;
+
+        if (dmax == 0) {
+            std::cerr << "ERROR 1 in GCM.\n";
+        }
+
+        unsigned int u = uIter->second;
+
+        if (method2) {
+            // Take the highest degree out of the vector
+            degrees.erase(uIter);
+
+            // Now assign randomly to the remaining vertices
+            // Since its shuffled, we can pick the first 'dmax' ones
+            auto vIter = degrees.begin();
+            while (dmax--) {
+                if (vIter->first == 0)
+                    std::cerr << "ERROR in GCM2.\n";
+                if (!g.addEdge({u, vIter->second}))
+                    std::cerr << "ERROR. Could not add edge in GCM2.\n";
+                vIter->first--;
+                if (vIter->first == 0)
+                    vIter = degrees.erase(vIter);
+                else
+                    vIter++;
+            }
+        } else {
+            // Pair with a random vertex that is not u itself and to which
+            // u has not paired yet!
+            available.clear();
+            for (auto iter = degrees.begin(); iter != degrees.end(); ++iter) {
+                if (iter->second != u && !g.hasEdge({u, iter->second}))
+                    available.push_back(iter);
+            }
+            if (available.empty())
+                return false;
+            std::uniform_int_distribution<> distr(0, available.size() - 1);
+            auto vIter = available[distr(rng)];
+            // pair u to v
+            if (vIter->first == 0)
+                std::cerr << "ERROR 2 in GCM1.\n";
+            if (!g.addEdge({u, vIter->second}))
+                std::cerr << "ERROR. Could not add edge in GCM1.\n";
+            // Purge anything with degree zero
+            // Be careful with invalidating the other iterator!
+            // Degree of u is always greater or equal to the degree of v
+            if (dmax == 1) {
+                // Remove both
+                // Erasure invalidates all iterators AFTER the erased one
+                if (vIter > uIter) {
+                    degrees.erase(vIter);
+                    degrees.erase(uIter);
+                } else {
+                    degrees.erase(uIter);
+                    degrees.erase(vIter);
+                }
+            } else {
+                // Remove only v if it reaches zero
+                uIter->first--;
+                vIter->first--;
+                if (vIter->first == 0)
+                    degrees.erase(vIter);
+            }
+            //degrees.erase(std::remove_if(degrees.begin(), degrees.end(),
+            //                             [](auto x) { return x.first == 0; }));
+        }
+    }
+    return true;
+}
+
+