math: Use atanhf from CORE-MATH

The CORE-MATH implementation is correctly rounded (for any rounding mode)
and shows slight better performance to the generic atanhf.

The code was adapted to glibc style and to use the definition of
math_config.h (to handle errno, overflow, and underflow).

Benchtest on x64_64 (Ryzen 9 5900X, gcc 14.2.1), aarch64 (Neoverse-N1,
gcc 13.3.1), and powerpc (POWER10, gcc 13.2.1):

Latency                      master        patched   improvement
x86_64                      59.4930        45.8568        22.92%
x86_64v2                    59.5705        45.5804        23.48%
x86_64v3                    53.1838        37.7155        29.08%
i686                        169.354       133.5940        21.12%
aarch64 (Neoverse)          26.0781        16.9829        34.88%
power10                     15.6591        10.7623        31.27%

reciprocal-throughput        master        patched   improvement
x86_64                      23.5903        18.5766        21.25%
x86_64v2                    22.6489        18.2683        19.34%
x86_64v3                    19.0401        13.9474        26.75%
i686                        97.6034       107.3260        -9.96%
aarch64 (Neoverse)          15.3664        9.57846        37.67%
power10                      6.8877        4.6242         32.86%

Signed-off-by: Alexei Sibidanov <sibid@uvic.ca>
Signed-off-by: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Signed-off-by: Adhemerval Zanella <adhemerval.zanella@linaro.org>
Reviewed-by: DJ Delorie <dj@redhat.com>

27 files changed, 157 insertions, 252 deletions

diff --git a/sysdeps/aarch64/libm-test-ulps b/sysdeps/aarch64/libm-test-ulps
index 931fac6..09f7f1c 100644
--- a/sysdeps/aarch64/libm-test-ulps
+++ b/sysdeps/aarch64/libm-test-ulps
@@ -206,7 +206,6 @@ ldouble: 2
 
 Function: "atanh":
 double: 2
-float: 2
 ldouble: 4
 
 Function: "atanh_advsimd":
@@ -215,7 +214,6 @@ float: 1
 
 Function: "atanh_downward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "atanh_sve":
@@ -224,12 +222,10 @@ float: 1
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "atanpi":
diff --git a/sysdeps/alpha/fpu/libm-test-ulps b/sysdeps/alpha/fpu/libm-test-ulps
index d098eaf..af11a87 100644
--- a/sysdeps/alpha/fpu/libm-test-ulps
+++ b/sysdeps/alpha/fpu/libm-test-ulps
@@ -98,22 +98,18 @@ ldouble: 2
 
 Function: "atanh":
 double: 2
-float: 2
 ldouble: 4
 
 Function: "atanh_downward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "cabs":
diff --git a/sysdeps/arc/fpu/libm-test-ulps b/sysdeps/arc/fpu/libm-test-ulps
index 5b89ee2..ef93b0b 100644
--- a/sysdeps/arc/fpu/libm-test-ulps
+++ b/sysdeps/arc/fpu/libm-test-ulps
@@ -75,19 +75,15 @@ double: 2
 
 Function: "atanh":
 double: 2
-float: 2
 
 Function: "atanh_downward":
 double: 3
-float: 3
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 
 Function: "cabs":
 double: 1
diff --git a/sysdeps/arc/nofpu/libm-test-ulps b/sysdeps/arc/nofpu/libm-test-ulps
index ac43ab8..0d2e660 100644
--- a/sysdeps/arc/nofpu/libm-test-ulps
+++ b/sysdeps/arc/nofpu/libm-test-ulps
@@ -20,7 +20,6 @@ Function: "atan2":
 
 Function: "atanh":
 double: 2
-float: 2
 
 Function: "cabs":
 double: 1
diff --git a/sysdeps/arm/libm-test-ulps b/sysdeps/arm/libm-test-ulps
index 234609c..bf16818 100644
--- a/sysdeps/arm/libm-test-ulps
+++ b/sysdeps/arm/libm-test-ulps
@@ -122,19 +122,15 @@ double: 1
 
 Function: "atanh":
 double: 2
-float: 2
 
 Function: "atanh_downward":
 double: 3
