[glibc] math: Update tanh from CORE-MATH (revision cf237fa0)
Adhemerval Zanella via Glibc-cvs <[email protected]> Sat, 18 Jul 2026 16:24:38 +0000 (GMT)
| Newsgroups | gmane.comp.lib.glibc.cvs |
|---|---|
| Message-ID | <[email protected]> |
https://sourceware.org/git/gitweb.cgi?p=glibc.git;h=6b29e74c731e2ac1bc8048a389fdaf029bd96072 commit 6b29e74c731e2ac1bc8048a389fdaf029bd96072 Author: Adhemerval Zanella <[email protected]> Date: Fri Jul 17 17:14:44 2026 +0000 math: Update tanh from CORE-MATH (revision cf237fa0) Sync the double-precision tanh implementation with CORE-MATH revision cf237fa0. Checked on aarch64-linux-gnu. Diff: --- 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 3ab7d41d98..cf2a19c015 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 73737235de..ca670768a2 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 04ecb56846..2f2a9166c0 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 6a31a005c3..28a910592e 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 0ebd8277df..5429a866a0 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 4be829b46e..6c632d2683 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 32aba5997f..2029de8fa5 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);