nd4j-common/src/main/java/org/nd4j/linalg/util/BigDecimalMath.java
static public BigDecimal cosh(final BigDecimal x) {
        if (x.compareTo(BigDecimal.ZERO) < 0) {
            return cos(x.negate());
        } else if (x.compareTo(BigDecimal.ZERO) == 0) {
            return BigDecimal.ONE;
        } else {
            if (x.doubleValue() > 1.5) {
                /* cosh^2(x) = 1+ sinh^2(x).
                 */
                return hypot(1, sinh(x));


            } else {
                BigDecimal xhighpr = scalePrec(x, 2);
                /* Simple Taylor expansion, sum_{0=1..infinity} x^(2i)/(2i)! */
                BigDecimal resul = BigDecimal.ONE;
                /* x^i */
                BigDecimal xpowi = BigDecimal.ONE;
                /* 2i factorial */
                BigInteger ifac = BigInteger.ONE;
                /* The absolute error in the result is the error in x^2/2 which is x times the error in x.
                 */


                double xUlpDbl = 0.5 * x.ulp().doubleValue() * x.doubleValue();
                /* The error in the result is set by the error in x^2/2 itself, xUlpDbl.
                 * We need at most k terms to push x^(2k)/(2k)! below this value.
                 * x^(2k) < xUlpDbl; (2k)*log(x) < log(xUlpDbl);
                 */


                int k = (int) (Math.log(xUlpDbl) / Math.log(x.doubleValue())) / 2;
                /* The individual terms are all smaller than 1, so an estimate of 1.0 for
                 * the absolute value will give a safe relative error estimate for the indivdual terms
                 */
                MathContext mcTay = new MathContext(err2prec(1., xUlpDbl / k));


                for (int i = 1;; i++) {
                    /* TBD: at which precision will 2*i-1 or 2*i overflow?
                     */
                    ifac = ifac.multiply(BigInteger.valueOf(2 * i - 1));
                    ifac = ifac.multiply(BigInteger.valueOf(2 * i));
                    xpowi = xpowi.multiply(xhighpr).multiply(xhighpr);
                    BigDecimal corr = xpowi.divide(new BigDecimal(ifac), mcTay);
                    resul = resul.add(corr);


                    if (corr.abs().doubleValue() < 0.5 * xUlpDbl) {
                        break;
                    }


                } /* The error in the result is governed by the error in x itself.
                  */
                MathContext mc = new MathContext(err2prec(resul.doubleValue(), xUlpDbl));


                return resul.round(mc);


            }
        }
    }
static public BigDecimal cosh(final BigDecimal x) {
        if (x.compareTo(BigDecimal.ZERO) < 0) {
            return cos(x.negate());
        } else if (x.compareTo(BigDecimal.ZERO) == 0) {
            return BigDecimal.ONE;
        } else {
            if (x.doubleValue() > 1.5) {
                /* cosh^2(x) = 1+ sinh^2(x).
                 */
                return hypot(1, sinh(x));


            } else {
                BigDecimal xhighpr = scalePrec(x, 2);
                /* Simple Taylor expansion, sum_{0=1..infinity} x^(2i)/(2i)! */
                BigDecimal resul = BigDecimal.ONE;
                /* x^i */
                BigDecimal xpowi = BigDecimal.ONE;
                /* 2i factorial */
                BigInteger ifac = BigInteger.ONE;
                /* The absolute error in the result is the error in x^2/2 which is x times the error in x.
                 */


                double xUlpDbl = 0.5 * x.ulp().doubleValue() * x.doubleValue();
                /* The error in the result is set by the error in x^2/2 itself, xUlpDbl.
                 * We need at most k terms to push x^(2k)/(2k)! below this value.
                 * x^(2k) < xUlpDbl; (2k)*log(x) < log(xUlpDbl);
                 */


                int k = (int) (Math.log(xUlpDbl) / Math.log(x.doubleValue())) / 2;
                /* The individual terms are all smaller than 1, so an estimate of 1.0 for
                 * the absolute value will give a safe relative error estimate for the indivdual terms
                 */
                MathContext mcTay = new MathContext(err2prec(1., xUlpDbl / k));


                for (int i = 1;; i++) {
                    /* TBD: at which precision will 2*i-1 or 2*i overflow?
                     */
                    ifac = ifac.multiply(BigInteger.valueOf(2 * i - 1));
                    ifac = ifac.multiply(BigInteger.valueOf(2 * i));
                    xpowi = xpowi.multiply(xhighpr).multiply(xhighpr);
                    BigDecimal corr = xpowi.divide(new BigDecimal(ifac), mcTay);
                    resul = resul.add(corr);


                    if (corr.abs().doubleValue() < 0.5 * xUlpDbl) {
                        break;
                    }


                } /* The error in the result is governed by the error in x itself.
                  */
                MathContext mc = new MathContext(err2prec(resul.doubleValue(), xUlpDbl));


                return resul.round(mc);


            }
        }
    }
The hyperbolic cosine.

@param x The argument.
@return The cosh(x) = (exp(x)+exp(-x))/2 .