-float: 3
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 
 Function: "atanpi":
 double: 1
diff --git a/sysdeps/csky/fpu/libm-test-ulps b/sysdeps/csky/fpu/libm-test-ulps
index b8de7af..f1b62e3 100644
--- a/sysdeps/csky/fpu/libm-test-ulps
+++ b/sysdeps/csky/fpu/libm-test-ulps
@@ -71,19 +71,15 @@ double: 1
 
 Function: "atanh":
 double: 2
-float: 2
 
 Function: "atanh_downward":
 double: 3
-float: 3
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 
 Function: "cabs":
 double: 1
diff --git a/sysdeps/csky/nofpu/libm-test-ulps b/sysdeps/csky/nofpu/libm-test-ulps
index a7594b6..9c2bfc6 100644
--- a/sysdeps/csky/nofpu/libm-test-ulps
+++ b/sysdeps/csky/nofpu/libm-test-ulps
@@ -71,19 +71,15 @@ double: 1
 
 Function: "atanh":
 double: 2
-float: 2
 
 Function: "atanh_downward":
 double: 3
-float: 3
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 
 Function: "cabs":
 double: 1
diff --git a/sysdeps/hppa/fpu/libm-test-ulps b/sysdeps/hppa/fpu/libm-test-ulps
index 5d2074b..f20c59b 100644
--- a/sysdeps/hppa/fpu/libm-test-ulps
+++ b/sysdeps/hppa/fpu/libm-test-ulps
@@ -122,19 +122,15 @@ double: 1
 
 Function: "atanh":
 double: 2
-float: 2
 
 Function: "atanh_downward":
 double: 3
-float: 3
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 
 Function: "atanpi":
 double: 1
