nd4j-common/src/main/java/org/nd4j/linalg/util/BigDecimalMath.java
static public BigDecimal Gamma(final BigDecimal x) {
        /* reduce to interval near 1.0 with the functional relation, Abramowitz-Stegun 6.1.33
         */
        if (x.compareTo(BigDecimal.ZERO) < 0) {
            return divideRound(Gamma(x.add(BigDecimal.ONE)), x);
        } else if (x.doubleValue() > 1.5) {
            /* Gamma(x) = Gamma(xmin+n) = Gamma(xmin)*Pochhammer(xmin,n).
             */
            int n = (int) (x.doubleValue() - 0.5);
            BigDecimal xmin1 = x.subtract(new BigDecimal(n));


            return multiplyRound(Gamma(xmin1), pochhammer(xmin1, n));


        } else {
            /* apply Abramowitz-Stegun 6.1.33
             */
            BigDecimal z = x.subtract(BigDecimal.ONE);
            /* add intermediately 2 digits to the partial sum accumulation
             */
            z = scalePrec(z, 2);
            MathContext mcloc = new MathContext(z.precision());
            /* measure of the absolute error is the relative error in the first, logarithmic term
             */


            double eps = x.ulp().doubleValue() / x.doubleValue();
            BigDecimal resul = log(scalePrec(x, 2)).negate();


            if (x.compareTo(BigDecimal.ONE) != 0) {
                BigDecimal gammCompl = BigDecimal.ONE.subtract(gamma(mcloc));
                resul = resul.add(multiplyRound(z, gammCompl));


                for (int n = 2;; n++) {
                    /* multiplying z^n/n by zeta(n-1) means that the two relative errors add.
                     * so the requirement in the relative error of zeta(n)-1 is that this is somewhat
                     * smaller than the relative error in z^n/n (the absolute error of thelatter is the
                     * absolute error in z)
                     */
                    BigDecimal c = divideRound(z.pow(n, mcloc), n);
                    MathContext m = new MathContext(err2prec(n * z.ulp().doubleValue() / 2. / z.doubleValue()));
                    c = c.round(m);
                    /* At larger n, zeta(n)-1 is roughly 1/2^n. The product is c/2^n.
                     * The relative error in c is c.ulp/2/c . The error in the product should be small versus eps/10.
                     * Error from 1/2^n is c*err(sigma-1).
                     * We need a relative error of zeta-1 of the order of c.ulp/50/c. This is an absolute
                     * error in zeta-1 of c.ulp/50/c/2^n, and also the absolute error in zeta, because zeta is
                     * of the order of 1.
                     */


                    if (eps / 100. / c.doubleValue() < 0.01) {
                        m = new MathContext(err2prec(eps / 100. / c.doubleValue()));
                    } else {
                        m = new MathContext(2);
                    }
                    /* zeta(n) -1 */
                    BigDecimal zetm1 = zeta(n, m).subtract(BigDecimal.ONE);
                    c = multiplyRound(c, zetm1);


                    if (n % 2 == 0) {
                        resul = resul.add(c);
                    } else {
                        resul = resul.subtract(c);
                    }
                    /* alternating sum, so truncating as eps is reached suffices
                     */


                    if (Math.abs(c.doubleValue()) < eps) {
                        break;
                    }


                }
            }
            /* The relative error in the result is the absolute error in the
             * input variable times the digamma (psi) value at that point.
             */
            double psi = 0.5772156649;


            double zdbl = z.doubleValue();


            for (int n = 1; n < 5; n++) {
                psi += zdbl / n / (n + zdbl);
            }
            eps = psi * x.ulp().doubleValue() / 2.;
            mcloc = new MathContext(err2prec(eps));


            return exp(resul).round(mcloc);


        }
    }
static public BigDecimal Gamma(final BigDecimal x) {
        /* reduce to interval near 1.0 with the functional relation, Abramowitz-Stegun 6.1.33
         */
        if (x.compareTo(BigDecimal.ZERO) < 0) {
            return divideRound(Gamma(x.add(BigDecimal.ONE)), x);
        } else if (x.doubleValue() > 1.5) {
            /* Gamma(x) = Gamma(xmin+n) = Gamma(xmin)*Pochhammer(xmin,n).
             */
            int n = (int) (x.doubleValue() - 0.5);
            BigDecimal xmin1 = x.subtract(new BigDecimal(n));


            return multiplyRound(Gamma(xmin1), pochhammer(xmin1, n));


        } else {
            /* apply Abramowitz-Stegun 6.1.33
             */
            BigDecimal z = x.subtract(BigDecimal.ONE);
            /* add intermediately 2 digits to the partial sum accumulation
             */
            z = scalePrec(z, 2);
            MathContext mcloc = new MathContext(z.precision());
            /* measure of the absolute error is the relative error in the first, logarithmic term
             */


            double eps = x.ulp().doubleValue() / x.doubleValue();
            BigDecimal resul = log(scalePrec(x, 2)).negate();


            if (x.compareTo(BigDecimal.ONE) != 0) {
                BigDecimal gammCompl = BigDecimal.ONE.subtract(gamma(mcloc));
                resul = resul.add(multiplyRound(z, gammCompl));


                for (int n = 2;; n++) {
                    /* multiplying z^n/n by zeta(n-1) means that the two relative errors add.
                     * so the requirement in the relative error of zeta(n)-1 is that this is somewhat
                     * smaller than the relative error in z^n/n (the absolute error of thelatter is the
                     * absolute error in z)
                     */
                    BigDecimal c = divideRound(z.pow(n, mcloc), n);
                    MathContext m = new MathContext(err2prec(n * z.ulp().doubleValue() / 2. / z.doubleValue()));
                    c = c.round(m);
                    /* At larger n, zeta(n)-1 is roughly 1/2^n. The product is c/2^n.
                     * The relative error in c is c.ulp/2/c . The error in the product should be small versus eps/10.
                     * Error from 1/2^n is c*err(sigma-1).
                     * We need a relative error of zeta-1 of the order of c.ulp/50/c. This is an absolute
                     * error in zeta-1 of c.ulp/50/c/2^n, and also the absolute error in zeta, because zeta is
                     * of the order of 1.
                     */


                    if (eps / 100. / c.doubleValue() < 0.01) {
                        m = new MathContext(err2prec(eps / 100. / c.doubleValue()));
                    } else {
                        m = new MathContext(2);
                    }
                    /* zeta(n) -1 */
                    BigDecimal zetm1 = zeta(n, m).subtract(BigDecimal.ONE);
                    c = multiplyRound(c, zetm1);


                    if (n % 2 == 0) {
                        resul = resul.add(c);
                    } else {
                        resul = resul.subtract(c);
                    }
                    /* alternating sum, so truncating as eps is reached suffices
                     */


                    if (Math.abs(c.doubleValue()) < eps) {
                        break;
                    }


                }
            }
            /* The relative error in the result is the absolute error in the
             * input variable times the digamma (psi) value at that point.
             */
            double psi = 0.5772156649;


            double zdbl = z.doubleValue();


            for (int n = 1; n < 5; n++) {
                psi += zdbl / n / (n + zdbl);
            }
            eps = psi * x.ulp().doubleValue() / 2.;
            mcloc = new MathContext(err2prec(eps));


            return exp(resul).round(mcloc);


        }
    }
