specialfunction.java

java 作图的程序

   * Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier<BR>
   * Direct inquiries to 30 Frost Street, Cambridge, MA 02140<BR>
   */

   static public double erf(double x)
                       throws ArithmeticException { 
       double y, z;
       double T[] = {
                     9.60497373987051638749E0,
                     9.00260197203842689217E1,
                     2.23200534594684319226E3,
                     7.00332514112805075473E3,
                     5.55923013010394962768E4
                    };
       double U[] = {
                   //1.00000000000000000000E0,
                     3.35617141647503099647E1,
                     5.21357949780152679795E2,
                     4.59432382970980127987E3,
                     2.26290000613890934246E4,
                     4.92673942608635921086E4
                    };

       if( Math.abs(x) > 1.0 ) return( 1.0 - erfc(x) );
       z = x * x;
       y = x * polevl( z, T, 4 ) / p1evl( z, U, 5 );
       return y;
}





      static  private double polevl( double x, double coef[], int N )
                              throws ArithmeticException {

           double ans;

           ans = coef[0];

           for(int i=1; i<=N; i++) { ans = ans*x+coef[i]; }

           return ans;
       }

      static  private double p1evl( double x, double coef[], int N )
                              throws ArithmeticException {

          double ans;

          ans = x + coef[0];

          for(int i=1; i<N; i++) { ans = ans*x+coef[i]; }

          return ans;
     }
/*
 *
 *	Natural logarithm of gamma function
 *
 */
/*
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 lgamma(double x)
                              throws ArithmeticException {
         double p, q, w, z;

         double A[] = {
                       8.11614167470508450300E-4,
                       -5.95061904284301438324E-4,
                        7.93650340457716943945E-4,
                       -2.77777777730099687205E-3,
                        8.33333333333331927722E-2
                       };
         double B[] = {
                       -1.37825152569120859100E3,
                       -3.88016315134637840924E4,
                       -3.31612992738871184744E5,
                       -1.16237097492762307383E6,
                       -1.72173700820839662146E6,
                       -8.53555664245765465627E5
                       };
         double C[] = {
                       /* 1.00000000000000000000E0, */
                       -3.51815701436523470549E2,
                       -1.70642106651881159223E4,
                       -2.20528590553854454839E5,
                       -1.13933444367982507207E6,
                       -2.53252307177582951285E6,
                       -2.01889141433532773231E6
                      };

         if( x < -34.0 ) {
  	   q = -x;
	   w = lgamma(q);
	   p = Math.floor(q);
	   if( p == q ) throw new ArithmeticException("lgam: Overflow");
	   z = q - p;
	   if( z > 0.5 ) {
		p += 1.0;
		z = p - q;
 	   }
	   z = q * Math.sin( Math.PI * z );
	   if( z == 0.0 ) throw new 
                               ArithmeticException("lgamma: Overflow");
	   z = LOGPI - Math.log( z ) - w;
	   return z;
	 }

         if( x < 13.0 ) {
  	   z = 1.0;
	   while( x >= 3.0 ) {
		x -= 1.0;
		z *= x;
	   }
	   while( x < 2.0 ) {
		if( x == 0.0 ) throw new 
                                ArithmeticException("lgamma: Overflow");
		z /= x;
		x += 1.0;
	   }
	   if( z < 0.0 ) z = -z;
	   if( x == 2.0 ) return Math.log(z);
	   x -= 2.0;
	   p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
 	   return( Math.log(z) + p );
	 }

         if( x > 2.556348e305 ) throw new 
                          ArithmeticException("lgamma: Overflow");

         q = ( x - 0.5 ) * Math.log(x) - x + 0.91893853320467274178;
         if( x > 1.0e8 ) return( q );

         p = 1.0/(x*x);
         if( x >= 1000.0 )
	     q += ((   7.9365079365079365079365e-4 * p
		      - 2.7777777777777777777778e-3) *p
		     + 0.0833333333333333333333) / x;
         else
	     q += polevl( p, A, 4 ) / x;
         return q;
     }



  /**
   * @param aa double value
   * @param bb double value
   * @param xx double value
   * @return The Incomplete Beta Function evaluated from zero to xx.
   * <P>
   * <FONT size=2>
   * Converted to Java from<BR>
   * Cephes Math Library Release 2.3:  July, 1995<BR>
   * Copyright 1984, 1995 by Stephen L. Moshier<BR>
   * Direct inquiries to 30 Frost Street, Cambridge, MA 02140<BR>
   */

     public static double ibeta( double aa, double bb, double xx )
                              throws ArithmeticException {
        double a, b, t, x, xc, w, y;
        boolean flag;

        if( aa <= 0.0 || bb <= 0.0 ) throw new 
                          ArithmeticException("ibeta: Domain error!");

        if( (xx <= 0.0) || ( xx >= 1.0) ) {
  	       if( xx == 0.0 ) return 0.0;
   	       if( xx == 1.0 ) return 1.0;
           throw new ArithmeticException("ibeta: Domain error!");
	    }

        flag = false;
        if( (bb * xx) <= 1.0 && xx <= 0.95) {
	        t = pseries(aa, bb, xx);
		    return t;
	    }

        w = 1.0 - xx;

        /* Reverse a and b if x is greater than the mean. */
        if( xx > (aa/(aa+bb)) ) {
	       flag = true;
	       a = bb;
	       b = aa;
	       xc = xx;
	       x = w;
	    } else {
  	       a = aa;
	       b = bb;
	       xc = w;
	       x = xx;
	    }

        if( flag  && (b * x) <= 1.0 && x <= 0.95) {
 	       t = pseries(a, b, x);
	       if( t <= MACHEP ) 	t = 1.0 - MACHEP;
	       else  		        t = 1.0 - t;
           return t;
	    }

        /* Choose expansion for better convergence. */
        y = x * (a+b-2.0) - (a-1.0);
        if( y < 0.0 )
	                  w = incbcf( a, b, x );
        else
	                  w = incbd( a, b, x ) / xc;

        /* Multiply w by the factor
           a      b   _             _     _
          x  (1-x)   | (a+b) / ( a | (a) | (b) ) .   */

        y = a * Math.log(x);
        t = b * Math.log(xc);
        if( (a+b) < MAXGAM && Math.abs(y) < MAXLOG && Math.abs(t) < MAXLOG ) {
	        t = Math.pow(xc,b);
	        t *= Math.pow(x,a);
	        t /= a;
	        t *= w;
	        t *= gamma(a+b) / (gamma(a) * gamma(b));
            if( flag ) {
 	           if( t <= MACHEP ) 	t = 1.0 - MACHEP;
	           else  		        t = 1.0 - t;
	        }
            return t;
	    }
        /* Resort to logarithms.  */
        y += t + lgamma(a+b) - lgamma(a) - lgamma(b);
        y += Math.log(w/a);
        if( y < MINLOG )
	                    t = 0.0;
        else
	                    t = Math.exp(y);

        if( flag ) {
 	           if( t <= MACHEP ) 	t = 1.0 - MACHEP;
	           else  		        t = 1.0 - t;
	    }
        return t;
   }

