1 /* mpfr_sinh_cosh -- hyperbolic sine and cosine 2 3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramel projects, INRIA. 5 6 This file is part of the GNU MPFR Library. 7 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. 12 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. 17 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. */ 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 34 int 35 mpfr_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