Class BoostGamma

java.lang.Object
org.apache.commons.numbers.gamma.BoostGamma

final class BoostGamma extends Object
Implementation of the Regularized Gamma functions and Incomplete Gamma functions.

This code has been adapted from the Boost c++ implementation <boost/math/special_functions/gamma.hpp>. All work is copyright to the original authors and subject to the Boost Software License.

See Also:
  • Nested Class Summary

    Nested Classes
    Modifier and Type
    Class
    Description
    (package private) static final class 
    53-bit precision implementation of the Lanczos approximation.
  • Field Summary

    Fields
    Modifier and Type
    Field
    Description
    private static final double
    Default epsilon value for relative error.
    private static final double
    Euler's constant.
    private static final double[]
    All factorials that can be represented as a double.
    private static final int
    The threshold value for choosing the Lanczos approximation.
    private static final int
    Approximate value for ln(Double.MAX_VALUE).
    private static final int
    Approximate value for ln(Double.MIN_VALUE).
    private static final double
    ln(pi).
    private static final double
    ln(sqrt(2 pi)).
    private static final int
    The largest factorial that can be represented as a double.
    private static final int
    The largest integer value for gamma(z) that can be represented as a double.
    private static final double
    Value for the sqrt of the epsilon for relative error.
    private static final double
    2^53.
  • Constructor Summary

    Constructors
    Modifier
    Constructor
    Description
    private
    Private constructor.
  • Method Summary

    Modifier and Type
    Method
    Description
    private static double
    finiteGammaQ(double a, double x)
    Upper gamma fraction for integer a.
    private static double
    finiteHalfGammaQ(double a, double x)
    Upper gamma fraction for half integer a.
    (package private) static double
    fullIgammaPrefix(double a, double z)
    Calculate power term prefix (z^a)(e^-z) used in the non-normalised incomplete gammas.
    private static double
    gammaIncompleteImp(double a, double x, boolean normalised, boolean invert, Policy pol)
    Main incomplete gamma entry point, handles all four incomplete gammas.
    (package private) static double
    gammaP(double a, double x)
    Regularised lower incomplete gamma.
    (package private) static double
    gammaP(double a, double x, Policy policy)
    Regularised lower incomplete gamma.
    (package private) static double
    gammaPDerivative(double a, double x)
    Derivative of the regularised lower incomplete gamma.
    (package private) static double
    gammaQ(double a, double x)
    Regularised upper incomplete gamma.
    (package private) static double
    gammaQ(double a, double x, Policy policy)
    Regularised upper incomplete gamma.
    (package private) static double[]
    All factorials that are representable as a double.
    (package private) static double
    igammaTemmeLarge(double a, double x)
    Implements the asymptotic expansions of the incomplete gamma functions P(a, x) and Q(a, x), used when a is large and x ~ a.
    (package private) static double
    incompleteTgammaLargeX(double a, double x, Policy pol)
    Incomplete tgamma for large X.
    private static boolean
    isOdd(double v)
    Checks if the value is odd.
    (package private) static double
    lgamma(double z)
    Log Gamma function.
    (package private) static double
    lgamma(double z, int[] sign)
    Log Gamma function.
    private static double
    lgammaSmall(double z, double zm1, double zm2)
    Log Gamma function for small z.
    (package private) static double
    lowerGammaSeries(double a, double z, double initValue, Policy pol)
    Lower gamma series.
    private static boolean
    nonZeroLength(int[] array)
    Return true if the array is not null and has non-zero length.
    (package private) static double
    regularisedGammaPrefix(double a, double z)
    Compute (z^a)(e^-z)/tgamma(a).
    (package private) static double
    sinpx(double x)
    Ad hoc function calculates x * sin(pi * x), taking extra care near when x is near a whole number.
    (package private) static double
    tgamma(double z)
    Gamma function.
    (package private) static double
    tgamma(double a, double x)
    Full upper incomplete gamma.
    (package private) static double
    tgamma(double a, double x, Policy policy)
    Full upper incomplete gamma.
    (package private) static double
    tgamma1pm1(double dz)
    Calculates tgamma(1+dz)-1.
    (package private) static double
    tgammaDeltaRatio(double z, double delta)
    Ratio of gamma functions.
    private static double
    tgammaDeltaRatioImpLanczos(double z, double delta)
    Ratio of gamma functions using Lanczos support.
    (package private) static double
    tgammaLower(double a, double x)
    Full lower incomplete gamma.
    (package private) static double
    tgammaLower(double a, double x, Policy policy)
    Full lower incomplete gamma.
    (package private) static double
    tgammaRatio(double x, double y)
    Ratio of gamma functions.
    private static double
    tgammaSmallUpperPart(double a, double x, Policy pol, double[] pgam, boolean invert)
    Upper gamma fraction for very small a.
    (package private) static double
    Returns the factorial of n.
    (package private) static double
    upperGammaFraction(double a, double z, Policy pol)
    Upper gamma fraction.

    Methods inherited from class java.lang.Object

    clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
  • Field Details

    • EPSILON

      private static final double EPSILON
      Default epsilon value for relative error. This is equal to the Boost constant boost::math::tools::epsilon<double>().
      See Also:
    • ROOT_EPSILON

      private static final double ROOT_EPSILON
      Value for the sqrt of the epsilon for relative error. This is equal to the Boost constant boost::math::tools::root_epsilon<double>().
      See Also:
    • LOG_MAX_VALUE

      private static final int LOG_MAX_VALUE
      Approximate value for ln(Double.MAX_VALUE). This is equal to the Boost constant boost::math::tools::log_max_value<double>(). No term x should be used in exp(x) if x > LOG_MAX_VALUE to avoid overflow.
      See Also:
    • LOG_MIN_VALUE

      private static final int LOG_MIN_VALUE
      Approximate value for ln(Double.MIN_VALUE). This is equal to the Boost constant boost::math::tools::log_min_value<double>(). No term x should be used in exp(x) if x < LOG_MIN_VALUE to avoid underflow to sub-normal or zero.
      See Also:
    • MAX_FACTORIAL

      private static final int MAX_FACTORIAL
      The largest factorial that can be represented as a double. This is equal to the Boost constant boost::math::max_factorial<double>::value.
      See Also:
    • MAX_GAMMA_Z

      private static final int MAX_GAMMA_Z
      The largest integer value for gamma(z) that can be represented as a double.
      See Also:
    • LOG_ROOT_TWO_PI

      private static final double LOG_ROOT_TWO_PI
      ln(sqrt(2 pi)). Computed to 25-digits precision.
      See Also:
    • LOG_PI

      private static final double LOG_PI
      ln(pi). Computed to 25-digits precision.
      See Also:
    • EULER

      private static final double EULER
      Euler's constant.
      See Also:
    • LANCZOS_THRESHOLD

      private static final int LANCZOS_THRESHOLD
      The threshold value for choosing the Lanczos approximation.
      See Also:
    • TWO_POW_53

      private static final double TWO_POW_53
      2^53.
      See Also:
    • FACTORIAL

      private static final double[] FACTORIAL
      All factorials that can be represented as a double. Size = 171.
  • Constructor Details

    • BoostGamma

      private BoostGamma()
      Private constructor.
  • Method Details

    • getFactorials

      static double[] getFactorials()
      All factorials that are representable as a double. This data is exposed for testing.
      Returns:
      factorials
    • uncheckedFactorial

      static double uncheckedFactorial(int n)
      Returns the factorial of n. This is unchecked as an index out of bound exception will occur if the value n is not within the range [0, 170]. This function is exposed for use in BoostBeta.
      Parameters:
      n - Argument n (must be in [0, 170])
      Returns:
      n!
    • tgamma

      static double tgamma(double z)
      Gamma function.

      For small z this is based on the NSWC Library of Mathematics Subroutines double precision implementation, DGAMMA.

      For large z this is an implementation of the Boost C++ tgamma function with Lanczos support.

      Integers are handled using a look-up table of factorials.

      Note: The Boost C++ implementation uses the Lanczos sum for all z. When promotion of double to long double is not available this has larger errors than the double precision specific NSWC implementation. For larger z the Boost C++ Lanczos implementation incorporates the sqrt(2 pi) factor and has lower error than the implementation using the LanczosApproximation class.

      Parameters:
      z - Argument z
      Returns:
      gamma value
    • sinpx

      static double sinpx(double x)
      Ad hoc function calculates x * sin(pi * x), taking extra care near when x is near a whole number.
      Parameters:
      x - Value (assumed to be negative)
      Returns:
      x * sin(pi * x)
    • isOdd

      private static boolean isOdd(double v)
      Checks if the value is odd.
      Parameters:
      v - Value (assumed to be positive and an integer)
      Returns:
      true if odd
    • lgamma

      static double lgamma(double z)
      Log Gamma function. Defined as the natural logarithm of the absolute value of tgamma(z).
      Parameters:
      z - Argument z
      Returns:
      log gamma value
    • lgamma

      static double lgamma(double z, int[] sign)
      Log Gamma function. Defined as the natural logarithm of the absolute value of tgamma(z).
      Parameters:
      z - Argument z
      sign - If a non-zero length array the first index is set on output to the sign of tgamma(z)
      Returns:
      log gamma value
    • lgammaSmall

      private static double lgammaSmall(double z, double zm1, double zm2)
      Log Gamma function for small z.
      Parameters:
      z - Argument z
      zm1 - z - 1
      zm2 - z - 2
      Returns:
      log gamma value
    • tgamma1pm1

      static double tgamma1pm1(double dz)
      Calculates tgamma(1+dz)-1.
      Parameters:
      dz - Argument
      Returns:
      tgamma(1+dz)-1
    • tgamma

      static double tgamma(double a, double x)
      Full upper incomplete gamma.
      Parameters:
      a - Argument a
      x - Argument x
      Returns:
      upper gamma value
    • tgamma

      static double tgamma(double a, double x, Policy policy)
      Full upper incomplete gamma.
      Parameters:
      a - Argument a
      x - Argument x
      policy - Function evaluation policy
      Returns:
      upper gamma value
    • tgammaLower

      static double tgammaLower(double a, double x)
      Full lower incomplete gamma.
      Parameters:
      a - Argument a
      x - Argument x
      Returns:
      lower gamma value
    • tgammaLower

      static double tgammaLower(double a, double x, Policy policy)
      Full lower incomplete gamma.
      Parameters:
      a - Argument a
      x - Argument x
      policy - Function evaluation policy
      Returns:
      lower gamma value
    • gammaQ

      static double gammaQ(double a, double x)
      Regularised upper incomplete gamma.
      Parameters:
      a - Argument a
      x - Argument x
      Returns:
      q
    • gammaQ

      static double gammaQ(double a, double x, Policy policy)
      Regularised upper incomplete gamma.
      Parameters:
      a - Argument a
      x - Argument x
      policy - Function evaluation policy
      Returns:
      q
    • gammaP

      static double gammaP(double a, double x)
      Regularised lower incomplete gamma.
      Parameters:
      a - Argument a
      x - Argument x
      Returns:
      p
    • gammaP

      static double gammaP(double a, double x, Policy policy)
      Regularised lower incomplete gamma.
      Parameters:
      a - Argument a
      x - Argument x
      policy - Function evaluation policy
      Returns:
      p
    • gammaPDerivative

      static double gammaPDerivative(double a, double x)
      Derivative of the regularised lower incomplete gamma.

      \( \frac{e^{-x} x^{a-1}}{\Gamma(a)} \)

      Adapted from boost::math::detail::gamma_p_derivative_imp

      Parameters:
      a - Argument a
      x - Argument x
      Returns:
      p derivative
    • gammaIncompleteImp

      private static double gammaIncompleteImp(double a, double x, boolean normalised, boolean invert, Policy pol)
      Main incomplete gamma entry point, handles all four incomplete gammas. Adapted from boost::math::detail::gamma_incomplete_imp.

      The Boost code has a pointer p_derivative that can be set to the value of the derivative. This is used for the inverse incomplete gamma functions gamma_(p|q)_inv_imp. It is not required for the forward evaluation functions.

      Parameters:
      a - Argument a
      x - Argument x
      normalised - true to compute the regularised value
      invert - true to compute the upper value Q (default is lower value P)
      pol - Function evaluation policy
      Returns:
      gamma value
    • upperGammaFraction

      static double upperGammaFraction(double a, double z, Policy pol)
      Upper gamma fraction. Multiply result by z^a * e^-z to get the full upper incomplete integral. Divide by tgamma(z) to normalise.
      Parameters:
      a - Argument a
      z - Argument z
      pol - Function evaluation policy
      Returns:
      upper gamma fraction
    • finiteGammaQ

      private static double finiteGammaQ(double a, double x)
      Upper gamma fraction for integer a. Called when a invalid input: '<' 30 and -x > LOG_MIN_VALUE.
      Parameters:
      a - Argument a (assumed to be small)
      x - Argument x
      Returns:
      upper gamma fraction
    • finiteHalfGammaQ

      private static double finiteHalfGammaQ(double a, double x)
      Upper gamma fraction for half integer a. Called when a invalid input: '<' 30 and -x > LOG_MIN_VALUE.
      Parameters:
      a - Argument a (assumed to be small)
      x - Argument x
      Returns:
      upper gamma fraction
    • lowerGammaSeries

      static double lowerGammaSeries(double a, double z, double initValue, Policy pol)
      Lower gamma series. Multiply result by ((z^a) * (e^-z) / a) to get the full lower incomplete integral. Then divide by tgamma(a) to get the normalised value.
      Parameters:
      a - Argument a
      z - Argument z
      initValue - Initial value
      pol - Function evaluation policy
      Returns:
      lower gamma series
    • tgammaSmallUpperPart

      private static double tgammaSmallUpperPart(double a, double x, Policy pol, double[] pgam, boolean invert)
      Upper gamma fraction for very small a.
      Parameters:
      a - Argument a (assumed to be small)
      x - Argument x
      pol - Function evaluation policy
      pgam - set to value of gamma(a) on output
      invert - true to invert the result
      Returns:
      upper gamma fraction
    • fullIgammaPrefix

      static double fullIgammaPrefix(double a, double z)
      Calculate power term prefix (z^a)(e^-z) used in the non-normalised incomplete gammas.
      Parameters:
      a - Argument a
      z - Argument z
      Returns:
      incomplete gamma prefix
    • regularisedGammaPrefix

      static double regularisedGammaPrefix(double a, double z)
      Compute (z^a)(e^-z)/tgamma(a).

      Most of the error occurs in this function.

      Parameters:
      a - Argument a
      z - Argument z
      Returns:
      regularized gamma prefix
    • igammaTemmeLarge

      static double igammaTemmeLarge(double a, double x)
      Implements the asymptotic expansions of the incomplete gamma functions P(a, x) and Q(a, x), used when a is large and x ~ a.

      The primary reference is:

       "The Asymptotic Expansion of the Incomplete Gamma Functions"
       N. M. Temme.
       Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
       

      A different way of evaluating these expansions, plus a lot of very useful background information is in:

       "A Set of Algorithms For the Incomplete Gamma Functions."
       N. M. Temme.
       Probability in the Engineering and Informational Sciences,
       8, 1994, 291.
       

      An alternative implementation is in:

       "Computation of the Incomplete Gamma Function Ratios and their Inverse."
       A. R. Didonato and A. H. Morris.
       ACM TOMS, Vol 12, No 4, Dec 1986, p377.
       

      This is a port of the function accurate for 53-bit mantissas (IEEE double precision or 10^-17). To understand the code, refer to Didonato and Morris, from Eq 17 and 18 onwards.

      The coefficients used here are not taken from Didonato and Morris: the domain over which these expansions are used is slightly different to theirs, and their constants are not quite accurate enough for 128-bit long doubles. Instead the coefficients were calculated using the methods described by Temme p762 from Eq 3.8 onwards. The values obtained agree with those obtained by Didonato and Morris (at least to the first 30 digits that they provide). At double precision the degrees of polynomial required for full machine precision are close to those recommended to Didonato and Morris, but of course many more terms are needed for larger types.

      Adapted from boost/math/special_functions/detail/igamma_large.hpp.

      Parameters:
      a - the a
      x - the x
      Returns:
      the double
    • incompleteTgammaLargeX

      static double incompleteTgammaLargeX(double a, double x, Policy pol)
      Incomplete tgamma for large X.

      This summation is a source of error as the series starts at 1 and descends to zero. It can have thousands of iterations when a and z are large and close in value.

      Parameters:
      a - Argument a
      x - Argument x
      pol - Function evaluation policy
      Returns:
      incomplete tgamma
    • nonZeroLength

      private static boolean nonZeroLength(int[] array)
      Return true if the array is not null and has non-zero length.
      Parameters:
      array - Array
      Returns:
      true if a non-zero length array
    • tgammaDeltaRatio

      static double tgammaDeltaRatio(double z, double delta)
      Ratio of gamma functions.

      \[ tgamma_ratio(z, delta) = \frac{\Gamma(z)}{\Gamma(z + delta)} \]

      Adapted from tgamma_delta_ratio_imp. The use of max_factorial<double>::value == 170 has been replaced with MAX_GAMMA_Z == 171. This threshold is used when it is possible to call the gamma function without overflow.

      Parameters:
      z - Argument z
      delta - The difference
      Returns:
      gamma ratio
    • tgammaDeltaRatioImpLanczos

      private static double tgammaDeltaRatioImpLanczos(double z, double delta)
      Ratio of gamma functions using Lanczos support.

      \[ tgamma_delta_ratio(z, delta) = \frac{\Gamma(z)}{\Gamma(z + delta)}

      Adapted from tgamma_delta_ratio_imp_lanczos. The use of max_factorial<double>::value == 170 has been replaced with MAX_GAMMA_Z == 171. This threshold is used when it is possible to use the precomputed factorial table.

      Parameters:
      z - Argument z
      delta - The difference
      Returns:
      gamma ratio
    • tgammaRatio

      static double tgammaRatio(double x, double y)
      Ratio of gamma functions.

      \[ tgamma_ratio(x, y) = \frac{\Gamma(x)}{\Gamma(y)}

      Adapted from tgamma_ratio_imp. The use of max_factorial<double>::value == 170 has been replaced with MAX_GAMMA_Z == 171. This threshold is used when it is possible to call the gamma function without overflow.

      Parameters:
      x - Argument x (must be positive finite)
      y - Argument y (must be positive finite)
      Returns:
      gamma ratio (or nan)