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