1 | /* Copyright (C) 2011-2023 Free Software Foundation, Inc. |
2 | This file is part of the GNU C Library. |
3 | |
4 | The GNU C Library is free software; you can redistribute it and/or |
5 | modify it under the terms of the GNU Lesser General Public |
6 | License as published by the Free Software Foundation; either |
7 | version 2.1 of the License, or (at your option) any later version. |
8 | |
9 | The GNU C Library is distributed in the hope that it will be useful, |
10 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
11 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
12 | Lesser General Public License for more details. |
13 | |
14 | You should have received a copy of the GNU Lesser General Public |
15 | License along with the GNU C Library; if not, see |
16 | <https://www.gnu.org/licenses/>. */ |
17 | |
18 | |
19 | /* __ieee754_atanh(x) |
20 | Method : |
21 | 1.Reduced x to positive by atanh(-x) = -atanh(x) |
22 | 2.For x>=0.5 |
23 | 1 2x x |
24 | atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) |
25 | 2 1 - x 1 - x |
26 | |
27 | For x<0.5 |
28 | atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) |
29 | |
30 | Special cases: |
31 | atanh(x) is NaN if |x| > 1 with signal; |
32 | atanh(NaN) is that NaN with no signal; |
33 | atanh(+-1) is +-INF with signal. |
34 | |
35 | */ |
36 | |
37 | #include <float.h> |
38 | #include <inttypes.h> |
39 | #include <math.h> |
40 | #include <math-barriers.h> |
41 | #include <math_private.h> |
42 | #include <math-underflow.h> |
43 | #include <libm-alias-finite.h> |
44 | |
45 | static const float huge = 1e30; |
46 | |
47 | float |
48 | __ieee754_atanhf (float x) |
49 | { |
50 | float xa = fabsf (x); |
51 | float t; |
52 | if (isless (xa, 0.5f)) |
53 | { |
54 | if (__glibc_unlikely (xa < 0x1.0p-28f)) |
55 | { |
56 | math_force_eval (huge + x); |
57 | math_check_force_underflow (x); |
58 | return x; |
59 | } |
60 | |
61 | t = xa + xa; |
62 | t = 0.5f * __log1pf (t + t * xa / (1.0f - xa)); |
63 | } |
64 | else if (__glibc_likely (isless (xa, 1.0f))) |
65 | t = 0.5f * __log1pf ((xa + xa) / (1.0f - xa)); |
66 | else |
67 | { |
68 | if (isgreater (xa, 1.0f)) |
69 | return (x - x) / (x - x); |
70 | |
71 | return x / 0.0f; |
72 | } |
73 | |
74 | return copysignf (t, x); |
75 | } |
76 | libm_alias_finite (__ieee754_atanhf, __atanhf) |
77 | |