1/* mpfr_sinh -- hyperbolic sine 2 3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao 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 sinh is done by 27 sinh(x) = 1/2 [e^(x)-e^(-x)] */ 28 29int 30mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mpfr_rnd_t rnd_mode) 31{ 32 mpfr_t x; 33 int inexact; 34 35 MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", xt, xt, rnd_mode), 36 ("y[%#R]=%R inexact=%d", y, y, inexact)); 37 38 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt))) 39 { 40 if (MPFR_IS_NAN (xt)) 41 { 42 MPFR_SET_NAN (y); 43 MPFR_RET_NAN; 44 } 45 else if (MPFR_IS_INF (xt)) 46 { 47 MPFR_SET_INF (y); 48 MPFR_SET_SAME_SIGN (y, xt); 49 MPFR_RET (0); 50 } 51 else /* xt is zero */ 52 { 53 MPFR_ASSERTD (MPFR_IS_ZERO (xt)); 54 MPFR_SET_ZERO (y); /* sinh(0) = 0 */ 55 MPFR_SET_SAME_SIGN (y, xt); 56 MPFR_RET (0); 57 } 58 } 59 60 /* sinh(x) = x + x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */ 61 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2 * MPFR_GET_EXP(xt), 2, 1, 62 rnd_mode, {}); 63 64 MPFR_TMP_INIT_ABS (x, xt); 65 66 { 67 mpfr_t t, ti; 68 mpfr_exp_t d; 69 mpfr_prec_t Nt; /* Precision of the intermediary variable */ 70 long int err; /* Precision of error */ 71 MPFR_ZIV_DECL (loop); 72 MPFR_SAVE_EXPO_DECL (expo); 73 MPFR_GROUP_DECL (group); 74 75 MPFR_SAVE_EXPO_MARK (expo); 76 77 /* compute the precision of intermediary variable */ 78 Nt = MAX (MPFR_PREC (x), MPFR_PREC (y)); 79 /* the optimal number of bits : see algorithms.ps */ 80 Nt = Nt + MPFR_INT_CEIL_LOG2 (Nt) + 4; 81 /* If x is near 0, exp(x) - 1/exp(x) = 2*x+x^3/3+O(x^5) */ 82 if (MPFR_GET_EXP (x) < 0) 83 Nt -= 2*MPFR_GET_EXP (x); 84 85 /* initialise of intermediary variables */ 86 MPFR_GROUP_INIT_2 (group, Nt, t, ti); 87 88 /* First computation of sinh */ 89 MPFR_ZIV_INIT (loop, Nt); 90 for (;;) 91 { 92 MPFR_BLOCK_DECL (flags); 93 94 /* compute sinh */ 95 MPFR_BLOCK (flags, mpfr_exp (t, x, MPFR_RNDD)); 96 if (MPFR_OVERFLOW (flags)) 97 /* exp(x) does overflow */ 98 { 99 /* sinh(x) = 2 * sinh(x/2) * cosh(x/2) */ 100 mpfr_div_2ui (ti, x, 1, MPFR_RNDD); /* exact */ 101 102 /* t <- cosh(x/2): error(t) <= 1 ulp(t) */ 103 MPFR_BLOCK (flags, mpfr_cosh (t, ti, MPFR_RNDD)); 104 if (MPFR_OVERFLOW (flags)) 105 /* when x>1 we have |sinh(x)| >= cosh(x/2), so sinh(x) 106 overflows too */ 107 { 108 inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN (xt)); 109 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW); 110 break; 111 } 112 113 /* ti <- sinh(x/2): , error(ti) <= 1 ulp(ti) 114 cannot overflow because 0 < sinh(x) < cosh(x) when x > 0 */ 115 mpfr_sinh (ti, ti, MPFR_RNDD); 116 117 /* multiplication below, error(t) <= 5 ulp(t) */ 118 MPFR_BLOCK (flags, mpfr_mul (t, t, ti, MPFR_RNDD)); 119 if (MPFR_OVERFLOW (flags)) 120 { 121 inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN (xt)); 122 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW); 123 break; 124 } 125 126 /* doubling below, exact */ 127 MPFR_BLOCK (flags, mpfr_mul_2ui (t, t, 1, MPFR_RNDN)); 128 if (MPFR_OVERFLOW (flags)) 129 { 130 inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN (xt)); 131 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW); 132 break; 133 } 134 135 /* we have lost at most 3 bits of precision */ 136 err = Nt - 3; 137 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, MPFR_PREC (y), 138 rnd_mode))) 139 { 140 inexact = mpfr_set4 (y, t, rnd_mode, MPFR_SIGN (xt)); 141 break; 142 } 143 err = Nt; /* double the precision */ 144 } 145 else 146 { 147 d = MPFR_GET_EXP (t); 148 mpfr_ui_div (ti, 1, t, MPFR_RNDU); /* 1/exp(x) */ 149 mpfr_sub (t, t, ti, MPFR_RNDN); /* exp(x) - 1/exp(x) */ 150 mpfr_div_2ui (t, t, 1, MPFR_RNDN); /* 1/2(exp(x) - 1/exp(x)) */ 151 152 /* it may be that t is zero (in fact, it can only occur when te=1, 153 and thus ti=1 too) */ 154 if (MPFR_IS_ZERO (t)) 155 err = Nt; /* double the precision */ 156 else 157 { 158 /* calculation of the error */ 159 d = d - MPFR_GET_EXP (t) + 2; 160 /* error estimate: err = Nt-(__gmpfr_ceil_log2(1+pow(2,d)));*/ 161 err = Nt - (MAX (d, 0) + 1); 162 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, MPFR_PREC (y), 163 rnd_mode))) 164 { 165 inexact = mpfr_set4 (y, t, rnd_mode, MPFR_SIGN (xt)); 166 break; 167 } 168 } 169 } 170 171 /* actualisation of the precision */ 172 Nt += err; 173 MPFR_ZIV_NEXT (loop, Nt); 174 MPFR_GROUP_REPREC_2 (group, Nt, t, ti); 175 } 176 MPFR_ZIV_FREE (loop); 177 MPFR_GROUP_CLEAR (group); 178 MPFR_SAVE_EXPO_FREE (expo); 179 } 180 181 return mpfr_check_range (y, inexact, rnd_mode); 182} 183