diff --git a/sysdeps/i386/fpu/e_atanhf.S b/sysdeps/i386/fpu/e_atanhf.S
deleted file mode 100644
index c94dc7b..0000000
--- a/sysdeps/i386/fpu/e_atanhf.S
+++ /dev/null
@@ -1,110 +0,0 @@
-/* ix87 specific implementation of arctanh function.
-   Copyright (C) 1996-2024 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-
-   The GNU C Library is free software; you can redistribute it and/or
-   modify it under the terms of the GNU Lesser General Public
-   License as published by the Free Software Foundation; either
-   version 2.1 of the License, or (at your option) any later version.
-
-   The GNU C Library is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-   Lesser General Public License for more details.
-
-   You should have received a copy of the GNU Lesser General Public
-   License along with the GNU C Library; if not, see
-   <https://www.gnu.org/licenses/>.  */
-
-#include <machine/asm.h>
-#include <i386-math-asm.h>
-#include <libm-alias-finite.h>
-
-	.section .rodata
-
-	.align ALIGNARG(4)
-	.type half,@object
-half:	.double 0.5
-	ASM_SIZE_DIRECTIVE(half)
-	.type one,@object
-one:	.double 1.0
-	ASM_SIZE_DIRECTIVE(one)
-	.type limit,@object
-limit:	.double 0.29
-	ASM_SIZE_DIRECTIVE(limit)
-	.align ALIGNARG(4)
-	.type ln2_2,@object
-ln2_2:	.quad   0xb17217f7d1cf79ac  /* 0.3465735902799726547086160  */
-	.short  0x3ffd
-	ASM_SIZE_DIRECTIVE(ln2_2)
-
-DEFINE_FLT_MIN
-
-#ifdef PIC
-#define MO(op) op##@GOTOFF(%edx)
-#else
-#define MO(op) op
-#endif
-
-	.text
-ENTRY(__ieee754_atanhf)
-	movl	4(%esp), %ecx
-
-	movl	%ecx, %eax
-	andl	$0x7fffffff, %eax
-	cmpl	$0x7f800000, %eax
-	ja	5f
-
-#ifdef PIC
-	LOAD_PIC_REG (dx)
-#endif
-
-	andl	$0x80000000, %ecx // ECX == 0 iff X >= 0
-
-	fldt	MO(ln2_2)	// 0.5*ln2
-	xorl	%ecx, 4(%esp)
-	flds	4(%esp)		// |x| : 0.5*ln2
-	fcoml	MO(half)	// |x| : 0.5*ln2
-	fld	%st(0)		// |x| : |x| : 0.5*ln2
-	fnstsw			// |x| : |x| : 0.5*ln2
-	sahf
-	jae	2f
-	fadd	%st, %st(1)	// |x| : 2*|x| : 0.5*ln2
-	fld	%st		// |x| : |x| : 2*|x| : 0.5*ln2
-	fsubrl	MO(one)		// 1-|x| : |x| : 2*|x| : 0.5*ln2
-	fxch			// |x| : 1-|x| : 2*|x| : 0.5*ln2
-	fmul	%st(2)		// 2*|x|^2 : 1-|x| : 2*|x| : 0.5*ln2
-	fdivp			// (2*|x|^2)/(1-|x|) : 2*|x| : 0.5*ln2
-	faddp			// 2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2
-	fcoml	MO(limit)	// 2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2
-	fnstsw			// 2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2
-	sahf
-	jae	4f
-	fyl2xp1			// 0.5*ln2*ld(1+2*|x|+(2*|x|^2)/(1-|x|))
-	FLT_CHECK_FORCE_UFLOW_NONNEG
-	jecxz	3f
-	fchs			// 0.5*ln2*ld(1+2*x+(2*x^2)/(1-x))
-3:	ret
-
-	.align ALIGNARG(4)
-4:	faddl	MO(one)		// 1+2*|x|+(2*|x|^2)/(1-|x|) : 0.5*ln2
-	fyl2x			// 0.5*ln2*ld(1+2*|x|+(2*|x|^2)/(1-|x|))
-	jecxz	3f
-	fchs			// 0.5*ln2*ld(1+2*x+(2*x^2)/(1-x))
-3:	ret
-
-	.align ALIGNARG(4)
-2:	faddl	MO(one)		// 1+|x| : |x| : 0.5*ln2
-	fxch			// |x| : 1+|x| : 0.5*ln2
-	fsubrl	MO(one)		// 1-|x| : 1+|x| : 0.5*ln2
-	fdivrp			// (1+|x|)/(1-|x|) : 0.5*ln2
-	fyl2x			// 0.5*ln2*ld((1+|x|)/(1-|x|))
-	jecxz	3f
-	fchs			// 0.5*ln2*ld((1+x)/(1-x))
-3:	ret
-
-	// x == NaN
-5:	flds	4(%esp)
-	ret
-END(__ieee754_atanhf)
-libm_alias_finite (__ieee754_atanhf, __atanhf)
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 3c16d88..de9b783 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -198,19 +198,16 @@ ldouble: 3
 
 Function: "atanh_downward":
 double: 1
-float: 1
 float128: 4
 ldouble: 4
 
 Function: "atanh_towardzero":
 double: 1
-float: 1
 float128: 2
 ldouble: 3
 
 Function: "atanh_upward":
 double: 1
-float: 1
 float128: 4
 ldouble: 5
 
diff --git a/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps b/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps
index 302d5c6..253dcc5 100644
--- a/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps
+++ b/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps
@@ -198,19 +198,16 @@ ldouble: 3
 
 Function: "atanh_downward":
 double: 1
-float: 1
 float128: 4
 ldouble: 4
 
 Function: "atanh_towardzero":
 double: 1
-float: 1
 float128: 2
 ldouble: 3
 
 Function: "atanh_upward":
 double: 1
-float: 1
 float128: 4
 ldouble: 5
 
diff --git a/sysdeps/ieee754/flt-32/e_atanhf.c b/sysdeps/ieee754/flt-32/e_atanhf.c
index d36eb15..5138408 100644
--- a/sysdeps/ieee754/flt-32/e_atanhf.c
+++ b/sysdeps/ieee754/flt-32/e_atanhf.c
@@ -1,76 +1,178 @@
-/* Copyright (C) 2011-2024 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
+/* Correctly-rounded inverse hyperbolic tangent function for binary32 value.
 
-   The GNU C Library is free software; you can redistribute it and/or
-   modify it under the terms of the GNU Lesser General Public
-   License as published by the Free Software Foundation; either
-   version 2.1 of the License, or (at your option) any later version.
+Copyright (c) 2023-2024 Alexei Sibidanov.
 
-   The GNU C Library is distributed in the hope that it will be useful,
-   but WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-   Lesser General Public License for more details.
+The original version of this file was copied from the CORE-MATH
+project (file src/binary32/acosh/acoshf.c, revision bc385c2).
 
-   You should have received a copy of the GNU Lesser General Public
-   License along with the GNU C Library; if not, see
-   <https://www.gnu.org/licenses/>.  */
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
 
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
 
