1/* mpfr_cosh -- hyperbolic cosine 2 3Copyright 2001, 2002, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#define MPFR_NEED_LONGLONG_H 24#include "mpfr-impl.h" 25 26/* The computation of cosh is done by * 27 * cosh= 1/2[e^(x)+e^(-x)] */ 28 29int 30mpfr_cosh (mpfr_ptr y, mpfr_srcptr xt , mpfr_rnd_t rnd_mode) 31{ 32 mpfr_t x; 33 int inexact; 34 MPFR_SAVE_EXPO_DECL (expo); 35 36 MPFR_LOG_FUNC ( 37 ("x[%Pu]=%*.Rg rnd=%d", mpfr_get_prec (xt), mpfr_log_prec, xt, rnd_mode), 38 ("y[%Pu]=%*.Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, 39 inexact)); 40 41 if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(xt))) 42 { 43 if (MPFR_IS_NAN(xt)) 44 { 45 MPFR_SET_NAN(y); 46 MPFR_RET_NAN; 47 } 48 else if (MPFR_IS_INF(xt)) 49 { 50 MPFR_SET_INF(y); 51 MPFR_SET_POS(y); 52 MPFR_RET(0); 53 } 54 else 55 { 56 MPFR_ASSERTD(MPFR_IS_ZERO(xt)); 57 return mpfr_set_ui (y, 1, rnd_mode); /* cosh(0) = 1 */ 58 } 59 } 60 61 MPFR_SAVE_EXPO_MARK (expo); 62 63 /* cosh(x) = 1+x^2/2 + ... <= 1+x^2 for x <= 2.9828..., 64 thus the error < 2^(2*EXP(x)). If x >= 1, then EXP(x) >= 1, 65 thus the following will always fail. */ 66 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, __gmpfr_one, -2 * MPFR_GET_EXP (xt), 0, 67 1, rnd_mode, inexact = _inexact; goto end); 68 69 MPFR_TMP_INIT_ABS(x, xt); 70 /* General case */ 71 { 72 /* Declaration of the intermediary variable */ 73 mpfr_t t, te; 74 /* Declaration of the size variable */ 75 mpfr_prec_t Ny = MPFR_PREC(y); /* Precision of output variable */ 76 mpfr_prec_t Nt; /* Precision of the intermediary variable */ 77 long int err; /* Precision of error */ 78 MPFR_ZIV_DECL (loop); 79 MPFR_GROUP_DECL (group); 80 81 /* compute the precision of intermediary variable */ 82 /* The optimal number of bits : see algorithms.tex */ 83 Nt = Ny + 3 + MPFR_INT_CEIL_LOG2 (Ny); 84 85 /* initialise of intermediary variables */ 86 MPFR_GROUP_INIT_2 (group, Nt, t, te); 87 88 /* First computation of cosh */ 89 MPFR_ZIV_INIT (loop, Nt); 90 for (;;) 91 { 92 MPFR_BLOCK_DECL (flags); 93 94 /* Compute cosh */ 95 MPFR_BLOCK (flags, mpfr_exp (te, x, MPFR_RNDD)); /* exp(x) */ 96 /* exp can overflow (but not underflow since x>0) */ 97 if (MPFR_OVERFLOW (flags)) 98 /* cosh(x) > exp(x), cosh(x) underflows too */ 99 { 100 inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN_POS); 101 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW); 102 break; 103 } 104 mpfr_ui_div (t, 1, te, MPFR_RNDU); /* 1/exp(x) */ 105 mpfr_add (t, te, t, MPFR_RNDU); /* exp(x) + 1/exp(x)*/ 106 mpfr_div_2ui (t, t, 1, MPFR_RNDN); /* 1/2(exp(x) + 1/exp(x))*/ 107 108 /* Estimation of the error */ 109 err = Nt - 3; 110 /* Check if we can round */ 111 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) 112 { 113 inexact = mpfr_set (y, t, rnd_mode); 114 break; 115 } 116 117 /* Actualisation of the precision */ 118 MPFR_ZIV_NEXT (loop, Nt); 119 MPFR_GROUP_REPREC_2 (group, Nt, t, te); 120 } 121 MPFR_ZIV_FREE (loop); 122 MPFR_GROUP_CLEAR (group); 123 } 124 125 end: 126 MPFR_SAVE_EXPO_FREE (expo); 127 return mpfr_check_range (y, inexact, rnd_mode); 128} 129