1/* mpfr_atanh -- Inverse Hyperbolic Tangente 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 /* The computation of atanh is done by 27 atanh= 1/2*ln(x+1)-1/2*ln(1-x) */ 28 29int 30mpfr_atanh (mpfr_ptr y, mpfr_srcptr xt , mpfr_rnd_t rnd_mode) 31{ 32 int inexact; 33 mpfr_t x, t, te; 34 mpfr_prec_t Nx, Ny, Nt; 35 mpfr_exp_t err; 36 MPFR_ZIV_DECL (loop); 37 MPFR_SAVE_EXPO_DECL (expo); 38 39 MPFR_LOG_FUNC 40 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (xt), mpfr_log_prec, xt, rnd_mode), 41 ("y[%Pu]=%.*Rg inexact=%d", 42 mpfr_get_prec (y), mpfr_log_prec, y, inexact)); 43 44 /* Special cases */ 45 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt))) 46 { 47 /* atanh(NaN) = NaN, and atanh(+/-Inf) = NaN since tanh gives a result 48 between -1 and 1 */ 49 if (MPFR_IS_NAN (xt) || MPFR_IS_INF (xt)) 50 { 51 MPFR_SET_NAN (y); 52 MPFR_RET_NAN; 53 } 54 else /* necessarily xt is 0 */ 55 { 56 MPFR_ASSERTD (MPFR_IS_ZERO (xt)); 57 MPFR_SET_ZERO (y); /* atanh(0) = 0 */ 58 MPFR_SET_SAME_SIGN (y,xt); 59 MPFR_RET (0); 60 } 61 } 62 63 /* atanh (x) = NaN as soon as |x| > 1, and arctanh(+/-1) = +/-Inf */ 64 if (MPFR_UNLIKELY (MPFR_GET_EXP (xt) > 0)) 65 { 66 if (MPFR_GET_EXP (xt) == 1 && mpfr_powerof2_raw (xt)) 67 { 68 MPFR_SET_INF (y); 69 MPFR_SET_SAME_SIGN (y, xt); 70 mpfr_set_divby0 (); 71 MPFR_RET (0); 72 } 73 MPFR_SET_NAN (y); 74 MPFR_RET_NAN; 75 } 76 77 /* atanh(x) = x + x^3/3 + ... so the error is < 2^(3*EXP(x)-1) */ 78 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2 * MPFR_GET_EXP (xt), 1, 1, 79 rnd_mode, {}); 80 81 MPFR_SAVE_EXPO_MARK (expo); 82 83 /* Compute initial precision */ 84 Nx = MPFR_PREC (xt); 85 MPFR_TMP_INIT_ABS (x, xt); 86 Ny = MPFR_PREC (y); 87 Nt = MAX (Nx, Ny); 88 /* the optimal number of bits : see algorithms.ps */ 89 Nt = Nt + MPFR_INT_CEIL_LOG2 (Nt) + 4; 90 91 /* initialise of intermediary variable */ 92 mpfr_init2 (t, Nt); 93 mpfr_init2 (te, Nt); 94 95 /* First computation of cosh */ 96 MPFR_ZIV_INIT (loop, Nt); 97 for (;;) 98 { 99 /* compute atanh */ 100 mpfr_ui_sub (te, 1, x, MPFR_RNDU); /* (1-xt)*/ 101 mpfr_add_ui (t, x, 1, MPFR_RNDD); /* (xt+1)*/ 102 mpfr_div (t, t, te, MPFR_RNDN); /* (1+xt)/(1-xt)*/ 103 mpfr_log (t, t, MPFR_RNDN); /* ln((1+xt)/(1-xt))*/ 104 mpfr_div_2ui (t, t, 1, MPFR_RNDN); /* (1/2)*ln((1+xt)/(1-xt))*/ 105 106 /* error estimate: see algorithms.tex */ 107 /* FIXME: this does not correspond to the value in algorithms.tex!!! */ 108 /* err=Nt-__gmpfr_ceil_log2(1+5*pow(2,1-MPFR_EXP(t)));*/ 109 err = Nt - (MAX (4 - MPFR_GET_EXP (t), 0) + 1); 110 111 if (MPFR_LIKELY (MPFR_IS_ZERO (t) 112 || MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) 113 break; 114 115 /* reactualisation of the precision */ 116 MPFR_ZIV_NEXT (loop, Nt); 117 mpfr_set_prec (t, Nt); 118 mpfr_set_prec (te, Nt); 119 } 120 MPFR_ZIV_FREE (loop); 121 122 inexact = mpfr_set4 (y, t, rnd_mode, MPFR_SIGN (xt)); 123 124 mpfr_clear(t); 125 mpfr_clear(te); 126 127 MPFR_SAVE_EXPO_FREE (expo); 128 return mpfr_check_range (y, inexact, rnd_mode); 129} 130 131