-/* __ieee754_atanh(x)
-   Method :
-      1.Reduced x to positive by atanh(-x) = -atanh(x)
-      2.For x>=0.5
-		    1              2x                          x
-	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
-		    2             1 - x                      1 - x
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
 
-	For x<0.5
-	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
-
-   Special cases:
-	atanh(x) is NaN if |x| > 1 with signal;
-	atanh(NaN) is that NaN with no signal;
-	atanh(+-1) is +-INF with signal.
-
- */
-
-#include <float.h>
-#include <inttypes.h>
 #include <math.h>
-#include <math-barriers.h>
-#include <math_private.h>
-#include <math-underflow.h>
+#include <stdint.h>
 #include <libm-alias-finite.h>
+#include "math_config.h"
 
-static const float huge = 1e30;
+static __attribute__((noinline)) float
+as_special (float x)
+{
+  uint32_t ix = asuint (x);
+  uint32_t ax = ix << 1;
+  if (ax == 0x7f000000u)  /* +-1 */
+    return __math_divzerof (ix >> 31);
+  if (ax > 0xff000000u)
+    return x + x; /* nan */
+  return __math_invalidf (x);
+}
 
 float
 __ieee754_atanhf (float x)
 {
-  float xa = fabsf (x);
-  float t;
-  if (isless (xa, 0.5f))
+  /* Calculate atanh(x) using the difference of two logarithms -- atanh(x) =
+     (ln(1+x) - ln(1-x))/2  */
+  static const double tr[] =
+    {
+      0x1.fc07f02p-1, 0x1.f44659ep-1, 0x1.ecc07b3p-1, 0x1.e573ac9p-1,
+      0x1.de5d6e4p-1, 0x1.d77b655p-1, 0x1.d0cb58fp-1, 0x1.ca4b305p-1,
+      0x1.c3f8f02p-1, 0x1.bdd2b8ap-1, 0x1.b7d6c3ep-1, 0x1.b20364p-1,
+      0x1.ac5701bp-1, 0x1.a6d01a7p-1, 0x1.a16d3f9p-1, 0x1.9c2d14fp-1,
+      0x1.970e4f8p-1, 0x1.920fb4ap-1, 0x1.8d3018dp-1, 0x1.886e5f1p-1,
+      0x1.83c977bp-1, 0x1.7f405fdp-1, 0x1.7ad2209p-1, 0x1.767dce4p-1,
+      0x1.724287fp-1, 0x1.6e1f76bp-1, 0x1.6a13cd1p-1, 0x1.661ec6ap-1,
+      0x1.623fa77p-1, 0x1.5e75bb9p-1, 0x1.5ac056bp-1, 0x1.571ed3cp-1,
+      0x1.5390949p-1, 0x1.5015015p-1, 0x1.4cab887p-1, 0x1.49539e4p-1,
+      0x1.460cbc8p-1, 0x1.42d6626p-1, 0x1.3fb014p-1,  0x1.3c995a4p-1,
+      0x1.3991c2cp-1, 0x1.3698df4p-1, 0x1.33ae45bp-1, 0x1.30d1901p-1,
+      0x1.2e025cp-1,  0x1.2b404adp-1, 0x1.288b013p-1, 0x1.25e2271p-1,
+      0x1.2345679p-1, 0x1.20b470cp-1, 0x1.1e2ef3bp-1, 0x1.1bb4a4p-1,
+      0x1.1945381p-1, 0x1.16e0689p-1, 0x1.1485f0ep-1, 0x1.12358e7p-1,
+      0x1.0fef011p-1, 0x1.0db20a9p-1, 0x1.0b7e6ecp-1, 0x1.0953f39p-1,
+      0x1.073260ap-1, 0x1.05197f8p-1, 0x1.03091b5p-1, 0x1.010101p-1
+    };
+  static const double tl[] =
+    {
+      0x1.fe02a69106789p-9, 0x1.7b91b1155b11bp-7, 0x1.39e87ba1ebd6p-6,
+      0x1.b42dd713971bfp-6, 0x1.16536ee637ae1p-5, 0x1.51b073c96183fp-5,
+      0x1.8c345da019b21p-5, 0x1.c5e5492abc743p-5, 0x1.fec912fbbeabbp-5,
+      0x1.1b72ad33f67ap-4,  0x1.371fc1f6e8f74p-4, 0x1.526e5e5a1b438p-4,
+      0x1.6d60fe601d21dp-4, 0x1.87fa06438c911p-4, 0x1.a23bc223ab563p-4,
+      0x1.bc28673a58cd6p-4, 0x1.d5c216b8fbb91p-4, 0x1.ef0adcaec5936p-4,
+      0x1.040259530d041p-3, 0x1.1058bf8d24ad5p-3, 0x1.1c898c09d99fbp-3,
+      0x1.2895a13e286a3p-3, 0x1.347dd9a447d55p-3, 0x1.404308716a7e4p-3,
+      0x1.4be5f963b78a1p-3, 0x1.5767718015a6cp-3, 0x1.62c82f3a5c795p-3,
+      0x1.6e08eab13a1e4p-3, 0x1.792a55fe147a2p-3, 0x1.842d1d9928b17p-3,
+      0x1.8f11e873a62c7p-3, 0x1.99d958207e08bp-3, 0x1.a484090c1bb0ap-3,
+      0x1.af129324b786bp-3, 0x1.b9858970710fbp-3, 0x1.c3dd7a6ddad4dp-3,
+      0x1.ce1af0b65f3ebp-3, 0x1.d83e725022f3ep-3, 0x1.e2488197c6c26p-3,
+      0x1.ec399d3d68ccp-3,  0x1.f6123fac028acp-3, 0x1.ffd2e07e7f498p-3,
+      0x1.04bdf9e3b26d2p-2, 0x1.0986f4fa93521p-2, 0x1.0e4498651cc8cp-2,
+      0x1.12f719595efbcp-2, 0x1.179eabb0a99a1p-2, 0x1.1c3b81e933c25p-2,
+      0x1.20cdcd0e0ab6ep-2, 0x1.2555bcf50f7cbp-2, 0x1.29d37ff34b08bp-2,
+      0x1.2e47437640268p-2, 0x1.32b1338401d71p-2, 0x1.37117b5c147b6p-2,
+      0x1.3b6844a13fc23p-2, 0x1.3fb5b857f6f42p-2, 0x1.43f9fe2f7ce67p-2,
+      0x1.48353d11488dfp-2, 0x1.4c679b014ee3ap-2, 0x1.50913cc03686bp-2,
+      0x1.54b2468259498p-2, 0x1.58cadb57d7989p-2, 0x1.5cdb1dcaa1765p-2,
+      0x1.60e32f46788d9p-2
+    };
+  static const double ln2n[] =
+    {
+      0x1.62e42fedb2a44p-2, 0x1.62e42feeab21ap-1, 0x1.0a2b23f33e789p+0,
+      0x1.62e42fef27604p+0, 0x1.bb9d3beb1048p+0,  0x1.0a2b23f37c97ep+1,
+      0x1.3687a9f1710bcp+1, 0x1.62e42fef657fap+1, 0x1.8f40b5ed59f38p+1,
+      0x1.bb9d3beb4e676p+1, 0x1.e7f9c1e942db4p+1, 0x1.0a2b23f39ba79p+2,
+      0x1.205966f295e18p+2, 0x1.3687a9f1901b7p+2, 0x1.4cb5ecf08a556p+2,
+      0x1.62e42fef848f5p+2, 0x1.791272ee7ec93p+2, 0x1.8f40b5ed79032p+2,
+      0x1.a56ef8ec733d1p+2, 0x1.bb9d3beb6d77p+2,  0x1.d1cb7eea67b0fp+2,
+      0x1.e7f9c1e961eaep+2, 0x1.fe2804e85c24dp+2, 0x1.0a2b23f3ab2f6p+3
+    };
+  static const double b[] =
     {
-      if (__glibc_unlikely (xa < 0x1.0p-28f))
+      0x1.fffffffce5a6ap-2, -0x1.0001f81ec0ab8p-2, 0x1.555a0f53d79a5p-3
+    };
+  static const double s[] = { 1, -1 };
+
+  uint32_t ux = asuint (x);
+  uint32_t ax = ux << 1;
+  if (__glibc_unlikely (ax < 0x7a300000u || ax >= 0x7f000000u))
+    {
+      if (__glibc_unlikely (ax >= 0x7f000000u))
+	return as_special (x);
+      if (__glibc_unlikely (ax < 0x73713744u))
 	{
-	  math_force_eval (huge + x);
-	  math_check_force_underflow (x);
-	  return x;
+	  if (ax == 0)
+	    return x; /* x = +-0 */
+	  return fmaf (x, 0x1p-25f, x); /* |x| < 0.000352112(0x1.713744p-12) */
+	}
+      else
+	{ /* |x| < 0x1.3p-5 */
+	  static const double c[] =
+	    {
+	      0x1.5555555555527p-2, 0x1.9999999ba4ee8p-3,
+	      0x1.24922c280990ap-3, 0x1.c8236aae809c6p-4
+	    };
+	  double z = x;
+	  double z2 = z * z;
+	  double z4 = z2 * z2;
+	  double r = c[0] + z2 * c[1] + z4 * (c[2] + z2 * c[3]);
+	  return z + (z * z2) * r;
 	}
-
-      t = xa + xa;
-      t = 0.5f * __log1pf (t + t * xa / (1.0f - xa));
     }
-  else if (__glibc_likely (isless (xa, 1.0f)))
-    t = 0.5f * __log1pf ((xa + xa) / (1.0f - xa));
-  else
+  double sgn = s[ux >> 31];
+  unsigned int e = ax >> 24;
+  unsigned int md = ((ux << 8) | 1 << 31) >> (126 - e);
+  unsigned int mn = -md;
+  int nz = __builtin_clz (mn) + 1;
+  mn <<= nz;
+  unsigned int jn = mn >> 26;
+  unsigned int jd = md >> 26;
+  double tn = asdouble (((int64_t) mn << 20) | (INT64_C(1023) << 52));
+  double td = asdouble (((int64_t) md << 20) | (INT64_C(1023) << 52));
+  double zn = tn * tr[jn] - 1;
+  double zd = td * tr[jd] - 1;
+  double zn2 = zn * zn;
+  double zd2 = zd * zd;
+  double rn = ((tl[jn] - ln2n[nz - 1]) + zn * b[0]) + zn2 * (b[1] + zn * b[2]);
+  double rd = (tl[jd] + zd * b[0]) + zd2 * (b[1] + zd * b[2]);
+  double r = sgn * (rd - rn);
+  float ub = r;
+  float lb = r + sgn * 0.226e-9;
+  if (__glibc_unlikely (ub != lb))
     {
-      if (isgreater (xa, 1.0f))
-	return (x - x) / (x - x);
-
-      return x / 0.0f;
+      static const double c[] =
+	{
+	   0x1p-1,               -0x1.000000000001bp-2,  0x1.55555555555bap-3,
+	  -0x1.fffffff26d72ep-4,  0x1.99999989035p-4,   -0x1.555c39cb9ee8p-4,
+	   0x1.24992d8b014a1p-4
+	};
+      double zn4 = zn2 * zn2;
+      double zd4 = zd2 * zd2;
+      double fn = zn * (((c[0] + zn * c[1]) + zn2 * (c[2] + zn * c[3]))
+		     + zn4 * ((c[4] + zn * c[5]) + zn2 * c[6]));
+      fn += 0x1.0ca86c3898dp-50 * nz;
+      fn += tl[jn];
+      double en = nz * 0x1.62e42fefa3ap-2;
+      double fd = zd * (((c[0] + zd * c[1]) + zd2 * (c[2] + zd * c[3]))
+		     + zd4 * ((c[4] + zd * c[5]) + zd2 * c[6]));
+      fd += tl[jd];
+      r = fd - fn + en;
+      ub = sgn * r;
     }
-
-  return copysignf (t, x);
+  return ub;
 }
 libm_alias_finite (__ieee754_atanhf, __atanhf)
