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;

    }

    // 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;

    unsigned int xmin, xmax;

};