diff --git a/cpp/powerlaw.hpp b/cpp/powerlaw.hpp
index 5dd0194ce349f77bf3f82f642f43d0e0bd2656eb..b325a2ab534a2dbeb29cde138d8d310b5e9a9817 100644
--- a/cpp/powerlaw.hpp
+++ b/cpp/powerlaw.hpp
@@ -1,10 +1,23 @@
 #pragma once
 #include <cassert>
 #include <cmath>
+#include <iostream>
 #include <random>
 
 // Discrete powerlaw distribution
 // obtained by rounding a continuous powerlaw distribution
+// Let m = minimum and M = maximum and a = alpha
+// Then for continuous version we have
+// PDF: f(x) = C x^{-a}
+// with C = (1-a)/(M^{1-a} - m^{1-a})
+// This gives CDF: F(x) = (x^{1-a} - m^{1-a}) / (M^{1-a} - m^{1-a})
+// with inversion F^{-1}(y) = ( y M^{1-a} + (1-y) m^{1-a} )^{1/(1-a)}
+// i.e. linear interpolate between M^{1-a} < m^{1-a}
+// For M = infty this becomes
+//  F'^{-1}(y) = m (1-y)^{1/(1-a)}
+//             = m (y')^{-1/(a-1)}
+// where higher y' means lower x
+
 class powerlaw_distribution {
   public:
     // xmin and xmax are inclusive
@@ -29,6 +42,30 @@ class powerlaw_distribution {
         return x;
     }
 
+    // Generate non-random ``canonical'' powerlaw distribution
+    // Same as above operator but where the random number between
+    // 0 and 1 is replaced by stepping from 0 to 1 in fixed steps
+    template <class FwdIterator>
+    void getFixedDistribution(int n, FwdIterator output) const {
+        // go in linear steps over interval
+        // [ M^{1-a} , m^{1-a} ]
+        double minVal = std::pow(double(xmax), 1.0f - alpha);
+        double maxVal = std::pow(double(xmin), 1.0f - alpha);
+        unsigned int x;
+        // Now it outputs in increasing order
+        for (int i = 0; i < n; ++i) {
+            double r1 = double(i) / double(n);
+            double r2 = r1 * minVal + (1.0 - r1) * maxVal;
+            x = std::round(std::pow(r2, _exponent));
+            if (x > xmax || x < xmin) {
+                std::cerr << "ERROR: x not in [xmin,xmax] : " << x
+                          << " not in [" << xmin << ',' << xmax
+                          << "]; i = " << i << std::endl;
+            }
+            *output++ = x;
+        }
+    }
+
   private:
     float alpha;
     double _exponent;