/* Continued fraction expansion #1
 * for incomplete beta integral
 */

    private static double incbcf( double a, double b, double x )
                              throws ArithmeticException {
       double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
       double k1, k2, k3, k4, k5, k6, k7, k8;
       double r, t, ans, thresh;
       int n;
       double big = 4.503599627370496e15;
       double biginv =  2.22044604925031308085e-16;

       k1 = a;
       k2 = a + b;
       k3 = a;
       k4 = a + 1.0;
       k5 = 1.0;
       k6 = b - 1.0;
       k7 = k4;
       k8 = a + 2.0;

       pkm2 = 0.0;
       qkm2 = 1.0;
       pkm1 = 1.0;
       qkm1 = 1.0;
       ans = 1.0;
       r = 1.0;
       n = 0;
       thresh = 3.0 * MACHEP;
       do {
	      xk = -( x * k1 * k2 )/( k3 * k4 );
	      pk = pkm1 +  pkm2 * xk;
	      qk = qkm1 +  qkm2 * xk;
	      pkm2 = pkm1;
	      pkm1 = pk;
	      qkm2 = qkm1;
	      qkm1 = qk;

	      xk = ( x * k5 * k6 )/( k7 * k8 );
	      pk = pkm1 +  pkm2 * xk;
	      qk = qkm1 +  qkm2 * xk;
	      pkm2 = pkm1;
	      pkm1 = pk;
	      qkm2 = qkm1;
	      qkm1 = qk;

	      if( qk != 0 )		r = pk/qk;
	      if( r != 0 ) {
		       t = Math.abs( (ans - r)/r );
		       ans = r;
		  }	else
		       t = 1.0;

	      if( t < thresh ) return ans;

	      k1 += 1.0;
          k2 += 1.0;
      	  k3 += 2.0;
      	  k4 += 2.0;
      	  k5 += 1.0;
      	  k6 -= 1.0;
      	  k7 += 2.0;
      	  k8 += 2.0;

      	  if( (Math.abs(qk) + Math.abs(pk)) > big ) {
      		pkm2 *= biginv;
      		pkm1 *= biginv;
      		qkm2 *= biginv;
      		qkm1 *= biginv;
		  }
      	  if( (Math.abs(qk) < biginv) || (Math.abs(pk) < biginv) ) {
      		pkm2 *= big;
      		pkm1 *= big;
      		qkm2 *= big;
      		qkm1 *= big;
		  }
	   } while( ++n < 300 );

    return ans;
   }
