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 .