diff --git a/sysdeps/loongarch/lp64/libm-test-ulps b/sysdeps/loongarch/lp64/libm-test-ulps
index a345397..1e1a289 100644
--- a/sysdeps/loongarch/lp64/libm-test-ulps
+++ b/sysdeps/loongarch/lp64/libm-test-ulps
@@ -98,22 +98,18 @@ ldouble: 2
 
 Function: "atanh":
 double: 2
-float: 2
 ldouble: 4
 
 Function: "atanh_downward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "cabs":
diff --git a/sysdeps/m68k/coldfire/fpu/libm-test-ulps b/sysdeps/m68k/coldfire/fpu/libm-test-ulps
index 989015f..998088f 100644
--- a/sysdeps/m68k/coldfire/fpu/libm-test-ulps
+++ b/sysdeps/m68k/coldfire/fpu/libm-test-ulps
@@ -4,7 +4,6 @@
 Function: "atan2":
 
 Function: "atanh":
-float: 1
 
 Function: Imaginary part of "cacosh":
 float: 1
diff --git a/sysdeps/microblaze/libm-test-ulps b/sysdeps/microblaze/libm-test-ulps
index 4fa242d..77018f4 100644
--- a/sysdeps/microblaze/libm-test-ulps
+++ b/sysdeps/microblaze/libm-test-ulps
@@ -17,7 +17,6 @@ Function: "atan2":
 
 Function: "atanh":
 double: 2
