specialfunction.java

java 作图的程序

   * @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, 1992
Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
Direct 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>