#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;
};