// Minimum value is min

    // Maximum value is min * n^{1/(alpha-1)}

        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 {