diff --git a/cpp/powerlaw.hpp b/cpp/powerlaw.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5dd0194ce349f77bf3f82f642f43d0e0bd2656eb
--- /dev/null
+++ b/cpp/powerlaw.hpp
@@ -0,0 +1,36 @@
+#pragma once
+#include <cassert>
+#include <cmath>
+#include <random>
+
+// Discrete powerlaw distribution
+// obtained by rounding a continuous powerlaw distribution
+class powerlaw_distribution {
+  public:
+    // xmin and xmax are inclusive
+    powerlaw_distribution(float a, unsigned int _xmin, unsigned int _xmax)
+        : alpha(a), xmin(_xmin), xmax(_xmax) {
+        assert(xmin > 0);
+        assert(xmin <= xmax);
+        assert(alpha > 1.0f);
+        _exponent = -1.0f / (alpha - 1.0f);
+    }
+    ~powerlaw_distribution() {}
+
+    template <class Generator>
+    unsigned int operator()(Generator& rng) const {
+        unsigned int x = xmax + 1;
+        while (x > xmax) {
+            // Generate uniform value in [0,1)
+            double r = std::generate_canonical<
+                double, std::numeric_limits<double>::digits>(rng);
+            x = std::round(std::pow(r, _exponent) * xmin);
+        }
+        return x;
+    }
+
+  private:
+    float alpha;
+    double _exponent;
+    unsigned int xmin, xmax;
+};