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);
        }
    }
Euler-Mascheroni constant.

@param mc The required precision of the result.
@return 0.577...