nd4j-common/src/main/java/org/nd4j/linalg/util/BigDecimalMath.java
static public BigDecimal gamma(MathContext mc) {
        /* look it up if possible */
        if (mc.getPrecision() < GAMMA.precision()) {
            return GAMMA.round(mc);
        } else {
            double eps = prec2err(0.577, mc.getPrecision());
            /* Euler-Stieltjes as shown in Dilcher, Aequat Math 48 (1) (1994) 55-85
            14
             */
            MathContext mcloc = new MathContext(2 + mc.getPrecision());
            BigDecimal resul = BigDecimal.ONE;
            resul = resul.add(log(2, mcloc));
            resul = resul.subtract(log(3, mcloc));
            /* how many terms: zeta-1 falls as 1/2^(2n+1), so the
             * terms drop faster than 1/2^(4n+2). Set 1/2^(4kmax+2) < eps.
             * Leading term zeta(3)/(4^1*3) is 0.017. Leading zeta(3) is 1.2. Log(2) is 0.7
             */
            int kmax = (int) ((Math.log(eps / 0.7) - 2.) / 4.);
            mcloc = new MathContext(1 + err2prec(1.2, eps / kmax));
            for (int n = 1;; n++) {
                /* zeta is close to 1. Division of zeta-1 through
                 * 4^n*(2n+1) means divion through roughly 2^(2n+1)
                 */
                BigDecimal c = zeta(2 * n + 1, mcloc).subtract(BigDecimal.ONE);
                BigInteger fourn = BigInteger.valueOf(2 * n + 1);
                fourn = fourn.shiftLeft(2 * n);
                c = divideRound(c, fourn);
                resul = resul.subtract(c);
                if (c.doubleValue() < 0.1 * eps) {
                    break;
                }
            }
            return resul.round(mc);
        }
    }
static public BigDecimal gamma(MathContext mc) {
        /* look it up if possible */
        if (mc.getPrecision() < GAMMA.precision()) {
            return GAMMA.round(mc);
        } else {
            double eps = prec2err(0.577, mc.getPrecision());
            /* Euler-Stieltjes as shown in Dilcher, Aequat Math 48 (1) (1994) 55-85
            14
             */
            MathContext mcloc = new MathContext(2 + mc.getPrecision());
            BigDecimal resul = BigDecimal.ONE;
            resul = resul.add(log(2, mcloc));
            resul = resul.subtract(log(3, mcloc));
            /* how many terms: zeta-1 falls as 1/2^(2n+1), so the
             * terms drop faster than 1/2^(4n+2). Set 1/2^(4kmax+2) < eps.
             * Leading term zeta(3)/(4^1*3) is 0.017. Leading zeta(3) is 1.2. Log(2) is 0.7
             */
            int kmax = (int) ((Math.log(eps / 0.7) - 2.) / 4.);
            mcloc = new MathContext(1 + err2prec(1.2, eps / kmax));
            for (int n = 1;; n++) {
                /* zeta is close to 1. Division of zeta-1 through
                 * 4^n*(2n+1) means divion through roughly 2^(2n+1)
                 */
                BigDecimal c = zeta(2 * n + 1, mcloc).subtract(BigDecimal.ONE);
                BigInteger fourn = BigInteger.valueOf(2 * n + 1);
                fourn = fourn.shiftLeft(2 * n);
                c = divideRound(c, fourn);
                resul = resul.subtract(c);
                if (c.doubleValue() < 0.1 * eps) {
                    break;
                }
            }
            return resul.round(mc);
        }
    }
static public BigDecimal gamma ( MathContext mc ) { /* look it up if possible */ if ( mc . getPrecision ( ) < GAMMA . precision ( ) ) { return GAMMA . round ( mc ) ; } else { double eps = prec2err ( 0.577 , mc . getPrecision ( ) ) ; /* Euler-Stieltjes as shown in Dilcher, Aequat Math 48 (1) (1994) 55-85
            14
             */ MathContext mcloc = new MathContext ( 2 + mc . getPrecision ( ) ) ; BigDecimal resul = BigDecimal . ONE ; resul = resul . add ( log ( 2 , mcloc ) ) ; resul = resul . subtract ( log ( 3 , mcloc ) ) ; /* how many terms: zeta-1 falls as 1/2^(2n+1), so the
             * terms drop faster than 1/2^(4n+2). Set 1/2^(4kmax+2) < eps.
             * Leading term zeta(3)/(4^1*3) is 0.017. Leading zeta(3) is 1.2. Log(2) is 0.7
             */ int kmax = ( int ) ( ( Math . log ( eps / 0.7 ) - 2. ) / 4. ) ; mcloc = new MathContext ( 1 + err2prec ( 1.2 , eps / kmax ) ) ; for ( int n = 1 ; ; n ++ ) { /* zeta is close to 1. Division of zeta-1 through
                 * 4^n*(2n+1) means divion through roughly 2^(2n+1)
                 */ BigDecimal c = zeta ( 2 * n + 1 , mcloc ) . subtract ( BigDecimal . ONE ) ; BigInteger fourn = BigInteger . valueOf ( 2 * n + 1 ) ; fourn = fourn . shiftLeft ( 2 * n ) ; c = divideRound ( c , fourn ) ; resul = resul . subtract ( c ) ; if ( c . doubleValue ( ) < 0.1 * eps ) { break ; } } return resul . round ( mc ) ; } }
Euler-Mascheroni constant.

@param mc The required precision of the result.
@return 0.577...
Euler - Mascheroni constant .