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


        }
    }
The Gamma function.

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