nd4j-common/src/main/java/org/nd4j/linalg/util/BigDecimalMath.java
static public BigDecimal log(int n, final MathContext mc) {
        /* the value is undefined if x is negative.
         */
        if (n <= 0) {
            throw new ArithmeticException("Cannot take log of negative " + n);
        } else if (n == 1) {
            return BigDecimal.ZERO;
        } else if (n == 2) {
            if (mc.getPrecision() < LOG2.precision()) {
                return LOG2.round(mc);
            } else {
                /* Broadhurst \protect\vrule width0pt\protect\href{http://arxiv.org/abs/math/9803067}{arXiv:math/9803067}
                 * Error propagation: the error in log(2) is twice the error in S(2,-5,...).
                 */
                int[] a = {2, -5, -2, -7, -2, -5, 2, -3};
                BigDecimal S = broadhurstBBP(2, 1, a, new MathContext(1 + mc.getPrecision()));
                S = S.multiply(new BigDecimal(8));
                S = sqrt(divideRound(S, 3));
                return S.round(mc);
            }
        } else if (n == 3) {
            /* summation of a series roughly proportional to (7/500)^k. Estimate count
             * of terms to estimate the precision (drop the favorable additional
             * 1/k here): 0.013^k <= 10^(-precision), so k*log10(0.013) <= -precision
             * so k>= precision/1.87.
             */
            int kmax = (int) (mc.getPrecision() / 1.87);
            MathContext mcloc = new MathContext(mc.getPrecision() + 1 + (int) (Math.log10(kmax * 0.693 / 1.098)));
            BigDecimal log3 = multiplyRound(log(2, mcloc), 19);
            /* log3 is roughly 1, so absolute and relative error are the same. The
             * result will be divided by 12, so a conservative error is the one
             * already found in mc
             */
            double eps = prec2err(1.098, mc.getPrecision()) / kmax;
            Rational r = new Rational(7153, 524288);
            Rational pk = new Rational(7153, 524288);
            for (int k = 1;; k++) {
                Rational tmp = pk.divide(k);
                if (tmp.doubleValue() < eps) {
                    break;
                }
                /* how many digits of tmp do we need in the sum?
                 */
                mcloc = new MathContext(err2prec(tmp.doubleValue(), eps));
                BigDecimal c = pk.divide(k).BigDecimalValue(mcloc);
                if (k % 2 != 0) {
                    log3 = log3.add(c);
                } else {
                    log3 = log3.subtract(c);
                }
                pk = pk.multiply(r);
            }
            log3 = divideRound(log3, 12);
            return log3.round(mc);
        } else if (n == 5) {
            /* summation of a series roughly proportional to (7/160)^k. Estimate count
             * of terms to estimate the precision (drop the favorable additional
             * 1/k here): 0.046^k <= 10^(-precision), so k*log10(0.046) <= -precision
             * so k>= precision/1.33.
             */
            int kmax = (int) (mc.getPrecision() / 1.33);
            MathContext mcloc = new MathContext(mc.getPrecision() + 1 + (int) (Math.log10(kmax * 0.693 / 1.609)));
            BigDecimal log5 = multiplyRound(log(2, mcloc), 14);
            /* log5 is roughly 1.6, so absolute and relative error are the same. The
             * result will be divided by 6, so a conservative error is the one
             * already found in mc
             */
            double eps = prec2err(1.6, mc.getPrecision()) / kmax;
            Rational r = new Rational(759, 16384);
            Rational pk = new Rational(759, 16384);
            for (int k = 1;; k++) {
                Rational tmp = pk.divide(k);
                if (tmp.doubleValue() < eps) {
                    break;
                }
                /* how many digits of tmp do we need in the sum?
                 */
                mcloc = new MathContext(err2prec(tmp.doubleValue(), eps));
                BigDecimal c = pk.divide(k).BigDecimalValue(mcloc);
                log5 = log5.subtract(c);
                pk = pk.multiply(r);
            }
            log5 = divideRound(log5, 6);
            return log5.round(mc);
        } else if (n == 7) {
            /* summation of a series roughly proportional to (1/8)^k. Estimate count
             * of terms to estimate the precision (drop the favorable additional
             * 1/k here): 0.125^k <= 10^(-precision), so k*log10(0.125) <= -precision
             * so k>= precision/0.903.
             */
            int kmax = (int) (mc.getPrecision() / 0.903);
            MathContext mcloc = new MathContext(mc.getPrecision() + 1 + (int) (Math.log10(kmax * 3 * 0.693 / 1.098)));
            BigDecimal log7 = multiplyRound(log(2, mcloc), 3);
            /* log7 is roughly 1.9, so absolute and relative error are the same.
             */
            double eps = prec2err(1.9, mc.getPrecision()) / kmax;
            Rational r = new Rational(1, 8);
            Rational pk = new Rational(1, 8);
            for (int k = 1;; k++) {
                Rational tmp = pk.divide(k);
                if (tmp.doubleValue() < eps) {
                    break;
                }
                /* how many digits of tmp do we need in the sum?
                 */
                mcloc = new MathContext(err2prec(tmp.doubleValue(), eps));
                BigDecimal c = pk.divide(k).BigDecimalValue(mcloc);
                log7 = log7.subtract(c);
                pk = pk.multiply(r);
            }
            return log7.round(mc);
        } else {
            /* At this point one could either forward to the log(BigDecimal) signature (implemented)
             * or decompose n into Ifactors and use an implemenation of all the prime bases.
             * Estimate of the result; convert the mc argument to an absolute error eps
             * log(n+errn) = log(n)+errn/n = log(n)+eps
             */
            double res = Math.log((double) n);
            double eps = prec2err(res, mc.getPrecision());
            /* errn = eps*n, convert absolute error in result to requirement on absolute error in input
             */
            eps *= n;

            /* Convert this absolute requirement of error in n to a relative error in n
             */
            final MathContext mcloc = new MathContext(1 + err2prec((double) n, eps));
            /* Padd n with a number of zeros to trigger the required accuracy in
             * the standard signature method
             */
            BigDecimal nb = scalePrec(new BigDecimal(n), mcloc);
            return log(nb);
        }
    }
