* @param x an integer value   * @return the factorial of the argument   */     static public int fac(int j) throws ArithmeticException {        int i = j;        int d = 1;        if(j < 0) i = Math.abs(j);                while( i > 1) { d *= i--; }        if(j < 0) return -d;        else      return d;     }  /**   * @param x a double value   * @return the Gamma function of the value.   * <P>   * <FONT size=2>   * Converted to Java from<BR>   * Cephes Math Library Release 2.2:  July, 1992<BR>   * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier<BR>   * Direct inquiries to 30 Frost Street, Cambridge, MA 02140<BR>   **/    static public double gamma(double x) throws ArithmeticException {     double P[] = {                   1.60119522476751861407E-4,                   1.19135147006586384913E-3,                   1.04213797561761569935E-2,                   4.76367800457137231464E-2,                   2.07448227648435975150E-1,                   4.94214826801497100753E-1,                   9.99999999999999996796E-1                  };     double Q[] = {                   -2.31581873324120129819E-5,                    5.39605580493303397842E-4,                   -4.45641913851797240494E-3,                    1.18139785222060435552E-2,                    3.58236398605498653373E-2,                   -2.34591795718243348568E-1,                    7.14304917030273074085E-2,                    1.00000000000000000320E0                   };     double MAXGAM = 171.624376956302725;     double LOGPI  = 1.14472988584940017414;     double p, z;     int i;     double q = Math.abs(x);     if( q > 33.0 ) {       if( x < 0.0 ) {            p = Math.floor(q);	    if( p == q ) throw new ArithmeticException("gamma: overflow");	    i = (int)p;	    z = q - p;	    if( z > 0.5 ) {			p += 1.0;			z = q - p;	    }	    z = q * Math.sin( Math.PI * z );	    if( z == 0.0 ) throw new ArithmeticException("gamma: overflow");	    z = Math.abs(z);	    z = Math.PI/(z * stirf(q) );            return -z;       } else {	    return stirf(x);       }     }     z = 1.0;       while( x >= 3.0 ) {  	     x -= 1.0;	     z *= x;       }       while( x < 0.0 ) {	     if( x == 0.0 ) {                throw new ArithmeticException("gamma: singular");             } else	     if( x > -1.E-9 ) {                 return( z/((1.0 + 0.5772156649015329 * x) * x) );             }	     z /= x;	     x += 1.0;       }       while( x < 2.0 ) {	     if( x == 0.0 ) {                throw new ArithmeticException("gamma: singular");             } else	     if( x < 1.e-9 ) {  	        return( z/((1.0 + 0.5772156649015329 * x) * x) );             }	     z /= x;	     x += 1.0;	}        if( (x == 2.0) || (x == 3.0) ) 	return z;        x -= 2.0;        p = polevl( x, P, 6 );        q = polevl( x, Q, 7 );        return  z * p / q;    }/* Gamma function computed by Stirling's formula. * The polynomial STIR is valid for 33 <= x <= 172.Cephes Math Library Release 2.2:  July, 1992Copyright 1984, 1987, 1989, 1992 by Stephen L. MoshierDirect inquiries to 30 Frost Street, Cambridge, MA 02140*/     static private double stirf(double x) throws ArithmeticException {        double STIR[] = {                     7.87311395793093628397E-4,                    -2.29549961613378126380E-4,                    -2.68132617805781232825E-3,                     3.47222221605458667310E-3,                     8.33333333333482257126E-2,                    };        double MAXSTIR = 143.01608;        double w = 1.0/x;        double  y = Math.exp(x);        w = 1.0 + w * polevl( w, STIR, 4 );        if( x > MAXSTIR ) {	       /* Avoid overflow in Math.pow() */	       double v = Math.pow( x, 0.5 * x - 0.25 );	       y = v * (v / y);	} else {               y = Math.pow( x, x - 0.5 ) / y;	}        y = SQTPI * y * w;        return y;     }  /**   * @param a double value   * @param x double value   * @return the Complemented Incomplete Gamma function.   * <P>   * <FONT size=2>   * Converted to Java from<BR>   * Cephes Math Library Release 2.2:  July, 1992<BR>   * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier<BR>   * Direct inquiries to 30 Frost Street, Cambridge, MA 02140<BR>   **/     static public double igamc( double a, double x )                         throws ArithmeticException {        double big    = 4.503599627370496e15;        double biginv =  2.22044604925031308085e-16;        double ans, ax, c, yc, r, t, y, z;        double pk, pkm1, pkm2, qk, qkm1, qkm2;        if( x <= 0 || a <= 0 ) return 1.0;        if( x < 1.0 || x < a ) return 1.0 - igam(a,x);        ax = a * Math.log(x) - x - lgamma(a);        if( ax < -MAXLOG ) return 0.0;        ax = Math.exp(ax);        /* continued fraction */        y = 1.0 - a;        z = x + y + 1.0;        c = 0.0;        pkm2 = 1.0;        qkm2 = x;        pkm1 = x + 1.0;        qkm1 = z * x;        ans = pkm1/qkm1;        do {  	    c += 1.0;	    y += 1.0;	    z += 2.0;	    yc = y * c;	    pk = pkm1 * z  -  pkm2 * yc;	    qk = qkm1 * z  -  qkm2 * yc;	    if( qk != 0 ) {		r = pk/qk;		t = Math.abs( (ans - r)/r );		ans = r;	    } else		t = 1.0;	    pkm2 = pkm1;	    pkm1 = pk;	    qkm2 = qkm1;	    qkm1 = qk;	    if( Math.abs(pk) > big ) {		pkm2 *= biginv;		pkm1 *= biginv;		qkm2 *= biginv;		qkm1 *= biginv;	    }	} while( t > MACHEP );        return ans * ax;     }  /**   * @param a double value   * @param x double value   * @return the Incomplete Gamma function.   * <P>   * <FONT size=2>   * Converted to Java from<BR>   * Cephes Math Library Release 2.2:  July, 1992<BR>   * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier<BR>   * Direct inquiries to 30 Frost Street, Cambridge, MA 02140<BR>   **/     static public double igam(double a, double x)                          throws ArithmeticException {        double ans, ax, c, r;        if( x <= 0 || a <= 0 ) return 0.0;        if( x > 1.0 && x > a ) return 1.0 - igamc(a,x);       /* Compute  x**a * exp(-x) / gamma(a)  */        ax = a * Math.log(x) - x - lgamma(a);        if( ax < -MAXLOG ) return( 0.0 );        ax = Math.exp(ax);        /* power series */        r = a;        c = 1.0;        ans = 1.0;        do {  	    r += 1.0;	    c *= x/r;	    ans += c;	}        while( c/ans > MACHEP );        return( ans * ax/a );     }  /**   * Returns the area under the left hand tail (from 0 to x)   * of the Chi square probability density function with   * v degrees of freedom.   *   * @param df degrees of freedom   * @param x double value   * @return the Chi-Square function.   **/   static public double chisq(double df, double x)                        throws ArithmeticException {         if( x < 0.0 || df < 1.0 ) return 0.0;        return igam( df/2.0, x/2.0 );   }  /**   * Returns the area under the right hand tail (from x to   * infinity) of the Chi square probability density function   * with v degrees of freedom:   *   * @param df degrees of freedom   * @param x double value   * @return the Chi-Square function.   * <P>   **/   static public double chisqc(double df, double x)                        throws ArithmeticException {         if( x < 0.0 || df < 1.0 ) return 0.0;        return igamc( df/2.0, x/2.0 );   }  /**   * Returns the sum of the first k terms of the Poisson   * distribution.   * @param k number of terms   * @param x double value   */   static public double poisson(int k, double x)                        throws ArithmeticException {        if( k < 0 || x < 0 ) return 0.0;    return igamc((double)(k+1) ,x);   }  /**    * Returns the sum of the terms k+1 to infinity of the Poisson   * distribution.   * @param k start   * @param x double value   */   static public double poissonc(int k, double x)                        throws ArithmeticException {        if( k < 0 || x < 0 ) return 0.0;    return igam((double)(k+1),x);   }  /**   * @param a double value   * @return The area under the Gaussian probability density   * function, integrated from minus infinity to x:   */   static public double normal( double a)                       throws ArithmeticException {       double x, y, z;      x = a * SQRTH;      z = Math.abs(x);      if( z < SQRTH )   y = 0.5 + 0.5 * erf(x);      else {                        y = 0.5 * erfc(z);                        if( x > 0 )  y = 1.0 - y;      }      return y;   }  /**   * @param a double value   * @return The complementary Error function   * <P>   * <FONT size=2>   * Converted to Java from<BR>   * Cephes Math Library Release 2.2:  July, 1992<BR>   * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier<BR>   * Direct inquiries to 30 Frost Street, Cambridge, MA 02140<BR>   */     static public double erfc(double a)                       throws ArithmeticException {        double x,y,z,p,q;       double P[] = {                     2.46196981473530512524E-10,                     5.64189564831068821977E-1,                     7.46321056442269912687E0,                     4.86371970985681366614E1,                     1.96520832956077098242E2,                     5.26445194995477358631E2,                     9.34528527171957607540E2,                     1.02755188689515710272E3,                     5.57535335369399327526E2                    };       double Q[] = {                    //1.0                      1.32281951154744992508E1,                      8.67072140885989742329E1,                      3.54937778887819891062E2,                      9.75708501743205489753E2,                      1.82390916687909736289E3,                      2.24633760818710981792E3,                      1.65666309194161350182E3,                      5.57535340817727675546E2                     };       double R[] = {                      5.64189583547755073984E-1,                      1.27536670759978104416E0,                      5.01905042251180477414E0,                      6.16021097993053585195E0,                      7.40974269950448939160E0,                      2.97886665372100240670E0                     };       double S[] = {                    //1.00000000000000000000E0,                       2.26052863220117276590E0,                      9.39603524938001434673E0,                      1.20489539808096656605E1,                      1.70814450747565897222E1,                      9.60896809063285878198E0,                      3.36907645100081516050E0                     };        if( a < 0.0 )   x = -a;        else            x = a;        if( x < 1.0 )   return 1.0 - erf(a);        z = -a * a;        if( z < -MAXLOG ) {             if( a < 0 )  return( 2.0 );             else         return( 0.0 );        }        z = Math.exp(z);        if( x < 8.0 ) {          p = polevl( x, P, 8 );          q = p1evl( x, Q, 8 );        } else {          p = polevl( x, R, 5 );          q = p1evl( x, S, 6 );        }        y = (z * p)/q;        if( a < 0 ) y = 2.0 - y;        if( y == 0.0 ) {                if( a < 0 ) return 2.0;                else        return( 0.0 );         }        return y;   }  /**   * @param a double value   * @return The Error function   * <P>   * <FONT size=2>   * Converted to Java from<BR>   * Cephes Math Library Release 2.2:  July, 1992<BR>