1 /* mpfr_atanh -- Inverse Hyperbolic Tangente
3 Copyright 2001-2015 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* The computation of atanh is done by
27 atanh= 1/2*ln(x+1)-1/2*ln(1-x) */
30 mpfr_atanh (mpfr_ptr y
, mpfr_srcptr xt
, mpfr_rnd_t rnd_mode
)
34 mpfr_prec_t Nx
, Ny
, Nt
;
37 MPFR_SAVE_EXPO_DECL (expo
);
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
));
45 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt
)))
47 /* atanh(NaN) = NaN, and atanh(+/-Inf) = NaN since tanh gives a result
49 if (MPFR_IS_NAN (xt
) || MPFR_IS_INF (xt
))
54 else /* necessarily xt is 0 */
56 MPFR_ASSERTD (MPFR_IS_ZERO (xt
));
57 MPFR_SET_ZERO (y
); /* atanh(0) = 0 */
58 MPFR_SET_SAME_SIGN (y
,xt
);
63 /* atanh (x) = NaN as soon as |x| > 1, and arctanh(+/-1) = +/-Inf */
64 if (MPFR_UNLIKELY (MPFR_GET_EXP (xt
) > 0))
66 if (MPFR_GET_EXP (xt
) == 1 && mpfr_powerof2_raw (xt
))
69 MPFR_SET_SAME_SIGN (y
, xt
);
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,
81 MPFR_SAVE_EXPO_MARK (expo
);
83 /* Compute initial precision */
85 MPFR_TMP_INIT_ABS (x
, xt
);
88 /* the optimal number of bits : see algorithms.ps */
89 Nt
= Nt
+ MPFR_INT_CEIL_LOG2 (Nt
) + 4;
91 /* initialise of intermediary variable */
95 /* First computation of cosh */
96 MPFR_ZIV_INIT (loop
, Nt
);
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))*/
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);
111 if (MPFR_LIKELY (MPFR_IS_ZERO (t
)
112 || MPFR_CAN_ROUND (t
, err
, Ny
, rnd_mode
)))
115 /* reactualisation of the precision */
116 MPFR_ZIV_NEXT (loop
, Nt
);
117 mpfr_set_prec (t
, Nt
);
118 mpfr_set_prec (te
, Nt
);
120 MPFR_ZIV_FREE (loop
);
122 inexact
= mpfr_set4 (y
, t
, rnd_mode
, MPFR_SIGN (xt
));
127 MPFR_SAVE_EXPO_FREE (expo
);
128 return mpfr_check_range (y
, inexact
, rnd_mode
);