static public BigDecimal log(int n, final MathContext mc) {
        /* the value is undefined if x is negative.
         */
        if (n <= 0) {
            throw new ArithmeticException("Cannot take log of negative " + n);
        } else if (n == 1) {
            return BigDecimal.ZERO;
        } else if (n == 2) {
            if (mc.getPrecision() < LOG2.precision()) {
                return LOG2.round(mc);
            } else {
                /* Broadhurst \protect\vrule width0pt\protect\href{http://arxiv.org/abs/math/9803067}{arXiv:math/9803067}
                 * Error propagation: the error in log(2) is twice the error in S(2,-5,...).
                 */
                int[] a = {2, -5, -2, -7, -2, -5, 2, -3};
                BigDecimal S = broadhurstBBP(2, 1, a, new MathContext(1 + mc.getPrecision()));
                S = S.multiply(new BigDecimal(8));
                S = sqrt(divideRound(S, 3));
                return S.round(mc);
            }
        } else if (n == 3) {
            /* summation of a series roughly proportional to (7/500)^k. Estimate count
             * of terms to estimate the precision (drop the favorable additional
             * 1/k here): 0.013^k <= 10^(-precision), so k*log10(0.013) <= -precision
             * so k>= precision/1.87.
             */
            int kmax = (int) (mc.getPrecision() / 1.87);
            MathContext mcloc = new MathContext(mc.getPrecision() + 1 + (int) (Math.log10(kmax * 0.693 / 1.098)));
            BigDecimal log3 = multiplyRound(log(2, mcloc), 19);
            /* log3 is roughly 1, so absolute and relative error are the same. The
             * result will be divided by 12, so a conservative error is the one
             * already found in mc
             */
            double eps = prec2err(1.098, mc.getPrecision()) / kmax;
            Rational r = new Rational(7153, 524288);
            Rational pk = new Rational(7153, 524288);
            for (int k = 1;; k++) {
                Rational tmp = pk.divide(k);
                if (tmp.doubleValue() < eps) {
                    break;
                }
                /* how many digits of tmp do we need in the sum?
                 */
                mcloc = new MathContext(err2prec(tmp.doubleValue(), eps));
                BigDecimal c = pk.divide(k).BigDecimalValue(mcloc);
                if (k % 2 != 0) {
                    log3 = log3.add(c);
                } else {
                    log3 = log3.subtract(c);
                }
                pk = pk.multiply(r);
            }
            log3 = divideRound(log3, 12);
            return log3.round(mc);
        } else if (n == 5) {
            /* summation of a series roughly proportional to (7/160)^k. Estimate count
             * of terms to estimate the precision (drop the favorable additional
             * 1/k here): 0.046^k <= 10^(-precision), so k*log10(0.046) <= -precision
             * so k>= precision/1.33.
             */
            int kmax = (int) (mc.getPrecision() / 1.33);
            MathContext mcloc = new MathContext(mc.getPrecision() + 1 + (int) (Math.log10(kmax * 0.693 / 1.609)));
            BigDecimal log5 = multiplyRound(log(2, mcloc), 14);
            /* log5 is roughly 1.6, so absolute and relative error are the same. The
             * result will be divided by 6, so a conservative error is the one
             * already found in mc
             */
            double eps = prec2err(1.6, mc.getPrecision()) / kmax;
            Rational r = new Rational(759, 16384);
            Rational pk = new Rational(759, 16384);
            for (int k = 1;; k++) {
                Rational tmp = pk.divide(k);
                if (tmp.doubleValue() < eps) {
                    break;
                }
                /* how many digits of tmp do we need in the sum?
                 */
                mcloc = new MathContext(err2prec(tmp.doubleValue(), eps));
                BigDecimal c = pk.divide(k).BigDecimalValue(mcloc);
                log5 = log5.subtract(c);
                pk = pk.multiply(r);
            }
            log5 = divideRound(log5, 6);
            return log5.round(mc);
        } else if (n == 7) {
            /* summation of a series roughly proportional to (1/8)^k. Estimate count
             * of terms to estimate the precision (drop the favorable additional
             * 1/k here): 0.125^k <= 10^(-precision), so k*log10(0.125) <= -precision
             * so k>= precision/0.903.
             */
            int kmax = (int) (mc.getPrecision() / 0.903);
            MathContext mcloc = new MathContext(mc.getPrecision() + 1 + (int) (Math.log10(kmax * 3 * 0.693 / 1.098)));
            BigDecimal log7 = multiplyRound(log(2, mcloc), 3);
            /* log7 is roughly 1.9, so absolute and relative error are the same.
             */
            double eps = prec2err(1.9, mc.getPrecision()) / kmax;
            Rational r = new Rational(1, 8);
            Rational pk = new Rational(1, 8);
            for (int k = 1;; k++) {
                Rational tmp = pk.divide(k);
                if (tmp.doubleValue() < eps) {
                    break;
                }
                /* how many digits of tmp do we need in the sum?
                 */
                mcloc = new MathContext(err2prec(tmp.doubleValue(), eps));
                BigDecimal c = pk.divide(k).BigDecimalValue(mcloc);
                log7 = log7.subtract(c);
                pk = pk.multiply(r);
            }
            return log7.round(mc);
        } else {
            /* At this point one could either forward to the log(BigDecimal) signature (implemented)
             * or decompose n into Ifactors and use an implemenation of all the prime bases.
             * Estimate of the result; convert the mc argument to an absolute error eps
             * log(n+errn) = log(n)+errn/n = log(n)+eps
             */
            double res = Math.log((double) n);
            double eps = prec2err(res, mc.getPrecision());
            /* errn = eps*n, convert absolute error in result to requirement on absolute error in input
             */
            eps *= n;

            /* Convert this absolute requirement of error in n to a relative error in n
             */
            final MathContext mcloc = new MathContext(1 + err2prec((double) n, eps));
            /* Padd n with a number of zeros to trigger the required accuracy in
             * the standard signature method
             */
            BigDecimal nb = scalePrec(new BigDecimal(n), mcloc);
            return log(nb);
        }
    }
The natural logarithm.

@param n  The main argument, a strictly positive integer.
@param mc The requirements on the precision.
@return ln(n).