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;
+}
+
+
diff --git a/cpp/switchchain.cpp b/cpp/switchchain.cpp
index 30e4389eed3e925e9bbe5b4d4dd5b790503d056a..36b2c464badc33cb3a262ab9d03504fe4be0f5ec 100644
--- a/cpp/switchchain.cpp
+++ b/cpp/switchchain.cpp
@@ -1,8 +1,9 @@
+#include "switchchain.hpp"
 #include "exports.hpp"
 #include "graph.hpp"
-#include "powerlaw.hpp"
+#include "graph_gcm.hpp"
 #include "graph_spectrum.hpp"
-#include "switchchain.hpp"
+#include "powerlaw.hpp"
 #include <algorithm>
 #include <array>
 #include <fstream>
@@ -32,116 +33,6 @@ void getTriangleDegrees(const Graph& g) {
     return;
 }
 
-//
-// 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;
-}
-
 int main(int argc, char* argv[]) {
     // Generate a random degree sequence
     std::mt19937 rng(std::random_device{}());
diff --git a/cpp/switchchain_initialtris.cpp b/cpp/switchchain_initialtris.cpp
index d84d426b9f59c86cb73b6f2aba7341b905da2811..3acda57dd7e8e959c1c51deb9c4a02f6151e8dc0 100644
--- a/cpp/switchchain_initialtris.cpp
+++ b/cpp/switchchain_initialtris.cpp
@@ -1,5 +1,6 @@
 #include "exports.hpp"
 #include "graph.hpp"
+#include "graph_gcm.hpp"
 #include "powerlaw.hpp"
 #include "switchchain.hpp"
 #include <algorithm>
@@ -9,116 +10,6 @@
 #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;
-}
-
 int main() {
     // Generate a random degree sequence
     std::mt19937 rng(std::random_device{}());