/***************************************************************************
 * Copyright (C) 1995, 1996, Jun Sun, jsun@junsun.net.                     *
 *                                                                         *
 ***************************************************************************
 */
/* implementatin of Gamma distribution */

#include <assert.h>
#include "random_variates.h"
#include <math.h>

// both gamma1 and gamma2 return gamma function values for variables in [0,1]
// gamma1 has error |e(x)| <= 5 x 10^-5, 
// and gamma2 has error |e(x)| <= 3 x 10^-7
double gamma1(double x)
{
  assert((x>0.0) && (x<=1.0));

  static const double a1 = -.5748646;
  static const double a2 = 0.9512363;
  static const double a3 = -.6998588;
  static const double a4 = 0.4245549;
  static const double a5 = -.1010678;

  double g;

  // gamma(x+1) = g, 0<=x<=1.0
  g = 1.0 + a1*x + a2*x*x + a3*x*x*x + a4*x*x*x*x + a5*x*x*x*x*x;
  
  // gamma(x+1) = x*gamma(x) => gamma(x) = gamma(x+1) / x, 0<x<=1.0
  g = g / x;

  return g;
}


double gamma2(double x)
{
  assert((x>0.0) && (x<=1.0));

  static const double a1 = -.577191652;
  static const double a2 = 0.988205891;
  static const double a3 = -.897056937;
  static const double a4 = 0.918206857;
  static const double a5 = -.756704078;
  static const double a6 = 0.482199394;
  static const double a7 = -.193527818;
  static const double a8 = 0.035868343;

  double x2 = x*x;
  double x4 = x2*x2;
  double x8 = x4*x4;

  double g;

  // gamma(x+1) = g, 0<=x<=1.0
  g = 1.0 + a1*x + a2*x2 + a3*x2*x + a4*x4 + a5*x4*x + a6*x4*x2 
      + a7*x4*x2*x + a8*x8;

  // gamma(x+1) = x*gamma(x) => gamma(x) = gamma(x+1) / x, 0<x<=1.0
  g = g / x;

  return g;
}


// gamma function 
double gamma(double x)		
{
  assert(x>0.0);

  if (x<=1.0) 
    return gamma2(x);
  else
    return (x-1.0)*gamma(x-1.0);
}

Gamma::Gamma(double alpha, double beta)
{ _alpha = alpha;
  _beta  = beta;
  expo   = new Exponential(1.0);
}

Gamma::Gamma(double *data, int num_data)
{
// need to be implemented
  mythrow("Gamma fitting data constructor is not implemented yet");
}

Gamma::~Gamma() { delete expo; }

double Gamma::density(double x)
{
  if (x<=0.0) 
    return 0.0;
  else
    return pow(_beta,-_alpha)*pow(x,_alpha-1.0)*exp(-x/_beta)/gamma(_alpha);
}

double Gamma::distribution(double x)
{
// to be implemented
mythrow("Gamma:: distribution is not implemented");
return -1;
}

const double theta = 4.5;
const double d     = 1.0 + log(theta);
const double ln4   = log(4.0);

double Gamma::sample()
{
  double Y=0;		// temp var for sample

  if (_alpha == 1.0)
    Y = expo->sample();
  else if (_alpha <1.0)
  { 
    double U1,U2,P;
    double b;

    b = (M_E + _alpha) / M_E;
    for(;;)
    {
      U1 = randp->sample();
      P  = b*U1;
      if (P<=1.0)
      { Y = pow(P,1.0/_alpha);
        U2 = randp->sample();
        if (U2 <= exp(-Y)) break;
      }
      else
      { Y = -log((b-P)/_alpha);
        U2 = randp->sample();
        if (U2 <= pow(Y,_alpha-1))  break;
      }
    }
  }
  else if (_alpha > 1.0)
  {
    double U1,U2,a,b,q,W,Z,V;

    a = 1.0 / sqrt(2.0*_alpha-1.0);
    b = _alpha - ln4;
    q = _alpha + 1.0/a;

    for(;;)
    {
      U1 = randp->sample();
      U2 = randp->sample();

      V = a*log(U1/(1.0-U1));
      Y = _alpha*exp(V);
      Z = U1*U1*U2;
      W = b + q*V - Y;
     
      if ((W+d-theta*Z) >= 0.0) break;
      if (W >= log(Z)) break;
    }
  }

  return _beta*Y;
}  


double Gamma::probability(double x1, double x2)
{
// to be implemented
mythrow("Gamma::probability is not implemented");
return -1;
}


double Gamma::inverse_distribution(double)
{
// to be implemented
mythrow("Gamma::inverse_distribution is not implemented");
return -1;
}