Real AddedRandom::Next() { return rv1->Next() + rv2->Next() ; }

ExtReal AddedRandom::Mean() const { return rv1->Mean() + rv2->Mean() ; }

ExtReal AddedRandom::Variance() const
   { return rv1->Variance() + rv2->Variance() ; }

Real SubtractedRandom::Next() { return rv1->Next() - rv2->Next() ; }

ExtReal SubtractedRandom::Mean() const { return rv1->Mean() - rv2->Mean() ; }

ExtReal SubtractedRandom::Variance() const
   { return rv1->Variance() + rv2->Variance() ; }

Real MultipliedRandom::Next() { return rv1->Next() * rv2->Next() ; }

ExtReal MultipliedRandom::Mean() const { return rv1->Mean() * rv2->Mean() ; }

ExtReal MultipliedRandom::Variance() const
{
   ExtReal e1 = square(rv1->Mean()); ExtReal e2 = square(rv2->Mean());
   ExtReal v1 = rv1->Variance(); ExtReal v2 = rv2->Variance();
   ExtReal r=v1*v2+v1*e2+e1*v2;
   return r;
}

Real DividedRandom::Next() { return rv1->Next() / rv2->Next() ; }

void Random::load(int*,Real*,Random**)
   { Throw(Logic_error("Newran: illegal combination")); }

void SelectedRandom::load(int* i, Real* probs, Random** rvx)
{
   probs[*i]=p; rvx[*i]=rv; (*i)++;
   delete this;
}

Real SelectedRandom::Next()
   { Throw(Logic_error("Newran: Next not defined")); return 0.0; }

Real AddedSelectedRandom::Next()
   { Throw(Logic_error("Newran: Next not defined")); return 0.0; }

Real RepeatedRandom::Next()
   { Real sum=0.0; for (int i=0; i<n; i++) sum += rv->Next(); return sum; }

MixedRandom::MixedRandom(int nx, Real* probs, Random** rvx)
{
   n = nx;
   rv = new Random*[n]; if (!rv) ErrorNoSpace();
   for (int i=0; i<n; i++) rv[i]=rvx[i];
   Build(probs);
}

MixedRandom::MixedRandom(AddedSelectedRandom& sr)
{
   n = sr.nelems();                       // number of terms;
   Real* probs = new Real[n]; rv = new Random*[n];
   if (!probs || !rv) ErrorNoSpace();
   int i=0; sr.load(&i,probs,rv);
   Build(probs); delete [] probs;
}

void MixedRandom::Build(Real* probs)
{
   int i;
   dg=new DiscreteGen(n,probs);
   if (!dg) ErrorNoSpace();
   mean=0.0; var=0.0;
   for (i=0; i<n; i++) mean = mean + (rv[i])->Mean()*probs[i];
   for (i=0; i<n; i++)
   {
      ExtReal sigsq=(rv[i])->Variance();
      ExtReal mudif=(rv[i])->Mean()-mean;
      var = var + (sigsq+square(mudif))*probs[i];
   }

}

MixedRandom::~MixedRandom()
{
   for (int i=0; i<n; i++) rv[i]->tDelete();
   delete dg; delete [] rv;
}

Real MixedRandom::Next()
   { int i = (int)(dg->Next()); return (rv[i])->Next(); }

int AddedSelectedRandom::nelems() const
   { return rv1->nelems() + rv2->nelems(); }

void AddedSelectedRandom::load(int* i, Real* probs, Random** rvx)
{
   rv1->load(i, probs, rvx); rv2->load(i, probs, rvx);
   delete this;
}

AddedSelectedRandom& operator+(SelectedRandom& rv1,
   SelectedRandom& rv2)
{
   AddedSelectedRandom* r = new AddedSelectedRandom(rv1, rv2);
   if (!r) ErrorNoSpace(); return *r;
}

AddedSelectedRandom& operator+(AddedSelectedRandom& rv1,
   SelectedRandom& rv2)
{
   AddedSelectedRandom* r = new AddedSelectedRandom(rv1, rv2);
   if (!r) ErrorNoSpace(); return *r;
}

AddedSelectedRandom& operator+(SelectedRandom& rv1,
   AddedSelectedRandom& rv2)
{
   AddedSelectedRandom* r = new AddedSelectedRandom(rv1, rv2);
   if (!r) ErrorNoSpace(); return *r;
}

