Logo
Distributed Genetic Programming Framework
print print

File org.sfc.math.Mathematics.java

Here you can find all the information about the file org.sfc.math.Mathematics.java. You may explore it here or download it onto your local disk.
/*
 * Copyright (c) 2005 Thomas Weise
 * 
 * E-Mail           : tweise@gmx.de
 * Creation Date    : 2005-06-17 08:05:49
 * Original Filename: org.sfc.math.Mathematics.java
 * Version          : 1.1.1
 * Last modification: 2006-03-07
 *                by: Thomas Weise
 * 
 * License          : GNU LESSER GENERAL PUBLIC LICENSE
 *                    Version 2.1, February 1999
 *                    You should have received a copy of this license along
 *                    with this library; if not, write to the Free Software
 *                    Foundation, Inc., 59 Temple Place, Suite 330, Boston,
 *                    MA 02111-1307, USA or download the license under
 *                    http://www.gnu.org/copyleft/lesser.html.
 *                    
 * Warranty         : This software is provided "as is" without any
 *                    warranty; without even the implied warranty of
 *                    merchantability or fitness for a particular purpose.
 *                    See the Gnu Lesser General Public License for more
 *                    details.
 */

 
package org.sfc.math;

import org.sfc.utils.Typesafe;



/**
 * This static class provides additional methods to the java
 * <code>Math</code> class.
 * 
 * TODO: Check "refine" for normal cdf
 * @see java.lang.Math
 * @author Thomas Weise 
 */

