nd4j-common/src/main/java/org/nd4j/linalg/util/BigDecimalMath.java
static public BigDecimal sinh(final BigDecimal x) {
        if (x.compareTo(BigDecimal.ZERO) < 0) {
            return sinh(x.negate()).negate();
        } else if (x.compareTo(BigDecimal.ZERO) == 0) {
            return BigDecimal.ZERO;
        } else {
            if (x.doubleValue() > 2.4) {
                /* Move closer to zero with sinh(2x)= 2*sinh(x)*cosh(x).
                 */
                BigDecimal two = new BigDecimal(2);
                BigDecimal xhalf = x.divide(two);


                BigDecimal resul = sinh(xhalf).multiply(cosh(xhalf)).multiply(two);
                /* The error in the result is set by the error in x itself.
                 * The first derivative of sinh(x) is cosh(x), so the absolute error
                 * in the result is cosh(x)*errx, and the relative error is coth(x)*errx = errx/tanh(x)
                 */


                double eps = Math.tanh(x.doubleValue());
                MathContext mc = new MathContext(err2prec(0.5 * x.ulp().doubleValue() / eps));


                return resul.round(mc);


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


                double xUlpDbl = x.ulp().doubleValue();
                /* The error in the result is set by the error in x itself.
                 * We need at most k terms to squeeze x^(2k+1)/(2k+1)! below this value.
                 * x^(2k+1) < x.ulp; (2k+1)*log10(x) < -x.precision; 2k*log10(x)< -x.precision;
                 * 2k*(-log10(x)) > x.precision; 2k*log10(1/x) > x.precision
                 */


                int k = (int) (x.precision() / Math.log10(1.0 / xhighpr.doubleValue())) / 2;
                MathContext mcTay = new MathContext(err2prec(x.doubleValue(), xUlpDbl / k));


                for (int i = 1;; i++) {
                    /* TBD: at which precision will 2*i or 2*i+1 overflow?
                     */
                    ifac = ifac.multiply(BigInteger.valueOf(2 * i));
                    ifac = ifac.multiply(BigInteger.valueOf(2 * i + 1));
                    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 set by the error in x itself.
                  */
                MathContext mc = new MathContext(x.precision());


                return resul.round(mc);


            }
        }
    }
static public BigDecimal sinh(final BigDecimal x) {
        if (x.compareTo(BigDecimal.ZERO) < 0) {
            return sinh(x.negate()).negate();
        } else if (x.compareTo(BigDecimal.ZERO) == 0) {
            return BigDecimal.ZERO;
        } else {
            if (x.doubleValue() > 2.4) {
                /* Move closer to zero with sinh(2x)= 2*sinh(x)*cosh(x).
                 */
                BigDecimal two = new BigDecimal(2);
                BigDecimal xhalf = x.divide(two);


                BigDecimal resul = sinh(xhalf).multiply(cosh(xhalf)).multiply(two);
                /* The error in the result is set by the error in x itself.
                 * The first derivative of sinh(x) is cosh(x), so the absolute error
                 * in the result is cosh(x)*errx, and the relative error is coth(x)*errx = errx/tanh(x)
                 */


                double eps = Math.tanh(x.doubleValue());
                MathContext mc = new MathContext(err2prec(0.5 * x.ulp().doubleValue() / eps));


                return resul.round(mc);


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


                double xUlpDbl = x.ulp().doubleValue();
                /* The error in the result is set by the error in x itself.
                 * We need at most k terms to squeeze x^(2k+1)/(2k+1)! below this value.
                 * x^(2k+1) < x.ulp; (2k+1)*log10(x) < -x.precision; 2k*log10(x)< -x.precision;
                 * 2k*(-log10(x)) > x.precision; 2k*log10(1/x) > x.precision
                 */


                int k = (int) (x.precision() / Math.log10(1.0 / xhighpr.doubleValue())) / 2;
                MathContext mcTay = new MathContext(err2prec(x.doubleValue(), xUlpDbl / k));


                for (int i = 1;; i++) {
                    /* TBD: at which precision will 2*i or 2*i+1 overflow?
                     */
                    ifac = ifac.multiply(BigInteger.valueOf(2 * i));
                    ifac = ifac.multiply(BigInteger.valueOf(2 * i + 1));
                    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 set by the error in x itself.
                  */
                MathContext mc = new MathContext(x.precision());


                return resul.round(mc);


            }
        }
    }
The hyperbolic sine.

@param x the argument.
@return the sinh(x) = (exp(x)-exp(-x))/2 .