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


            }
        }
    }
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 .
The hyperbolic cosine .