#ifndef Rnd_IS_INCLUDED
#define Rnd_IS_INCLUDED

//Included from Gnu Scientific Library:
#include <gsl/gsl_rng.h>     //Provides uniform random number generators.
#include <gsl/gsl_randist.h> //Extends uniform random number generator to
			     //produce random numbers distrbutions
			     //(i.e. gaussian).

using namespace std;

struct Rnd {
    gsl_rng* r;
    /* gsl-specific content:
     * struct gsl_rng; (random number generator data structure)
     */

    Rnd() {}
    Rnd(unsigned long int s) {construct(s);}
    void construct(unsigned long int s) {
	r = gsl_rng_alloc(gsl_rng_default);	
	/* gsl-specific content:
	 * gsl_rng* gsl_rng_alloc(const gsl_rng_type*);
	 * (data structure allocator for the random number generator)
	 * (the parameter gsl_rng_default is equavalent to gsl_rng_mt19937
	 *  which initializes the random number generator to use the Mersenne
	 *  Twister algorithm)
	 */

	gsl_rng_set(r, s);
	/* gsl-specific content:
	 * void gsl_rng_set (const gsl_rng*, unsigned long int);
	 * (sets the seed of the random number generator *r to s)
	 */

    } // end constructor

    inline double getUniform() {
	return gsl_rng_uniform_pos(r);
	/* gsl-specific content:
	 * double gsl_rng_uniform_pos(const gsl_rng*)
	 * (samples a random number from the random number generator *r,
	 * distributed uniformly on the interval (0,1), excluding both
	 * endpoints)
	 */
    }

    inline double getGaussian(double sigma) {
	return gsl_ran_gaussian_ziggurat(r, sigma);
	/* gsl-specific content:
	 * double gsl_ran_gaussian_ziggurat(const gsl_rng*, double);
	 * (samples a random number from the random number generator *r,
	 *  distributed normally with mean 0 and standard deviation sigma,
	 *  using the Marsaglia-Tsang ziggurat method)
	 */
    }

    inline double getGaussian() {
	return gsl_ran_gaussian_ziggurat(r, 1);
	/* gsl-specific content:
	 * double gsl_ran_gaussian_ziggurat(const gsl_rng*, double);
	 * (samples a random number from the random number generator *r,
	 *  distributed normally with mean 0 and standard deviation set to 1,
	 *  using the Marsaglia-Tsang ziggurat method)
	 */
    }

}; // end struct Rnd

#endif