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


/**
 * This class implements inverse gaussian random variate generators using
 * the many-to-one transformation method of Michael, Schucany and Haas (MHS).
 * 
 */
public class InverseGaussianMSHGen extends InverseGaussianGen  {
   private NormalGen genN;


   /**
    * Creates an inverse gaussian random variate generator with
    *  parameters μ = mu and λ = lambda,
    *  using streams s and sn.
    * 
    */
   public InverseGaussianMSHGen (RandomStream s, NormalGen sn,
                                 double mu, double lambda)  {
      super (s, null);
      setParams (mu, lambda, sn);
   }

 
   /**
    * Creates a new generator for the distribution dist
    *  using streams s and sn.
    * 
    */
   public InverseGaussianMSHGen (RandomStream s, NormalGen sn,
                                 InverseGaussianDist dist)  {
      super (s, dist);
      if (dist != null)
         setParams (dist.getMu(), dist.getLambda(), sn);
   }


   /**
    * Generates a new variate from the inverse gaussian
    * distribution with parameters μ = mu and λ =
    * lambda, using streams s and sn.
    * 
    */
   public static double nextDouble (RandomStream s, NormalGen sn,
                                    double mu, double lambda)  {
      return mhs(s, sn, mu, lambda);
   }
 

   public double nextDouble() {
      return mhs(stream, genN, mu, lambda);
   }


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

   private static double mhs (RandomStream stream, NormalGen genN,
                              double mu, double lambda) {
      if (lambda <= 0.0)
         throw new IllegalArgumentException ("lambda <= 0");
      if (mu <= 0.0)
         throw new IllegalArgumentException ("mu <= 0");

      double z = genN.nextDouble ();
      double ymu = mu*z*z;
      double x1 = mu + 0.5*mu*ymu/lambda - 0.5*mu/lambda *
                  Math.sqrt(4.0*ymu*lambda + ymu*ymu);
      double u = stream.nextDouble();
      if (u <= mu/(mu + x1))
         return x1;
      return mu*mu/x1;
   }


   protected void setParams (double mu, double lambda, NormalGen sn) {
      setParams (mu, lambda);
      this.genN = sn;
   }

}