1 /* mpfr_cosh -- hyperbolic cosine
3 Copyright 2001-2002, 2004-2016 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
6 This file is part of the GNU MPFR Library.
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* The computation of cosh is done by *
27 * cosh= 1/2[e^(x)+e^(-x)] */
30 mpfr_cosh (mpfr_ptr y
, mpfr_srcptr xt
, mpfr_rnd_t rnd_mode
)
34 MPFR_SAVE_EXPO_DECL (expo
);
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
,
41 if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(xt
)))
48 else if (MPFR_IS_INF(xt
))
56 MPFR_ASSERTD(MPFR_IS_ZERO(xt
));
57 return mpfr_set_ui (y
, 1, rnd_mode
); /* cosh(0) = 1 */
61 MPFR_SAVE_EXPO_MARK (expo
);
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
);
69 MPFR_TMP_INIT_ABS(x
, xt
);
72 /* Declaration of the intermediary variable */
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 */
79 MPFR_GROUP_DECL (group
);
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
);
85 /* initialise of intermediary variables */
86 MPFR_GROUP_INIT_2 (group
, Nt
, t
, te
);
88 /* First computation of cosh */
89 MPFR_ZIV_INIT (loop
, Nt
);
92 MPFR_BLOCK_DECL (flags
);
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 */
100 inexact
= mpfr_overflow (y
, rnd_mode
, MPFR_SIGN_POS
);
101 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, MPFR_FLAGS_OVERFLOW
);
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))*/
108 /* Estimation of the error */
110 /* Check if we can round */
111 if (MPFR_LIKELY (MPFR_CAN_ROUND (t
, err
, Ny
, rnd_mode
)))
113 inexact
= mpfr_set (y
, t
, rnd_mode
);
117 /* Actualisation of the precision */
118 MPFR_ZIV_NEXT (loop
, Nt
);
119 MPFR_GROUP_REPREC_2 (group
, Nt
, t
, te
);
121 MPFR_ZIV_FREE (loop
);
122 MPFR_GROUP_CLEAR (group
);
126 MPFR_SAVE_EXPO_FREE (expo
);
127 return mpfr_check_range (y
, inexact
, rnd_mode
);