static public BigDecimal Gamma ( final BigDecimal x ) { /* reduce to interval near 1.0 with the functional relation, Abramowitz-Stegun 6.1.33
         */ if ( x . compareTo ( BigDecimal . ZERO ) < 0 ) { return divideRound ( Gamma ( x . add ( BigDecimal . ONE ) ) , x ) ; } else if ( x . doubleValue ( ) > 1.5 ) { /* Gamma(x) = Gamma(xmin+n) = Gamma(xmin)*Pochhammer(xmin,n).
             */ int n = ( int ) ( x . doubleValue ( ) - 0.5 ) ; BigDecimal xmin1 = x . subtract ( new BigDecimal ( n ) ) ; return multiplyRound ( Gamma ( xmin1 ) , pochhammer ( xmin1 , n ) ) ; } else { /* apply Abramowitz-Stegun 6.1.33
             */ BigDecimal z = x . subtract ( BigDecimal . ONE ) ; /* add intermediately 2 digits to the partial sum accumulation
             */ z = scalePrec ( z , 2 ) ; MathContext mcloc = new MathContext ( z . precision ( ) ) ; /* measure of the absolute error is the relative error in the first, logarithmic term
             */ double eps = x . ulp ( ) . doubleValue ( ) / x . doubleValue ( ) ; BigDecimal resul = log ( scalePrec ( x , 2 ) ) . negate ( ) ; if ( x . compareTo ( BigDecimal . ONE ) != 0 ) { BigDecimal gammCompl = BigDecimal . ONE . subtract ( gamma ( mcloc ) ) ; resul = resul . add ( multiplyRound ( z , gammCompl ) ) ; for ( int n = 2 ; ; n ++ ) { /* multiplying z^n/n by zeta(n-1) means that the two relative errors add.
                     * so the requirement in the relative error of zeta(n)-1 is that this is somewhat
                     * smaller than the relative error in z^n/n (the absolute error of thelatter is the
                     * absolute error in z)
                     */ BigDecimal c = divideRound ( z . pow ( n , mcloc ) , n ) ; MathContext m = new MathContext ( err2prec ( n * z . ulp ( ) . doubleValue ( ) / 2. / z . doubleValue ( ) ) ) ; c = c . round ( m ) ; /* At larger n, zeta(n)-1 is roughly 1/2^n. The product is c/2^n.
                     * The relative error in c is c.ulp/2/c . The error in the product should be small versus eps/10.
                     * Error from 1/2^n is c*err(sigma-1).
                     * We need a relative error of zeta-1 of the order of c.ulp/50/c. This is an absolute
                     * error in zeta-1 of c.ulp/50/c/2^n, and also the absolute error in zeta, because zeta is
                     * of the order of 1.
                     */ if ( eps / 100. / c . doubleValue ( ) < 0.01 ) { m = new MathContext ( err2prec ( eps / 100. / c . doubleValue ( ) ) ) ; } else { m = new MathContext ( 2 ) ; } /* zeta(n) -1 */ BigDecimal zetm1 = zeta ( n , m ) . subtract ( BigDecimal . ONE ) ; c = multiplyRound ( c , zetm1 ) ; if ( n % 2 == 0 ) { resul = resul . add ( c ) ; } else { resul = resul . subtract ( c ) ; } /* alternating sum, so truncating as eps is reached suffices
                     */ if ( Math . abs ( c . doubleValue ( ) ) < eps ) { break ; } } } /* The relative error in the result is the absolute error in the
             * input variable times the digamma (psi) value at that point.
             */ double psi = 0.5772156649 ; double zdbl = z . doubleValue ( ) ; for ( int n = 1 ; n < 5 ; n ++ ) { psi += zdbl / n / ( n + zdbl ) ; } eps = psi * x . ulp ( ) . doubleValue ( ) / 2. ; mcloc = new MathContext ( err2prec ( eps ) ) ; return exp ( resul ) . round ( mcloc ) ; } }
The Gamma function.

@param x The argument.
@return Gamma(x).
The Gamma function .