AddedSelectedRandom& operator+(AddedSelectedRandom& rv1,
   AddedSelectedRandom& rv2)
{
   AddedSelectedRandom* r = new AddedSelectedRandom(rv1, rv2);
   if (!r) ErrorNoSpace(); return *r;
}

SelectedRandom& Random::operator()(double p)
{
   SelectedRandom* r = new SelectedRandom(*this, p);
   if (!r) ErrorNoSpace(); return *r;
}




// Identification routines for each class - may not work on all compilers?

char* Random::Name()            { return "Random";           }
char* Uniform::Name()           { return "Uniform";          }
char* Constant::Name()          { return "Constant";         }
char* PosGen::Name()            { return "PosGen";           }
char* SymGen::Name()            { return "SymGen";           }
char* AsymGen::Name()           { return "AsymGen";          }
char* PosGenX::Name()           { return "PosGenX";          }
char* SymGenX::Name()           { return "SymGenX";          }
char* AsymGenX::Name()          { return "AsymGenX";         }
char* Normal::Name()            { return "Normal";           }
char* ChiSq::Name()             { return "ChiSq";            }
char* Cauchy::Name()            { return "Cauchy";           }
char* Exponential::Name()       { return "Exponential";      }
char* Poisson::Name()           { return "Poisson";          }
char* Binomial::Name()          { return "Binomial";         }
char* NegativeBinomial::Name()  { return "NegativeBinomial"; }
char* Gamma::Name()             { return "Gamma";            }
char* Pareto::Name()            { return "Pareto";           }
char* DiscreteGen::Name()       { return "DiscreteGen";      }
char* SumRandom::Name()         { return "SumRandom";        }
char* MixedRandom::Name()       { return "MixedRandom";      }


// ********************** permutation generator ***************************

void RandomPermutation::Next(int N, int M, int p[], int start)
{
   // N = size of urn; M = number of draws
   if (N < M) Throw(Logic_error("Newran: N < M in RandomPermutation"));
   int i;
   int* q = new int [N];
   if (!q) ErrorNoSpace();
   for (i = 0; i < N; i++) q[i] = start + i;
   for (i = 0; i < M; i++)
   {
      int k = i + (int)(U.Next() * (N - i));       // uniform on i ... N-1
      p[i] = q[k]; q[k] = q[i];                    // swap i and k terms
                                                   // but put i-th term into p
   }
   delete [] q;
}

// from Sedgewick
void ShellSortAscending(int* a, int N)
{
   int h;
   for (h = 1; h <= N / 9; h = 3*h + 1) {}
   for (; h > 0; h /= 3) for (int i = h; i < N; ++i)
   {
      int v = a[i]; int j = i;
      while (j>=h && a[j-h] > v) { a[j] = a[j-h]; j-= h; }
      a[j] = v;
   }
}

void RandomCombination::SortAscending(int n, int gm[])
{
   ShellSortAscending(gm, n);
}

//***************** Generators with variable parameters ********************

VariPoisson::VariPoisson() : P100(100.0), P200(200.0) {}

int VariPoisson::iNext_very_small(Real mu)
{
   // mu is expected value of Poisson random number
   // Efficient when mu is small
   // generate a Poisson variable with mean mu
   if (mu == 0) return 0;
   Real t = exp(-mu); int k = 0; Real u = U.Next();
   for (Real s = t; s < u; s += t) { ++k; t *= mu / k; }
   return k;
}

int VariPoisson::iNext_small(Real mu)
{
   // mu is expected value of Poisson random number
   // Efficient when mu is reasonably small
   // generate a Poisson variable with mean mu
   // start iteration at mode
   if (mu < 10) return iNext_very_small(mu);
   int k_start = (int)mu;  Real u = U.Next();
   Real t1 = exp(k_start * log(mu) - mu - ln_gamma(k_start+1));
   if (t1 > u) return k_start;
   int k1 = k_start; int k2 = k_start; Real t2 = t1; Real s = t1;
   for(;;)
   {
      ++k1; t1 *= mu / k1; s += t1; if (s > u) return k1;
      if (k2 > 0) { t2 *= k2 / mu; --k2; s += t2; if (s > u) return k2; }
   }
}

