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


/**
 * This class implements symmetrical beta random variate generators using
 * Ulrich's polar method. The method generates two uniform
 * random variables 
 * x∈[0, 1] and 
 * y∈[- 1, 1] until
 *  
 * x2 + y2 <= 1. Then it returns
 * 

 * 
 * 

 * (1/2) + (xy(1 - S^2/(2α- 1))1/2)/S
 * 

 * where 
 * S = x2 + y2, and α is the shape parameter of the beta
 * distribution. The method is valid only when 
 * α > 1/2.
 * 
 */
public class BetaSymmetricalPolarGen extends BetaSymmetricalGen  {
   private double afactor;      // = 2/(2*alpha - 1)
   private RandomStream stream2;



   /**
    * Creates a symmetrical beta random variate generator with
    *  parameter α = alpha,  using stream s1 to generate x
    *   and stream s2 to generate y, as in above. Restriction: 
    * α > 1/2.
    * 
    */
   public BetaSymmetricalPolarGen (RandomStream s1, RandomStream s2,
                                   double alpha)  {
      super (s1, null);
      stream2 = s2;
      if (alpha <= 0.5)
         throw new IllegalArgumentException ("  must have alpha > 1/2");
      afactor = 2.0/(2.0*alpha - 1.0);
      setParams (alpha, alpha, 0.0, 1.0);
   }


   /**
    * Creates a symmetrical beta random variate generator with
    *  parameter α = alpha,  using stream s1 to generate x
    *  and y, as in above.  Restriction: 
    * α > 1/2.
    * 
    */
   public BetaSymmetricalPolarGen (RandomStream s1, double alpha)  {
      this (s1, s1, alpha);
   }


   /**
    * Creates a new generator for the distribution dist, using
    *    stream s1 to generate x and stream s2 to generate y,
    *    as in above.
    *     Restriction: dist must have 
    * α > 1/2.
    * 
    */
   public BetaSymmetricalPolarGen (RandomStream s1, RandomStream s2,
                                   BetaSymmetricalDist dist)  {
      super (s1, dist);
      stream2 = s2;
      double alp = dist.getAlpha();
      if (alp <= 0.5)
         throw new IllegalArgumentException ("  must have alpha > 1/2");
      afactor = 2.0/(2.0*dist.getAlpha() - 1.0);
      setParams (alp, alp, 0.0, 1.0);
   }


   /**
    * Creates a new generator for the distribution dist,
    *      using only one stream s1.
    *     Restriction: dist must have 
    * α > 1/2.
    * 
    */
   public BetaSymmetricalPolarGen (RandomStream s1,
                                   BetaSymmetricalDist dist)  {
      this (s1, s1, dist);
   }


   /**
    * Generates a random number using Ulrich's polar method. Stream
    *   s1 generates x and stream s2  generates y
    *    [see eq.].
    *   Restriction:  
    * α > 1/2.
    * 
    */
   public static double nextDouble (RandomStream s1, RandomStream s2,
                                    double alpha)   {
      double u, v, S;
      do {
         u = s1.nextDouble();
         v = -1.0 + 2.0*s2.nextDouble();
         S = u*u + v*v;
      } while (S > 1.0);
      return 0.5 + u*v/S* Math.sqrt(1.0 - Math.pow(S, 2.0/(2.0*alpha - 1.0)));
   }


   /**
    * Generates a random number by Ulrich's polar method using
    *   stream s.  Restriction:  
    * α > 1/2.
    * 
    */
   public static double nextDouble (RandomStream s, double alpha)   {
      return nextDouble (s, s, alpha);
   }
 

   public double nextDouble() {
      // Generates a random number using Ulrich's polar method.
      double u, v, S;
      do {
         u = stream.nextDouble();
         v = -1.0 + 2.0*stream2.nextDouble();
         S = u*u + v*v;
      } while (S > 1.0);
      return 0.5 + u*v/S* Math.sqrt(1.0 - Math.pow(S, afactor));
  }

   /**
    * Returns stream s2 associated with this object.
    * 
    */
   public RandomStream getStream2() {
      return stream2;
   }

}