// newran.cpp -----------------------------------------------------------

// NEWRAN02B - 22 July 2002

#define WANT_STREAM
#define WANT_MATH
#include "include.h"

#include "newran.h"
//#include "mother.h"

#ifdef use_namespace
namespace NEWRAN {
#endif

//********* classes which are used internally only **********************

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

class ChiSq1 : public Normal              // generate non-central chi-square
					  // rv with 1 degree of freedom
{
   Real deltasq;                          // non-centrality parameter
   Real delta;

public:
   ChiSq1(Real);                          // non-centrality parameter
   ExtReal Mean() const { return 1.0+deltasq; }
   ExtReal Variance() const { return 2.0+4.0*deltasq; }
   Real Next();
};

//*********** Poisson random number generator, larger mu ****************

class Poisson1 : public AsymGen           // generate poisson rv, large mu
{
   Real mu, ln_mu;
public:
   Poisson1(Real);                        // constructor (Real=mean)
   Real Density(Real) const;              // Poisson density function
   Real Next() { return floor(AsymGen::Next()); }
   ExtReal Mean() const { return mu; }
   ExtReal Variance() const { return mu; }
};

//*********** Poisson random number generator, mu <= 10 ****************

class Poisson2 : public Random            // generate poisson rv, large mu
{
   DiscreteGen* dg;
public:
   Poisson2(Real);                        // constructor (Real=mean)
   ~Poisson2();
   Real Next() { return dg->Next(); }
   ExtReal Mean() const { return dg->Mean(); }
   ExtReal Variance() const { return dg->Variance(); }
};

//********** Gamma random number generator, alpha <= 1 *****************

class Gamma1 : public PosGen              // generate gamma rv
					  // shape parameter <= 1
{
   Real ln_gam, ralpha, alpha;
public:
   Gamma1(Real);                          // constructor (Real=shape)
   Real Density(Real) const;              // gamma density function
   Real Next();                           // carries out power transform
   ExtReal Mean() const { return alpha; }
   ExtReal Variance() const { return alpha; }
};

//********** Gamma random number generator, alpha > 1 ******************

class Gamma2 : public AsymGen             // generate gamma rv
					  // shape parameter > 1
{
   Real alpha, ln_gam;
public:
   Gamma2(Real);                          // constructor (Real=shape)
   Real Density(Real) const;              // gamma density function
   ExtReal Mean() const { return alpha; }
   ExtReal Variance() const { return alpha; }
};

//*********** Binomial random number generator, n > 40 *****************

class Binomial1 : public AsymGen           // generate binomial rv, larger n
{
   Real p, q, ln_p, ln_q, ln_n_fac; int n;
public:
   Binomial1(int nx, Real px);
   Real Density(Real) const;
   Real Next() { return floor(AsymGen::Next()); }
   ExtReal Mean() const { return p * n; }
   ExtReal Variance() const { return p * q * n; }
};

//******* Binomial random number generator, n < 40 or n*p <= 8 *************

class Binomial2 : public Random            // generate binomial rv, smaller n
{
   DiscreteGen* dg;
public:
   Binomial2(int nx, Real px);
   ~Binomial2();
   Real Next() { return dg->Next(); }
   ExtReal Mean() const { return dg->Mean(); }
   ExtReal Variance() const { return dg->Variance(); }
};

//************************ static variables ***************************

double Random::seed;
//unsigned long Random::iseed;                // for Mother
Real Random::Buffer[128];
Real Normal::Nxi;
Real* Normal::Nsx;
Real* Normal::Nsfx;
long Normal::count=0;

//**************************** utilities ******************************

inline Real square(Real x) { return x*x; }
inline ExtReal square(const ExtReal& x) { return x*x; }

static void ErrorNoSpace() { Throw(Bad_alloc("Newran: out of space")); }

//************************* end of definitions ************************


Real Random::Raw()                           // get new uniform random number
{
   // m = 2147483647 = 2^31 - 1; a = 16807;
   // 127773 = m div a; 2836 = m mod a
   long iseed = (long)seed;
   long hi = iseed / 127773L;                 // integer division
   long lo = iseed - hi * 127773L;            // modulo
   iseed = 16807 * lo - 2836 * hi;
   if (iseed <= 0) iseed += 2147483647L;
   seed = (double)iseed; return seed*4.656612875e-10;
}

Real Random::Density(Real) const
{ Throw(Logic_error("density function not defined")); return 0.0; }

#ifdef _MSC_VER
static void DoNothing(int) {}
#endif

Real Random::Next()                          // get new mixed random number
{
   if (!seed)
      Throw(Logic_error("Random number generator not initialised"));
   int i = (int)(Raw()*128);               // 0 <= i < 128
#ifdef _MSC_VER
   DoNothing(i); DoNothing(i);
#endif
   Real f = Buffer[i]; Buffer[i] = Raw();  // Microsoft release gets this wrong
   return f;

   // return Mother(&iseed);
}

double Random::Get()                  // get random number seed
{ return seed/2147483648L; }

void Random::Set(double s)            // set random number seed
                                      // s must be between 0 and 1
{
   if (s>=1.0 || s<=0.0)
      Throw(Logic_error("Newran: seed out of range"));
   //iseed = 2147483648L * s;         // for Mother
   seed = (long)(s*2147483648L);
   for (int i = 0; i<128; i++) Buffer[i] = Raw();
}


//Real Constant::Next() { return value; }


PosGen::PosGen()                             // Constructor
{
   #ifdef MONITOR
      cout << "constructing PosGen\n";
   #endif
   NotReady=true;
}

PosGen::~PosGen()
{
   if (!NotReady)
   {
      #ifdef MONITOR
	 cout << "freeing PosGen arrays\n";
      #endif
      delete [] sx; delete [] sfx;
   }
   #ifdef MONITOR
      cout << "destructing PosGen\n";
   #endif
}

void PosGen::Build(bool sym)                 // set up arrays
{
   #ifdef MONITOR
      cout << "building PosGen arrays\n";
   #endif
   int i;
   NotReady=false;
   sx=new Real[60]; sfx=new Real[60];
   if (!sx || !sfx) ErrorNoSpace();
   Real sxi=0.0; Real inc = sym ? 0.01 : 0.02;
   for (i=0; i<60; i++)
   {
      sx[i]=sxi; Real f1=Density(sxi); sfx[i]=f1;
      if (f1<=0.0) goto L20;
      sxi+=inc/f1;
   }
   Throw(Runtime_error("Newran: area too large"));
L20:
   if (i<50) Throw(Runtime_error("Newran: area too small"));
   xi = sym ? 2*i : i;
   return;
}

Real PosGen::Next()
{
   Real ak,y; int ir;
   if (NotReady) Build(false);
   do
   {
      Real r1=Random::Next();
      ir = (int)(r1*xi); Real sxi=sx[ir];
      ak=sxi+(sx[ir+1]-sxi)*Random::Next();
      y=sfx[ir]*Random::Next();
   }
   while ( y>=sfx[ir+1] && y>=Density(ak) );
   return ak;
}

Real SymGen::Next()
{
   Real s,ak,y; int ir;
   if (NotReady) Build(true);
   do
   {
      s=1.0;
      Real r1=Random::Next();
      if (r1>0.5) { s=-1.0; r1=1.0-r1; }
      ir = (int)(r1*xi); Real sxi=sx[ir];
      ak=sxi+(sx[ir+1]-sxi)*Random::Next();
      y=sfx[ir]*Random::Next();
   }
   while ( y>=sfx[ir+1] && y>=Density(ak) );
   return s*ak;
}

AsymGen::AsymGen(Real modex)                 // Constructor
{
   #ifdef MONITOR
      cout << "constructing AsymGen\n";
   #endif
   mode=modex; NotReady=true;
}

void AsymGen::Build()                        // set up arrays
{
   #ifdef MONITOR
      cout << "building AsymGen arrays\n";
   #endif
   int i;
   NotReady=false;
   sx=new Real[121]; sfx=new Real[121];
   if (!sx || !sfx)  ErrorNoSpace();
   Real sxi=mode;
   for (i=0; i<120; i++)
   {
      sx[i]=sxi; Real f1=Density(sxi); sfx[i]=f1;
      if (f1<=0.0) goto L20;
      sxi+=0.01/f1;
   }
   Throw(Runtime_error("Newran: area too large (a)"));
L20:
   ic=i-1; sx[120]=sxi; sfx[120]=0.0;
   sxi=mode;
   for (; i<120; i++)
   {
      sx[i]=sxi; Real f1=Density(sxi); sfx[i]=f1;
      if (f1<=0.0) goto L30;
      sxi-=0.01/f1;
   }
   Throw(Runtime_error("Newran: area too large (b)"));
L30:
   if (i<100)  Throw(Runtime_error("Newran: area too small"));
   xi=i;
   return;
}

Real AsymGen::Next()
{
   Real ak,y; int ir1;
   if (NotReady) Build();
   do
   {
      Real r1=Random::Next();
      int ir=(int)(r1*xi); Real sxi=sx[ir];
      ir1 = (ir==ic) ? 120 : ir+1;
      ak=sxi+(sx[ir1]-sxi)*Random::Next();
      y=sfx[ir]*Random::Next();
   }
   while ( y>=sfx[ir1] && y>=Density(ak) );
   return ak;
}

AsymGen::~AsymGen()
{
   if (!NotReady)
   {
      #ifdef MONITOR
	 cout << "freeing AsymGen arrays\n";
      #endif
      delete [] sx; delete [] sfx;
   }
   #ifdef MONITOR
      cout << "destructing AsymGen\n";
   #endif
}

PosGenX::PosGenX(PDF fx) { f=fx; }

SymGenX::SymGenX(PDF fx) { f=fx; }

AsymGenX::AsymGenX(PDF fx, Real mx) : AsymGen(mx) { f=fx; }


Normal::Normal()
{
   if (count) { NotReady=false; xi=Nxi; sx=Nsx; sfx=Nsfx; }
   else { Build(true); Nxi=xi; Nsx=sx; Nsfx=sfx; }
   count++;
}

Normal::~Normal()
{
   count--;
   if (count) NotReady=true;                     // disable freeing arrays
}

Real Normal::Density(Real x) const               // normal density
{ return (fabs(x)>8.0) ? 0 : 0.398942280 * exp(-x*x / 2); }

ChiSq1::ChiSq1(Real d) : Normal()                // chisquare with 1 df
{ deltasq=d; delta=sqrt(d); }

Real ChiSq1::Next()
{ Real z=Normal::Next()+delta; return z*z; }

ChiSq::ChiSq(int df, Real noncen)
{
  if (df<=0 || noncen<0.0) Throw(Logic_error("Newran: illegal parameters"));
  else if (df==1) { version=0; c1=new ChiSq1(noncen); }
  else if (noncen==0.0)
  {
     if (df==2) { version=1; c1=new Exponential(); }
     else { version=2; c1=new Gamma2(0.5*df); }
  }