-float: 2
 
 Function: "cabs":
 double: 1
diff --git a/sysdeps/mips/mips32/libm-test-ulps b/sysdeps/mips/mips32/libm-test-ulps
index 2f6f7fff..2191d57 100644
--- a/sysdeps/mips/mips32/libm-test-ulps
+++ b/sysdeps/mips/mips32/libm-test-ulps
@@ -74,19 +74,15 @@ double: 1
 
 Function: "atanh":
 double: 2
-float: 2
 
 Function: "atanh_downward":
 double: 3
-float: 3
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 
 Function: "cabs":
 double: 1
diff --git a/sysdeps/mips/mips64/libm-test-ulps b/sysdeps/mips/mips64/libm-test-ulps
index 037fb88..7d789b9 100644
--- a/sysdeps/mips/mips64/libm-test-ulps
+++ b/sysdeps/mips/mips64/libm-test-ulps
@@ -98,22 +98,18 @@ ldouble: 2
 
 Function: "atanh":
 double: 2
-float: 2
 ldouble: 4
 
 Function: "atanh_downward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "cabs":
diff --git a/sysdeps/or1k/fpu/libm-test-ulps b/sysdeps/or1k/fpu/libm-test-ulps
index 3fee2b5..357d0a6 100644
--- a/sysdeps/or1k/fpu/libm-test-ulps
+++ b/sysdeps/or1k/fpu/libm-test-ulps
@@ -74,19 +74,15 @@ double: 1
 
 Function: "atanh":
 double: 2
