1/* mpfr_sinh_cosh -- hyperbolic sine and cosine 2 3Copyright 2001, 2002, 2003, 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#define INEXPOS(y) ((y) == 0 ? 0 : (((y) > 0) ? 1 : 2)) 27#define INEX(y,z) (INEXPOS(y) | (INEXPOS(z) << 2)) 28 29 /* The computations are done by 30 cosh(x) = 1/2 [e^(x)+e^(-x)] 31 sinh(x) = 1/2 [e^(x)-e^(-x)] 32 Adapted from mpfr_sinh.c */ 33 34int 35mpfr_sinh_cosh (mpfr_ptr sh, mpfr_ptr ch, mpfr_srcptr xt, mpfr_rnd_t rnd_mode) 36{ 37 mpfr_t x; 38 int inexact_sh, inexact_ch; 39 40 MPFR_ASSERTN (sh != ch); 41 42 MPFR_LOG_FUNC 43 (("x[%Pu]=%.*Rg rnd=%d", 44 mpfr_get_prec (xt), mpfr_log_prec, xt, rnd_mode), 45 ("sh[%Pu]=%.*Rg ch[%Pu]=%.*Rg", 46 mpfr_get_prec (sh), mpfr_log_prec, sh, 47 mpfr_get_prec (ch), mpfr_log_prec, ch)); 48 49 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt))) 50 { 51 if (MPFR_IS_NAN (xt)) 52 { 53 MPFR_SET_NAN (ch); 54 MPFR_SET_NAN (sh); 55 MPFR_RET_NAN; 56 } 57 else if (MPFR_IS_INF (xt)) 58 { 59 MPFR_SET_INF (sh); 60 MPFR_SET_SAME_SIGN (sh, xt); 61 MPFR_SET_INF (ch); 62 MPFR_SET_POS (ch); 63 MPFR_RET (0); 64 } 65 else /* xt is zero */ 66 { 67 MPFR_ASSERTD (MPFR_IS_ZERO (xt)); 68 MPFR_SET_ZERO (sh); /* sinh(0) = 0 */ 69 MPFR_SET_SAME_SIGN (sh, xt); 70 inexact_sh = 0; 71 inexact_ch = mpfr_set_ui (ch, 1, rnd_mode); /* cosh(0) = 1 */ 72 return INEX(inexact_sh,inexact_ch); 73 } 74 } 75 76 /* Warning: if we use MPFR_FAST_COMPUTE_IF_SMALL_INPUT here, make sure 77 that the code also works in case of overlap (see sin_cos.c) */ 78 79 MPFR_TMP_INIT_ABS (x, xt); 80 81 { 82 mpfr_t s, c, ti; 83 mpfr_exp_t d; 84 mpfr_prec_t N; /* Precision of the intermediary variables */ 85 long int err; /* Precision of error */ 86 MPFR_ZIV_DECL (loop); 87 MPFR_SAVE_EXPO_DECL (expo); 88 MPFR_GROUP_DECL (group); 89 90 MPFR_SAVE_EXPO_MARK (expo); 91 92 /* compute the precision of intermediary variable */ 93 N = MPFR_PREC (ch); 94 N = MAX (N, MPFR_PREC (sh)); 95 /* the optimal number of bits : see algorithms.ps */ 96 N = N + MPFR_INT_CEIL_LOG2 (N) + 4; 97 98 /* initialise of intermediary variables */ 99 MPFR_GROUP_INIT_3 (group, N, s, c, ti); 100 101 /* First computation of sinh_cosh */ 102 MPFR_ZIV_INIT (loop, N); 103 for (;;) 104 { 105 MPFR_BLOCK_DECL (flags); 106 107 /* compute sinh_cosh */ 108 MPFR_BLOCK (flags, mpfr_exp (s, x, MPFR_RNDD)); 109 if (MPFR_OVERFLOW (flags)) 110 /* exp(x) does overflow */ 111 { 112 /* since cosh(x) >= exp(x), cosh(x) overflows too */ 113 inexact_ch = mpfr_overflow (ch, rnd_mode, MPFR_SIGN_POS); 114 /* sinh(x) may be representable */ 115 inexact_sh = mpfr_sinh (sh, xt, rnd_mode); 116 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW); 117 break; 118 } 119 d = MPFR_GET_EXP (s); 120 mpfr_ui_div (ti, 1, s, MPFR_RNDU); /* 1/exp(x) */ 121 mpfr_add (c, s, ti, MPFR_RNDU); /* exp(x) + 1/exp(x) */ 122 mpfr_sub (s, s, ti, MPFR_RNDN); /* exp(x) - 1/exp(x) */ 123 mpfr_div_2ui (c, c, 1, MPFR_RNDN); /* 1/2(exp(x) + 1/exp(x)) */ 124 mpfr_div_2ui (s, s, 1, MPFR_RNDN); /* 1/2(exp(x) - 1/exp(x)) */ 125 126 /* it may be that s is zero (in fact, it can only occur when exp(x)=1, 127 and thus ti=1 too) */ 128 if (MPFR_IS_ZERO (s)) 129 err = N; /* double the precision */ 130 else 131 { 132 /* calculation of the error */ 133 d = d - MPFR_GET_EXP (s) + 2; 134 /* error estimate: err = N-(__gmpfr_ceil_log2(1+pow(2,d)));*/ 135 err = N - (MAX (d, 0) + 1); 136 if (MPFR_LIKELY (MPFR_CAN_ROUND (s, err, MPFR_PREC (sh), 137 rnd_mode) && \ 138 MPFR_CAN_ROUND (c, err, MPFR_PREC (ch), 139 rnd_mode))) 140 { 141 inexact_sh = mpfr_set4 (sh, s, rnd_mode, MPFR_SIGN (xt)); 142 inexact_ch = mpfr_set (ch, c, rnd_mode); 143 break; 144 } 145 } 146 /* actualisation of the precision */ 147 N += err; 148 MPFR_ZIV_NEXT (loop, N); 149 MPFR_GROUP_REPREC_3 (group, N, s, c, ti); 150 } 151 MPFR_ZIV_FREE (loop); 152 MPFR_GROUP_CLEAR (group); 153 MPFR_SAVE_EXPO_FREE (expo); 154 } 155 156 /* now, let's raise the flags if needed */ 157 inexact_sh = mpfr_check_range (sh, inexact_sh, rnd_mode); 158 inexact_ch = mpfr_check_range (ch, inexact_ch, rnd_mode); 159 160 return INEX(inexact_sh,inexact_ch); 161} 162