// std_random.cpp


#include <cmath>

#include "std_random.hxx"


namespace blogs
{


knuth_b Std_Random::__Engine__;


uniform_real_distribution<double> Std_Random::__Distri__(0, 1);


void
Std_Random::set_seed(long seed)
{
    __Engine__.seed(seed);
}


double
Std_Random::random()
{
    return __Distri__(__Engine__);
}


int
Std_Random::uniform(int N)
{
    return (int)(random() * N);
}


int
Std_Random::uniform(int lo, int hi)
{
    return lo + (int)(random()*(hi - lo));
}


long
Std_Random::uniform(long N)
{
    return (long)(random() * N);
}


long
Std_Random::uniform(long lo, long hi)
{
    return lo + (long)(random()*(hi - lo));
}


double
Std_Random::uniform(double lo, double hi)
{
    return lo + random() * (hi - lo);
}


bool
Std_Random::bernoulli(double p)
{
    return random() < p;
}


double
Std_Random::gaussian(double mu, double sgm)
{
    double r, x, y;
    do {
        x = uniform(-1.0, +1.0);
        y = uniform(-1.0, +1.0);
        r = x*x + y*y;
    } while (r >= 1 || r == 0);
    return mu + sgm * x * std::sqrt(-2 * std::log(r) / r);
}


int
Std_Random::poisson(double lmd)
{
    int k = 0;
    double p = 1.0;
    double exp = std::exp(-lmd);
    do {
        ++k;
        p *= random();
    } while (p >= exp);
    return k-1;
}


int
Std_Random::geometric(double p)
{
    return (int)std::ceil(std::log(random()) / std::log(1.0 - p));
}


int
Std_Random::discrete(Array<double> const& probs)
{
    double r, sum;
    while (true) {
        r = random();
        sum = 0.0;
        for (int i = 0; i < probs.size(); ++i) {
            sum += probs[i];
            if (sum > r) return i;
        }
    }
}


} // end of namespace blogs