-float: 2
 
 Function: "atanh_downward":
 double: 3
-float: 3
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 
 Function: "cabs":
 double: 1
diff --git a/sysdeps/or1k/nofpu/libm-test-ulps b/sysdeps/or1k/nofpu/libm-test-ulps
index a380709..ac488ab 100644
--- a/sysdeps/or1k/nofpu/libm-test-ulps
+++ b/sysdeps/or1k/nofpu/libm-test-ulps
@@ -86,19 +86,15 @@ double: 1
 
 Function: "atanh":
 double: 2
-float: 2
 
 Function: "atanh_downward":
 double: 3
-float: 3
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 
 Function: "atanpi":
 double: 1
diff --git a/sysdeps/powerpc/fpu/libm-test-ulps b/sysdeps/powerpc/fpu/libm-test-ulps
index dd166aa..4165161 100644
--- a/sysdeps/powerpc/fpu/libm-test-ulps
+++ b/sysdeps/powerpc/fpu/libm-test-ulps
@@ -198,25 +198,21 @@ ldouble: 2
 
 Function: "atanh":
 double: 2
-float: 2
 float128: 4
 ldouble: 2
 
 Function: "atanh_downward":
 double: 3
-float: 3
 float128: 4
 ldouble: 3
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 float128: 2
 ldouble: 4
 
 Function: "atanh_upward":
 double: 3