int VariPoisson::iNext_large(Real mu)
{
   // mu is expected value of Poisson random number
   // reasonably accurate when mu is large, but should try to improve
   // generate a Poisson variable with mean mu
   // method is to start with normal variable X and use X with prob 1/3
   // and X**2 with prob 2/3. In each case X has mean and variance to
   // agree with Poisson rv after adjusting with Sheppard's correction

   const Real sc = 1.0 / 12.0;        // Sheppard correction
   int k;
   if (U.Next() > 1.0 / 3.0)
   {
      Real musc = 0.5 * (mu - sc); Real meansq = sqrt(mu * mu - musc);
      Real sigma = sqrt(musc / (mu + meansq));
      k = (int)floor(square(N.Next() * sigma + sqrt(meansq)) + 0.5);
   }
   else k = (int)floor(N.Next() * sqrt(mu - sc) + mu + 0.5);
   return k < 0 ? 0 : k;
}

int VariPoisson::iNext(Real mu)
{
   if (mu <= 0)
   {
      if (mu == 0) return 0;
      else Throw(Logic_error("Newran: illegal parameters"));
   }
   if (mu < 10) return iNext_very_small(mu);
   else if (mu < 100) return iNext_small(mu);
   else if (mu < 200)
   {
      // do in two statements so order of evaluation is preserved
      int i = (int)P100.Next();
      return i + iNext_small(mu - 100);
   }
   else if (mu < 300)
   {
      // do in two statements so order of evaluation is preserved
      int i = (int)P200.Next();
      return i + iNext_small(mu - 200);
   }
   else return iNext_large(mu);
}

VariBinomial::VariBinomial() {}

int VariBinomial::iNext_very_small(int n, Real p)
{
   // Efficient when n * p is small
   // generate a Binomial variable with parameters n and p

   Real q = 1 - p; if (q == 0) return n;
   Real r = p / q; Real t = pow(q, n); Real s = t;
   Real u = U.Next();
   for (int k = 0; k <= n; ++k)
   {
      if (s >= u) return k;
      t *= r * (Real)(n - k) / (Real)(k + 1);
      s += t;
   }
   return 0;    // can happen only if we have round-off error
}

int VariBinomial::iNext_small(int n, Real p)
{
   // Efficient when sqrt(n) * p is small
   // Assume p <= 1/2
   // same as small but start at mode
   // generate a Binomial variable with parameters n and p

   Real q = 1 - p; Real r = p / q; Real u = U.Next();
   int k_start = (int)( (n * r) / (r + 1) );
   Real t1 = exp( ln_gamma(n+1) - ln_gamma(k_start+1) - ln_gamma(n-k_start+1)
      + k_start * log(p) + (n-k_start) * log(q) );
   if (t1 >= u) return k_start;
   Real t2 = t1; Real s = t1;
   int k1 = k_start; int k2 = k_start;
   for (;;)
   {
      t1 *= r * (Real)(n - k1) / (Real)(k1 + 1); ++k1; s += t1;
      if (s >= u) return k1;
      if (k2)
      {
         --k2; t2 *= (Real)(k2 + 1) / (Real)(n - k2) / r; s += t2;
         if (s >= u) return k2;
      }
      else if (k1 == n) return 0;  // can happen only if we have round-off error
   }
}

int VariBinomial::iNext(int n, Real p)
{
   if (p > 0.5) return n - iNext(n, 1.0 - p);

   if (n <= 0 || p <= 0)
   {
      if (n < 0 || p < 0) Throw(Logic_error("Newran: illegal parameters"));
      else return 0;
   }

   Real mean = n * p;
   if (mean <= 10) return iNext_very_small(n, p);
   else if (mean <= 200) return iNext_small(n, p);
   else return iNext_large(n,p);

}

// probably not sufficiently accurate
int VariBinomial::iNext_large(int n, Real p)
{
   const Real sc = 1.0 / 12.0;        // Sheppard correction
   Real mean = n * p;
   Real sd = sqrt(mean * (1.0 - p) - sc);
   return (int)floor(N.Next() * sd + mean + 0.5);
}

Real VariLogNormal::Next(Real mean, Real sd)
{
   // should have special version of log for small sd/mean
   Real n_var = log(1 + square(sd / mean));
   return mean * exp(N.Next() * sqrt(n_var) - 0.5 * n_var);
}






#ifdef use_namespace
}
#endif