1 /* mpfr_atanh -- Inverse Hyperbolic Tangente
2
3 Copyright 2001-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 https://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 /* Put in y an approximation of atanh(x) for x small.
27 We assume x <= 1/2, in which case:
28 x <= y ~ atanh(x) = x + x^3/3 + x^5/5 + x^7/7 + ... <= 2*x.
29 Return k such that the error is bounded by 2^k*ulp(y).
30 */
31 static int
32 mpfr_atanh_small (mpfr_ptr y, mpfr_srcptr x)
33 {
34 mpfr_prec_t p = MPFR_PREC(y), err;
35 mpfr_t x2, t, u;
36 unsigned long i;
37 int k;
38
39 MPFR_ASSERTD(MPFR_GET_EXP (x) <= -1);
40
41 /* in the following, theta represents a value with |theta| <= 2^(1-p)
42 (might be a different value each time) */
43
44 mpfr_init2 (t, p);
45 mpfr_init2 (u, p);
46 mpfr_init2 (x2, p);
47 mpfr_set (t, x, MPFR_RNDF); /* t = x * (1 + theta) */
48 mpfr_set (y, t, MPFR_RNDF); /* exact */
49 mpfr_sqr (x2, x, MPFR_RNDF); /* x2 = x^2 * (1 + theta) */
50 for (i = 3; ; i += 2)
51 {
52 mpfr_mul (t, t, x2, MPFR_RNDF); /* t = x^i * (1 + theta)^i */
53 mpfr_div_ui (u, t, i, MPFR_RNDF); /* u = x^i/i * (1 + theta)^(i+1) */
54 if (MPFR_GET_EXP (u) <= MPFR_GET_EXP (y) - p) /* |u| < ulp(y) */
55 break;
56 mpfr_add (y, y, u, MPFR_RNDF); /* error <= ulp(y) */
57 }
58 /* We assume |(1 + theta)^(i+1)| <= 2.
59 The neglected part is at most |u| + |u|/4 + |u|/16 + ... <= 4/3*|u|,
60 which has to be multiplied by |(1 + theta)^(i+1)| <= 2, thus at most
61 3 ulp(y).
62 The rounding error on y is bounded by:
63 * for the (i-3)/2 add/sub, each error is bounded by ulp(y_i),
64 where y_i is the current value of y, which is bounded by ulp(y)
65 for y the final value (since it increases in absolute value),
66 this yields (i-3)/2*ulp(y)
67 * from Lemma 3.1 from [Higham02] (see algorithms.tex),
68 the relative error on u at step i is bounded by:
69 (i+1)*epsilon/(1-(i+1)*epsilon) where epsilon = 2^(1-p).
70 If (i+1)*epsilon <= 1/2, then the relative error on u at
71 step i is bounded by 2*(i+1)*epsilon, and since |u| <= 1/2^(i+1)
72 at step i, this gives an absolute error bound of;
73 2*epsilon*x*(4/2^4 + 6/2^6 + 8/2^8 + ...) = 2*2^(1-p)*x*(7/18) =
74 14/9*2^(-p)*x <= 2*ulp(x).
75
76 If (i+1)*epsilon <= 1/2, then the relative error on u at step i
77 is bounded by (i+1)*epsilon/(1-(i+1)*epsilon) <= 1, thus it follows
78 |(1 + theta)^(i+1)| <= 2.
79
80 Finally the total error is bounded by 3*ulp(y) + (i-3)/2*ulp(y) +2*ulp(x).
81 Since x <= 2*y, we have ulp(x) <= 2*ulp(y), thus the error is bounded by:
82 (i+7)/2*ulp(y).
83 */
84 err = (i + 8) / 2; /* ceil((i+7)/2) */
85 k = __gmpfr_int_ceil_log2 (err);
86 MPFR_ASSERTN(k + 2 < p);
87 /* if k + 2 < p, since k = ceil(log2(err)), we have err <= 2^k <= 2^(p-3),
88 thus i+7 <= 2*err <= 2^(p-2), thus (i+7)*epsilon <= 1/2, which implies
89 our assumption (i+1)*epsilon <= 1/2. */
90 mpfr_clear (t);
91 mpfr_clear (u);
92 mpfr_clear (x2);
93 return k;
94 }
95
96 /* The computation of atanh is done by:
97 atanh = ln((1+x)/(1-x)) / 2
98 except when x is very small, in which case atanh = x + tiny error,
99 and when x is small, where we use directly the Taylor expansion.
100 */
101
102 int
103 mpfr_atanh (mpfr_ptr y, mpfr_srcptr xt, mpfr_rnd_t rnd_mode)
104 {
105 int inexact;
106 mpfr_t x, t, te;
107 mpfr_prec_t Nx, Ny, Nt;
108 mpfr_exp_t err;
109 MPFR_ZIV_DECL (loop);
110 MPFR_SAVE_EXPO_DECL (expo);
111
112 MPFR_LOG_FUNC
113 (("x[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (xt), mpfr_log_prec, xt, rnd_mode),
114 ("y[%Pd]=%.*Rg inexact=%d",
115 mpfr_get_prec (y), mpfr_log_prec, y, inexact));
116
117 /* Special cases */
118 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
119 {
120 /* atanh(NaN) = NaN, and atanh(+/-Inf) = NaN since tanh gives a result
121 between -1 and 1 */
122 if (MPFR_IS_NAN (xt) || MPFR_IS_INF (xt))
123 {
124 MPFR_SET_NAN (y);
125 MPFR_RET_NAN;
126 }
127 else /* necessarily xt is 0 */
128 {
129 MPFR_ASSERTD (MPFR_IS_ZERO (xt));
130 MPFR_SET_ZERO (y); /* atanh(0) = 0 */
131 MPFR_SET_SAME_SIGN (y,xt);
132 MPFR_RET (0);
133 }
134 }
135
136 /* atanh (x) = NaN as soon as |x| > 1, and arctanh(+/-1) = +/-Inf */
137 if (MPFR_UNLIKELY (MPFR_GET_EXP (xt) > 0))
138 {
139 if (MPFR_GET_EXP (xt) == 1 && mpfr_powerof2_raw (xt))
140 {
141 MPFR_SET_INF (y);
142 MPFR_SET_SAME_SIGN (y, xt);
143 MPFR_SET_DIVBY0 ();
144 MPFR_RET (0);
145 }
146 MPFR_SET_NAN (y);
147 MPFR_RET_NAN;
148 }
149
150 /* atanh(x) = x + x^3/3 + ... so the error is < 2^(3*EXP(x)-1) */
151 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2 * MPFR_GET_EXP (xt), 1, 1,
152 rnd_mode, {});
153
154 MPFR_SAVE_EXPO_MARK (expo);
155
156 /* Compute initial precision */
157 Nx = MPFR_PREC (xt);
158 MPFR_TMP_INIT_ABS (x, xt);
159 Ny = MPFR_PREC (y);
160 Nt = MAX (Nx, Ny);
161 Nt = Nt + MPFR_INT_CEIL_LOG2 (Nt) + 4;
162
163 /* initialize of intermediary variable */
164 mpfr_init2 (t, Nt);
165 mpfr_init2 (te, Nt);
166
167 MPFR_ZIV_INIT (loop, Nt);
168 for (;;)
169 {
170 int k;
171
172 /* small case: assuming the AGM algorithm used by mpfr_log uses
173 log2(p) steps for a precision of p bits, we try the special
174 variant whenever EXP(x) <= -p/log2(p). */
175 k = 1 + __gmpfr_int_ceil_log2 (Ny); /* the +1 avoids a division by 0
176 when Ny=1 */
177 if (MPFR_GET_EXP (x) <= - 1 - (mpfr_exp_t) (Ny / k))
178 /* this implies EXP(x) <= -1 thus x < 1/2 */
179 {
180 err = Nt - mpfr_atanh_small (t, x);
181 goto round;
182 }
183
184 /* compute atanh */
185 mpfr_ui_sub (te, 1, x, MPFR_RNDU); /* (1-x) with x = |xt| */
186 mpfr_add_ui (t, x, 1, MPFR_RNDD); /* (1+x) */
187 mpfr_div (t, t, te, MPFR_RNDN); /* (1+x)/(1-x) */
188 mpfr_log (t, t, MPFR_RNDN); /* ln((1+x)/(1-x)) */
189 mpfr_div_2ui (t, t, 1, MPFR_RNDN); /* ln((1+x)/(1-x)) / 2 */
190
191 /* error estimate: see algorithms.tex */
192 /* FIXME: this does not correspond to the value in algorithms.tex!!! */
193 /* err = Nt - __gmpfr_ceil_log2(1+5*pow(2,1-MPFR_EXP(t))); */
194 err = Nt - (MAX (4 - MPFR_GET_EXP (t), 0) + 1);
195
196 round:
197 if (MPFR_LIKELY (MPFR_IS_ZERO (t)
198 || MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
199 break;
200
201 /* reactualisation of the precision */
202 MPFR_ZIV_NEXT (loop, Nt);
203 mpfr_set_prec (t, Nt);
204 mpfr_set_prec (te, Nt);
205 }
206 MPFR_ZIV_FREE (loop);
207
208 inexact = mpfr_set4 (y, t, rnd_mode, MPFR_SIGN (xt));
209
210 mpfr_clear (t);
211 mpfr_clear (te);
212
213 MPFR_SAVE_EXPO_FREE (expo);
214 return mpfr_check_range (y, inexact, rnd_mode);
215 }