[PATCH] math: Update tanh from CORE-MATH (revision cf237fa0)
Adhemerval Zanella <[email protected]>
| Newsgroups | gmane.comp.lib.glibc.alpha |
|---|---|
| Message-ID | <[email protected]> |
Sync the double-precision tanh implementation with CORE-MATH revision
cf237fa0.
Checked on aarch64-linux-gnu.
---
SHARED-FILES | 2 +-
math/auto-libm-test-in | 2 +
math/auto-libm-test-out-tanh | 138 +++++++++++++++++++++++++++
sysdeps/ieee754/dbl-64/ddcoremath.h | 10 ++
sysdeps/ieee754/dbl-64/e_tanh_data.c | 29 +++---
sysdeps/ieee754/dbl-64/e_tanh_data.h | 6 +-
sysdeps/ieee754/dbl-64/s_tanh.c | 111 ++++++++++++---------
7 files changed, 232 insertions(+), 66 deletions(-)
diff --git a/SHARED-FILES b/SHARED-FILES
index 3ab7d41d989..cf2a19c015b 100644
--- a/SHARED-FILES
+++ b/SHARED-FILES
@@ -280,7 +280,7 @@ core-math:
sysdeps/ieee754/dbl-64/s_erf.c
# src/binary64/erfc/erfc.c, revision 55e9869e
sysdeps/ieee754/dbl-64/s_erfc.c
- # src/binary64/tanh/tanh.c, revision 8ea8ea35
+ # src/binary64/tanh/tanh.c, revision cf237fa0
sysdeps/ieee754/dbl-64/s_tanh.c
# src/binary32/acos/acosf.c, revision 56dd347
sysdeps/ieee754/flt-32/e_acosf.c
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 73737235de2..ca670768a29 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -10109,6 +10109,8 @@ tanh 0x1p-10000
tanh 0x0.fffffffffffffp-1022
# the next value generates larger error bounds on x86_64 (ldbl-96)
tanh 0x3.b9979a543d0fbfa8p-4
+tanh 0x1.a0112a16e9318p+1
+tanh -0x1.a0112a16e9318p+1
tanh max
tanh -max
tanh min
diff --git a/math/auto-libm-test-out-tanh b/math/auto-libm-test-out-tanh
index 04ecb568463..2f2a9166c0f 100644
--- a/math/auto-libm-test-out-tanh
+++ b/math/auto-libm-test-out-tanh
@@ -2463,6 +2463,144 @@ tanh 0x3.b9979a543d0fbfa8p-4
= tanh tonearest ibm128 0x3.b9979a543d0fbfa8p-4 : 0x3.a8baae38037e30ab26bdc37e74p-4 : inexact-ok
= tanh towardzero ibm128 0x3.b9979a543d0fbfa8p-4 : 0x3.a8baae38037e30ab26bdc37e74p-4 : inexact-ok
= tanh upward ibm128 0x3.b9979a543d0fbfa8p-4 : 0x3.a8baae38037e30ab26bdc37e75p-4 : inexact-ok
+tanh 0x1.a0112a16e9318p+1
+= tanh downward binary32 0x3.402258p+0 : 0xf.f3b71p-4 : inexact-ok
+= tanh tonearest binary32 0x3.402258p+0 : 0xf.f3b71p-4 : inexact-ok
+= tanh towardzero binary32 0x3.402258p+0 : 0xf.f3b71p-4 : inexact-ok
+= tanh upward binary32 0x3.402258p+0 : 0xf.f3b72p-4 : inexact-ok
+= tanh downward binary64 0x3.402258p+0 : 0xf.f3b715e6c4268p-4 : inexact-ok
+= tanh tonearest binary64 0x3.402258p+0 : 0xf.f3b715e6c427p-4 : inexact-ok
+= tanh towardzero binary64 0x3.402258p+0 : 0xf.f3b715e6c4268p-4 : inexact-ok
+= tanh upward binary64 0x3.402258p+0 : 0xf.f3b715e6c427p-4 : inexact-ok
+= tanh downward intel96 0x3.402258p+0 : 0xf.f3b715e6c426c87p-4 : inexact-ok
+= tanh tonearest intel96 0x3.402258p+0 : 0xf.f3b715e6c426c88p-4 : inexact-ok
+= tanh towardzero intel96 0x3.402258p+0 : 0xf.f3b715e6c426c87p-4 : inexact-ok
+= tanh upward intel96 0x3.402258p+0 : 0xf.f3b715e6c426c88p-4 : inexact-ok
+= tanh downward m68k96 0x3.402258p+0 : 0xf.f3b715e6c426c87p-4 : inexact-ok
+= tanh tonearest m68k96 0x3.402258p+0 : 0xf.f3b715e6c426c88p-4 : inexact-ok
+= tanh towardzero m68k96 0x3.402258p+0 : 0xf.f3b715e6c426c87p-4 : inexact-ok
+= tanh upward m68k96 0x3.402258p+0 : 0xf.f3b715e6c426c88p-4 : inexact-ok
+= tanh downward binary128 0x3.402258p+0 : 0xf.f3b715e6c426c87d0ef0472d18ap-4 : inexact-ok
+= tanh tonearest binary128 0x3.402258p+0 : 0xf.f3b715e6c426c87d0ef0472d18a8p-4 : inexact-ok
+= tanh towardzero binary128 0x3.402258p+0 : 0xf.f3b715e6c426c87d0ef0472d18ap-4 : inexact-ok
+= tanh upward binary128 0x3.402258p+0 : 0xf.f3b715e6c426c87d0ef0472d18a8p-4 : inexact-ok
+= tanh downward ibm128 0x3.402258p+0 : 0xf.f3b715e6c426c87d0ef0472d18p-4 : inexact-ok
+= tanh tonearest ibm128 0x3.402258p+0 : 0xf.f3b715e6c426c87d0ef0472d18p-4 : inexact-ok
+= tanh towardzero ibm128 0x3.402258p+0 : 0xf.f3b715e6c426c87d0ef0472d18p-4 : inexact-ok
+= tanh upward ibm128 0x3.402258p+0 : 0xf.f3b715e6c426c87d0ef0472d1cp-4 : inexact-ok
+= tanh downward binary32 0x3.402254p+0 : 0xf.f3b71p-4 : inexact-ok
+= tanh tonearest binary32 0x3.402254p+0 : 0xf.f3b71p-4 : inexact-ok
+= tanh towardzero binary32 0x3.402254p+0 : 0xf.f3b71p-4 : inexact-ok
+= tanh upward binary32 0x3.402254p+0 : 0xf.f3b72p-4 : inexact-ok
+= tanh downward binary64 0x3.402254p+0 : 0xf.f3b71584a28fp-4 : inexact-ok
+= tanh tonearest binary64 0x3.402254p+0 : 0xf.f3b71584a28fp-4 : inexact-ok
+= tanh towardzero binary64 0x3.402254p+0 : 0xf.f3b71584a28fp-4 : inexact-ok
+= tanh upward binary64 0x3.402254p+0 : 0xf.f3b71584a28f8p-4 : inexact-ok
+= tanh downward intel96 0x3.402254p+0 : 0xf.f3b71584a28f24fp-4 : inexact-ok
+= tanh tonearest intel96 0x3.402254p+0 : 0xf.f3b71584a28f25p-4 : inexact-ok
+= tanh towardzero intel96 0x3.402254p+0 : 0xf.f3b71584a28f24fp-4 : inexact-ok
+= tanh upward intel96 0x3.402254p+0 : 0xf.f3b71584a28f25p-4 : inexact-ok
+= tanh downward m68k96 0x3.402254p+0 : 0xf.f3b71584a28f24fp-4 : inexact-ok
+= tanh tonearest m68k96 0x3.402254p+0 : 0xf.f3b71584a28f25p-4 : inexact-ok
+= tanh towardzero m68k96 0x3.402254p+0 : 0xf.f3b71584a28f24fp-4 : inexact-ok
+= tanh upward m68k96 0x3.402254p+0 : 0xf.f3b71584a28f25p-4 : inexact-ok
+= tanh downward binary128 0x3.402254p+0 : 0xf.f3b71584a28f24fdc24f0ecd71c8p-4 : inexact-ok
+= tanh tonearest binary128 0x3.402254p+0 : 0xf.f3b71584a28f24fdc24f0ecd71dp-4 : inexact-ok
+= tanh towardzero binary128 0x3.402254p+0 : 0xf.f3b71584a28f24fdc24f0ecd71c8p-4 : inexact-ok
+= tanh upward binary128 0x3.402254p+0 : 0xf.f3b71584a28f24fdc24f0ecd71dp-4 : inexact-ok
+= tanh downward ibm128 0x3.402254p+0 : 0xf.f3b71584a28f24fdc24f0ecd7p-4 : inexact-ok
+= tanh tonearest ibm128 0x3.402254p+0 : 0xf.f3b71584a28f24fdc24f0ecd7p-4 : inexact-ok
+= tanh towardzero ibm128 0x3.402254p+0 : 0xf.f3b71584a28f24fdc24f0ecd7p-4 : inexact-ok
+= tanh upward ibm128 0x3.402254p+0 : 0xf.f3b71584a28f24fdc24f0ecd74p-4 : inexact-ok
+= tanh downward binary64 0x3.4022542dd263p+0 : 0xf.f3b7158906b28p-4 : inexact-ok
+= tanh tonearest binary64 0x3.4022542dd263p+0 : 0xf.f3b7158906b28p-4 : inexact-ok
+= tanh towardzero binary64 0x3.4022542dd263p+0 : 0xf.f3b7158906b28p-4 : inexact-ok
+= tanh upward binary64 0x3.4022542dd263p+0 : 0xf.f3b7158906b3p-4 : inexact-ok
+= tanh downward intel96 0x3.4022542dd263p+0 : 0xf.f3b7158906b28p-4 : inexact-ok
+= tanh tonearest intel96 0x3.4022542dd263p+0 : 0xf.f3b7158906b28p-4 : inexact-ok
+= tanh towardzero intel96 0x3.4022542dd263p+0 : 0xf.f3b7158906b28p-4 : inexact-ok
+= tanh upward intel96 0x3.4022542dd263p+0 : 0xf.f3b7158906b2801p-4 : inexact-ok
+= tanh downward m68k96 0x3.4022542dd263p+0 : 0xf.f3b7158906b28p-4 : inexact-ok
+= tanh tonearest m68k96 0x3.4022542dd263p+0 : 0xf.f3b7158906b28p-4 : inexact-ok
+= tanh towardzero m68k96 0x3.4022542dd263p+0 : 0xf.f3b7158906b28p-4 : inexact-ok
+= tanh upward m68k96 0x3.4022542dd263p+0 : 0xf.f3b7158906b2801p-4 : inexact-ok
+= tanh downward binary128 0x3.4022542dd263p+0 : 0xf.f3b7158906b2800000506e94ef6p-4 : inexact-ok
+= tanh tonearest binary128 0x3.4022542dd263p+0 : 0xf.f3b7158906b2800000506e94ef68p-4 : inexact-ok
+= tanh towardzero binary128 0x3.4022542dd263p+0 : 0xf.f3b7158906b2800000506e94ef6p-4 : inexact-ok
+= tanh upward binary128 0x3.4022542dd263p+0 : 0xf.f3b7158906b2800000506e94ef68p-4 : inexact-ok
+= tanh downward ibm128 0x3.4022542dd263p+0 : 0xf.f3b7158906b2800000506e94ecp-4 : inexact-ok
+= tanh tonearest ibm128 0x3.4022542dd263p+0 : 0xf.f3b7158906b2800000506e94fp-4 : inexact-ok
+= tanh towardzero ibm128 0x3.4022542dd263p+0 : 0xf.f3b7158906b2800000506e94ecp-4 : inexact-ok
+= tanh upward ibm128 0x3.4022542dd263p+0 : 0xf.f3b7158906b2800000506e94fp-4 : inexact-ok
+tanh -0x1.a0112a16e9318p+1
+= tanh downward binary32 -0x3.402254p+0 : -0xf.f3b72p-4 : inexact-ok
+= tanh tonearest binary32 -0x3.402254p+0 : -0xf.f3b71p-4 : inexact-ok
+= tanh towardzero binary32 -0x3.402254p+0 : -0xf.f3b71p-4 : inexact-ok
+= tanh upward binary32 -0x3.402254p+0 : -0xf.f3b71p-4 : inexact-ok
+= tanh downward binary64 -0x3.402254p+0 : -0xf.f3b71584a28f8p-4 : inexact-ok
+= tanh tonearest binary64 -0x3.402254p+0 : -0xf.f3b71584a28fp-4 : inexact-ok
+= tanh towardzero binary64 -0x3.402254p+0 : -0xf.f3b71584a28fp-4 : inexact-ok
+= tanh upward binary64 -0x3.402254p+0 : -0xf.f3b71584a28fp-4 : inexact-ok
+= tanh downward intel96 -0x3.402254p+0 : -0xf.f3b71584a28f25p-4 : inexact-ok
+= tanh tonearest intel96 -0x3.402254p+0 : -0xf.f3b71584a28f25p-4 : inexact-ok
+= tanh towardzero intel96 -0x3.402254p+0 : -0xf.f3b71584a28f24fp-4 : inexact-ok
+= tanh upward intel96 -0x3.402254p+0 : -0xf.f3b71584a28f24fp-4 : inexact-ok
+= tanh downward m68k96 -0x3.402254p+0 : -0xf.f3b71584a28f25p-4 : inexact-ok
+= tanh tonearest m68k96 -0x3.402254p+0 : -0xf.f3b71584a28f25p-4 : inexact-ok
+= tanh towardzero m68k96 -0x3.402254p+0 : -0xf.f3b71584a28f24fp-4 : inexact-ok
+= tanh upward m68k96 -0x3.402254p+0 : -0xf.f3b71584a28f24fp-4 : inexact-ok
+= tanh downward binary128 -0x3.402254p+0 : -0xf.f3b71584a28f24fdc24f0ecd71dp-4 : inexact-ok
+= tanh tonearest binary128 -0x3.402254p+0 : -0xf.f3b71584a28f24fdc24f0ecd71dp-4 : inexact-ok
+= tanh towardzero binary128 -0x3.402254p+0 : -0xf.f3b71584a28f24fdc24f0ecd71c8p-4 : inexact-ok
+= tanh upward binary128 -0x3.402254p+0 : -0xf.f3b71584a28f24fdc24f0ecd71c8p-4 : inexact-ok
+= tanh downward ibm128 -0x3.402254p+0 : -0xf.f3b71584a28f24fdc24f0ecd74p-4 : inexact-ok
+= tanh tonearest ibm128 -0x3.402254p+0 : -0xf.f3b71584a28f24fdc24f0ecd7p-4 : inexact-ok
+= tanh towardzero ibm128 -0x3.402254p+0 : -0xf.f3b71584a28f24fdc24f0ecd7p-4 : inexact-ok
+= tanh upward ibm128 -0x3.402254p+0 : -0xf.f3b71584a28f24fdc24f0ecd7p-4 : inexact-ok
+= tanh downward binary32 -0x3.402258p+0 : -0xf.f3b72p-4 : inexact-ok
+= tanh tonearest binary32 -0x3.402258p+0 : -0xf.f3b71p-4 : inexact-ok
+= tanh towardzero binary32 -0x3.402258p+0 : -0xf.f3b71p-4 : inexact-ok
+= tanh upward binary32 -0x3.402258p+0 : -0xf.f3b71p-4 : inexact-ok
+= tanh downward binary64 -0x3.402258p+0 : -0xf.f3b715e6c427p-4 : inexact-ok
+= tanh tonearest binary64 -0x3.402258p+0 : -0xf.f3b715e6c427p-4 : inexact-ok
+= tanh towardzero binary64 -0x3.402258p+0 : -0xf.f3b715e6c4268p-4 : inexact-ok
+= tanh upward binary64 -0x3.402258p+0 : -0xf.f3b715e6c4268p-4 : inexact-ok
+= tanh downward intel96 -0x3.402258p+0 : -0xf.f3b715e6c426c88p-4 : inexact-ok
+= tanh tonearest intel96 -0x3.402258p+0 : -0xf.f3b715e6c426c88p-4 : inexact-ok
+= tanh towardzero intel96 -0x3.402258p+0 : -0xf.f3b715e6c426c87p-4 : inexact-ok
+= tanh upward intel96 -0x3.402258p+0 : -0xf.f3b715e6c426c87p-4 : inexact-ok
+= tanh downward m68k96 -0x3.402258p+0 : -0xf.f3b715e6c426c88p-4 : inexact-ok
+= tanh tonearest m68k96 -0x3.402258p+0 : -0xf.f3b715e6c426c88p-4 : inexact-ok
+= tanh towardzero m68k96 -0x3.402258p+0 : -0xf.f3b715e6c426c87p-4 : inexact-ok
+= tanh upward m68k96 -0x3.402258p+0 : -0xf.f3b715e6c426c87p-4 : inexact-ok
+= tanh downward binary128 -0x3.402258p+0 : -0xf.f3b715e6c426c87d0ef0472d18a8p-4 : inexact-ok
+= tanh tonearest binary128 -0x3.402258p+0 : -0xf.f3b715e6c426c87d0ef0472d18a8p-4 : inexact-ok
+= tanh towardzero binary128 -0x3.402258p+0 : -0xf.f3b715e6c426c87d0ef0472d18ap-4 : inexact-ok
+= tanh upward binary128 -0x3.402258p+0 : -0xf.f3b715e6c426c87d0ef0472d18ap-4 : inexact-ok
+= tanh downward ibm128 -0x3.402258p+0 : -0xf.f3b715e6c426c87d0ef0472d1cp-4 : inexact-ok
+= tanh tonearest ibm128 -0x3.402258p+0 : -0xf.f3b715e6c426c87d0ef0472d18p-4 : inexact-ok
+= tanh towardzero ibm128 -0x3.402258p+0 : -0xf.f3b715e6c426c87d0ef0472d18p-4 : inexact-ok
+= tanh upward ibm128 -0x3.402258p+0 : -0xf.f3b715e6c426c87d0ef0472d18p-4 : inexact-ok
+= tanh downward binary64 -0x3.4022542dd263p+0 : -0xf.f3b7158906b3p-4 : inexact-ok
+= tanh tonearest binary64 -0x3.4022542dd263p+0 : -0xf.f3b7158906b28p-4 : inexact-ok
+= tanh towardzero binary64 -0x3.4022542dd263p+0 : -0xf.f3b7158906b28p-4 : inexact-ok
+= tanh upward binary64 -0x3.4022542dd263p+0 : -0xf.f3b7158906b28p-4 : inexact-ok
+= tanh downward intel96 -0x3.4022542dd263p+0 : -0xf.f3b7158906b2801p-4 : inexact-ok
+= tanh tonearest intel96 -0x3.4022542dd263p+0 : -0xf.f3b7158906b28p-4 : inexact-ok
+= tanh towardzero intel96 -0x3.4022542dd263p+0 : -0xf.f3b7158906b28p-4 : inexact-ok
+= tanh upward intel96 -0x3.4022542dd263p+0 : -0xf.f3b7158906b28p-4 : inexact-ok
+= tanh downward m68k96 -0x3.4022542dd263p+0 : -0xf.f3b7158906b2801p-4 : inexact-ok
+= tanh tonearest m68k96 -0x3.4022542dd263p+0 : -0xf.f3b7158906b28p-4 : inexact-ok
+= tanh towardzero m68k96 -0x3.4022542dd263p+0 : -0xf.f3b7158906b28p-4 : inexact-ok
+= tanh upward m68k96 -0x3.4022542dd263p+0 : -0xf.f3b7158906b28p-4 : inexact-ok
+= tanh downward binary128 -0x3.4022542dd263p+0 : -0xf.f3b7158906b2800000506e94ef68p-4 : inexact-ok
+= tanh tonearest binary128 -0x3.4022542dd263p+0 : -0xf.f3b7158906b2800000506e94ef68p-4 : inexact-ok
+= tanh towardzero binary128 -0x3.4022542dd263p+0 : -0xf.f3b7158906b2800000506e94ef6p-4 : inexact-ok
+= tanh upward binary128 -0x3.4022542dd263p+0 : -0xf.f3b7158906b2800000506e94ef6p-4 : inexact-ok
+= tanh downward ibm128 -0x3.4022542dd263p+0 : -0xf.f3b7158906b2800000506e94fp-4 : inexact-ok
+= tanh tonearest ibm128 -0x3.4022542dd263p+0 : -0xf.f3b7158906b2800000506e94fp-4 : inexact-ok
+= tanh towardzero ibm128 -0x3.4022542dd263p+0 : -0xf.f3b7158906b2800000506e94ecp-4 : inexact-ok
+= tanh upward ibm128 -0x3.4022542dd263p+0 : -0xf.f3b7158906b2800000506e94ecp-4 : inexact-ok
tanh max
= tanh downward binary32 0xf.fffffp+124 : 0xf.fffffp-4 : inexact-ok
= tanh tonearest binary32 0xf.fffffp+124 : 0x1p+0 : inexact-ok
diff --git a/sysdeps/ieee754/dbl-64/ddcoremath.h b/sysdeps/ieee754/dbl-64/ddcoremath.h
index 6a31a005c3a..28a910592ed 100644
--- a/sysdeps/ieee754/dbl-64/ddcoremath.h
+++ b/sysdeps/ieee754/dbl-64/ddcoremath.h
@@ -176,6 +176,16 @@ mulddd (double xh, double xl, double ch, double *l)
return ch;
}
+/* Same as mulddd, except the final normalization uses the fasttwosum
+ variant discussed in the muldd_acc2 comment above. */
+static inline double
+mulddd_acc (double xh, double xl, double ch, double *l)
+{
+ double ahlh = ch * xl, ahhh = ch * xh, ahhl = fma (ch, xh, -ahhh);
+ ahhl += ahlh;
+ return fasttwosum (ahhh, ahhl, l);
+}
+
static inline double
mulddd2 (double x, double ch, double cl, double *l)
{
diff --git a/sysdeps/ieee754/dbl-64/e_tanh_data.c b/sysdeps/ieee754/dbl-64/e_tanh_data.c
index 0ebd8277df2..5429a866a04 100644
--- a/sysdeps/ieee754/dbl-64/e_tanh_data.c
+++ b/sysdeps/ieee754/dbl-64/e_tanh_data.c
@@ -1,6 +1,6 @@
/* Data definitions used by tanh implementation.
-Copyright (c) 2023-2026 Alexei Sibidanov.
+Copyright (c) 2023-2026 Alexei Sibidanov, Cyprien Peignier, Paul Zimmermann.
This file is part of the CORE-MATH project
(https://core-math.gitlabpages.inria.fr/).
@@ -26,24 +26,23 @@ SOFTWARE.
#include "e_tanh_data.h"
-const double __tanh_data_ch[][2] =
+const double __tanh_data_ch[9][2] =
{
- { -0x1.5555555555555p-2, -0x1.5555555555555p-56 },
- { 0x1.1111111111111p-3, 0x1.1111111110916p-59 },
- { -0x1.ba1ba1ba1ba1cp-5, 0x1.7917917a46f2cp-59 },
- { 0x1.664f4882c10fap-6, -0x1.9a52a06f1e599p-63 },
- { -0x1.226e355e6c23dp-7, 0x1.c297394c24e38p-61 },
- { 0x1.d6d3d0e157dep-9, -0x1.311087e5b1526p-63 },
- { -0x1.7da36452b75e1p-10, -0x1.2868cde54ea0cp-65 },
- { 0x1.355824803667bp-11, 0x1.2cd8fc406c3f7p-66 },
- { -0x1.f57d7734c821dp-13, 0x1.da22861b4ca8p-70 },
- { 0x1.967e18ad3facfp-14, -0x1.0831108273a74p-68 }
+ { -0x1.5555555555555p-2, -0x1.5555555554cc4p-56 },
+ { 0x1.1111111111111p-3, 0x1.111110f8c0178p-59 },
+ { -0x1.ba1ba1ba1ba1cp-5, 0x1.7917c1d676ff5p-59 },
+ { 0x1.664f4882c10fap-6, -0x1.9d5cb27c0af28p-63 },
+ { -0x1.226e355e6c23cp-7, -0x1.c9674586913f3p-61 },
+ { 0x1.d6d3d0e157db3p-9, -0x1.71376fa06ce94p-65 },
+ { -0x1.7da36452b5e46p-10, -0x1.aba8d51bd9cp-65 },
+ { 0x1.3558247faa32dp-11, -0x1.e0cfb423aedfdp-65 },
+ { -0x1.f57d76ea30928p-13, -0x1.c30601213cae9p-67 }
};
-const double __tanh_data_cl[] =
+const double __tanh_data_cl[5] =
{
- -0x1.497d8e6462927p-15, 0x1.0b1318c243bd7p-16, -0x1.b0f2935e9a12p-18,
- 0x1.5e9444536e654p-19, -0x1.174ff2a31908cp-20, 0x1.749698c8d338dp-22
+ 0x1.967e0a63ca836p-14, -0x1.497b99d2a77d1p-15, 0x1.0ae346258cbdep-16,
+ -0x1.aade68fb2f076p-18, 0x1.22e609bf8671fp-19
};
const double __tanh_data_db[12][3] =
diff --git a/sysdeps/ieee754/dbl-64/e_tanh_data.h b/sysdeps/ieee754/dbl-64/e_tanh_data.h
index 4be829b46e4..6c632d26832 100644
--- a/sysdeps/ieee754/dbl-64/e_tanh_data.h
+++ b/sysdeps/ieee754/dbl-64/e_tanh_data.h
@@ -1,6 +1,6 @@
/* Correctly rounded hyperbolic tangent function for binary64 values.
-Copyright (c) 2023-2026 Alexei Sibidanov.
+Copyright (c) 2023-2026 Alexei Sibidanov, Cyprien Peignier, Paul Zimmermann.
This file is part of the CORE-MATH project
(https://core-math.gitlabpages.inria.fr/).
@@ -27,9 +27,9 @@ SOFTWARE.
#ifndef _E_TANH_DATA_H
#define _E_TANH_DATA_H
-extern const double __tanh_data_ch[][2] attribute_hidden;
+extern const double __tanh_data_ch[9][2] attribute_hidden;
#define CH __tanh_data_ch
-extern const double __tanh_data_cl[] attribute_hidden;
+extern const double __tanh_data_cl[5] attribute_hidden;
#define CL __tanh_data_cl
extern const double __tanh_data_db[12][3] attribute_hidden;
diff --git a/sysdeps/ieee754/dbl-64/s_tanh.c b/sysdeps/ieee754/dbl-64/s_tanh.c
index 32aba5997f4..2029de8fa5e 100644
--- a/sysdeps/ieee754/dbl-64/s_tanh.c
+++ b/sysdeps/ieee754/dbl-64/s_tanh.c
@@ -1,9 +1,13 @@
/* Correctly rounded hyperbolic tangent function for binary64 values.
-Copyright (c) 2023-2026 Alexei Sibidanov.
+Copyright (c) 2023-2026 Alexei Sibidanov, Cyprien Peignier, Paul Zimmermann.
+
+Alexei Sibidanov designed the original algorithm, while Cyprien Peignier and
+Paul Zimmermann extended the fma formula for |x0| <= 0x1.d12ed0af1a27fp-27,
+and improved the minimax polynomial for x0 <= |x| < 0.25.
The original version of this file was copied from the CORE-MATH
-project (file src/binary64/tanh/tanh.c, revision 80f5e8b0).
+project (file src/binary64/tanh/tanh.c, revision cf237fa0).
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
@@ -35,17 +39,22 @@ SOFTWARE.
static __attribute__((noinline)) double as_tanh_database(double, double);
+/* At input, *l is the approximation of the upper part of the polynomial
+ (evaluated with double arithmetic only). */
static inline double
polydd_tanh (double xh, double xl, int n, const double c[][2], double *l)
{
int i = n - 1;
- double ch = c[i][0] + *l, cl = ((c[i][0] - ch) + *l) + c[i][1];
+ double t;
+ double ch = fasttwosum (c[i][0], *l, &t);
+ double cl = t + c[i][1];
+ /* ch + cl ~= c[i][0] + c[i][1] + *l */
while (--i >= 0)
{
- ch = muldd_acc (xh, xl, ch, cl, &cl);
- double th = ch + c[i][0], tl = (c[i][0] - th) + ch;
- ch = th;
- cl += tl + c[i][1];
+ double tl;
+ ch = muldd_acc2 (xh, xl, ch, cl, &cl);
+ ch = fasttwosum (c[i][0], ch, &tl);
+ cl = (cl + c[i][1]) + tl;
}
*l = cl;
return ch;
@@ -69,9 +78,9 @@ as_exp_accurate (double x, double t, double th, double tl, double *l)
double fl = dxh
* (0x1.5555555555555p-5
+ dxh * (0x1.11111113e93e9p-7 + dxh * 0x1.6c16c169400a7p-10));
- double fh = polydd_tanh (dxh, dxl, 3, ch, &fl);
- fh = muldd_acc (dxh, dxl, fh, fl, &fl);
- fh = muldd_acc (th, tl, fh, fl, &fl);
+ double fh = polydd_tanh (dxh, dxl, array_length (ch), ch, &fl);
+ fh = muldd_acc2 (dxh, dxl, fh, fl, &fl);
+ fh = muldd_acc2 (th, tl, fh, fl, &fl);
double zh = th + fh, zl = (th - zh) + fh;
double uh = zh + tl, ul = ((zh - uh) + tl) + zl;
double vh = uh + fl, vl = ((uh - vh) + fl) + ul;
@@ -79,21 +88,16 @@ as_exp_accurate (double x, double t, double th, double tl, double *l)
return vh;
}
+/* Assumes |x| < 0.25. */
static double __attribute__ ((noinline))
as_tanh_zero (double x)
-{ // |x|<0.25
+{
double x2 = x * x, x2l = fma (x, x, -x2);
double y2
- = x2
- * (CL[0]
- + x2
- * (CL[1]
- + x2
- * (CL[2]
- + x2 * (CL[3] + x2 * (CL[4] + x2 * (CL[5]))))));
- double y1 = polydd_tanh (x2, x2l, 10, CH, &y2);
- y1 = mulddd (y1, y2, x, &y2);
- y1 = muldd_acc (y1, y2, x2, x2l, &y2);
+ = x2 * (CL[0] + x2 * (CL[1] + x2 * (CL[2] + x2 * (CL[3] + x2 * CL[4]))));
+ double y1 = polydd_tanh (x2, x2l, array_length (CH), CH, &y2);
+ y1 = mulddd_acc (y1, y2, x, &y2);
+ y1 = muldd_acc2 (y1, y2, x2, x2l, &y2);
double y0 = fasttwosum (x, y1, &y1);
y1 = fasttwosum (y1, y2, &y2);
uint64_t t = asuint64 (y1);
@@ -105,7 +109,7 @@ as_tanh_zero (double x)
else
t++;
y1 = asdouble (t);
- if (__glibc_unlikely (y2 == 0x0p+0))
+ if (__glibc_unlikely (fabs (x) == 0x1.ac343b179fec4p-3))
return as_tanh_database (x, y0 + y1);
}
return y0 + y1;
@@ -148,7 +152,7 @@ __tanh (double x)
if (__glibc_unlikely (aix >= UINT64_C(0x40330fc1931f09ca)))
{
if (aix > UINT64_C(0x7ff0000000000000))
- return x + x; // nan
+ return x + x; /* nan */
double f = copysign (1.0, x);
if (aix == UINT64_C(0x7ff0000000000000))
return f;
@@ -164,34 +168,28 @@ __tanh (double x)
double t = asdouble (v) - 0x1.8p25;
int64_t i1 = (jt >> 27) & 0x3f, i0 = (jt >> 33) & 0x3f,
ie = (int64_t) (jt << 13) >> MANTISSA_WIDTH;
- const double sp = asdouble ((1023 + ie) << MANTISSA_WIDTH);
+ const double sp = asdouble ((uint64_t) (1023 + ie) << MANTISSA_WIDTH);
static const double ch[]
= { 0x1p+1, 0x1p+1, 0x1.55555557e54ffp+0, 0x1.55555553a12f4p-1 };
double t0h = T0[i0][1], t1h = T1[i1][1], th = t0h * t1h, tl;
if (aix < UINT64_C(0x400d76c8b4395810))
- { // |x| ~< 3.683
+ { /* |x| ~< 3.683 */
if (__glibc_unlikely (aix < UINT64_C(0x3fd0000000000000)))
- { // |x| < 0x1p-2
- if (__glibc_unlikely (aix < UINT64_C(0x3e10000000000000)))
- { // |x| < 0x1p-30
- if (__glibc_unlikely (aix < UINT64_C(0x3df0000000000000)))
- { // |x| < 0x1p-32
- if (__glibc_unlikely (!aix))
- return x;
- /* We have underflow when 0 < |x| < 2^-1022 or when |x| =
- 2^-1022 and rounding towards zero. */
- double res = fma (x, -0x1p-55, x);
- return res;
- }
- double x3 = x * x * x;
- return x - x3 / 3;
+ { /* |x| < 0x1p-2 */
+ if (__glibc_unlikely (aix <= UINT64_C(0x3e4d12ed0af1a27f)))
+ { /* |x| <= 0x1.d12ed0af1a27fp-27 */
+ if (__glibc_unlikely (!aix))
+ return x;
+ /* We have underflow when 0 < |x| < 2^-1022 or when |x| =
+ 2^-1022 and rounding towards zero. */
+ return fma (x, -0x1p-55, x);
}
static const double c[] =
{
- -0x1.5555555555554p-2, 0x1.1111111110d61p-3,
- -0x1.ba1ba1b983d8bp-5, 0x1.664f4820e99fp-6,
- -0x1.226e11e4ac7cfp-7, 0x1.d6c4ab70668b6p-9,
- -0x1.7bbecb57ce996p-10, 0x1.1451443697dd8p-11
+ -0x1.5555555555555p-2, 0x1.1111111110f33p-3,
+ -0x1.ba1ba1b9b8ea6p-5, 0x1.664f4838e0a43p-6,
+ -0x1.226e17d1bc09bp-7, 0x1.d6c64dfba2565p-9,
+ -0x1.7bdd094d327afp-10, 0x1.1535ad0c31d0ep-11
};
double x2 = x * x, x3 = x2 * x, x4 = x2 * x2, x8 = x4 * x4;
double p1 = (c[4] + x2 * c[5]) + x4 * (c[6] + x2 * c[7]);
@@ -199,6 +197,16 @@ __tanh (double x)
p0 += x8 * p1;
p0 *= x3;
double rl, rh = fasttwosum (x, p0, &rl);
+ /* The branch 0x1.d12ed0af1a27fp-27 <= x < 0x1p-26 was checked
+ exhaustively (with and without fma contraction) with revision
+ 1820535, with the error bound e = x3*0x1.4dp-52. It fails with
+ 0x1.4cp-52 and x=0x1.27a0e7f47f0fap-4 (rndz, no fma contraction).
+ The interval [0x1p-3, 0x1.00cp-3] was checked exhaustively with
+ rndz and without fma contraction, with error bound
+ e = x3*0x1.80p-52: no failure. The interval
+ [0x1.015891c9eaef8p-3, 0x1.019891c9eaef8p-3] was checked
+ exhaustively with rndz and without fma contraction, with error
+ bound e = x3*0x1.80p-52: no failure. */
double e = x3 * 0x1.ap-52, lb = rh + (rl - e), ub = rh + (rl + e);
if (lb == ub)
return lb;
@@ -221,9 +229,13 @@ __tanh (double x)
ql += qd;
double rqh = 1 / qh, rql = (ql * rqh + fma (rqh, qh, -1)) * -rqh;
- ph = muldd_acc (ph, pl, rqh, rql, &pl);
+ ph = muldd_acc2 (ph, pl, rqh, rql, &pl);
- double e = rh * 0x1p-62;
+ /* This branch was tested exhaustively with/without fma contraction.
+ During this search, a failure was found with the original error
+ bound (e = rh*0x1p-62) and x=0x1.a0112a16e9318p+1 (rndu, no fma
+ contraction). */
+ double e = rh * 0x1.0bp-62;
rh = fasttwosub (0.5, ph, &rl);
rl -= pl;
rh *= copysign (2, x);
@@ -233,12 +245,17 @@ __tanh (double x)
return lb;
}
else
- {
+ { /* 3.683 ~< |x| < 0x1.30fc1931f09cap+4 */
static const double l2 = -0x1.62e42fefa39efp-14;
double dx = fma (l2, t, -ax), dx2 = dx * dx;
double p = dx * ((ch[0] + dx * ch[1]) + dx2 * (ch[2] + dx * ch[3]));
double rh = th * sp;
rh += (p + ((2 * 0x1.3p-55) * ax)) * rh;
+ /* This branch was tested exhaustively with/without fma contraction.
+ During this search, the largest 9-bit value of e for which it fails
+ was found to be e = rh*0x1.fap-50 with x=0x1.09cc2de69e78cp+2
+ (rndu, with/without fma contraction). Thus the bound below can be
+ reduced to rh*0x1.fbp-50. */
double e = rh * 0x1.1p-49;
rh = (2 * rh) / (1 + rh);
double one = copysign (1, x);
@@ -256,7 +273,7 @@ __tanh (double x)
double qd, qh = fasttwosum (1, rh, &qd), ql = rl + qd;
qh = fasttwosum (qh, ql, &ql);
double rqh = 1 / qh, rql = (ql * rqh + fma (rqh, qh, -1)) * -rqh;
- double pl, ph = muldd_acc (rh, rl, rqh, rql, &pl);
+ double pl, ph = muldd_acc2 (rh, rl, rqh, rql, &pl);
rh = fasttwosub (0.5, ph, &rl);
rl -= pl;
rh = fasttwosum (rh, rl, &rl);
--
2.43.0