diff --git a/cpp/powerlaw.hpp b/cpp/powerlaw.hpp
index b325a2ab534a2dbeb29cde138d8d310b5e9a9817..80ab891b09bb2ced7064964f1049530180e67076 100644
--- a/cpp/powerlaw.hpp
+++ b/cpp/powerlaw.hpp
@@ -45,8 +45,31 @@ class powerlaw_distribution {
     // 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
+    // Minimum value is min
+    // Maximum value is min * n^{1/(alpha-1)}
     template <class FwdIterator>
     void getFixedDistribution(int n, FwdIterator output) const {
+        unsigned int x;
+        // Now it outputs in increasing order
+        for (int i = 0; i < n; ++i) {
+            double r = 1.0 - double(i) / double(n);
+            x = std::round(std::pow(r, _exponent) * xmin);
+            if (x > xmax || x < xmin) {
+                std::cerr << "ERROR: x not in [xmin,xmax] : " << x
+                          << " not in [" << xmin << ',' << xmax
+                          << "]; i = " << i << std::endl;
+            }
+            *output++ = 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
+    // Minimum value is min
+    // Maximum value is (M^{-1} + M^{-(a-1)} - M^{-a})^{-1/(alpha-1)}
+    template <class FwdIterator>
+    void getFixedDistribution2(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);