// newran.h ------------------------------------------------------------

// NEWRAN02

#ifndef NEWRAN_LIB
#define NEWRAN_LIB 0

//******************* utilities and definitions *************************

#include "include.h"
#include "boolean.h"
#include "myexcept.h"
#include "extreal.h"

#ifdef use_namespace
namespace NEWRAN { using namespace RBD_COMMON; }
namespace RBD_LIBRARIES { using namespace NEWRAN; }
namespace NEWRAN {
#endif

typedef Real (*PDF)(Real);                // probability density function

Real ln_gamma(Real);                      // log gamma function

//**************** uniform random number generator **********************

class RepeatedRandom;
class SelectedRandom;

class Random                              // uniform random number generator
{
   static double seed;                    // seed
   //static unsigned long iseed;          // for Mother
   static Real Buffer[128];               // for mixing random numbers
   static Real Raw();                     // unmixed random numbers
   void operator=(const Random&) {}       // private so can't access

public:
   static void Set(double s);             // set seed (0 < seed < 1)
   static double Get();                   // get seed
   virtual Real Next();                   // get new value
   virtual char* Name();                  // identification
   virtual Real Density(Real) const;      // used by PosGen & Asymgen
   Random() {}                            // do nothing
   virtual ~Random() {}                   // make destructors virtual
   SelectedRandom& operator()(double);    // used by MixedRandom
   RepeatedRandom& operator()(int);       // for making combinations
   virtual ExtReal Mean() const { return 0.5; }
                                          // mean of distribution
   virtual ExtReal Variance() const { return 1.0/12.0; }
					  // variance of distribution
   virtual void tDelete() {}              // delete components of sum
   virtual int nelems() const { return 1; }
                                          // used by MixedRandom
   virtual void load(int*,Real*,Random**);
   friend class RandomPermutation;
};


//****************** uniform random number generator *********************

class Uniform : public Random
{
   void operator=(const Uniform&) {}      // private so can't access

public:
   char* Name();                          // identification
   Uniform() {}                           // set value
   Real Next() { return Random::Next(); }
   ExtReal Mean() const { return 0.5; }
   ExtReal Variance() const { return 1.0/12.0; }
   Real Density(Real x) const { return (x < 0.0 || x > 1.0) ? 0 : 1.0; }
};


//************************* return constant ******************************

class Constant : public Random
{
   void operator=(const Constant&) {}     // private so can't access
   Real value;                            // value to be returned

public:
   char* Name();                          // identification
   Constant(Real v) { value=v; }          // set value
   Real Next() { return value; }
   ExtReal Mean() const { return value; }
   ExtReal Variance() const { return 0.0; }
};

//**************** positive random number generator **********************

class PosGen : public Random              // generate positive rv
{
   void operator=(const PosGen&) {}       // private so can't access

protected:
   Real xi, *sx, *sfx;
   bool NotReady;
   void Build(bool);                      // called on first call to Next

public:
   char* Name();                          // identification
   PosGen();                              // constructor
   ~PosGen();                             // destructor
   Real Next();                           // to get a single new value
   ExtReal Mean() const { return (ExtReal)Missing; }
   ExtReal Variance() const { return (ExtReal)Missing; }
};

//**************** symmetric random number generator **********************

class SymGen : public PosGen              // generate symmetric rv
{
   void operator=(const SymGen&) {}       // private so can't access

public:
   char* Name();                          // identification
   Real Next();                           // to get a single new value
};

//**************** normal random number generator **********************

class Normal : public SymGen              // generate standard normal rv
{
   void operator=(const Normal&) {}       // private so can't access
   static Real Nxi, *Nsx, *Nsfx;          // so we need initialise only once
   static long count;                     // assume initialised to 0

public:
   char* Name();                          // identification
   Normal();
   ~Normal();
   Real Density(Real) const;              // normal density function
   ExtReal Mean() const { return 0.0; }
   ExtReal Variance() const { return 1.0; }
};

//*********** chi-square random number generator **********************

class ChiSq : public Random               // generate non-central chi-sq rv
{
   void operator=(const ChiSq&) {}        // private so can't access
   Random* c1;                            // pointers to generators
   Random* c2;                            // pointers to generators
   int version;                           // indicates method of generation
   Real mean, var;

public:
   char* Name();                          // identification
   ChiSq(int, Real=0.0);                  // df and non-centrality parameter
   ~ChiSq();
   ExtReal Mean() const { return mean; }
   ExtReal Variance() const { return var; }
   Real Next();
};

//**************** Cauchy random number generator **********************

class Cauchy : public SymGen              // generate standard cauchy rv
{
   void operator=(const Cauchy&) {}       // private so can't access

public:
   char* Name();                          // identification
   Real Density(Real) const;              // Cauchy density function
   ExtReal Mean() const { return Indefinite; }
   ExtReal Variance() const { return PlusInfinity; }
};

//**************** Exponential random number generator **********************

class Exponential : public PosGen         // generate standard exponential rv
{
   void operator=(const Exponential&) {}  // private so can't access

public:
   char* Name();                          // identification
   Real Density(Real) const;              // Exponential density function
   ExtReal Mean() const { return 1.0; }
   ExtReal Variance() const { return 1.0; }
};

//**************** asymmetric random number generator **********************

class AsymGen : public Random             // generate asymmetric rv
{
   void operator=(const AsymGen&) {}      // private so can't access
   Real xi, *sx, *sfx; int ic;
   bool NotReady;
   void Build();                          // called on first call to Next

protected:
   Real mode;

public:
   char* Name();                          // identification
   AsymGen(Real);                         // constructor (Real=mode)
   ~AsymGen();                            // destructor
   Real Next();                           // to get a single new value
   ExtReal Mean() const { return (ExtReal)Missing; }
   ExtReal Variance() const { return (ExtReal)Missing; }
};

//**************** Gamma random number generator **********************

class Gamma : public Random               // generate gamma rv
{
   void operator=(const Gamma&) {}        // private so can't access
   Random* method;

public:
   char* Name();                          // identification
   Gamma(Real);                           // constructor (Real=shape)
   ~Gamma() { delete method; }
   Real Next() { return method->Next(); }
   ExtReal Mean() const { return method->Mean(); }
   ExtReal Variance() const { return method->Variance(); }
};

//**************** Generators with pointers to pdf ***********************

class PosGenX : public PosGen
{
   void operator=(const PosGenX&) {}      // private so can't access
   PDF f;

public:
   char* Name();                          // identification
   PosGenX(PDF fx);
   Real Density(Real x) const { return (*f)(x); }
};

class SymGenX : public SymGen
{
   void operator=(const SymGenX&) {}      // private so can't access
   PDF f;

public:
   char* Name();                          // identification
   SymGenX(PDF fx);
   Real Density(Real x) const { return (*f)(x); }
};

class AsymGenX : public AsymGen
{
   void operator=(const AsymGenX&) {}     // private so can't access
   PDF f;

public:
   char* Name();                          // identification
   AsymGenX(PDF fx, Real mx);
   Real Density(Real x) const { return (*f)(x); }
};

//***************** Pareto random number generator ************************

class Pareto : public Random
// Use definition of Kotz and Johnson: "Continuous univariate distributions 1",
// chapter 19 with k = 1.
{
   void operator=(const Pareto&) {}       // private so can't access
   Real Shape, RS;

public:
   char* Name();                          // identification
   Pareto(Real shape);                    // constructor (Real=shape)
   ~Pareto() {}
   Real Next();
   ExtReal Mean() const;
   ExtReal Variance() const;
};


//**************** discrete random number generator **********************

class DiscreteGen : public Random         // discrete random generator
{
   void operator=(const DiscreteGen&) {}  // private so can't access
   Real *p; int *ialt; int n; Real *val;
   void Gen(int, Real*);                  // called by constructors
   Real mean, var;                        // calculated by the constructor

public:
   char* Name();                          // identification
   DiscreteGen(int,Real*);                // constructor
   DiscreteGen(int,Real*,Real*);          // constructor
   ~DiscreteGen();                        // destructor
   Real Next();                           // new single value
   ExtReal Mean() const { return mean; }
   ExtReal Variance() const { return var; }
};

//***************** Poisson random number generator *******************