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 } 216
Indexes created Tue Oct 20 00:27:02 UTC 2026