Class BoostGamma

    • Nested Class Summary

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

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

      Constructors 
      Modifier Constructor Description
      private BoostGamma()
      Private constructor.
    • Method Summary

      All Methods Static Methods Concrete Methods 
      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[] getFactorials()
      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 uncheckedFactorial​(int n)
      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 Detail

      • 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:
        Constant Field Values
      • 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:
        Constant Field Values
      • 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:
        Constant Field Values
      • 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:
        Constant Field Values
      • 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:
        Constant Field Values
      • 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:
        Constant Field Values
      • LOG_ROOT_TWO_PI

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

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

        private static final int LANCZOS_THRESHOLD
        The threshold value for choosing the Lanczos approximation.
        See Also:
        Constant Field Values
      • FACTORIAL

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

      • BoostGamma

        private BoostGamma()
        Private constructor.
    • Method Detail

      • 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 < 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 < 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)