/* Continued fraction expansion #2
 * for incomplete beta integral
 */

      static private double incbd( double a, double b, double x )
                              throws ArithmeticException {
         double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
         double k1, k2, k3, k4, k5, k6, k7, k8;
         double r, t, ans, z, thresh;
         int n;
         double big = 4.503599627370496e15;
         double biginv =  2.22044604925031308085e-16;

         k1 = a;
         k2 = b - 1.0;
         k3 = a;
         k4 = a + 1.0;
         k5 = 1.0;
         k6 = a + b;
         k7 = a + 1.0;;
         k8 = a + 2.0;

         pkm2 = 0.0;
         qkm2 = 1.0;
         pkm1 = 1.0;
         qkm1 = 1.0;
         z = x / (1.0-x);
         ans = 1.0;
         r = 1.0;
         n = 0;
         thresh = 3.0 * MACHEP;
         do {
	         xk = -( z * k1 * k2 )/( k3 * k4 );
	         pk = pkm1 +  pkm2 * xk;
	         qk = qkm1 +  qkm2 * xk;
	         pkm2 = pkm1;
	         pkm1 = pk;
	         qkm2 = qkm1;
	         qkm1 = qk;

	         xk = ( z * k5 * k6 )/( k7 * k8 );
	         pk = pkm1 +  pkm2 * xk;
	         qk = qkm1 +  qkm2 * xk;
	         pkm2 = pkm1;
	         pkm1 = pk;
	         qkm2 = qkm1;
	         qkm1 = qk;

	         if( qk != 0 )  r = pk/qk;
	         if( r != 0 ) {
		         t = Math.abs( (ans - r)/r );
		         ans = r;
		     } else
		         t = 1.0;

	         if( t < thresh ) return ans;

	         k1 += 1.0;
	         k2 -= 1.0;
	         k3 += 2.0;
	         k4 += 2.0;
	         k5 += 1.0;
	         k6 += 1.0;
	         k7 += 2.0;
	         k8 += 2.0;

	         if( (Math.abs(qk) + Math.abs(pk)) > big ) {
		        pkm2 *= biginv;
		        pkm1 *= biginv;
		        qkm2 *= biginv;
		        qkm1 *= biginv;
		     }
	         if( (Math.abs(qk) < biginv) || (Math.abs(pk) < biginv) ) {
		        pkm2 *= big;
		        pkm1 *= big;
		        qkm2 *= big;
		        qkm1 *= big;
		     }
	    } while( ++n < 300 );

        return ans;
     }
/* Power series for incomplete beta integral.
   Use when b*x is small and x not too close to 1.  */

     static private  double pseries( double a, double b, double x )
                              throws ArithmeticException {
        double s, t, u, v, n, t1, z, ai;

        ai = 1.0 / a;
        u = (1.0 - b) * x;
        v = u / (a + 1.0);
        t1 = v;
        t = u;
        n = 2.0;
        s = 0.0;
        z = MACHEP * ai;
        while( Math.abs(v) > z ) {
	       u = (n - b) * x / n;
	       t *= u;
	       v = t / (a + n);
	       s += v; 
	       n += 1.0;
	    }
        s += t1;
        s += ai;

        u = a * Math.log(x);
        if( (a+b) < MAXGAM && Math.abs(u) < MAXLOG ) {
	        t = gamma(a+b)/(gamma(a)*gamma(b));
	        s = s * t * Math.pow(x,a);
	    } else {
	       t = lgamma(a+b) - lgamma(a) - lgamma(b) + u + Math.log(s);
	       if( t < MINLOG ) 	s = 0.0;
	       else  	            s = Math.exp(t);
	    }
        return s;
     }

}