public final class Mathematics
  {
/**
 * The sqare root of the maximum integer.
 */

  private static  final int     SQRT_MAX_INT  =
                              ((int)(Math.sqrt(Integer.MAX_VALUE)));

/**
 * The sqare root of the maximum long.
 */

  private static  final long    SQRT_MAX_LONG =
                              ((long)(Math.sqrt(Long.MAX_VALUE)));
  
/**
 * The natural logarithm of 10.
 */

  public  static  final double  LN_10       = Math.log(10.0d);
  
/**
 * The natural logarithm of 2.
 */

  public  static  final double  LN_2        = Math.log(2.0d);
  
/**
 *  The smallest difference of two numbers, if they differ smaller,
 *  they are considered as equal.
 */

  public  static  final double  EPS         ;
  
/**
 * The natural logarithm of EP.
 */

  public  static  final double  LN_EPS      ;
  
/**
 * 1.0 / EPS
 * @see #EPS
 */

  public  static  final double  INV_EPS     ;
  
/**
 * The square root of pi.
 */

  public  static  final double  SQRT_PI     = Math.sqrt(Math.PI);
  
/**
 * The natural logarithm of the square root of pi.
 */

  public  static  final double  LN_SQRT_PI  = Math.log(SQRT_PI);

/**
 * 1.0 / square root of pi.
 */

  public  static  final double  INV_SQRT_PI = 1.0d / SQRT_PI;
  
/**
 * The square root of 2.0.
 */

  public  static  final double  SQRT_2      = Math.sqrt(2.0d);
/**
 * The square root of 2*pi.
 */

  public  static  final double  SQRT_2_PI   = SQRT_2 * SQRT_PI;
    
/**
 * Euler's (Mascheroni's) constant.
 */

  public  static  final double  EULER_CONSTANT =
    0.57721566490153286060651209008240243104215933593992359880576723488486772677766467093694706329174674951463144724980708248096050401448654283622417399764492353625350033374293733773767394279259525824709491600873520394816567d; 

/**
 * Exactly 1/3
 */

  public  static  final double  ONE_THIRD = (1.0d / 3.0d);

/**
 * Exactly 2/3
 */

  public  static  final double  TWO_THIRD = (1.0d - ONE_THIRD);
  
/**
 * The natural logarithm of Double.MAX_VALUE, or, in other words, the
 * biggest value where <code>Math.exp(x)</code> produces an exact result.
 * @see Double#MAX_VALUE
 */

  public  static  final double  LN_MAX_DOUBLE = Math.log(Double.MAX_VALUE);
/**
 * The natural logarithm of Double.MIN_VALUE, or, in other words, the
 * smallest value where <code>Math.exp(x)</code> produces an exact result.
 * @see Double#MIN_VALUE
 */

  public  static  final double  LN_MIN_DOUBLE = Math.log(Double.MIN_VALUE);
  
///
/// Constants used to calculate the gamma function.
///  
/**
 * The calculus constant GAMMA_A.
 */

  private static final double GAMMA_A =
    (Mathematics.SQRT_2_PI *   1.000000000190015d);
/**
 * The calculus constant GAMMA_B.
 */

  private static final double GAMMA_B =
    (Mathematics.SQRT_2_PI *  76.18009172947146d);
/**
 * The calculus constant GAMMA_C.
 */
 
  private static final double GAMMA_C =
    (Mathematics.SQRT_2_PI * -86.50532032941677d);
/**
 * The calculus constant GAMMA_D.
 */

  private static final double GAMMA_D =
    (Mathematics.SQRT_2_PI *  24.01409824083091d);
/**
 * The calculus constant GAMMA_E.
 */

  private static final double GAMMA_E =
    (Mathematics.SQRT_2_PI *  -1.231739572450155d);
/**
 * The calculus constant GAMMA_F.
 */

  private static final double GAMMA_F =
    (Mathematics.SQRT_2_PI *   0.1208650973866179E-2d);
/**
 * The calculus constant GAMMA_G.
 */

  private static final double GAMMA_G =
    (Mathematics.SQRT_2_PI *  -0.5395239384953E-5d);
/**
 * Another precision factor.
 */
  
  private static final double GAMMA_MIN_DP                   = 1e-30d;
  
///
/// Constants for the computation of the Inverse Normal Cummulative
/// Distribution.
///
  
/**
 * A helper value for the computation of the Inverse Normal Cummulative
 * Distribution.
 */

  private static final double ICDF_P_LOW  = 0.02425D;
/**
 * A helper value for the computation of the Inverse Normal Cummulative
 * Distribution.
 */

  private static final double ICDF_P_HIGH = 1.0D - ICDF_P_LOW;

/**
 * Coefficients A in rational approximations of the Inverse Normal
 * Cummulative Distribution.
 */

  private static final double ICDF_A[] =
  { -3.969683028665376e+01,  2.209460984245205e+02,
    -2.759285104469687e+02,  1.383577518672690e+02,
    -3.066479806614716e+01,  2.506628277459239e+00 };

/**
 * Coefficients B in rational approximations of the Inverse Normal
 * Cummulative Distribution.
 */

  private static final double ICDF_B[] =
  { -5.447609879822406e+01,  1.615858368580409e+02,
    -1.556989798598866e+02,  6.680131188771972e+01,
    -1.328068155288572e+01 };


/**
 * Coefficients C in rational approximations of the Inverse Normal
 * Cummulative Distribution.
 */

  private static final double ICDF_C[] =
  { -7.784894002430293e-03, -3.223964580411365e-01,
    -2.400758277161838e+00, -2.549732539343734e+00,
    4.374664141464968e+00,  2.938163982698783e+00 };



/**
 * Coefficients D in rational approximations of the Inverse Normal
 * Cummulative Distribution.
 */

  private static final double ICDF_D[] =
  { 7.784695709041462e-03,  3.224671290700398e-01,
    2.445134137142996e+00,  3.754408661907416e+00 };
  
///
/// Constants used for the erfc-functions
///
/**
 *  Coefficients for approximation to  erf  in first interval
*/

  private static final double ERF_A[] =
  { 3.16112374387056560E00, 1.13864154151050156E02,
    3.77485237685302021E02, 3.20937758913846947E03,
    1.85777706184603153E-1 };

/**
 *  Coefficients for approximation to  erf  in first interval
*/

  private static final double ERF_B[] =
  { 2.36012909523441209E01, 2.44024637934444173E02,
    1.28261652607737228E03, 2.84423683343917062E03 };

/**
 *  Coefficients for approximation to  erfc  in second interval
 */

  private static final double ERF_C[] =
  { 5.64188496988670089E-1, 8.88314979438837594E0,
    6.61191906371416295E01, 2.98635138197400131E02,
    8.81952221241769090E02, 1.71204761263407058E03,
    2.05107837782607147E03, 1.23033935479799725E03,
    2.15311535474403846E-8 };

/**
 *  Coefficients for approximation to  erfc  in second interval
 */

  private static final double ERF_D[] =
  { 1.57449261107098347E01,1.17693950891312499E02,
    5.37181101862009858E02,1.62138957456669019E03,
    3.29079923573345963E03,4.36261909014324716E03,
    3.43936767414372164E03,1.23033935480374942E03 };

/**
 *  Coefficients for approximation to  erfc  in third interval
 */

  private static final double ERF_P[] =
  { 3.05326634961232344E-1,3.60344899949804439E-1,
    1.25781726111229246E-1,1.60837851487422766E-2,
    6.58749161529837803E-4,1.63153871373020978E-2 };


/**
 *  Coefficients for approximation to  erfc  in third interval
 */

  private static final double ERF_Q[] =
  { 2.56852019228982242E00,1.87295284992346047E00,
    5.27905102951428412E-1,6.05183413124413191E-2,
    2.33520497626869185E-3 };

//
// Constants for the digamma function
//
/**
 * The helper constant for the digamme (psi) function.
 */

  private static  final double[]  DIGAMMA_KNOCE = {
     0.30459198558715155634315638246624251d,
     0.72037977439182833573548891941219706d,
    -0.12454959243861367729528855995001087d,
     0.27769457331927827002810119567456810e-1d,
    -0.67762371439822456447373550186163070e-2d,
     0.17238755142247705209823876688592170e-2d,
    -0.44817699064252933515310345718960928e-3d,
     0.11793660000155572716272710617753373e-3d,
    -0.31253894280980134452125172274246963e-4d,
     0.83173997012173283398932708991137488e-5d,
    -0.22191427643780045431149221890172210e-5d,
     0.59302266729329346291029599913617915e-6d,
    -0.15863051191470655433559920279603632e-6d,
     0.42459203983193603241777510648681429e-7d,
    -0.11369129616951114238848106591780146e-7d,
     0.304502217295931698401459168423403510e-8d,
    -0.81568455080753152802915013641723686e-9d,
     0.21852324749975455125936715817306383e-9d,
    -0.58546491441689515680751900276454407e-10d,
     0.15686348450871204869813586459513648e-10d,
    -0.42029496273143231373796179302482033e-11d,
     0.11261435719264907097227520956710754e-11d,
    -0.30174353636860279765375177200637590e-12d,
     0.80850955256389526647406571868193768e-13d,
    -0.21663779809421233144009565199997351e-13d,
     0.58047634271339391495076374966835526e-14d,
    -0.15553767189204733561108869588173845e-14d,
     0.41676108598040807753707828039353330e-15d,
    -0.11167065064221317094734023242188463e-15d };

/**
 * The erf threshold
 */

  private static final double ERF_THRESHOLD = 0.46875D;

/**
 * Hardware dependant constant for erf.
 */

  private static final double ERF_X_NEG = -9.38241396824444;
/**
 * Hardware dependant constant for erf.
 */

  private static final double ERF_X_BIG = 9.194E0;
/**
 * Hardware dependant constant for erf.
 */

  private static final double ERF_X_HUGE ;
/**
 * Hardware dependant constant for erf.
 */

  private static final double ERF_X_MAX =
    Math.min(Double.MAX_VALUE, (1.0d / (Math.sqrt(Math.PI) * Double.MIN_VALUE)));
  
///
/// Factorial table.
///
/**
 * The table with the factorials.
 */

  private static  final double[]  FACTORIALS      = new double[40];
/**
 * The table with the natural logarithms of the factorials.
 */

  private static  final double[]  FACTORIALS_LN   = new double[100];
  
    
///
/// Bessel constants.
///
/**
 * The bessel acc value.
 */

  private static  final double    BESSEL_ACC    = 40.0d;
/**
 * The bessel big value.
 */

  private static  final double    BESSEL_BIGNO  = 1.0e+10d;
/**
 * The bessek small value.
 */

  private static  final double    BESSEL_BIGNI  = 1.0e-10d;
  
  
///
/// Constants for the incomplete complemented gamma function.
///
/**
 * First helper.
 */

  private static  final double    ICG_BIG     ;
/**
 * Second helper.
 */

  private static  final double    ICG_BIGINV  ;
  
  
/**
 * The maximum value where gamma(x) returns an exact result.
 */

  public static  final double  MAX_GAMMA     ;
  
  
  static
    {
    double  l_d, l_s, l_m;
    int     l_i;   
    
//
// Calculate constants derived from machine accuracy.
//
    l_d = 1.0d;
    l_s = 0.5d;
    l_m = 0.0d;
    
    while( (l_d - (l_s * l_d)) != l_d)
      {
      l_m  = l_s;
      l_s *= 0.5d;
      }
    
    EPS         = l_m;
    LN_EPS      = Math.log(l_m);
    INV_EPS     = (1.0d / EPS);
    ICG_BIG     = (2.0d * EPS);
    ICG_BIGINV  = (1.0d / ICG_BIG);
    
//
// Fill the factorial table.
//
    l_d = 1.0d;
    for(l_i = 1; l_i < FACTORIALS.length; l_i++)
      {
      l_d = Math.rint(l_d * l_i);
      FACTORIALS[l_i]    = l_d;
      FACTORIALS_LN[l_i] = Math.log(l_i);
      }
    
    ERF_X_HUGE = (1.0d / (2.0d * Math.sqrt(EPS)));
    
    for(; l_i < FACTORIALS_LN.length; l_i++)
      {
      FACTORIALS_LN[l_i] = gamma_ln(l_i+1);
      }
    
    FACTORIALS_LN[0] = 1.0d;
    
    
//
// Determine the maximum value for the gamma function.
//
    l_d = 1.0;
    l_s = 1E3;
    
    while(l_d != l_s)
      {
      l_m = (l_d + l_s) * 0.5;
      if(Double.isInfinite(gamma(l_m)))
        {
        l_s = l_m - (EPS * l_m);
        }
      else
        {
        l_d = l_m + (EPS * l_m);
        }
      }
    
    MAX_GAMMA = l_d;
    }
  
/**
 * Prevent anyone from creatin an instance of <code>Mathematics</code>.
 */

  private  Mathematics ()
    {
    Typesafe.do_not_call();
    }
  
/**
 * Returns wether or not two double values are approximately equal.
 * @param p_p1    The first double to compare.
 * @param p_p2    The second double to compare.
 * @return  <code>true</code> if and only if the two values are
 *          approximately equal.
 */

  public  static  final boolean equals  (final  double  p_p1,
                                         final  double  p_p2)
    {
    if(Double.compare(p_p1, p_p2) == 0) return true;
    
    return ( Math.abs(p_p1 - p_p2) < Math.ulp(Math.max(p_p1, p_p2)) );
    }
  
/**
 * Divides p_a by p_b, but increases the result if
 * p_b doesn't fit exactly n times into p_a.
 * 
 * @param p_a   The number to divide by p_b
 * @param p_b   The number to divide p_a by
 * @return      (p_a / p_b) + (((p_a % p_b) == 0) ? 0 : 1)
 * @exception ArithmeticException Pf p_b == 0.
 */

  public  static  final long  ceil_div  (long p_a,
                                         long p_b)
          throws ArithmeticException
    {
    long  l_r;
    
    l_r = (p_a / p_b);
    if((p_a % p_b) != 0)
      {
      if(l_r < 0) return (l_r - 1);
      return (l_r + 1);
      }
    
    return l_r;
    }

/**
 * Divides p_a by p_b, but increases the result if
 * p_b doesn't fit exactly n times into p_a.
 * 
 * @param p_a   The number to divide by p_b
 * @param p_b   The number to divide p_a by
 * @return      (p_a / p_b) + (((p_a % p_b) == 0) ? 0 : 1)
 * @exception ArithmeticException If p_b == 0.
 */

  public  static  final int   ceil_div  (int  p_a,
                                         int  p_b)
          throws ArithmeticException
    {
    int l_r;
    
    l_r = (p_a / p_b);
    if((p_a % p_b) != 0)
      {
      if(l_r < 0) return (l_r - 1);
      return (l_r + 1);
      }
    
    return l_r;
    }
  
  
/**
 * Returns the highest bit that's set in p_d.
 * 
 * @param   p_d   The number to find the highest bit in.
 * @return        The index of the mvb, or -1 if p_d == 0.
 */

  public  static  final int     ld (int p_d)
    {
    int l_r;
    
    l_r = -1;
    while(p_d != 0)
      {
      l_r++;
      p_d >>= 1;
      }
    
    return l_r;
    }
  

  
/**
 * Returns the highest bit that's set in p_d.
 * 
 * @param   p_d   The number to find the highest bit in.
 * @return        The index of the mvb, or -1 if p_d == 0.
 */

  public  static  final int     ld (long p_d)
    {
    int l_r;
    
    l_r = -1;
    while(p_d != 0)
      {
      l_r++;
      p_d >>= 1;
      }
    
    return l_r;
    }
  
/**
 * Returns the count of digits of <code>p_d</code>.
 * 
 * @param   p_d   The number to find the count of digits of.
 * @return        The count of digits, or -1 if p_d == 0.
 */

  public  static  final int     lg (long p_d)
    {
    int l_r;
    
    l_r = 1;
    do
      {
      l_r++;
      p_d /= 10;
      } while(p_d != 0);
    
    return l_r;
    }
  

/**
 * Returns the count of digits of <code>p_d</code>.
 * 
 * @param   p_d   The number to find the count of digits of.
 * @return        The count of digits, or -1 if p_d == 0.
 */

  public  static  final int     lg (int p_d)
    {
    int l_r;
    
    l_r = 1;
    do
      {
      l_r++;
      p_d /= 10;
      } while(p_d != 0);
    
    return l_r;
    }
  
/**
 * Calculate the closest integer <= sqare root <code>p_d</code>.
 * @param p_d The input value.
 * @return The integer closest to <= sqare root <code>p_d</code>.
 */

  public  static  final int sqrt  (int p_d)
    {
    int l_i, l_m, l_l, l_j;
    
    if(p_d < 0) throw new ArithmeticException();
    
    l_m = SQRT_MAX_INT;
    l_l = 0;
    
    while(l_m > l_l)
      {
      l_i = ((l_m + l_l) >> 1);
      l_j = l_i * l_i;
            if(l_j > p_d) l_m = (l_i-1);
      else  if(l_j < p_d) l_l = (l_i+1);
      else  return l_i;      
      }
    
    return (((l_l*l_l) > p_d) ? (l_l-1) : l_l);
    }
  
/**
 * Calculate the closest long <= sqare root <code>p_d</code>.
 * @param p_d The input value.
 * @return The long closest to <= sqare root <code>p_d</code>.
 */

  public  static  final long sqrt  (long p_d)
    {
    long l_i, l_m, l_l, l_j;
    
    if(p_d < 0) throw new ArithmeticException();
    
    l_m = SQRT_MAX_LONG;
    l_l = 0;
    
    while(l_m > l_l)
      {
      l_i = ((l_m + l_l) >> 1);
      l_j = l_i * l_i;
            if(l_j > p_d) l_m = (l_i-1);
      else  if(l_j < p_d) l_l = (l_i+1);
      else  return l_i;      
      }
    
    return (((l_l*l_l) > p_d) ? (l_l-1) : l_l);
    }
  
/**
 * Returns the 10s logarithm of p_d.
 * 
 * @param   p_d   The number to take the logarithm of.
 * @return        log(p_d) / log(10)
 */

  public  static  final double  lg    (double p_d)
    {
    return Math.log10(p_d);
    }

/**
 * Returns 10<sup>p_d</sup>.
 * @param   p_d   The number to 10^p_d.
 * @return        exp(p_d * log(10))
 */

  public  static  final double  exp_10    (double p_d)
    {
    return  Math.exp(p_d * LN_10);
    }

/**
 * Returns the dual logarithm of p_d.
 * 
 * @param   p_d   The number to take the dual logarithm of.
 * @return        log(p_d) / log(2)
 */

  public  static  final double  ld   (final double p_d)
    {
    return Math.log(p_d) / LN_2;
    }

  
/**
 * Returns the 2^p_d.
 * 
 * @param   p_d   The number for 2^p_d.
 * @return        exp(p_d * log(2))
 */

  public  static  final double  exp_2   (double p_d)
    {
    return  Math.exp(p_d * LN_2);
    }

/**
 * Rounds p_r to p_decimals decimals.
 * 
 * @param   p_r         The number to round.
 * @param   p_decimals  The count of decimals to round p_r to.
 *                      If negative, maybe -1, p_r would be rounded to
 *                      full decades.
 * @return              p_r rounded to p_decimals decimals.
 */

  public  static  final double  round   (double p_r,
                                         int    p_decimals)
    {
    double  l_z;
    
    if(p_r == 0.0d) return p_r;

    l_z                = Math.rint(Math.log10(p_r));

    if(l_z < 0.0d) l_z = (p_decimals - l_z - 1);
    else           l_z = p_decimals;

    if(l_z < 0.0) l_z = 1.0d / Math.rint(exp_10(Math.rint(-l_z)));
    else          l_z = Math.rint(exp_10(Math.rint(l_z)));

    return (Math.rint(p_r * l_z) / l_z);
    }
  
/**
 * Calculates the sigmoid function.
 * 
 * @param p_d   The number to calculate the sigmoid function of.
 * 
 * @return  sigmoid(p_d) = 1 / (1 + exp(p_d). 
 */

  public  static  final double  sigmoid (double p_d)
    {
    return (1.0d / (1.0d + Math.exp(-p_d)));
    }
  
/**
 * This routine calculates sqrt(a^2 + b^2) 
 * without under- or overflow.
 * 
 * @param   p_a   first real number
 * @param   p_b   second real number
 * @return        sqrt(a<sup>2</sup> + b<sup>2</sup>)
 */

  public  static  final double  hypot     (double p_a,
                                           double p_b)
    {
    return Math.hypot(p_a, p_b);    
    }
  
/**
 * Returns the cosinus hyperbolicus of p_d.
 * 
 * @param   p_d   The number to take the cosh of.
 * @return        0.5 (e<sup>x</sup> + x<sup>-x</sup>)
 */

  public  static  final double  cosh   (double p_d)
    {
    return Math.cosh(p_d);
    }
  
/**
 * Returns the tangens hyperbolicus of p_d.
 * 
 * @param   p_d   The number to take the tanh of.
 * @return        (e<sup>x</sup> - x<sup>-x</sup>) / (e<sup>x</sup> + x<sup>-x</sup>)
 */

  public  static  final double  tanh   (double p_d)
    {
    return Math.tan(p_d);
    }
  
/**
 * Compute the arcus cosinus hyperbolicus.
 * @param p_x a double value
 * @return The hyperbolic arc cosine of the argument.
 */


  public static  final double acosh(double p_x)
    {
    return Math.log(p_x + Math.sqrt((p_x*p_x)-1.0));
    }

/**
 * Compute the arcus sinus hyperbolicus.
 * @param p_x a double value
 * @return The hyperbolic arc sine of the argument.
 */

  public static final double asinh(double p_x)
    {
    double l_x;
    int    l_sign;
      
    if(p_x == 0.0d) return p_x;
    
    if(p_x < 0.0d)
      {
      l_sign = -1;
      l_x    = -p_x;
      }
    else
      {
      l_sign = 1;
      l_x    = p_x;
      }
     
    return l_sign*Math.log(l_x + Math.sqrt((p_x*p_x)+1));
    }

/**
 * Compute the arcus tangens hyperbolicus.
 * @param p_x a double value
 * @return The hyperbolic arc tangent of the argument.
 */

  public static final double atanh(double p_x)
    {
    return 0.5d * Math.log( (1.0d+p_x)/(1.0d-p_x) );
    }
  
/**
 * Returns the signum of p_d.
 * 
 * @param   p_d   The number to get the sign of.
 * @return        -1.0 if p_d < 0.0
 *                 1.0 if p_d > 0.0
 *                 0.0 else                
 */

  public  static  final double  signum   (double p_d)
    {
    return Math.signum(p_d);
    }
  
/**
 * Returns the gamma function at <code>p_z</code>.
 * @param p_z a real number
 * @return <code> &Gamma;(p_z) = &int;<sub>0</sub><sup>&infin;</sup> t<sup>p_z-1</sup>e<sup>-t</sup>dt</code>
 */

  public static final double gamma(double p_z)
    {
    double  l_d;
    boolean l_b;
    
    l_b = (p_z <= 0.5d);    
    l_d = (l_b ? (1.0d - p_z) : p_z);    
    l_d = Math.exp(Math.log(l_d + 5.5d) * (l_d + 0.5d) - (l_d + 5.5d)) *
                  ((GAMMA_A + (GAMMA_B/(l_d+1.0d)) +
                   (GAMMA_C/(l_d+2.0d)) +
                   (GAMMA_D/(l_d+3.0d)) + 
                   (GAMMA_E/(l_d+4.0d)) +
                   (GAMMA_F/(l_d+5.0d)) +
                   (GAMMA_G/(l_d+6.0d))) /
                        l_d);
    if(l_b)
      {
      l_d = ((Math.PI / l_d) / Math.sin(Math.PI * p_z)); 
      }
    
    if(Math.abs(p_z - Math.rint(p_z)) <= Mathematics.EPS)
      {
      return Math.rint(l_d);
      }
    return l_d;
    }

  
/**
 * Returns the logarithm of the gamma function at <code>z</code>.
 * @param p_z a real number
 * @return <code>log( &Gamma;(p_z) ) = log ( &int;<sub>0</sub><sup>&infin;</sup> t<sup>p_z-1</sup>e<sup>-t</sup>dt )</code>
 */

  public static final double gamma_ln(double p_z )
    {
    double  l_d;
    boolean l_b;
    
    l_b = (p_z <= 0.5d);    
    l_d = (l_b ? (1.0d - p_z) : p_z);    
    
    l_d = ((Math.log(l_d + 5.5d) * (l_d + 0.5d) - (l_d + 5.5d)) +     
            Math.log( (GAMMA_A + (GAMMA_B/(l_d+1.0d)) +
                (GAMMA_C/(l_d+2.0d)) + (GAMMA_D/(l_d+3.0d)) + 
                (GAMMA_E/(l_d+4.0d)) + (GAMMA_F/(l_d+5.0d)) +
                (GAMMA_G/(l_d+6.0d))  ) / l_d ));
    if(l_b)
      {
      return (Math.log(Math.PI / Math.sin(Math.PI * p_z) ) - l_d);
      }
    return l_d;
    }

/**
 * Returns the lower incomplete gamma function at <code>p_a, p_x</code>.
 * @param p_a a positive real number
 * @param p_x a real number greater than <code>p_a</code>
 * @return The lower incomplete gamma function at <code>p_a, p_x</code>.
 */

  public final  static double gamma_lower_incomplete( double p_a,
                                                      double p_x )
    {
    return gamma(p_a) - gamma_upper_incomplete(p_a, p_x);
    }

/**
 * Returns the upper incomplete gamma function at <code>p_a, p_x</code>.
 * @param p_a a positive real number
 * @param p_x a real number greater than <code>p_a</code>
 * @return <code>&Gamma;( p_a, p_x ) = &int;<sub>0</sub><sup>&infin;</sup> t<sup>p_z-1</sup>e<sup>-t</sup>dt</code>
 */

  public final  static double gamma_upper_incomplete( double p_a,
                                                      double p_x )
    {
    double  l_b, l_delta, l_od, l_c, l_d, l_h, l_an, l_i;
   
    if(p_x < (p_a + 1.0d))
      {
      l_d           = gamma(p_a);
      l_b           = p_a;
      l_delta       = 1.0d / p_a;
      l_od          = -l_delta;
      l_c           = l_delta;  
  
      while(l_delta != l_od) // (Math.abs(l_delta) > Math.abs(l_sum) * EPS) )
        {
        l_od     = l_delta;
        l_b      = (1.0d + l_b);
        l_delta *= (p_x/l_b);
        l_c     += l_delta;      
        }
  
      return l_d  - (l_c * Math.exp(p_a*Math.log(p_x) - p_x)); 
      }
    
    l_b     = ((p_x + 1.0d) - p_a);
    l_c     = (1.0d/GAMMA_MIN_DP);
    l_d     = (1.0d/l_b);
    l_h     = l_d;
    l_delta = 1.0d; 
    l_i     = 1.0d;
    
    do
      {
      l_od = l_h;
      l_an = (-l_i * (l_i-p_a));
      l_i  = (l_i + 1.0d);
      l_b += 2.0d;
      l_d  = ((l_an*l_d) + l_b);
      l_c  = ((l_an/l_c) + l_b);

      if(Math.abs(l_c) < GAMMA_MIN_DP ) l_c = GAMMA_MIN_DP;    
      if(Math.abs(l_d) < GAMMA_MIN_DP ) l_d = GAMMA_MIN_DP;

      l_d     = (1.0d / l_d);
      l_delta = (l_c * l_d);
      l_h    *= l_delta;
      
      } while(l_od != l_h);
    
    return l_h * Math.exp(p_a * Math.log(p_x) - p_x);
    }
  
/**
 * Compute the regularized q gamma function.
 * @param p_a double value
 * @param p_x double value
 * @return The regularized q gamma function. 
 */

  public  static  final double gamma_regularized_q(double   p_a,
                                                   double   p_x)
    {
    double l_ans, l_ax, l_c, l_yc, l_r, l_t, l_y, l_z, l_pk, l_pkm1,
           l_pkm2, l_qk, l_qkm1, l_qkm2, l_ot;

    if((p_x <= 0.0d) || (p_a <= 0.0d) )
      {
      return 1.0d;
      }

    if((p_x < 1.0d) || (p_x < p_a))
      {
      return 1.0d - gamma_regularized_p(p_a, p_x);
      }

    l_ax = p_a * Math.log(p_x) - p_x - gamma_ln(p_a);
    if(l_ax < -LN_MAX_DOUBLE ) return 0.0d;

    l_ax    = Math.exp(l_ax);
    l_y     = (1.0d - p_a);
    l_z     = p_x + l_y + 1.0d;
    l_c     = 0.0d;
    l_pkm2  = 1.0d;
    l_qkm2  = p_x;
    l_pkm1  = p_x + 1.0d;
    l_qkm1  = l_z * p_x;
    l_ans   = l_pkm1/l_qkm1;
    l_t     = Double.NaN;
    
    do
      {
      l_ot = l_t;
      l_c += 1.0;
      l_y += 1.0;
      l_z += 2.0;
      l_yc = l_y * l_c;
      l_pk = (l_pkm1 * l_z) - (l_pkm2 * l_yc);
      l_qk = (l_qkm1 * l_z) - (l_qkm2 * l_yc);
      
      if(l_qk != 0.0d)
        {
        l_r   = l_pk/l_qk;
        l_t   = Math.abs( (l_ans - l_r)/l_r );
        l_ans = l_r;
        }
      else
        {
        l_t = 1.0d;
        }

      l_pkm2 = l_pkm1;
      l_pkm1 = l_pk;
      l_qkm2 = l_qkm1;
      l_qkm1 = l_qk;
      
      if( Math.abs(l_pk) > ICG_BIG )
        {
        l_pkm2 *= ICG_BIGINV;
        l_pkm1 *= ICG_BIGINV;
        l_qkm2 *= ICG_BIGINV;
        l_qkm1 *= ICG_BIGINV;
        }
      } while( l_ot != l_t );

     return l_ans * l_ax;
     }
  
/**
 * Computes the regularized p gamma function.
 * @param p_a double value
 * @param p_x double value
 * @return The regularized p gamma function.
 */

  public final  static double gamma_regularized_p (double   p_a,
                                                   double   p_x)
    {
    double l_ans, l_ax, l_c, l_r;

    if((p_x <= 0.0d) || (p_a <= 0.0d))
      {
      return 0.0d;
      }

    if((p_x > 1.0d) && (p_x > p_a))
      {
      return 1.0d - Mathematics.gamma_regularized_q(p_a, p_x);
      }

    l_ax = p_a * Math.log(p_x) - p_x - gamma_ln(p_a);
    if(l_ax < -LN_MAX_DOUBLE) return 0.0d;

    l_ax = Math.exp(l_ax);

    l_r   = p_a;
    l_c   = 1.0d;
    l_ans = 1.0d;

    do
      {
      l_r   += 1.0d;
      l_c   *= p_x/l_r;
      l_ans += l_c;
      } while( l_c/l_ans > 1E-10 );

    return (l_ans * (l_ax/p_a));
    }

  
/**
 * Compute the Inverse Normal Cummulative Distribution. 
 * <p>
 * Original algorythm and Perl implementation can
 * be found at: http://www.math.uio.no/~jacklam/notes/invnorm/index.html
 * Author:
 *  Peter J. Acklam
 *  jacklam@math.uio.no
 * </p>
 * 
 * @param p_d The value to compute the Inverse Normal Cummulative
 *            Distribution for.
 * @return  The value of the inverse normal cummulative for the
 *          gaussian distribution.
 */

  public static final double normal_inv_cdf(double   p_d)
    {
    double l_z, l_r;
    
    l_z = 0;

         if(p_d == 0.0d) l_z = Double.NEGATIVE_INFINITY;
    else if(p_d == 1.0d) l_z = Double.POSITIVE_INFINITY;
    else if((p_d < 0) || (p_d > 1) || Double.isNaN(p_d)) l_z = Double.NaN;
    else if(p_d < ICDF_P_LOW)
      {
      l_z = Math.sqrt(-2.0d * Math.log(p_d));
      l_z = (((((ICDF_C[0]*l_z+ICDF_C[1])*l_z+ICDF_C[2])*l_z+
                 ICDF_C[3])*l_z+ICDF_C[4])*l_z+ICDF_C[5]) /
            ((((ICDF_D[0]*l_z+ICDF_D[1])*l_z+ICDF_D[2])*l_z+
                ICDF_D[3])*l_z+1.0);
      }
    else if(ICDF_P_HIGH < p_d)
      {
      l_z = Math.sqrt(-2.0 * Math.log(1-p_d));
      l_z = -(((((ICDF_C[0]*l_z+ICDF_C[1])*l_z+ICDF_C[2])*l_z+
                  ICDF_C[3])*l_z+ICDF_C[4])*l_z+ICDF_C[5]) /
             ((((ICDF_D[0]*l_z+ICDF_D[1])*l_z+ICDF_D[2])*l_z+
                 ICDF_D[3])*l_z+1.0d);
      }
    else
      {
      l_z = (p_d - 0.5d);
      l_r = (l_z * l_z);
      l_z = (((((ICDF_A[0]*l_r+ICDF_A[1])*l_r+ICDF_A[2])*l_r+
                 ICDF_A[3])*l_r+ICDF_A[4])*l_r+ICDF_A[5])*l_z /
            (((((ICDF_B[0]*l_r+ICDF_B[1])*l_r+ICDF_B[2])*l_r+
                 ICDF_B[3])*l_r+ICDF_B[4])*l_r+1.0d);
      }
         
         /*
          * TODO
    if((p_d > 0) && (p_d < 1))      
      {
      l_r = 0.5D * erfc(-l_z/Math.sqrt(2.0D)) - p_d;
      p_d = l_z * Math.sqrt(2.0D*Math.PI) * Math.exp((l_z*l_z)/2.0D);
      l_z = l_z - p_d/(1.0D + l_z*p_d/2.0D);
      }
*/

    return l_z;
    }


/**
 * Computes the error function
 * 
 * The code was taken from FORTRAN and translated to Java.
 * The next paragraph contains information from the original
 * source. 
 * 
 * [
 *  The main computation evaluates near-minimax approximations
 *  from "Rational Chebyshev approximations for the error function"
 *  by W. J. Cody, Math. Comp., 1969, PP. 631-638.  This
 *  transportable program uses rational functions that theoretically
 *  approximate  erf(x)  and  erfc(x)  to at least 18 significant
 *  decimal digits.  The accuracy achieved depends on the arithmetic
 *  system, the compiler, the intrinsic functions, and proper
 *  selection of the machine-dependent constants.
 * 
 *  Author: W. J. Cody
 *  Mathematics and Computer Science Division
 *  Argonne National Laboratory
 *  Argonne, IL 60439  
 * ]
 * 
 * @param   p_x     a real number
 * @param   p_flag  integer indicator
 * @return          see below
 * 
 * <code>p_flag==0</code>:
 *  computes the error function:
 *                    y = 2/pi * integral(0, x, exp(-t*t) dt)
 * 
 * <code>p_flag==1</code>:
 *  computes the complementary error function:
 *                    y = 2/pi * integral(x, infinite, exp(-t*t) dt)
 * 
 * <code>p_flag==2</code>:
 *  computes the scaled complementary error function for large x:
 *                    y = exp(x*x) * erfc(x) * 1/(x*sqrt(pi))
 */

  private static double calerf(double p_x,
                               int    p_flag)
    {
    double  l_result, l_y, l_y_sq, l_x_num, l_x_den, l_del;
    int     l_i;
    
    l_result = 0.0d;
    l_y      = Math.abs(p_x);
    
    if(l_y <= ERF_THRESHOLD)
      {
      l_y_sq = 0.0D;
      if(l_y > EPS) l_y_sq = (l_y * l_y);
      
      l_x_num = (ERF_A[4] * l_y_sq);
      l_x_den = l_y_sq;
      
      for(l_i = 0; l_i < 3; l_i++)
        {
        l_x_num = (l_x_num + ERF_A[l_i]) * l_y_sq;
        l_x_den = (l_x_den + ERF_B[l_i]) * l_y_sq;
        }
      
      l_result = p_x * (l_x_num + ERF_A[3]) / (l_x_den + ERF_B[3]);
      if(p_flag != 0) l_result = (1.0 - l_result);
      if(p_flag == 2) l_result = Math.exp(l_y_sq) * l_result;
      return l_result;
      }
    else if(l_y <= 4.0d)
      {
      l_x_num = (ERF_C[8] * l_y);
      l_x_den = l_y;
      
      for(l_i=0; l_i < 7; l_i++)
        {
        l_x_num = (l_x_num + ERF_C[l_i]) * l_y;
        l_x_den = (l_x_den + ERF_D[l_i]) * l_y;
        }
      
      l_result = (l_x_num + ERF_C[7]) / (l_x_den + ERF_D[7]);
      if(p_flag != 2)
        {
        l_y_sq   = Math.round(l_y*16.0D)/16.0D;
        l_del    = (l_y - l_y_sq) * (l_y + l_y_sq);
        l_result = Math.exp(-l_y_sq*l_y_sq) * Math.exp(-l_del) * l_result;
        }
      }
    else
      {
      l_result = 0.0d;
      
      
      if((l_y < ERF_X_BIG) || ((p_flag == 2) && (l_y < ERF_X_MAX)))
        {
        if((l_y >= ERF_X_BIG) && (l_y >= ERF_X_HUGE))
          {
          l_result = SQRT_PI / l_y;
          }
        else
          {
          l_y_sq  = 1.0D / (l_y * l_y);
          l_x_num = ERF_P[5] * l_y_sq;
          l_x_den = l_y_sq;
          
          for(l_i = 0; l_i<4; l_i++)
            {
            l_x_num = (l_x_num + ERF_P[l_i]) * l_y_sq;
            l_x_den = (l_x_den + ERF_Q[l_i]) * l_y_sq;
            }
          
          l_result = l_y_sq * (l_x_num + ERF_P[4]) / (l_x_den + ERF_Q[4]);
          l_result = (SQRT_PI - l_result) / l_y;
          
          if(p_flag != 2)
            {
            l_y_sq   = Math.round(l_y*16.0D)/16.0D;
            l_del    = (l_y - l_y_sq) * (l_y + l_y_sq);
            l_result = Math.exp(-l_y_sq * l_y_sq) * Math.exp(-l_del) * l_result;
            }
          }
        }
      }

    if(p_flag == 0)
      {
      l_result = (0.5d - l_result) + 0.5D;
      if(p_x < 0) l_result = -l_result;
      }
    else if(p_flag == 1)
      {
      if(p_x < 0) l_result = 2.0D - l_result;
      }
    else
      {
      if(p_x < 0)
        {
        if(p_x < ERF_X_NEG) l_result = Double.MAX_VALUE;
        else
          {
          l_y_sq   = Math.round(p_x*16.0D)/16.0D;
          l_del    = (p_x - l_y_sq) * (p_x + l_y_sq);
          l_y      = Math.exp(l_y_sq * l_y_sq) * Math.exp(l_del);
          l_result = (l_y + l_y) - l_result;
          }
        }
      }
    return l_result;
    }


/**
 * Compute the error function for the value <code>p_d</code>.
 * @param p_d   The parameter of the error function.
 * @return  The value of the error function at
 *          <code>p_d</code>.
 */

  public static final double erf    (double p_d)
    {
    return calerf(p_d, 0);
    }
  
/**
 * Compute the complementary error function for the value <code>p_d</code>.
 * @param p_d   The parameter of the complementary error function.
 * @return  The value of the complementary error function at
 *          <code>p_d</code>.
 */

  public static final double erfc   (double p_d)
    {
    return calerf(p_d, 1);
    }
  
/**
 * Compute the scaled complementary error function for the value
 * <code>p_d</code>.
 * @param p_d   The parameter of the scaled complementary error function.
 * @return  The value of the scaled complementary error function at
 *          <code>p_d</code>.
 * TODO: CHECK
 */

  public static final double erfcx  (double p_d)
    {
    return calerf(p_d, 2);
    }
  
  
/**
 * Computes the cumulative of the normal distribution which is the
 * probability of p_d. 
 * @param   p_d   the x-value
 * @return        1/sqrt(2*pi) * integral(-infinit, p_d, exp(p_d*p_d/2))
 */

  public  static  final double  normal_cdf (double p_d)
    {
    return  0.5d + (0.5d*erf(p_d / Mathematics.SQRT_2));
    }
  
/**
 * Returns the probability density function of the normal distribution.
 * @param p_d The position to get the probability density of.
 * @return  The value of the probability density at position p_d.
 */

  public  static  final double  normal_pdf  (double p_d)
    {
    return (1.0d / SQRT_2_PI) * Math.exp(-0.5d * p_d*p_d);
    }
  
  
/**
 * Computes the cumulative of the chi square distribution which is the
 * probability of p_d. 
 * @param   p_x   the x-value
 * @param   p_dof the degrees of freedom
 * @return        the probability of p_x
 * TODO: CHECK
 */
  
  public  static  final  double  chi_square_cdf(double p_x,
                                                int    p_dof)
    {
    double  l_a, l_y, l_s, l_e, l_c, l_z;
    boolean l_even;
  
     if((p_x <= 0.0d) || (p_dof < 1))
       {
       return 1.0d;
       }
        
     l_a    = 0.5 * p_x;
     l_even = (p_dof & 1) == 0;
     
     if(p_dof > 1)
       {
       l_y = Math.exp(-l_a);
       }
     else l_y = 0.0;
     
     l_s = (l_even ? l_y : (2.0 * normal_cdf(-Math.sqrt(p_x))));
      
     if(p_dof > 2)
       {
       p_x = 0.5 * (p_dof - 1.0);
       l_z = (l_even ? 1.0 : 0.5);
       
       if(l_a > (-Mathematics.LN_EPS))
         {
         l_e = (l_even ? 0.0 : Mathematics.LN_SQRT_PI);
         l_c = Math.log(l_a);
         
         while(l_z <= p_x)
           {
           l_e += Math.log(l_z);
           l_s += Math.exp( (l_c * l_z) - l_a - l_e);
           l_z += 1.0;
           }
         
         return l_s;
         }
       
       l_e = (l_even ? 1.0 : (Mathematics.INV_SQRT_PI / Math.sqrt(l_a)));
       l_c = 0.0;
         
       while(l_z <= p_x)
         {
         l_e *= (l_a / l_z);
         l_c += l_e;
         l_z += 1.0;
         }
       return (l_c * l_y) + l_s;       
       }
     
    return l_s;
    }
  
/**
 * Computes the cumulative of the chi square distribution which is the
 * probability of p_d. 
 * @param   p_x   the x-value
 * @param   p_dof the degrees of freedom
 * @return        the probability of p_x
 * TODO: CHECK
 */
  
  public  static  final double  chi_square_compl_cdf(double p_x,
                                                     int    p_dof)
    {
    if( (p_x < 0.0d) || (p_dof < 1)) return 0.0d;
    
    return gamma_regularized_p( p_dof >> 1, p_x * 0.5d );
    }
  
/**
 * Calculates the chi-square-quantil for the probability p_d with p_dof
 * degrees of freedom.
 * @param   p_d     the probability (often called alpha)
 * @param   p_dof   the degrees of freedom
 * @return          the quantile chi^2 (p_d, p_dof)
 * TODO: CHECK
 */
    
  public  static  final double  chi_quantil (double   p_d,
                                             int      p_dof)
    {
    double  l_max, l_min, l_x;
        
    if(p_d <= 0.0) return Double.POSITIVE_INFINITY;
    if(p_d >= 1.0) return 0.0;
    
    l_max = 1E4d;
    l_min = 0.0;

    for(;;)
      {
      l_x = 0.5 * (l_max + l_min);
            
      if(chi_square_cdf(l_x, p_dof) <= p_d)  
        {
        if(l_max == l_x) return l_x;
        l_max = l_x;        
        }
      else
        {
        if(l_min == l_x) return l_x;
        l_min = l_x;        
        }
      }    
    }
  
  
/**
 * Returns the probability with that a random distribution with
 * the calculated chi-square-value p_d and p_dof degrees of freedom passes
 * the chi-square org.sfc.ztest. 
 * @param     p_d     chi-square value of a measurement series
 * @param     p_dof   degrees of freedom of this measurement series
 * @return    the probability with that the series of measurements passes
 *            the chi-square org.sfc.ztest
 * TODO: CHECK
 */

  public  static  final double  chi_square_pass (double p_d,
                                                 int    p_dof)
    {
    double l_max, l_min, l_x;
  
    l_max = 1.0;
    l_min = 0.0;
    for(;;)
      {
      l_x = 0.5 * (l_max + l_min);
      if(chi_quantil(l_x, p_dof) < p_d)
        {
        if(l_max == l_x) return l_x;
        l_max = l_x;
        }
      else
        {
        if(l_min == l_x) return l_x;
        l_min = l_x;
        }
      }
    }
  
/**
 * Returns the factorial of <code>p_d</code>.
 * @param p_d The number to calculate the factorial of.
 * @return  p_d!
 */

  public  static  final double  factorial (double p_d)
    {
    if(p_d < FACTORIALS.length)
      {
      return FACTORIALS[(int)p_d];
      }
    if(p_d < FACTORIALS_LN.length)
      {
      return Math.rint(FACTORIALS_LN[(int)p_d]);
      }
    return Math.rint(Math.exp(gamma_ln(p_d+1.0d)));
    }
  

/**
 * Returns the natural logarithm of the factorial of <code>p_d</code>.
 * @param p_d The number to calculate the natural logarithm of the
 *            factorial of.
 * @return  ln(p_d!)
 */

  private  static  final double  factorial_ln (double p_d)
    {
    if(p_d < FACTORIALS_LN.length) return FACTORIALS_LN[(int)p_d];
    return gamma_ln(p_d+1.0d);
    }
  
/**
 * Calculate the binominal coeffcient <code>p_n over p_k</code>.
 * @param p_n The number to put above <code>p_k</code>.
 * @param p_k The number to put below <code>p_n</code>.
 * @return  The binominal coefficient <code>p_n over p_k</code>.
 */

  public  static  final double  binominal (double p_n,
                                           double p_k)
    {
    return Math.rint(Math.exp(factorial_ln(p_n) - factorial_ln(p_k) - 
                              factorial_ln(p_n-p_k)));
    }
  
/**
 * Calculate the digamma-function. The digamma-function is often also
 * refered as "psi"-function.
 * @param p_x The argument.
 * @return  The digamma-function's value of <code>p_x</code>.
 */

  public  static  final double digamma  (double   p_x)
    {
    double  l_tn1, l_tn, l_resul, l_tnx;
    int     l_n;
    
    
    /* force into the interval 1..3 */
    if(p_x < 0.0d)
      {
      return digamma(1.0d - p_x) +
             Math.PI + Math.tan(Math.PI*(1.0d - p_x));
      }
    
    if( p_x < 1.0d )
      {
      return digamma(1.0d+p_x)-(1.0d/p_x) ;
      }
    
    if(Math.abs(p_x - 1.0d) <= EPS) return -EULER_CONSTANT ;
    if(Math.abs(p_x - 2.0d) <= EPS) return 1.0d - EULER_CONSTANT;
    if(Math.abs(p_x - 3.0d) <= EPS) return 1.5d - EULER_CONSTANT;
    
    if(p_x > 3.0d)
      {
      return ((0.5d*(digamma(p_x*0.5d)+digamma((p_x+1.0d)*0.5d)))+LN_2);
      }
  
   
    l_tn1   = 1.0d;
    l_tn    = (p_x - 2.0d) ;
    l_resul = DIGAMMA_KNOCE[0] + (DIGAMMA_KNOCE[1]*l_tn) ;

    p_x -= 2.0d ;

    for(l_n = 2 ; l_n < DIGAMMA_KNOCE.length; l_n++)
      {
      l_tnx = 2.0d * p_x * l_tn - l_tn1 ;  
      l_resul += DIGAMMA_KNOCE[l_n]*l_tnx;
      l_tn1 = l_tn ;
      l_tn = l_tnx ;
      }
    
    return l_resul ;
    }
  
  
/**
 * Compute the Bessel function of order 0 of the argument.
 * @param p_x a double value
 * @return The Bessel function of order 0 of the argument.
 * TODO: CHECK
 */


  public static final double j0(double p_x)
    {
    double l_ax, l_y, l_z, l_xx;

    if( (l_ax=Math.abs(p_x)) < 8.0d)
      {
      l_y = (p_x*p_x);
           
      return (57568490574.0d+l_y*(-13362590354.0d+l_y*(651619640.7d
               +l_y*(-11214424.18+l_y*(77392.33017+l_y*(-184.9052456))))))/
             (57568490411.0d+l_y*(1029532985.0d+l_y*(9494680.718d
               +l_y*(59272.64853d+l_y*(267.8532712d+l_y*1.0)))));
      }
    
    l_z     = 8.0d/l_ax;
    l_y     = l_z*l_z;
    l_xx    = l_ax-0.785398164d;
         
    return Math.sqrt(0.636619772d/l_ax)*
           (Math.cos(l_xx)*
               (1.0d+l_y*(-0.1098628627e-2d+l_y*(0.2734510407e-4d
                  +l_y*(-0.2073370639e-5d+l_y*0.2093887211e-6d))))
               -
               l_z*Math.sin(l_xx)*
                 (-0.1562499995e-1d+l_y*(0.1430488765e-3d
              +l_y*(-0.6911147651e-5d+l_y*(0.7621095161e-6d
              -l_y*0.934935152e-7d)))));
    }
  
/**
 * Compute the Bessel function of order 1 of the argument.
 * @param p_x a double value
 * @return The Bessel function of order 1 of the argument.
 * TODO: CHECK
 */

  public static final double j1(double p_x)
    {
    double  l_ax, l_y, l_z, l_xx;

    if((l_ax = Math.abs(p_x)) < 8.0d)
      {
      l_y = (p_x*p_x);
      
      return (p_x*(72362614232.0d+l_y*(-7895059235.0d+l_y*(242396853.1d
             +l_y*(-2972611.439d+l_y*(15704.48260d+l_y*(-30.16036606d)))))))/
              (144725228442.0d+l_y*(2300535178.0d+l_y*(18583304.74d
               +l_y*(99447.43394d+l_y*(376.9991397d+l_y*1.0d)))));
      } 
         
    l_z   = 8.0d/l_ax;
    l_xx  = l_ax-2.356194491d;
    l_y   = l_z*l_z;

    
    l_y    = Math.sqrt(0.636619772/l_ax)*
                   (Math.cos(l_xx)*
                       ((1.0d+l_y*(0.183105e-2d+l_y*(-0.3516396496e-4d
              +l_y*(0.2457520174e-5d+l_y*(-0.240337019e-6d))))))
                       -l_z*Math.sin(l_xx)*
                       (0.04687499995d+l_y*(-0.2002690873e-3d
              +l_y*(0.8449199096e-5d+l_y*(-0.88228987e-6d
              +l_y*0.105787412e-6d)))));
    
    if(p_x < 0.0d) return -l_y;
    return l_y;
    }

/**
 * Compute the Bessel function of order n of the argument.
 * @param p_n integer order
 * @param p_x a double value
 * @return The Bessel function of order n of the argument.
 * TODO: CHECK
 */

  public static final double jn(int     p_n,
                                double  p_x)
    {
    int       l_j,  l_m;
    double    l_ax,  l_bj, l_bjm,  l_bjp,  l_sum,  l_tox,  l_ans;
    boolean   l_jsum;

    if(p_n == 0) return j0(p_x);
    if(p_n == 1) return j1(p_x);

    l_ax = Math.abs(p_x);
    
    if(l_ax == 0.0) return 0.0d;
      
    l_tox = (2.0 / l_ax);
    
    if(l_ax > p_n)
      {
      l_bjm = j0(l_ax);         
      l_bj  = j1(l_ax);
      
      for(l_j = 1; l_j < p_n; l_j++)
        {
        l_bjp = l_j*l_tox*l_bj-l_bjm;
        l_bjm = l_bj;
        l_bj  = l_bjp;
        }
      l_ans = l_bj;
      }
    else
      {
      l_m    = ((p_n + ((int)(Math.sqrt(BESSEL_ACC*p_n))) >> 1)) << 1;
      l_jsum = false;
      l_bjp  = l_ans = l_sum = 0.0d;
      l_bj   = 1.0d;
        
      for(l_j = l_m; l_j > 0; l_j--)
        {
        l_bjm = l_j*l_tox*l_bj - l_bjp;
        l_bjp = l_bj;
        l_bj  = l_bjm;
            
        if(Math.abs(l_bj) > BESSEL_BIGNO)
          {
          l_bj  *= BESSEL_BIGNI;
          l_bjp *= BESSEL_BIGNI;
          l_ans *= BESSEL_BIGNI;
          l_sum *= BESSEL_BIGNI;
          }
          
        if(l_jsum) l_sum += l_bj;
            
        l_jsum = !l_jsum;
        if(l_j == p_n) l_ans = l_bjp;
        }
      
      l_sum = (2.0d*l_sum) - l_bj;
      l_ans /= l_sum;
      }
    return  ((p_x < 0.0d) && ((p_n%2) == 1)) ? -l_ans : l_ans;
    }
  
/**
 * Compute the Bessel function of the second kind, of order 0 of the
 * argument.
 * @param p_x a double value
 * @return The Bessel function of the second kind, of order 0 of the
 *         argument.
 * TODO: CHECK
 */


  public final  static  double y0(double p_x)
    {
    double  l_y, l_z, l_xx;
    
    if(p_x < 8.0d)
      {
      l_y    = (p_x*p_x);
      return ((-2957821389.0d+l_y*(7062834065.0d+l_y*(-512359803.6d
               +l_y*(10879881.29d+l_y*(-86327.92757d+l_y*228.4622733)))))/
              (40076544269.0d+l_y*(745249964.8d+l_y*(7189466.438d
               +l_y*(47447.26470d+l_y*(226.1030244d+l_y*1.0d))))))
          +0.636619772d*j0(p_x)*Math.log(p_x);
      }
    
    l_z   = 8.0d/p_x;
    l_y   = l_z*l_z;
    l_xx  = p_x-0.785398164d;
      
    return Math.sqrt(0.636619772d/p_x)*
              (Math.sin(l_xx)*
                  (1.0d+l_y*(-0.1098628627e-2d+l_y*(0.2734510407e-4d
                   +l_y*(-0.2073370639e-5d+l_y*0.2093887211e-6d))))
                  +l_z*Math.cos(l_xx)*
                  (-0.1562499995e-1d+l_y*(0.1430488765e-3d
                    +l_y*(-0.6911147651e-5d+l_y*(0.7621095161e-6d
                    +l_y*(-0.934945152e-7d))))));      
    }

/**
 * Compute the Bessel function of the second kind, of order 1 of the
 * argument.
 * @param p_x a double value
 * @return The Bessel function of the second kind, of order 1 of the
 *         argument.
 * TODO: CHECK
 */

  public static final double y1(double p_x)
    {
    double  l_y, l_z, l_xx;
   
    if(p_x < 8.0d)
      {
      l_y = (p_x*p_x);      
      return ((p_x*(-0.4900604943e13d+l_y*(0.1275274390e13d
                     +l_y*(-0.5153438139e11d+l_y*(0.7349264551e9d
                     +l_y*(-0.4237922726e7d+l_y*0.8511937935e4d))))))/
                (0.2499580570e14d+l_y*(0.4244419664e12d
                     +l_y*(0.3733650367e10d+l_y*(0.2245904002e8d
                     +l_y*(0.1020426050e6d+l_y*(0.3549632885e3d+l_y)))))))
                     +0.636619772d*(j1(p_x)*Math.log(p_x)-1.0/p_x);
      }
    
    l_z   = 8.0d/p_x;
    l_y   = l_z*l_z;
    l_xx  = p_x-2.356194491d;
        
    
    return Math.sqrt(0.636619772d/p_x)*
          (Math.sin(l_xx)*
              (1.0d+l_y*(0.183105e-2d+l_y*(-0.3516396496e-4d
             +l_y*(0.2457520174e-5d+l_y*(-0.240337019e-6d)))))
              +l_z*Math.cos(l_xx)*
              (0.04687499995d+l_y*(-0.2002690873e-3d
                 +l_y*(0.8449199096e-5d+l_y*(-0.88228987e-6d
                 +l_y*0.105787412e-6d)))));
    }
  
  
/**
 * Compute the Bessel function of the second kind, of order n of the
 * argument.
 * 
 * @param p_n integer order
 * @param p_x a double value
 * @return The Bessel function of the second kind, of order n of the
 *         argument.
 * TODO: CHECK
 */

  public static final double yn (int    p_n,
                                 double p_x)
    {
    double l_by, l_bym, l_byp, l_tox, l_j;

    if(p_n == 0) return y0(p_x);
    if(p_n == 1) return y1(p_x);

    l_tox = (2.0d/p_x);
    l_by  = y1(p_x);
    l_bym = y0(p_x);
    
    for(l_j = 1; l_j < p_n; l_j++)
      {
      l_byp = l_j*l_tox*l_by-l_bym;
      l_bym = l_by;
      l_by  = l_byp;
      }
    
    return l_by;
    }


/**
 * Returns the sum of the first p_k terms of the Poisson distribution.
 * @param p_k number of terms
 * @param p_x double value
 * @return The sum of the first p_k terms of the Poisson distribution.
 * TODO: CHECK
 */


  public static final double poisson_cdf(int      p_k,
                                         double   p_x)
    {
    if( (p_k < 0) || (p_x < 0.0d)) return 0.0d;

    return Mathematics.gamma_regularized_q(p_k+1, p_x);
    }

/** 
 * Returns the sum of the terms p_k+1 to infinity of the Poisson
 * distribution.
 * @param p_k start
 * @param p_x double value
 * @return The sum of the terms p_k+1 to infinity of the Poisson
 *          distribution.
 * TODO: CHECK
 */


  public static final double poisson_compl_cdf(int    p_k,
                                               double p_x)
    {
    if( (p_k < 0) || (p_x < 0.0d) ) return 0.0d;

    return Mathematics.gamma_regularized_p(p_k+1, p_x);
    }


  
/**
 * Calculates the beta function.
 * @param p_z   The first (z) parameter of the beta function.
 * @param p_w   The second (w) parameter of the beta function.
 * @return  The value of the beta function <code>beta(p_z, p_w)</code>. 
 * TODO: CHECK
 */

  public  static  final double  beta      (double p_z,
                                           double p_w)
    {
    return Math.exp(gamma_ln(p_z) + gamma_ln(p_w) - gamma_ln(p_z + p_w));
    }
  
/**
 * Compute the Incomplete Beta Function evaluated from zero to p_xx.
 * @param p_aa double value
 * @param p_bb double value
 * @param p_xx double value
 * @return The Incomplete Beta Function evaluated from zero to p_xx.
 * TODO: CHECK
 */

  public static final double beta_incomplete(double   p_aa,
                                             double   p_bb,
                                             double   p_xx )
    {
    double  l_a, l_b, l_t, l_x, l_xc, l_w, l_y;
    boolean l_flag;

    l_flag = false;
    
    if( ((p_bb * p_xx) <= 1.0d) && (p_xx <= 0.95))
      {
      return pseries(p_aa, p_bb, p_xx);
      }

    l_w = (1.0d - p_xx);

    if( p_xx > (p_aa/(p_aa+p_bb)) )
      {
      l_flag  = true;
      l_a     = p_bb;
      l_b     = p_aa;
      l_xc    = p_xx;
      l_x     = l_w;
      }
    else
      {
      l_a   = p_aa;
      l_b   = p_bb;
      l_xc  = l_w;
      l_x   = p_xx;
      }

    if( l_flag  && ((l_b * l_x) <= 1.0) && (l_x <= 0.95))
      {
      l_t = pseries(l_a, l_b, l_x);
      
      if( l_t <= EPS )  l_t = 1.0d - EPS;
      else                 l_t = 1.0 - l_t;
           
      return l_t;
      }

    l_y = l_x * (l_a+l_b-2.0d) - (l_a-1.0);
    
    if(l_y < 0.0d) l_w = incbcf( l_a, l_b, l_x );
    else           l_w = incbd( l_a, l_b, l_x ) / l_xc;

    l_y = l_a * Math.log(l_x);
    l_t = l_b * Math.log(l_xc);
        
    if( ((l_a+l_b) < MAX_GAMMA) &&
        (Math.abs(l_y) < LN_MAX_DOUBLE) &&
        (Math.abs(l_t) < LN_MAX_DOUBLE))
      {
      l_t = Math.pow(l_xc, l_b);
      l_t *= Math.pow(l_x, l_a);
      l_t /= l_a;
      l_t *= l_w;
      l_t *= gamma(l_a+l_b) / (gamma(l_a) * gamma(l_b));
            
      if( l_flag )
        {
        if( l_t <= EPS )  l_t = 1.0d - EPS;
        else                 l_t = 1.0 - l_t;
        }
      return l_t;
      }
        
    l_y += l_t + gamma_ln(l_a+l_b) - gamma_ln(l_a) - gamma(l_b);
    l_y += Math.log(l_w/l_a);
        
    if( l_y < LN_MIN_DOUBLE ) l_t = 0.0;
    else               l_t = Math.exp(l_y);

    if(l_flag )
      {
      if( l_t <= EPS )  l_t = 1.0d - EPS;
      else                 l_t = 1.0d - l_t;
      }
        
    return l_t;
    }

/**
 * Continued fraction expansion #1 for incomplete beta integral.
 * @param p_a double value
 * @param p_b double value
 * @param p_x double value
 * @return The Incomplete Beta Function evaluated from zero to p_xx.
 */

  private static final  double incbcf( double p_a,
                                       double p_b,
                                       double p_x )
    {
    double  l_xk, l_pk, l_pkm1, l_pkm2, l_qk, l_qkm1, l_qkm2,
            l_k1, l_k2, l_k3, l_k4, l_k5, l_k6, l_k7, l_k8,
            l_r, l_t, l_ans, l_thresh;
       
    int     l_n;
    
    l_k1    = p_a;
    l_k2    = (p_a + p_b);
    l_k3    = p_a;
    l_k4    = (p_a + 1.0d);
    l_k5    = 1.0d;
    l_k6    = (p_b - 1.0d);
    l_k7    = l_k4;
    l_k8    = (p_a + 2.0);

    l_pkm2  = 0.0d;
    l_qkm2  = 1.0d;
    l_pkm1  = 1.0d;
    l_qkm1  = 1.0d;
    l_ans   = 1.0d;
    l_r     = 1.0;
    l_n     = 0;
       
    l_thresh = 3.0 * EPS;
       
    do  
      {
      l_xk    = -( p_x * l_k1 * l_k2 )/( l_k3 * l_k4 );
      l_pk    = l_pkm1 + l_pkm2 * l_xk;
      l_qk    = l_qkm1 + l_qkm2 * l_xk;
      l_pkm2  = l_pkm1;
      l_pkm1  = l_pk;
      l_qkm2  = l_qkm1;
      l_qkm1  = l_qk;

      l_xk    = ( p_x * l_k5 * l_k6 )/( l_k7 * l_k8 );
      l_pk    = l_pkm1 + l_pkm2 * l_xk;
      l_qk    = l_qkm1 + l_qkm2 * l_xk;
      l_pkm2  = l_pkm1;
      l_pkm1  = l_pk;
      l_qkm2  = l_qkm1;
      l_qkm1  = l_qk;

      if( l_qk != 0 )   l_r = l_pk/l_qk;
        
      if( l_r != 0 )
        {
        l_t   = Math.abs( (l_ans - l_r)/l_r );
        l_ans = l_r;
        }
      else
        {
        l_t = 1.0d;
        }

      if( l_t < l_thresh ) return l_ans;

      l_k1 += 1.0d;
      l_k2 += 1.0d;
      l_k3 += 2.0d;
      l_k4 += 2.0d;
      l_k5 += 1.0d;
      l_k6 -= 1.0d;
      l_k7 += 2.0d;
      l_k8 += 2.0d;

      if( (Math.abs(l_qk) + Math.abs(l_pk)) > ICG_BIG)
        {        
        l_pkm2 *= ICG_BIGINV;
        l_pkm1 *= ICG_BIGINV;
        l_qkm2 *= ICG_BIGINV;
        l_qkm1 *= ICG_BIGINV;
        }
          
      if( (Math.abs(l_qk) < ICG_BIG) ||
          (Math.abs(l_pk) < ICG_BIGINV) )
        {
        l_pkm2 *= ICG_BIG;
        l_pkm1 *= ICG_BIG;
        l_qkm2 *= ICG_BIG;
        l_qkm1 *= ICG_BIG;
        }
      } while( (++l_n) < 300 );

    return l_ans;
    }
  
/**
 * Continued fraction expansion #2 for incomplete beta integral
 * @param p_a double value
 * @param p_b double value
 * @param p_x double value
 * @return The Incomplete Beta Function evaluated from zero to p_xx.
 */

  private static  final double incbd( double    p_a,
                                      double    p_b,
                                      double    p_x )
    {
    double l_xk, l_pk, l_pkm1, l_pkm2, l_qk, l_qkm1, l_qkm2,
           l_k1, l_k2, l_k3, l_k4, l_k5, l_k6, l_k7, l_k8,
           l_r, l_t, l_ans, l_z, l_thresh;
    int   l_n;
    
    l_k1    = p_a;
    l_k2    = p_b - 1.0d;
    l_k3    = p_a;
    l_k4    = p_a + 1.0d;
    l_k5    = 1.0d;
    l_k6    = p_a + p_b;
    l_k7    = p_a + 1.0d;
    l_k8    = p_a + 2.0d;

    l_pkm2  = 0.0;
    l_qkm2  = 1.0;
    l_pkm1  = 1.0;
    l_qkm1  = 1.0;
    l_z     = p_x / (1.0d-p_x);
    l_ans   = 1.0;
    l_r     = 1.0;
    l_n     = 0;
    
    l_thresh = 3.0 * EPS;
      
    do
      {
      l_xk    = -( l_z * l_k1 * l_k2 )/( l_k3 * l_k4 );
      l_pk    = l_pkm1 + l_pkm2 * l_xk;
      l_qk    = l_qkm1 + l_qkm2 * l_xk;
      l_pkm2  = l_pkm1;
      l_pkm1  = l_pk;
      l_qkm2  = l_qkm1;
      l_qkm1  = l_qk;

      l_xk    = (l_z * l_k5 * l_k6 )/( l_k7 * l_k8 );
      l_pk    = l_pkm1 + l_pkm2 * l_xk;
      l_qk    = l_qkm1 + l_qkm2 * l_xk;
      l_pkm2  = l_pkm1;
      l_pkm1  = l_pk;
      l_qkm2  = l_qkm1;
      l_qkm1  = l_qk;

      if(l_qk != 0)  l_r = l_pk/l_qk;
      
      if( l_r != 0 )
        {
        l_t   = Math.abs( (l_ans - l_r)/l_r );
        l_ans = l_r;
        }
      else
        {
        l_t   = 1.0d;
        }
       

      if( l_t < l_thresh ) return l_ans;

      l_k1 += 1.0;
      l_k2 -= 1.0;
      l_k3 += 2.0;
      l_k4 += 2.0;
      l_k5 += 1.0;
      l_k6 += 1.0;
      l_k7 += 2.0;
      l_k8 += 2.0;

      if( (Math.abs(l_qk) + Math.abs(l_pk)) > ICG_BIG )
        {
        l_pkm2 *= ICG_BIGINV;
        l_pkm1 *= ICG_BIGINV;
        l_qkm2 *= ICG_BIGINV;
        l_qkm1 *= ICG_BIGINV;
        }
        
      if( (Math.abs(l_qk) < ICG_BIGINV) ||
          (Math.abs(l_pk) < ICG_BIGINV) )
        {
        l_pkm2 *= ICG_BIG;
        l_pkm1 *= ICG_BIG;
        l_qkm2 *= ICG_BIG;
        l_qkm1 *= ICG_BIG;
        }
      } while( ++l_n < 300 );

    return l_ans;
    }
  
/**
 * Power series for incomplete beta integral.
 * Use when b*x is small and x not too close to 1.
 * @param p_a double value
 * @param p_b double value
 * @param p_x double value
 * @return The Incomplete Beta Function evaluated from zero to p_xx.
 */

  private static  final double pseries( double  p_a,
                                        double  p_b,
                                        double  p_x )
    {
    double l_s, l_t, l_u, l_v, l_n, l_t1, l_z, l_ai;

    l_ai    = (1.0d / p_a);
    l_u     = (1.0d - p_b) * p_x;
    l_v     = l_u / (p_a + 1.0);
    l_t1    = l_v;
    l_t     = l_u;
    l_n     = 2.0d;
    l_s     = 0.0d;
    l_z     = EPS * l_ai;
        
    while( Math.abs(l_v) > l_z )
      {
      l_u    = (l_n - p_b) * p_x / l_n;
      l_t   *= l_u;
      l_v    = l_t / (p_a + l_n);
      l_s   += l_v; 
      l_n   += 1.0d;
      }
     
    l_s += l_t1;
    l_s += l_ai;

    l_u = p_a * Math.log(p_x);
    
    if( ((p_a+p_b) < MAX_GAMMA) && (Math.abs(l_u) < LN_MAX_DOUBLE ))
      {
      l_t = gamma(p_a+p_b)/(gamma(p_a)*gamma(p_b));
      l_s = l_s * l_t * Math.pow(p_x, p_a);
      }
    else
      {
      l_t = gamma_ln(p_a+p_b) - gamma_ln(p_a) - gamma_ln(p_b) + l_u +
            Math.log(l_s);
         
      if( l_t < LN_MIN_DOUBLE ) l_s = 0.0d;
      else                      l_s = Math.exp(l_t);
      }
        
    return l_s;
    }
  
  
/**
 * Calculate the <code>p_n</code>'th harmonic number.
 * @param p_n The index of the harmonic number to be created.
 * @return  The <code>p_n</code>'th harmonic number.
 * TODO: improve code - it's really weak
 */

  public  static  final double  harmonic  (double p_n)
    {
    double  l_sum;
    
    l_sum = 1.0d;
    while(p_n > 1.0d)
      {
      l_sum += (1.0d / p_n);
      p_n -= 1.0d;
      }
    
    return l_sum;
    }
  
/**
 * Convert a binary integer to gray code.
 * @param p_data  The int to be converted.
 * @return  The conversation result.
 */

  public  static  final int binary_to_gray (final int p_data)
    {
    int     l_i, l_shift;
    boolean l_b, l_c;
    
    l_shift = (1 << 31);
    l_b     = ((p_data & l_shift) != 0);    
    l_i     = (l_b ? l_shift : 0);
    
    do
      {          
      l_shift >>>= 1;
      l_c = ((p_data & l_shift) != 0);
      if(l_b ^ l_c)
        {
        l_i |= l_shift;
        }
      
      l_b = l_c;
      } while(l_shift != 1);
    
    
    return l_i;
    }
  
/**
 * Convert a gray-code integer to binary format.
 * @param p_data  The int to be converted.
 * @return  The conversation result.
 */

  public  static  final int gray_to_binary (final int p_data)
    {
    int     l_i, l_shift;
    boolean l_b;
    
    l_shift = (1 << 31);
    l_b     = ((p_data & l_shift) != 0);    
    l_i     = (l_b ? l_shift : 0);
    
    do
      {          
      l_shift >>>= 1;
      
      l_b ^= ((p_data & l_shift) != 0);
      if(l_b)
        {
        l_i |= l_shift;
        }
      
      } while(l_shift != 1);
    
    
    return l_i;
    }
  
/**
 * Checks wether the specified number represents a single numerical value,
 * (that means that it is neither infinite nor nan).
 * @param p_number  The number to be checked.
 * @return <code>true</code> if <code>p_number</code> is a normal number,
 *         <code>false</code> it it is either infinite or nan.
 */

  public  static  final boolean is_number  (final double p_number)
    {
    return ((p_number == p_number)                 &&
            (p_number != Double.POSITIVE_INFINITY) &&   
            (p_number != Double.NEGATIVE_INFINITY));
            
    }
  }

File Information:

file name:Mathematics.java
package:org.sfc.math
qualified name:org.sfc.math.Mathematics.java
file type:Java Source File
download location:download http://dgpf.sourceforge.net/source/org/sfc/math/Mathematics.java
size:61.777 KB (63260 B)
uploaded: 2015-07-22 04:11:12 GMT+0000
last update: 2006-03-07 04:20:02 GMT+0000
last access: 2018-01-22 03:53:05 GMT+0000

statistics online since 2006-01-02.   RSS Feed
Contact us by sending an email to tweise@gmx.de to receive further information, to report errors, or to join our project.
All content on this site (http://dgpf.sourceforge.net/) is LGPL-licensed.
http://dgpf.sourceforge.net/scripts/source/source.php last modified at 2015-07-22 04:10:53 GMT+0000 served at 2018-01-22 03:53:05 GMT+0000.
Valid CSS Valid XHTML 1.1
Valid RSS SourceForge.net Logo