How do you simulate an Erlang(k) variable using C++?

[K33p3r]

Good day,

I have a college project in which I need to create a program which simulates an Erlang(k) variable. I am not talking about the Erlang language, it's about a standard Gamma(0,1,k), k ϵ N+. I appreciate if someone can help me out or give me advice on how to approach this problem.

I've written the formulas and theorem that comes with this problem: Erlang distribution equations

All the best,
K33p3r

[m4ster_r0shi]

If you have to code this yourself, one way to do it would be to use the inverse transformation and acceptance-rejection methods.

I had to simulate normal distribution a while ago and that's how I did it:

#include <iostream>
#include <algorithm>
#include <vector>
#include <cstdlib>
#include <ctime>
#include <cmath>

using namespace std;

double random() { return rand() / double(RAND_MAX); }

struct Distribution
{
    virtual double Draw() { return 0.0; }
    virtual ~Distribution() {}
};

struct UniformDistribution : Distribution
{
    double a, b;

    UniformDistribution(double a_ = 0.0, double b_ = 1.0) :  a(a_), b(b_) {}
    ~UniformDistribution() {}

    // inverse transform sampling
    double Draw() { return a + (b - a) * random(); }
};

struct NormalDistribution : Distribution
{
    UniformDistribution ud1, ud2;

    double mean, dev, max;

    NormalDistribution(double mean_, double dev_) : mean(mean_), dev(dev_)
    {
        max = 1.0 / sqrt(2 * 3.14159265 * dev * dev);

        ud1 = UniformDistribution(mean - 4 * dev, mean + 4 * dev);
        ud2 = UniformDistribution(0, max);
    }
    ~NormalDistribution() {}

    double density(double x) { return max * exp(- (x - mean) * (x - mean) / (2 * dev * dev)); }

    // rejection sampling
    double Draw()
    {
        double x = ud1.Draw();
        double y = ud2.Draw();

        if (y <= density(x)) return x;
        else return Draw();
    }
};

int main()
{
    srand(time(0));

    NormalDistribution nd(10, 2);

    vector<double> random_numbers;

    for (size_t i = 0; i < 100; ++i) random_numbers.push_back(nd.Draw());

    sort(random_numbers.begin(), random_numbers.end());

    for (size_t i = 0; i < 100; ++i) cout << random_numbers[i] << endl;
}

If you are allowed to use ready-made stuff, take a look at Boost.Random.