-float: 3
 float128: 4
 ldouble: 4
 
diff --git a/sysdeps/powerpc/nofpu/libm-test-ulps b/sysdeps/powerpc/nofpu/libm-test-ulps
index 8b7ed72..b13e465 100644
--- a/sysdeps/powerpc/nofpu/libm-test-ulps
+++ b/sysdeps/powerpc/nofpu/libm-test-ulps
@@ -102,22 +102,18 @@ ldouble: 2
 
 Function: "atanh":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "atanh_downward":
 double: 3
-float: 3
 ldouble: 3
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 ldouble: 4
 
 Function: "atanh_upward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "cabs":
diff --git a/sysdeps/riscv/nofpu/libm-test-ulps b/sysdeps/riscv/nofpu/libm-test-ulps
index 19e3674..dc78136 100644
--- a/sysdeps/riscv/nofpu/libm-test-ulps
+++ b/sysdeps/riscv/nofpu/libm-test-ulps
@@ -98,22 +98,18 @@ ldouble: 2
 
 Function: "atanh":
 double: 2
-float: 2
 ldouble: 4
 
 Function: "atanh_downward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "cabs":
diff --git a/sysdeps/riscv/rvd/libm-test-ulps b/sysdeps/riscv/rvd/libm-test-ulps
index ad3242f..9477c4a 100644
--- a/sysdeps/riscv/rvd/libm-test-ulps
+++ b/sysdeps/riscv/rvd/libm-test-ulps
@@ -98,22 +98,18 @@ ldouble: 2
 
 Function: "atanh":
 double: 2
-float: 2
 ldouble: 4
 
 Function: "atanh_downward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "cabs":
diff --git a/sysdeps/s390/fpu/libm-test-ulps b/sysdeps/s390/fpu/libm-test-ulps
index 7b43058..d7be586 100644
--- a/sysdeps/s390/fpu/libm-test-ulps
+++ b/sysdeps/s390/fpu/libm-test-ulps
@@ -98,22 +98,18 @@ ldouble: 2
 
 Function: "atanh":
 double: 2
-float: 2
 ldouble: 4
 
 Function: "atanh_downward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "cabs":
diff --git a/sysdeps/sh/libm-test-ulps b/sysdeps/sh/libm-test-ulps
index 15655ad..4b308af 100644
--- a/sysdeps/sh/libm-test-ulps
+++ b/sysdeps/sh/libm-test-ulps
@@ -35,11 +35,9 @@ double: 1
 
 Function: "atanh":
 double: 2
-float: 2
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 
 Function: "cabs":
 double: 1
diff --git a/sysdeps/sparc/fpu/libm-test-ulps b/sysdeps/sparc/fpu/libm-test-ulps
index 0d86746..324fd21 100644
--- a/sysdeps/sparc/fpu/libm-test-ulps
+++ b/sysdeps/sparc/fpu/libm-test-ulps
@@ -98,22 +98,18 @@ ldouble: 2
 
 Function: "atanh":
 double: 2
-float: 2
 ldouble: 4
 
 Function: "atanh_downward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "atanh_upward":
 double: 3
-float: 3
 ldouble: 4
 
 Function: "cabs":
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index f09e128..8c1fafe 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -311,25 +311,21 @@ float: 1
 
 Function: "atanh":
 double: 2
-float: 2
 float128: 4
 ldouble: 3
 
 Function: "atanh_downward":
 double: 3
-float: 3
 float128: 4
 ldouble: 5
 
 Function: "atanh_towardzero":
 double: 2
-float: 2
 float128: 2
 ldouble: 4
 
 Function: "atanh_upward":
 double: 3
-float: 3
 float128: 4
 ldouble: 5