MATH-1416: Delete functionality now in "Commons Numbers".
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/44ab2569 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/44ab2569 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/44ab2569 Branch: refs/heads/master Commit: 44ab256961029c9e104e70fb804617fd582412cd Parents: 2ec4dea Author: Gilles <er...@apache.org> Authored: Sun May 14 23:30:00 2017 +0200 Committer: Gilles <er...@apache.org> Committed: Sun May 14 23:31:25 2017 +0200 ---------------------------------------------------------------------- .../org/apache/commons/math4/special/Beta.java | 484 --------- .../org/apache/commons/math4/special/Gamma.java | 746 -------------- .../apache/commons/math4/special/BetaTest.java | 979 ------------------ .../apache/commons/math4/special/GammaTest.java | 998 ------------------- 4 files changed, 3207 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/44ab2569/src/main/java/org/apache/commons/math4/special/Beta.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/special/Beta.java b/src/main/java/org/apache/commons/math4/special/Beta.java deleted file mode 100644 index 9355099..0000000 --- a/src/main/java/org/apache/commons/math4/special/Beta.java +++ /dev/null @@ -1,484 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math4.special; - -import org.apache.commons.numbers.gamma.Gamma; -import org.apache.commons.numbers.gamma.LogGamma; -import org.apache.commons.math4.exception.NumberIsTooSmallException; -import org.apache.commons.math4.exception.OutOfRangeException; -import org.apache.commons.math4.util.ContinuedFraction; -import org.apache.commons.math4.util.FastMath; - -/** - * <p> - * This is a utility class that provides computation methods related to the - * Beta family of functions. - * </p> - * <p> - * Implementation of {@link #logBeta(double, double)} is based on the - * algorithms described in - * <ul> - * <li><a href="http://dx.doi.org/10.1145/22721.23109";>Didonato and Morris - * (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios - * and their Inverse</em>, TOMS 12(4), 377-393,</li> - * <li><a href="http://dx.doi.org/10.1145/131766.131776";>Didonato and Morris - * (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the - * Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li> - * </ul> - * and implemented in the - * <a href="http://www.dtic.mil/docs/citations/ADA476840";>NSWC Library of Mathematical Functions</a>, - * available - * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html";>here</a>. - * This library is "approved for public release", and the - * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf";>Copyright guidance</a> - * indicates that unless otherwise stated in the code, all FORTRAN functions in - * this library are license free. Since no such notice appears in the code these - * functions can safely be ported to Commons-Math. - */ -public class Beta { - /** Maximum allowed numerical error. */ - private static final double DEFAULT_EPSILON = 1E-14; - - /** The constant value of ½log 2Ï. */ - private static final double HALF_LOG_TWO_PI = .9189385332046727; - - /** - * <p> - * The coefficients of the series expansion of the Î function. This function - * is defined as follows - * </p> - * <center>Î(x) = log Î(x) - (x - 0.5) log a + a - 0.5 log 2Ï,</center> - * <p> - * see equation (23) in Didonato and Morris (1992). The series expansion, - * which applies for x ⥠10, reads - * </p> - * <pre> - * 14 - * ==== - * 1 \ 2 n - * Î(x) = --- > d (10 / x) - * x / n - * ==== - * n = 0 - * <pre> - */ - private static final double[] DELTA = { - .833333333333333333333333333333E-01, - -.277777777777777777777777752282E-04, - .793650793650793650791732130419E-07, - -.595238095238095232389839236182E-09, - .841750841750832853294451671990E-11, - -.191752691751854612334149171243E-12, - .641025640510325475730918472625E-14, - -.295506514125338232839867823991E-15, - .179643716359402238723287696452E-16, - -.139228964661627791231203060395E-17, - .133802855014020915603275339093E-18, - -.154246009867966094273710216533E-19, - .197701992980957427278370133333E-20, - -.234065664793997056856992426667E-21, - .171348014966398575409015466667E-22 - }; - - /** - * Default constructor. Prohibit instantiation. - */ - private Beta() {} - - /** - * Returns the - * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html";> - * regularized beta function</a> I(x, a, b). - * - * @param x Value. - * @param a Parameter {@code a}. - * @param b Parameter {@code b}. - * @return the regularized beta function I(x, a, b). - * @throws org.apache.commons.math4.exception.MaxCountExceededException - * if the algorithm fails to converge. - */ - public static double regularizedBeta(double x, double a, double b) { - return regularizedBeta(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE); - } - - /** - * Returns the - * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html";> - * regularized beta function</a> I(x, a, b). - * - * @param x Value. - * @param a Parameter {@code a}. - * @param b Parameter {@code b}. - * @param epsilon When the absolute value of the nth item in the - * series is less than epsilon the approximation ceases to calculate - * further elements in the series. - * @return the regularized beta function I(x, a, b) - * @throws org.apache.commons.math4.exception.MaxCountExceededException - * if the algorithm fails to converge. - */ - public static double regularizedBeta(double x, - double a, double b, - double epsilon) { - return regularizedBeta(x, a, b, epsilon, Integer.MAX_VALUE); - } - - /** - * Returns the regularized beta function I(x, a, b). - * - * @param x the value. - * @param a Parameter {@code a}. - * @param b Parameter {@code b}. - * @param maxIterations Maximum number of "iterations" to complete. - * @return the regularized beta function I(x, a, b) - * @throws org.apache.commons.math4.exception.MaxCountExceededException - * if the algorithm fails to converge. - */ - public static double regularizedBeta(double x, - double a, double b, - int maxIterations) { - return regularizedBeta(x, a, b, DEFAULT_EPSILON, maxIterations); - } - - /** - * Returns the regularized beta function I(x, a, b). - * - * The implementation of this method is based on: - * <ul> - * <li> - * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html";> - * Regularized Beta Function</a>.</li> - * <li> - * <a href="http://functions.wolfram.com/06.21.10.0001.01";> - * Regularized Beta Function</a>.</li> - * </ul> - * - * @param x the value. - * @param a Parameter {@code a}. - * @param b Parameter {@code b}. - * @param epsilon When the absolute value of the nth item in the - * series is less than epsilon the approximation ceases to calculate - * further elements in the series. - * @param maxIterations Maximum number of "iterations" to complete. - * @return the regularized beta function I(x, a, b) - * @throws org.apache.commons.math4.exception.MaxCountExceededException - * if the algorithm fails to converge. - */ - public static double regularizedBeta(double x, - final double a, final double b, - double epsilon, int maxIterations) { - double ret; - - if (Double.isNaN(x) || - Double.isNaN(a) || - Double.isNaN(b) || - x < 0 || - x > 1 || - a <= 0 || - b <= 0) { - ret = Double.NaN; - } else if (x > (a + 1) / (2 + b + a) && - 1 - x <= (b + 1) / (2 + b + a)) { - ret = 1 - regularizedBeta(1 - x, b, a, epsilon, maxIterations); - } else { - ContinuedFraction fraction = new ContinuedFraction() { - - /** {@inheritDoc} */ - @Override - protected double getB(int n, double x) { - double ret; - double m; - if (n % 2 == 0) { // even - m = n / 2.0; - ret = (m * (b - m) * x) / - ((a + (2 * m) - 1) * (a + (2 * m))); - } else { - m = (n - 1.0) / 2.0; - ret = -((a + m) * (a + b + m) * x) / - ((a + (2 * m)) * (a + (2 * m) + 1.0)); - } - return ret; - } - - /** {@inheritDoc} */ - @Override - protected double getA(int n, double x) { - return 1.0; - } - }; - ret = FastMath.exp((a * FastMath.log(x)) + (b * FastMath.log1p(-x)) - - FastMath.log(a) - logBeta(a, b)) * - 1.0 / fraction.evaluate(x, epsilon, maxIterations); - } - - return ret; - } - - /** - * Returns the value of log Î(a + b) for 1 ⤠a, b ⤠2. Based on the - * <em>NSWC Library of Mathematics Subroutines</em> double precision - * implementation, {@code DGSMLN}. In {@code BetaTest.testLogGammaSum()}, - * this private method is accessed through reflection. - * - * @param a First argument. - * @param b Second argument. - * @return the value of {@code log(Gamma(a + b))}. - * @throws OutOfRangeException if {@code a} or {@code b} is lower than - * {@code 1.0} or greater than {@code 2.0}. - */ - private static double logGammaSum(final double a, final double b) - throws OutOfRangeException { - - if ((a < 1.0) || (a > 2.0)) { - throw new OutOfRangeException(a, 1.0, 2.0); - } - if ((b < 1.0) || (b > 2.0)) { - throw new OutOfRangeException(b, 1.0, 2.0); - } - - final double x = (a - 1.0) + (b - 1.0); - if (x <= 0.5) { - return org.apache.commons.math4.special.Gamma.logGamma1p(1.0 + x); - } else if (x <= 1.5) { - return org.apache.commons.math4.special.Gamma.logGamma1p(x) + FastMath.log1p(x); - } else { - return org.apache.commons.math4.special.Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x)); - } - } - - /** - * Returns the value of log[Î(b) / Î(a + b)] for a ⥠0 and b ⥠10. Based on - * the <em>NSWC Library of Mathematics Subroutines</em> double precision - * implementation, {@code DLGDIV}. In - * {@code BetaTest.testLogGammaMinusLogGammaSum()}, this private method is - * accessed through reflection. - * - * @param a First argument. - * @param b Second argument. - * @return the value of {@code log(Gamma(b) / Gamma(a + b))}. - * @throws NumberIsTooSmallException if {@code a < 0.0} or {@code b < 10.0}. - */ - private static double logGammaMinusLogGammaSum(final double a, - final double b) - throws NumberIsTooSmallException { - - if (a < 0.0) { - throw new NumberIsTooSmallException(a, 0.0, true); - } - if (b < 10.0) { - throw new NumberIsTooSmallException(b, 10.0, true); - } - - /* - * d = a + b - 0.5 - */ - final double d; - final double w; - if (a <= b) { - d = b + (a - 0.5); - w = deltaMinusDeltaSum(a, b); - } else { - d = a + (b - 0.5); - w = deltaMinusDeltaSum(b, a); - } - - final double u = d * FastMath.log1p(a / b); - final double v = a * (FastMath.log(b) - 1.0); - - return u <= v ? (w - u) - v : (w - v) - u; - } - - /** - * Returns the value of Î(b) - Î(a + b), with 0 ⤠a ⤠b and b ⥠10. Based - * on equations (26), (27) and (28) in Didonato and Morris (1992). - * - * @param a First argument. - * @param b Second argument. - * @return the value of {@code Delta(b) - Delta(a + b)} - * @throws OutOfRangeException if {@code a < 0} or {@code a > b} - * @throws NumberIsTooSmallException if {@code b < 10} - */ - private static double deltaMinusDeltaSum(final double a, - final double b) - throws OutOfRangeException, NumberIsTooSmallException { - - if ((a < 0) || (a > b)) { - throw new OutOfRangeException(a, 0, b); - } - if (b < 10) { - throw new NumberIsTooSmallException(b, 10, true); - } - - final double h = a / b; - final double p = h / (1.0 + h); - final double q = 1.0 / (1.0 + h); - final double q2 = q * q; - /* - * s[i] = 1 + q + ... - q**(2 * i) - */ - final double[] s = new double[DELTA.length]; - s[0] = 1.0; - for (int i = 1; i < s.length; i++) { - s[i] = 1.0 + (q + q2 * s[i - 1]); - } - /* - * w = Delta(b) - Delta(a + b) - */ - final double sqrtT = 10.0 / b; - final double t = sqrtT * sqrtT; - double w = DELTA[DELTA.length - 1] * s[s.length - 1]; - for (int i = DELTA.length - 2; i >= 0; i--) { - w = t * w + DELTA[i] * s[i]; - } - return w * p / b; - } - - /** - * Returns the value of Î(p) + Î(q) - Î(p + q), with p, q ⥠10. Based on - * the <em>NSWC Library of Mathematics Subroutines</em> double precision - * implementation, {@code DBCORR}. In - * {@code BetaTest.testSumDeltaMinusDeltaSum()}, this private method is - * accessed through reflection. - * - * @param p First argument. - * @param q Second argument. - * @return the value of {@code Delta(p) + Delta(q) - Delta(p + q)}. - * @throws NumberIsTooSmallException if {@code p < 10.0} or {@code q < 10.0}. - */ - private static double sumDeltaMinusDeltaSum(final double p, - final double q) { - - if (p < 10.0) { - throw new NumberIsTooSmallException(p, 10.0, true); - } - if (q < 10.0) { - throw new NumberIsTooSmallException(q, 10.0, true); - } - - final double a = FastMath.min(p, q); - final double b = FastMath.max(p, q); - final double sqrtT = 10.0 / a; - final double t = sqrtT * sqrtT; - double z = DELTA[DELTA.length - 1]; - for (int i = DELTA.length - 2; i >= 0; i--) { - z = t * z + DELTA[i]; - } - return z / a + deltaMinusDeltaSum(a, b); - } - - /** - * Returns the value of {@code log B(p, q)} for {@code 0 ⤠x ⤠1} and {@code p, q > 0}. Based on the - * <em>NSWC Library of Mathematics Subroutines</em> implementation, - * {@code DBETLN}. - * - * @param p First argument. - * @param q Second argument. - * @return the value of {@code log(Beta(p, q))}, {@code NaN} if - * {@code p <= 0} or {@code q <= 0}. - */ - public static double logBeta(final double p, final double q) { - if (Double.isNaN(p) || Double.isNaN(q) || (p <= 0.0) || (q <= 0.0)) { - return Double.NaN; - } - - final double a = FastMath.min(p, q); - final double b = FastMath.max(p, q); - if (a >= 10.0) { - final double w = sumDeltaMinusDeltaSum(a, b); - final double h = a / b; - final double c = h / (1.0 + h); - final double u = -(a - 0.5) * FastMath.log(c); - final double v = b * FastMath.log1p(h); - if (u <= v) { - return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - u) - v; - } else { - return (((-0.5 * FastMath.log(b) + HALF_LOG_TWO_PI) + w) - v) - u; - } - } else if (a > 2.0) { - if (b > 1000.0) { - final int n = (int) FastMath.floor(a - 1.0); - double prod = 1.0; - double ared = a; - for (int i = 0; i < n; i++) { - ared -= 1.0; - prod *= ared / (1.0 + ared / b); - } - return (FastMath.log(prod) - n * FastMath.log(b)) + - (LogGamma.value(ared) + - logGammaMinusLogGammaSum(ared, b)); - } else { - double prod1 = 1.0; - double ared = a; - while (ared > 2.0) { - ared -= 1.0; - final double h = ared / b; - prod1 *= h / (1.0 + h); - } - if (b < 10.0) { - double prod2 = 1.0; - double bred = b; - while (bred > 2.0) { - bred -= 1.0; - prod2 *= bred / (ared + bred); - } - return FastMath.log(prod1) + - FastMath.log(prod2) + - (LogGamma.value(ared) + - (LogGamma.value(bred) - - logGammaSum(ared, bred))); - } else { - return FastMath.log(prod1) + - LogGamma.value(ared) + - logGammaMinusLogGammaSum(ared, b); - } - } - } else if (a >= 1.0) { - if (b > 2.0) { - if (b < 10.0) { - double prod = 1.0; - double bred = b; - while (bred > 2.0) { - bred -= 1.0; - prod *= bred / (a + bred); - } - return FastMath.log(prod) + - (LogGamma.value(a) + - (LogGamma.value(bred) - - logGammaSum(a, bred))); - } else { - return LogGamma.value(a) + - logGammaMinusLogGammaSum(a, b); - } - } else { - return LogGamma.value(a) + - LogGamma.value(b) - - logGammaSum(a, b); - } - } else { - if (b >= 10.0) { - return LogGamma.value(a) + - logGammaMinusLogGammaSum(a, b); - } else { - // The following command is the original NSWC implementation. - // return LogGamma.value(a) + - // (LogGamma.value(b) - LogGamma.value(a + b)); - // The following command turns out to be more accurate. - return FastMath.log(Gamma.value(a) * Gamma.value(b) / - Gamma.value(a + b)); - } - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/44ab2569/src/main/java/org/apache/commons/math4/special/Gamma.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/special/Gamma.java b/src/main/java/org/apache/commons/math4/special/Gamma.java deleted file mode 100644 index 3c4ff97..0000000 --- a/src/main/java/org/apache/commons/math4/special/Gamma.java +++ /dev/null @@ -1,746 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math4.special; - -import org.apache.commons.math4.exception.MaxCountExceededException; -import org.apache.commons.math4.exception.NumberIsTooLargeException; -import org.apache.commons.math4.exception.NumberIsTooSmallException; -import org.apache.commons.math4.util.ContinuedFraction; -import org.apache.commons.math4.util.FastMath; - -/** - * <p> - * This is a utility class that provides computation methods related to the - * Γ (Gamma) family of functions. - * </p> - * <p> - * Implementation of {@link #invGamma1pm1(double)} and - * {@link #logGamma1p(double)} is based on the algorithms described in - * <ul> - * <li><a href="http://dx.doi.org/10.1145/22721.23109";>Didonato and Morris - * (1986)</a>, <em>Computation of the Incomplete Gamma Function Ratios and - * their Inverse</em>, TOMS 12(4), 377-393,</li> - * <li><a href="http://dx.doi.org/10.1145/131766.131776";>Didonato and Morris - * (1992)</a>, <em>Algorithm 708: Significant Digit Computation of the - * Incomplete Beta Function Ratios</em>, TOMS 18(3), 360-373,</li> - * </ul> - * and implemented in the - * <a href="http://www.dtic.mil/docs/citations/ADA476840";>NSWC Library of Mathematical Functions</a>, - * available - * <a href="http://www.ualberta.ca/CNS/RESEARCH/Software/NumericalNSWC/site.html";>here</a>. - * This library is "approved for public release", and the - * <a href="http://www.dtic.mil/dtic/pdf/announcements/CopyrightGuidance.pdf";>Copyright guidance</a> - * indicates that unless otherwise stated in the code, all FORTRAN functions in - * this library are license free. Since no such notice appears in the code these - * functions can safely be ported to Commons-Math. - * - */ -public class Gamma { - /** - * <a href="http://en.wikipedia.org/wiki/Euler-Mascheroni_constant";>Euler-Mascheroni constant</a> - * - * @since 2.0 - */ - public static final double GAMMA = 0.577215664901532860606512090082; - - /** - * Constant \( g = \frac{607}{128} \) in the {@link #lanczos(double) Lanczos approximation}. - * - * @since 3.1 - */ - public static final double LANCZOS_G = 607.0 / 128.0; - - /** Maximum allowed numerical error. */ - private static final double DEFAULT_EPSILON = 10e-15; - - /** Lanczos coefficients */ - private static final double[] LANCZOS = { - 0.99999999999999709182, - 57.156235665862923517, - -59.597960355475491248, - 14.136097974741747174, - -0.49191381609762019978, - .33994649984811888699e-4, - .46523628927048575665e-4, - -.98374475304879564677e-4, - .15808870322491248884e-3, - -.21026444172410488319e-3, - .21743961811521264320e-3, - -.16431810653676389022e-3, - .84418223983852743293e-4, - -.26190838401581408670e-4, - .36899182659531622704e-5, - }; - - /** Avoid repeated computation of log of 2 PI in logGamma */ - private static final double HALF_LOG_2_PI = 0.5 * FastMath.log(2.0 * FastMath.PI); - - /** The constant value of √(2π). */ - private static final double SQRT_TWO_PI = 2.506628274631000502; - - // limits for switching algorithm in digamma - /** C limit. */ - private static final double C_LIMIT = 49; - - /** S limit. */ - private static final double S_LIMIT = 1e-5; - - /* - * Constants for the computation of double invGamma1pm1(double). - * Copied from DGAM1 in the NSWC library. - */ - - /** The constant {@code A0} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_A0 = .611609510448141581788E-08; - - /** The constant {@code A1} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_A1 = .624730830116465516210E-08; - - /** The constant {@code B1} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_B1 = .203610414066806987300E+00; - - /** The constant {@code B2} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_B2 = .266205348428949217746E-01; - - /** The constant {@code B3} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_B3 = .493944979382446875238E-03; - - /** The constant {@code B4} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_B4 = -.851419432440314906588E-05; - - /** The constant {@code B5} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_B5 = -.643045481779353022248E-05; - - /** The constant {@code B6} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_B6 = .992641840672773722196E-06; - - /** The constant {@code B7} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_B7 = -.607761895722825260739E-07; - - /** The constant {@code B8} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_B8 = .195755836614639731882E-09; - - /** The constant {@code P0} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_P0 = .6116095104481415817861E-08; - - /** The constant {@code P1} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_P1 = .6871674113067198736152E-08; - - /** The constant {@code P2} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_P2 = .6820161668496170657918E-09; - - /** The constant {@code P3} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_P3 = .4686843322948848031080E-10; - - /** The constant {@code P4} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_P4 = .1572833027710446286995E-11; - - /** The constant {@code P5} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_P5 = -.1249441572276366213222E-12; - - /** The constant {@code P6} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_P6 = .4343529937408594255178E-14; - - /** The constant {@code Q1} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_Q1 = .3056961078365221025009E+00; - - /** The constant {@code Q2} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_Q2 = .5464213086042296536016E-01; - - /** The constant {@code Q3} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_Q3 = .4956830093825887312020E-02; - - /** The constant {@code Q4} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_Q4 = .2692369466186361192876E-03; - - /** The constant {@code C} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C = -.422784335098467139393487909917598E+00; - - /** The constant {@code C0} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C0 = .577215664901532860606512090082402E+00; - - /** The constant {@code C1} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C1 = -.655878071520253881077019515145390E+00; - - /** The constant {@code C2} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C2 = -.420026350340952355290039348754298E-01; - - /** The constant {@code C3} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C3 = .166538611382291489501700795102105E+00; - - /** The constant {@code C4} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C4 = -.421977345555443367482083012891874E-01; - - /** The constant {@code C5} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C5 = -.962197152787697356211492167234820E-02; - - /** The constant {@code C6} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C6 = .721894324666309954239501034044657E-02; - - /** The constant {@code C7} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C7 = -.116516759185906511211397108401839E-02; - - /** The constant {@code C8} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C8 = -.215241674114950972815729963053648E-03; - - /** The constant {@code C9} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C9 = .128050282388116186153198626328164E-03; - - /** The constant {@code C10} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C10 = -.201348547807882386556893914210218E-04; - - /** The constant {@code C11} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C11 = -.125049348214267065734535947383309E-05; - - /** The constant {@code C12} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C12 = .113302723198169588237412962033074E-05; - - /** The constant {@code C13} defined in {@code DGAM1}. */ - private static final double INV_GAMMA1P_M1_C13 = -.205633841697760710345015413002057E-06; - - /** - * Class contains only static methods. - */ - private Gamma() {} - - /** - * Computes the function \( \ln \Gamma(x) \) for \( x \gt 0 \). - * - * <p> - * For \( x \leq 8 \), the implementation is based on the double precision - * implementation in the <em>NSWC Library of Mathematics Subroutines</em>, - * {@code DGAMLN}. For \( x \geq 8 \), the implementation is based on - * </p> - * - * <ul> - * <li><a href="http://mathworld.wolfram.com/GammaFunction.html";>Gamma - * Function</a>, equation (28).</li> - * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html";> - * Lanczos Approximation</a>, equations (1) through (5).</li> - * <li><a href="http://my.fit.edu/~gabdo/gamma.txt";>Paul Godfrey, A note on - * the computation of the convergent Lanczos complex Gamma - * approximation</a></li> - * </ul> - * - * @param x Argument. - * @return \( \ln \Gamma(x) \), or {@code NaN} if {@code x <= 0}. - */ - public static double logGamma(double x) { - double ret; - - if (Double.isNaN(x) || (x <= 0.0)) { - ret = Double.NaN; - } else if (x < 0.5) { - return logGamma1p(x) - FastMath.log(x); - } else if (x <= 2.5) { - return logGamma1p((x - 0.5) - 0.5); - } else if (x <= 8.0) { - final int n = (int) FastMath.floor(x - 1.5); - double prod = 1.0; - for (int i = 1; i <= n; i++) { - prod *= x - i; - } - return logGamma1p(x - (n + 1)) + FastMath.log(prod); - } else { - double sum = lanczos(x); - double tmp = x + LANCZOS_G + .5; - ret = ((x + .5) * FastMath.log(tmp)) - tmp + - HALF_LOG_2_PI + FastMath.log(sum / x); - } - - return ret; - } - - /** - * Computes the regularized gamma function \( P(a, x) \). - * - * @param a Parameter \( a \). - * @param x Value. - * @return \( P(a, x) \) - * @throws MaxCountExceededException if the algorithm fails to converge. - */ - public static double regularizedGammaP(double a, double x) { - return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE); - } - - /** - * Computes the regularized gamma function \( P(a, x) \). - * - * The implementation of this method is based on: - * <ul> - * <li> - * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html";> - * Regularized Gamma Function</a>, equation (1) - * </li> - * <li> - * <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html";> - * Incomplete Gamma Function</a>, equation (4). - * </li> - * <li> - * <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html";> - * Confluent Hypergeometric Function of the First Kind</a>, equation (1). - * </li> - * </ul> - * - * @param a Parameter \( a \). - * @param x Argument. - * @param epsilon When the absolute value of the nth item in the - * series is less than epsilon the approximation ceases to calculate - * further elements in the series. - * @param maxIterations Maximum number of "iterations" to complete. - * @return \( P(a, x) \) - * @throws MaxCountExceededException if the algorithm fails to converge. - */ - public static double regularizedGammaP(double a, - double x, - double epsilon, - int maxIterations) { - double ret; - - if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) { - ret = Double.NaN; - } else if (x == 0.0) { - ret = 0.0; - } else if (x >= a + 1) { - // use regularizedGammaQ because it should converge faster in this - // case. - ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations); - } else { - // calculate series - double n = 0.0; // current element index - double an = 1.0 / a; // n-th element in the series - double sum = an; // partial sum - while (FastMath.abs(an/sum) > epsilon && - n < maxIterations && - sum < Double.POSITIVE_INFINITY) { - // compute next element in the series - n += 1.0; - an *= x / (a + n); - - // update partial sum - sum += an; - } - if (n >= maxIterations) { - throw new MaxCountExceededException(maxIterations); - } else if (Double.isInfinite(sum)) { - ret = 1.0; - } else { - ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * sum; - } - } - - return ret; - } - - /** - * Computes the regularized gamma function \( Q(a, x) = 1 - P(a, x) \). - * - * @param a Parameter \( a \). - * @param x Argument. - * @return \( Q(a, x) \) - * @throws MaxCountExceededException if the algorithm fails to converge. - */ - public static double regularizedGammaQ(double a, double x) { - return regularizedGammaQ(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE); - } - - /** - * Computes the regularized gamma function \( Q(a, x) = 1 - P(a, x) \). - * - * The implementation of this method is based on: - * <ul> - * <li> - * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html";> - * Regularized Gamma Function</a>, equation (1). - * </li> - * <li> - * <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/";> - * Regularized incomplete gamma function: Continued fraction representations - * (formula 06.08.10.0003)</a> - * </li> - * </ul> - * - * @param a Parameter \( a \). - * @param x Argument. - * @param epsilon When the absolute value of the nth item in the - * series is less than epsilon the approximation ceases to calculate - * further elements in the series. - * @param maxIterations Maximum number of "iterations" to complete. - * @return \( Q(a, x) \) - * @throws MaxCountExceededException if the algorithm fails to converge. - */ - public static double regularizedGammaQ(final double a, - double x, - double epsilon, - int maxIterations) { - double ret; - - if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) { - ret = Double.NaN; - } else if (x == 0.0) { - ret = 1.0; - } else if (x < a + 1.0) { - // use regularizedGammaP because it should converge faster in this - // case. - ret = 1.0 - regularizedGammaP(a, x, epsilon, maxIterations); - } else { - // create continued fraction - ContinuedFraction cf = new ContinuedFraction() { - - /** {@inheritDoc} */ - @Override - protected double getA(int n, double x) { - return ((2.0 * n) + 1.0) - a + x; - } - - /** {@inheritDoc} */ - @Override - protected double getB(int n, double x) { - return n * (a - n); - } - }; - - ret = 1.0 / cf.evaluate(x, epsilon, maxIterations); - ret = FastMath.exp(-x + (a * FastMath.log(x)) - logGamma(a)) * ret; - } - - return ret; - } - - - /** - * Computes the digamma function, defined as the logarithmic derivative - * of the \( \Gamma \) function: - * \( \frac{d}{dx}(\ln \Gamma(x)) = \frac{\Gamma^\prime(x)}{\Gamma(x)} \). - * - * <p>This is an independently written implementation of the algorithm described in - * Jose Bernardo, Algorithm AS 103: Psi (Digamma) Function, Applied Statistics, 1976. - * A <a href="https://en.wikipedia.org/wiki/Digamma_function#Reflection_formula";> - * reflection formula</a> is incorporated to improve performance on negative values.</p> - * - * <p>Some of the constants have been changed to increase accuracy at the moderate - * expense of run-time. The result should be accurate to within \( 10^{-8} \) - * relative tolerance for \( 0 \le x \le 10^{-5} \) and within \( 10^{-8} \) absolute - * tolerance otherwise.</p> - * - * @param x Argument. - * @return digamma(x) to within \( 10^{-8} \) relative or absolute error whichever is larger. - * - * @see <a href="http://en.wikipedia.org/wiki/Digamma_function";>Digamma</a> - * @see <a href="http://www.uv.es/~bernardo/1976AppStatist.pdf";>Bernardo's original article</a> - * - * @since 2.0 - */ - public static double digamma(double x) { - if (Double.isNaN(x) || Double.isInfinite(x)) { - return x; - } - - double digamma = 0.0; - if (x < 0) { - // use reflection formula to fall back into positive values - digamma -= FastMath.PI / FastMath.tan(FastMath.PI * x); - x = 1 - x; - } - - if (x > 0 && x <= S_LIMIT) { - // use method 5 from Bernardo AS103 - // accurate to O(x) - return digamma -GAMMA - 1 / x; - } - - while (x < C_LIMIT) { - digamma -= 1.0 / x; - x += 1.0; - } - - // use method 4 (accurate to O(1/x^8) - double inv = 1 / (x * x); - // 1 1 1 1 - // log(x) - --- - ------ + ------- - ------- - // 2 x 12 x^2 120 x^4 252 x^6 - digamma += FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252)); - - return digamma; - } - - /** - * Computes the trigamma function \( \psi_1(x) = \frac{d^2}{dx^2} (\ln \Gamma(x)) \). - * <p> - * This function is the derivative of the {@link #digamma(double) digamma function}. - * </p> - * - * @param x Argument. - * @return \( \psi_1(x) \) to within \( 10^{-8} \) relative or absolute - * error whichever is smaller - * - * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function";>Trigamma</a> - * @see #digamma(double) - * - * @since 2.0 - */ - public static double trigamma(double x) { - if (Double.isNaN(x) || Double.isInfinite(x)) { - return x; - } - - if (x > 0 && x <= S_LIMIT) { - return 1 / (x * x); - } - - if (x >= C_LIMIT) { - double inv = 1 / (x * x); - // 1 1 1 1 1 - // - + ---- + ---- - ----- + ----- - // x 2 3 5 7 - // 2 x 6 x 30 x 42 x - return 1 / x + inv / 2 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv / 42)); - } - - return trigamma(x + 1) + 1 / (x * x); - } - - /** - * Computes the Lanczos approximation used to compute the gamma function. - * - * <p> - * The Lanczos approximation is related to the Gamma function by the - * following equation - * \[ - * \Gamma(x) = \sqrt{2\pi} \, \frac{(g + x + \frac{1}{2})^{x + \frac{1}{2}} \, e^{-(g + x + \frac{1}{2})} \, \mathrm{lanczos}(x)} - * {x} - * \] - * where \(g\) is the Lanczos constant. - * </p> - * - * @param x Argument. - * @return The Lanczos approximation. - * - * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html";>Lanczos Approximation</a> - * equations (1) through (5), and Paul Godfrey's - * <a href="http://my.fit.edu/~gabdo/gamma.txt";>Note on the computation - * of the convergent Lanczos complex Gamma approximation</a> - * - * @since 3.1 - */ - public static double lanczos(final double x) { - double sum = 0.0; - for (int i = LANCZOS.length - 1; i > 0; --i) { - sum += LANCZOS[i] / (x + i); - } - return sum + LANCZOS[0]; - } - - /** - * Computes the function \( \frac{1}{\Gamma(1 + x)} - 1 \) for \( -0.5 \leq x \leq 1.5 \). - * <p> - * This implementation is based on the double precision implementation in - * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGAM1}. - * </p> - * - * @param x Argument. - * @return \( \frac{1}{\Gamma(1 + x)} - 1 \) - * @throws NumberIsTooSmallException if {@code x < -0.5} - * @throws NumberIsTooLargeException if {@code x > 1.5} - * - * @since 3.1 - */ - public static double invGamma1pm1(final double x) { - - if (x < -0.5) { - throw new NumberIsTooSmallException(x, -0.5, true); - } - if (x > 1.5) { - throw new NumberIsTooLargeException(x, 1.5, true); - } - - final double ret; - final double t = x <= 0.5 ? x : (x - 0.5) - 0.5; - if (t < 0.0) { - final double a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1; - double b = INV_GAMMA1P_M1_B8; - b = INV_GAMMA1P_M1_B7 + t * b; - b = INV_GAMMA1P_M1_B6 + t * b; - b = INV_GAMMA1P_M1_B5 + t * b; - b = INV_GAMMA1P_M1_B4 + t * b; - b = INV_GAMMA1P_M1_B3 + t * b; - b = INV_GAMMA1P_M1_B2 + t * b; - b = INV_GAMMA1P_M1_B1 + t * b; - b = 1.0 + t * b; - - double c = INV_GAMMA1P_M1_C13 + t * (a / b); - c = INV_GAMMA1P_M1_C12 + t * c; - c = INV_GAMMA1P_M1_C11 + t * c; - c = INV_GAMMA1P_M1_C10 + t * c; - c = INV_GAMMA1P_M1_C9 + t * c; - c = INV_GAMMA1P_M1_C8 + t * c; - c = INV_GAMMA1P_M1_C7 + t * c; - c = INV_GAMMA1P_M1_C6 + t * c; - c = INV_GAMMA1P_M1_C5 + t * c; - c = INV_GAMMA1P_M1_C4 + t * c; - c = INV_GAMMA1P_M1_C3 + t * c; - c = INV_GAMMA1P_M1_C2 + t * c; - c = INV_GAMMA1P_M1_C1 + t * c; - c = INV_GAMMA1P_M1_C + t * c; - if (x > 0.5) { - ret = t * c / x; - } else { - ret = x * ((c + 0.5) + 0.5); - } - } else { - double p = INV_GAMMA1P_M1_P6; - p = INV_GAMMA1P_M1_P5 + t * p; - p = INV_GAMMA1P_M1_P4 + t * p; - p = INV_GAMMA1P_M1_P3 + t * p; - p = INV_GAMMA1P_M1_P2 + t * p; - p = INV_GAMMA1P_M1_P1 + t * p; - p = INV_GAMMA1P_M1_P0 + t * p; - - double q = INV_GAMMA1P_M1_Q4; - q = INV_GAMMA1P_M1_Q3 + t * q; - q = INV_GAMMA1P_M1_Q2 + t * q; - q = INV_GAMMA1P_M1_Q1 + t * q; - q = 1.0 + t * q; - - double c = INV_GAMMA1P_M1_C13 + (p / q) * t; - c = INV_GAMMA1P_M1_C12 + t * c; - c = INV_GAMMA1P_M1_C11 + t * c; - c = INV_GAMMA1P_M1_C10 + t * c; - c = INV_GAMMA1P_M1_C9 + t * c; - c = INV_GAMMA1P_M1_C8 + t * c; - c = INV_GAMMA1P_M1_C7 + t * c; - c = INV_GAMMA1P_M1_C6 + t * c; - c = INV_GAMMA1P_M1_C5 + t * c; - c = INV_GAMMA1P_M1_C4 + t * c; - c = INV_GAMMA1P_M1_C3 + t * c; - c = INV_GAMMA1P_M1_C2 + t * c; - c = INV_GAMMA1P_M1_C1 + t * c; - c = INV_GAMMA1P_M1_C0 + t * c; - - if (x > 0.5) { - ret = (t / x) * ((c - 0.5) - 0.5); - } else { - ret = x * c; - } - } - - return ret; - } - - /** - * Computes the function \( \ln \Gamma(1 + x) \) for \( -0.5 \leq x \leq 1.5 \). - * <p> - * This implementation is based on the double precision implementation in - * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}. - * </p> - * - * @param x Argument. - * @return \( \ln \Gamma(1 + x) \) - * @throws NumberIsTooSmallException if {@code x < -0.5}. - * @throws NumberIsTooLargeException if {@code x > 1.5}. - * @since 3.1 - */ - public static double logGamma1p(final double x) - throws NumberIsTooSmallException, NumberIsTooLargeException { - - if (x < -0.5) { - throw new NumberIsTooSmallException(x, -0.5, true); - } - if (x > 1.5) { - throw new NumberIsTooLargeException(x, 1.5, true); - } - - return -FastMath.log1p(invGamma1pm1(x)); - } - - - /** - * Computes the value of \( \Gamma(x) \). - * <p> - * Based on the <em>NSWC Library of Mathematics Subroutines</em> double - * precision implementation, {@code DGAMMA}. - * </p> - * - * @param x Argument. - * @return \( \Gamma(x) \) - * - * @since 3.1 - */ - public static double gamma(final double x) { - - if ((x == FastMath.rint(x)) && (x <= 0.0)) { - return Double.NaN; - } - - final double ret; - final double absX = FastMath.abs(x); - if (absX <= 20.0) { - if (x >= 1.0) { - /* - * From the recurrence relation - * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n), - * then - * Gamma(t) = 1 / [1 + invGamma1pm1(t - 1)], - * where t = x - n. This means that t must satisfy - * -0.5 <= t - 1 <= 1.5. - */ - double prod = 1.0; - double t = x; - while (t > 2.5) { - t -= 1.0; - prod *= t; - } - ret = prod / (1.0 + invGamma1pm1(t - 1.0)); - } else { - /* - * From the recurrence relation - * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)] - * then - * Gamma(x + n + 1) = 1 / [1 + invGamma1pm1(x + n)], - * which requires -0.5 <= x + n <= 1.5. - */ - double prod = x; - double t = x; - while (t < -0.5) { - t += 1.0; - prod *= t; - } - ret = 1.0 / (prod * (1.0 + invGamma1pm1(t))); - } - } else { - final double y = absX + LANCZOS_G + 0.5; - final double gammaAbs = SQRT_TWO_PI / absX * - FastMath.pow(y, absX + 0.5) * - FastMath.exp(-y) * lanczos(absX); - if (x > 0.0) { - ret = gammaAbs; - } else { - /* - * From the reflection formula - * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi, - * and the recurrence relation - * Gamma(1 - x) = -x * Gamma(-x), - * it is found - * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)]. - */ - ret = -FastMath.PI / - (x * FastMath.sin(FastMath.PI * x) * gammaAbs); - } - } - return ret; - } -}
