/* *
 * Title:        

 * Description:  

 * Copyright:    Copyright (c) 

 * Company:      

 * @author
 * @version 1.0
 */

package umontreal.iro.lecuyer.probdist;


/**
 * Provides a variant of {@link ChiSquareDist} with
 * faster but less accurate methods.
 * The non-static version of 
 * inverseF calls the static version.
 * This method is not very accurate for small n but becomes
 * better as n increases.
 * The other methods are the same as in {@link ChiSquareDist}.
 * 
 */
public class ChiSquareDistQuick extends ChiSquareDist {



   /**
    * Constructs a chi-square distribution with n degrees of freedom.
    * 
    */
   public ChiSquareDistQuick (int n) {
      super (n);
   }


   public double inverseF (double u) {
      return inverseF (n, u);
   }

   /**
    * Computes a quick-and-dirty approximation of F-1(u), 
    *   where F is the chi-square distribution with n degrees of freedom.
    *   Uses the approximation given in  Figure L.24 of  
    *   Bratley, Fox and Schrage (1987) over most of the range.
    *   For u < 0.02 or u > 0.98, it uses the approximation given in 
    *    Goldstein for n >= 10, and returns 
    *   2.0 *
    *   {@link GammaDist#inverseF(double,int,double) GammaDist.inverseF} (n/2, 6, u)
    *    for n < 10 in order to avoid
    *   the loss of precision of the above approximations. 
    *    When n >= 10 or 
    * 0.02 < u < 0.98,
    *   it is between 20 to 30 times faster than the same method in
    *   {@link ChiSquareDist} for n between 10 and 1000 and even faster
    *   for larger n. 
    * 
    * 

    * Note that the number d of decimal digits of precision
    *   generally increases with n. For n = 3, we only have d = 3 over
    *   most of the range. For n = 10, d = 5 except far in the tails where d = 3.
    *   For n = 100, one has more than d = 7 over most of the range and for
    *   n = 1000, at least d = 8.
    *   The cases n = 1 and n = 2 are exceptions, with precision of about d = 10.
    * 
    */
   public static double inverseF (int n, double u) {
      /*
       * Returns an approximation of the inverse of Chi square cdf
       * with n degrees of freedom.
       * As in Figure L.24 of P.Bratley, B.L.Fox, and L.E.Schrage.
       *         A Guide to Simulation Springer-Verlag,
       *         New York, second edition, 1987.
       */

      if (u < 0.0 || u > 1.0)
         throw new IllegalArgumentException ("u is not in [0,1]");
      if (u <= 0.0)
         return 0.0;
      if (u >= 1.0)
         return Double.POSITIVE_INFINITY;

      final double SQP5 = 0.70710678118654752440;
      final double DWARF = 0.1e-15;
      final double ULOW = 0.02;
      double Z, arg, v, ch, sqdf;

      if (n == 1) {
          Z = NormalDist.inverseF01 ((1.0 + u)/2.0);
          return Z*Z;

      } else if (n == 2) {
         arg = 1.0 - u;
         if (arg < DWARF)
            arg = DWARF;
         return -Math.log (arg)*2.0;

     } else if ((u > ULOW) && (u < 1.0 - ULOW)) {
        Z = NormalDist.inverseF01 (u);
        sqdf = Math.sqrt ((double)n);
        v = Z * Z;

        ch = -(((3753.0 * v + 4353.0) * v - 289517.0) * v -
           289717.0) * Z * SQP5 / 9185400;

        ch = ch / sqdf + (((12.0 * v - 243.0) * v - 923.0)
           * v + 1472.0) / 25515.0;

        ch = ch / sqdf + ((9.0 * v + 256.0) * v - 433.0)
           * Z * SQP5 / 4860;

        ch = ch / sqdf - ((6.0 * v + 14.0) * v - 32.0) / 405.0;
        ch = ch / sqdf + (v - 7.0) * Z * SQP5 / 9;
        ch = ch / sqdf + 2.0 * (v - 1.0) / 3.0;
        ch = ch / sqdf + Z / SQP5;
        return n * (ch / sqdf + 1.0);

     } else if (n >= 10) {
        Z = NormalDist.inverseF01 (u);
        v = Z * Z;
        double temp;
        temp = 1.0 / 3.0 + (-v + 3.0) / (162.0 * n) -
           (3.0 * v * v + 40.0 * v + 45.0) / (5832.0 * n * n) +
           (301.0 * v * v * v - 1519.0 * v * v - 32769.0 * v -
           79349.0) / (7873200.0 * n * n * n);
        temp *= Z * Math.sqrt (2.0 / n);

        ch = 1.0 - 2.0 / (9.0 * n) + (4.0 * v * v + 16.0 * v -
           28.0) / (1215.0 * n * n) + (8.0 * v * v * v + 720.0 * v * v +
           3216.0 * v + 2904.0) / (229635.0 * n * n * n) + temp;

        return n * ch * ch * ch;

    } else {
    // Note: this implementation is quite slow.
    // Since it is restricted to the tails, we could perhaps replace
    // this with some other approximation (series expansion ?)
        return 2.0*GammaDist.inverseF (n/2.0, 6, u);
    }
}

}