?? specialfunction.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, 1992Copyright 1984, 1987, 1989, 1992 by Stephen L. MoshierDirect 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; }}?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -