sfun.java

一个用于排队系统仿真的开源软件,有非常形象的图象仿真过程!

		}
		return ans;
	}

	/**
	 *	Returns the complementary error function of a double.
	 *	@param	x	A double value.
	 *	@return  The complementary error function of x.
	 */
	static public double erfc(double x) {
		double ans;
		double y = Math.abs(x);

		if (x <= -6.013687357) {
			// -6.013687357 = -Math.sqrt(-Math.getLog(1.77245385090551602729816748334 * EPSILON_SMALL))
			ans = 2;
		} else if (y < 1.49012e-08) {
			// 1.49012e-08 = Math.sqrt(2*EPSILON_SMALL)
			ans = 1 - 2 * x / 1.77245385090551602729816748334;
		} else {
			double ysq = y * y;
			if (y < 1) {
				ans = 1 - x * (1 + csevl(2 * ysq - 1, ERFC_COEF));
			} else if (y <= 4.0) {
				ans = Math.exp(-ysq) / y * (0.5 + csevl((8.0 / ysq - 5.0) / 3.0, ERFC2_COEF));
				if (x < 0) ans = 2.0 - ans;
				if (x < 0) ans = 2.0 - ans;
				if (x < 0) ans = 2.0 - ans;
			} else {
				ans = Math.exp(-ysq) / y * (0.5 + csevl(8.0 / ysq - 1, ERFCC_COEF));
				if (x < 0) ans = 2.0 - ans;
			}
		}
		return ans;
	}

	/**
	 *	Returns the factorial of an integer.
	 *	@param	n	An integer value.
	 *	@return  The factorial of n, n!.
	 *	If x is negative, the result is NaN.
	 */
	static public double fact(int n) {
		double ans = 1;

		if (Double.isNaN(n) || n < 0) {
			ans = Double.NaN;
		} else if (n > 170) {
			// The 171! is too large to fit in a double.
			ans = Double.POSITIVE_INFINITY;
		} else {
			for (int k = 2; k <= n; k++)
				ans *= k;
		}
		return ans;
	}

	/**
	 *	Returns the Gamma function of a double.
	 *	@param	x	A double value.
	 *	@return  The Gamma function of x.
	 *	If x is a negative integer, the result is NaN.
	 */
	static public double gamma(double x) {
		double ans;
		double y = Math.abs(x);

		if (y <= 10.0) {
			/*
			 * Compute gamma(x) for |x|<=10.
			 * First reduce the interval and  find gamma(1+y) for 0 <= y < 1.
			 */
			int n = (int) x;
			if (x < 0.0) n--;
			y = x - n;
			n--;
			ans = 0.9375 + csevl(2.0 * y - 1.0, GAMMA_COEF);
			if (n == 0) {
			} else if (n < 0) {
				// Compute gamma(x) for x < 1
				n = -n;
				if (x == 0.0) {
					ans = Double.NaN;
				} else if (y < 1.0 / Double.MAX_VALUE) {
					ans = Double.POSITIVE_INFINITY;
				} else {
					double xn = n - 2;
					if (x < 0.0 && x + xn == 0.0) {
						ans = Double.NaN;
					} else {
						for (int i = 0; i < n; i++) {
							ans /= x + i;
						}
					}
				}
			} else {	// gamma(x) for x >= 2.0
				for (int i = 1; i <= n; i++) {
					ans *= y + i;
				}
			}
		} else {  // gamma(x) for |x| > 10
			if (x > 171.614) {
				ans = Double.POSITIVE_INFINITY;
			} else if (x < -170.56) {
				ans = 0.0; // underflows
			} else {
				// 0.9189385332046727 = 0.5*getLog(2*PI)
				ans = Math.exp((y - 0.5) * Math.log(y) - y + 0.9189385332046727 + r9lgmc(y));
				if (x < 0.0) {
					double sinpiy = Math.sin(Math.PI * y);
					if (sinpiy == 0 || Math.round(y) == y) {
						ans = Double.NaN;
					} else {
						ans = -Math.PI / (y * sinpiy * ans);
					}
				}
			}
		}
		return ans;
	}

	/**
	 *	Returns the common (base 10) logarithm of a double.
	 *	@param	x	A double value.
	 *	@return  The common logarithm of x.
	 */
	static public double log10(double x) {
		//if (Double.isNaN(x)) return Double.NaN;
		return 0.43429448190325182765 * Math.log(x);
	}

	/**
	 *	Returns the logarithm of the Beta function.
	 *	@param	a	A double value.
	 *	@param	b	A double value.
	 *	@return  The natural logarithm of the Beta function.
	 */
	static public double logBeta(double a, double b) {
		double corr, ans;
		double p = Math.min(a, b);
		double q = Math.max(a, b);

		if (p <= 0.0) {
			ans = Double.NaN;
		} else if (p >= 10.0) {
			// P and Q are large;
			corr = r9lgmc(p) + r9lgmc(q) - r9lgmc(p + q);
			double temp = dlnrel(-p / (p + q));
			ans = -0.5 * Math.log(q) + 0.918938533204672741780329736406 + corr + (p - 0.5) * Math.log(p / (p + q)) + q * temp;
		} else if (q >= 10.0) {
			// P is small, but Q is large
			corr = Sfun.r9lgmc(q) - r9lgmc(p + q);
			//  Check from underflow from r9lgmc
			ans = logGamma(p) + corr + p - p * Math.log(p + q) + (q - 0.5) * dlnrel(-p / (p + q));
		} else {
			// P and Q are small;
			ans = Math.log(gamma(p) * (gamma(q) / gamma(p + q)));
		}
		return ans;
	}

	/**
	 *	Returns the logarithm of the Gamma function of a double.
	 *	@param	x	A double value.
	 *	@return  The natural logarithm of the Gamma function of x.
	 *	If x is a negative integer, the result is NaN.
	 */
	static public double logGamma(double x) {
		double ans, sinpiy, y;

		y = Math.abs(x);

		if (y <= 10) {
			ans = Math.log(Math.abs(gamma(x)));
		} else if (x > 0) {
			// A&S 6.1.40
			// 0.9189385332046727 = 0.5*getLog(2*PI)
			ans = 0.9189385332046727 + (x - 0.5) * Math.log(x) - x + r9lgmc(y);
		} else {
			sinpiy = Math.abs(Math.sin(Math.PI * y));
			if (sinpiy == 0 || Math.round(y) == y) {
				// The argument for the function can not be a negative integer.
				ans = Double.NaN;
			} else {
				ans = 0.22579135264472743236 + (x - 0.5) * Math.log(y) - x - Math.log(sinpiy) - r9lgmc(y);
			}
		}
		return ans;
	}

	/*
	 *	Returns the getLog gamma correction term for argument
	 *	values greater than or equal to 10.0.
	 */
	static double r9lgmc(double x) {
		double ans;

		if (x < 10.0) {
			ans = Double.NaN;
		} else if (x < 9.490626562e+07) {
			// 9.490626562e+07 = 1/Math.sqrt(EPSILON_SMALL)
			double y = 10.0 / x;
			ans = csevl(2.0 * y * y - 1.0, R9LGMC_COEF) / x;
		} else if (x < 1.39118e+11) {
			// 1.39118e+11 = exp(min(getLog(amach(2) / 12.0), -getLog(12.0 * amach(1))));
			// See A&S 6.1.41
			ans = 1.0 / (12.0 * x);
		} else {
			ans = 0.0; // underflows
		}
		return ans;
	}

	/*
	 *	Returns the value of x with the sign of y.
	 */
	static private double sign(double x, double y) {
		double abs_x = ((x < 0) ? -x : x);
		return (y < 0.0) ? -abs_x : abs_x;
	}

	/**
	 *	Returns the inverse (arc) hyperbolic sine of a double.
	 *	@param	x	A double value.
	 *	@return  The arc hyperbolic sine of x.
	 *	If x is NaN or less than one, the result is NaN.
	 */
	static public double sinh(double x) {
		double ans;
		double y = Math.abs(x);

		if (Double.isNaN(x)) {
			ans = Double.NaN;
		} else if (Double.isInfinite(y)) {
			return x;
		} else if (y < 2.58096e-08) {
			// 2.58096e-08 = Math.sqrt(6.0*EPSILON_SMALL)
			ans = x;
		} else if (y <= 1.0) {
			ans = x * (1.0 + csevl(2.0 * x * x - 1.0, SINH_COEF));
		} else {
			y = Math.exp(y);
			if (y >= 94906265.62) {
				// 94906265.62 = 1.0/Math.sqrt(EPSILON_SMALL)
				ans = sign(0.5 * y, x);
			} else {
				ans = sign(0.5 * (y - 1.0 / y), x);
			}
		}
		return ans;
	}

	/**
	 *	Returns the hyperbolic tangent of a double.
	 *	@param	x	A double value.
	 *	@return  The hyperbolic tangent of x.
	 */
	static public double tanh(double x) {
		double ans, y;
		y = Math.abs(x);

		if (Double.isNaN(x)) {
			ans = Double.NaN;
		} else if (y < 1.82501e-08) {
			// 1.82501e-08 = Math.sqrt(3.0*EPSILON_SMALL)
			ans = x;
		} else if (y <= 1.0) {
			ans = x * (1.0 + csevl(2.0 * x * x - 1.0, TANH_COEF));
		} else if (y < 7.977294885) {
			// 7.977294885 = -0.5*Math.getLog(EPSILON_SMALL)
			y = Math.exp(y);
			ans = sign((y - 1.0 / y) / (y + 1.0 / y), x);
		} else {
			ans = sign(1.0, x);
		}
		return ans;
	}
}