package umontreal.iro.lecuyer.randvar;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.util.Num;


/**
 * This class implements random variate generators for the (discrete)
 * logarithmic distribution. Its  mass function is
 * 

 * 
 * 

 * p(x) = -θx/x log(1-θ)         for x = 1, 2,…,
 * 

 * where 
 * 0 < θ < 1.
 * It uses inversion with the LS chop-down algorithm if 
 * θ < θ0
 * and the LK transformation algorithm if 
 * θ >= θ0,
 * as described in.
 * The threshold θ0 can be specified when invoking the constructor.
 * Its default value is 
 * θ0 = 0.96, as suggested in.
 * 
 * 

 * A local copy of the parameter θ is maintained in this class.
 * 
 */
public class LogarithmicGen extends RandomVariateGenInt  {
   private static final double default_theta_limit = 0.96;

   private double theta_limit = default_theta_limit;
   private double theta;
   private double t;      // = log (1.0-theta).
   private double h;      // = -theta/log (1.0-theta).



   /**
    * Creates a logarithmic random variate generator with parameters
    *   θ = theta and default value 
    * θ0 = 0.96,
    *   using stream s.
    * 
    */
   public LogarithmicGen (RandomStream s, double theta) {
      this (s, theta, default_theta_limit);
   }


   /**
    * Creates a logarithmic random variate generator with parameters
    *   θ = theta and 
    * θ0 = theta0,
    *   using stream s.
    * 
    */
   public LogarithmicGen (RandomStream s, double theta, double theta0) {
      super (s, null);
      this.theta = theta;
      theta_limit = theta0;
      init();
   }


   /**
    * Creates a new generator with distribution dist and
    *   stream s, with default value 
    * θ0 = 0.96.
    * 
    */
   public LogarithmicGen (RandomStream s, LogarithmicDist dist) {
      this (s, dist, default_theta_limit);
   }


   /**
    * Creates a new generator with distribution dist
    *    and stream s, with 
    * θ0 = theta0.
    * 
    */
   public LogarithmicGen (RandomStream s, LogarithmicDist dist,
                          double theta0) {
      super (s, dist);
      theta_limit = theta0;
      if (dist != null)
         theta = dist.getTheta();
      init();
   }


   private void init () {
      if (theta <= 0.0 || theta >= 1.0)
         throw new IllegalArgumentException ("theta not in (0, 1)");
      if (theta >= theta_limit)
         h = Math.log1p(-theta);
      else
         t = -theta / Math.log1p(-theta);
   }

   public int nextInt() {
      if (theta < theta_limit)
         return ls (stream, theta, t);
      else   // Transformation
         return lk (stream, theta, h);
   }

   /**
    * Uses stream s to generate
    *    a new variate from the logarithmic distribution with parameter
    *  θ = theta.
    * 
    */
   public static int nextInt (RandomStream s, double theta) {
      if (theta < default_theta_limit)
         return ls (s, theta, -theta/Math.log1p(-theta));
      else   // Transformation
         return lk (s, theta, Math.log1p(-theta));
   }



//>>>>>>>>>>>>>>>>>>>>  P R I V A T E    M E T H O D S   <<<<<<<<<<<<<<<<<<<<


   private static int ls (RandomStream stream, double theta, double t) {
      double u = stream.nextDouble();
      int x = 1;

      double p =  t;

      while (u > p) {
            u -= p;
            x++;
            p *= theta*((double) x - 1.0)/((double)x);
      }
      return x;
   }

   private static int lk (RandomStream stream, double theta, double h) {
      double u, v, p, q;
      int x;

      u = stream.nextDouble();
      if (u > theta)
            return 1;
      v = stream.nextDouble();
      q = 1.0 - Math.exp(v * h);
      if ( u <= q * q) {
           x = (int)(1. + (Math.log(u) / Math.log(q)));
           return x;
      }
      return ((u > q) ? 1 : 2);
   }

   /**
    * Returns the θ associated with this object.
    * 
    */
   public double getTheta() {
      return theta;
   }


   /**
    * Returns the θ0 associated with this object.
    * 
    */
   public double getTheta0() {
      return theta_limit;
   }

}