diff --git a/cpp/Makefile b/cpp/Makefile
index 569442ea869954ee83b0529705aaccb75384da60..0c4f0e46df1871d2d1630ec3d3dde19c2a8e542c 100644
--- a/cpp/Makefile
+++ b/cpp/Makefile
@@ -20,6 +20,7 @@ TARGETS += switchchain_mixingtime
 TARGETS += switchchain_spectrum
 TARGETS += switchchain_successrates
 TARGETS += switchchain_timeevol
+TARGETS += bruteforce_triangles
 
 all: $(TARGETS)
 
diff --git a/cpp/bruteforce_triangles.cpp b/cpp/bruteforce_triangles.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9b0eb2f4fcf045fffc28d85f07fc1685f8cf14d3
--- /dev/null
+++ b/cpp/bruteforce_triangles.cpp
@@ -0,0 +1,58 @@
+#include "exports.hpp"
+#include "graph.hpp"
+#include <iostream>
+#include <map>
+
+int main() {
+    // Bruteforce over all graphs of n nodes
+    constexpr int n = 8;
+    constexpr int E = n * (n - 1) / 2;
+
+    std::map<DegreeSequence, unsigned int> optimalTriangles;
+
+    std::cout << "Iterating over all graphs with " << n << " vertices." << std::endl;
+    Graph g;
+    for (int edgeMask = 0; edgeMask < (1 << E); ++edgeMask) {
+        g.reset(n);
+        int mask = edgeMask;
+        for (auto i = 0u; i < n; ++i) {
+            for (auto j = i + 1; j < n; ++j) {
+                if (mask & 1) {
+                    g.addEdge({i, j});
+                }
+                mask >>= 1;
+            }
+        }
+        auto triangles = g.countTriangles();
+
+        auto ds = g.getDegreeSequence();
+        std::sort(ds.begin(), ds.end());
+        auto iter = optimalTriangles.find(ds);
+        if (iter != optimalTriangles.end()) {
+            if (iter->second < triangles)
+                iter->second = triangles;
+        } else {
+            optimalTriangles[ds] = triangles;
+        }
+    }
+
+    std::cout << "Comparing with Erdos-Gallai triangles." << std::endl;
+    int total = 0;
+    int unequal = 0;
+    for (auto &keyvalue : optimalTriangles) {
+        auto ds = keyvalue.first;
+        auto optimalTris = keyvalue.second;
+        g.createFromDegreeSequence(ds);
+        auto EGtris = g.countTriangles();
+        if (EGtris < optimalTris) {
+            std::cout << "Erdos-Gallai triangles: " << EGtris
+                      << "; Optimal triangles: " << optimalTris
+                      << "; Difference: " << (optimalTris - EGtris)
+                      << "; ds = " << ds << std::endl;
+            unequal++;
+        }
+        total++;
+    }
+    std::cout << "Done." << std::endl;
+    std::cout << unequal << " out of " << total << " degree sequences gave non-optimal triangles.\n";
+}
diff --git a/cpp/graph.hpp b/cpp/graph.hpp
index 13f997a413f3e565161162a31b8808f43b32d99d..987d9a89ddf802d397f3af178816d475e21f92ec 100644
--- a/cpp/graph.hpp
+++ b/cpp/graph.hpp
@@ -79,7 +79,12 @@ class Graph {
         reset(n);
 
         while (!degrees.empty()) {
-            std::sort(degrees.begin(), degrees.end());
+            // Construction will find maximum triangles only if sort is stable
+            // and does NOT sort on vertex id
+            std::stable_sort(degrees.begin(), degrees.end(),
+                             [](const auto &p1, const auto &p2) {
+                                 return p1.first < p2.first;
+                             });
             // Highest degree is at back of the vector
             // Take it out
             unsigned int degree = degrees.back().first;