mirror of
git://sourceware.org/git/glibc.git
synced 2026-09-08 23:58:31 +08:00
math: Use sinh from CORE-MATH
The current implementation precision shows the following accuracy, on
three ranges ([-DBL_MAX,-10], [-10,10], [10,DBL_MAX]) with 10e9 uniform
randomly generated numbers for each range (first column is the
accuracy in ULP, with '0' being correctly rounded, second is the
number of samples with the corresponding precision):
* Range [-DBL_MAX, -10]
* FE_TONEAREST
0: 10000000000 100.00%
* FE_UPWARD
0: 10000000000 100.00%
* FE_DOWNWARD
0: 10000000000 100.00%
* FE_TOWARDZERO
0: 10000000000 100.00%
* Range [-10, -10]
* FE_TONEAREST
0: 3169388892 73.79%
1: 1125270674 26.20%
2: 307729 0.01%
* FE_UPWARD
0: 1450068660 33.76%
1: 2146926394 49.99%
2: 697404986 16.24%
3: 567255 0.01%
* FE_DOWNWARD
0: 1449727976 33.75%
1: 2146957381 49.99%
2: 697719649 16.25%
3: 562289 0.01%
* FE_TOWARDZERO
0: 2519351889 58.66%
1: 1773434502 41.29%
2: 2180904 0.05%
* Range [10, DBL_MAX]
* FE_TONEAREST
0: 10000000000 100.00%
* FE_UPWARD
0: 10000000000 100.00%
* FE_DOWNWARD
0: 10000000000 100.00%
* FE_TOWARDZERO
0: 10000000000 100.00%
The CORE-MATH implementation is correctly rounded for any rounding mode.
The code was adapted to glibc style and to use the definition of
math_config.h (to handle errno, overflow, and underflow).
Performance-wise, it shows:
latency master patched improvement
x86_64 101.0710 129.4710 -28.10%
x86_64v2 101.1810 127.6370 -26.15%
x86_64v3 96.0685 48.5911 49.42%
aarch64 41.4229 22.3971 45.93%
armhf-vpfv4 42.8620 25.6011 40.27%
powerpc64le 29.2630 13.1450 55.08%
reciprocal-throughput master patched improvement
x86_64 42.6895 105.7150 -147.64%
x86_64v2 42.7255 104.7480 -145.17%
x86_64v3 39.6949 25.9087 34.73%
aarch64 26.0104 19.2236 26.09%
armhf-vpfv4 29.4362 23.6350 19.71%
powerpc64le 12.9170 8.34582 35.39%
* x86_64: gcc version 15.2.1 20260112, Ryzen 9 5900X, --disable-multi-arch
* aarch64: gcc version 15.2.1 20251105, Neoverse-N1
* armv7a-vpfv4: gcc version 15.2.1 20251105, Neoverse-N1
* powerpc64le: gcc version 15.2.1 20260128, POWER10
Checked on x86_64-linux-gnu, aarch64-linux-gnu, and
powerpc64le-linux-gnu.
Reviewed-by: DJ Delorie <dj@redhat.com>
This commit is contained in:
Adhemerval Zanella
@@ -10,7 +10,7 @@ Version 2.44
|
||||
Major new features:
|
||||
|
||||
* Additional optimized and correctly rounded mathematical functions have
|
||||
been imported from the CORE-MATH project, in particular cosh.
|
||||
been imported from the CORE-MATH project, in particular cosh, and sinh.
|
||||
|
||||
* The SVID handling for cosh was moved to compat symbols, allowing
|
||||
improvements in performance.
|
||||
|
||||
@@ -254,6 +254,8 @@ core-math:
|
||||
sysdeps/ieee754/dbl-64/e_gamma_r.c
|
||||
# src/binary64/lgamma/lgamma.c, revision 0413bb7e
|
||||
sysdeps/ieee754/dbl-64/e_lgamma_r.c
|
||||
# src/binary64/sinh/sinh.c, revision 8127b7ac
|
||||
sysdeps/ieee754/dbl-64/e_sinh.c
|
||||
# src/binary64/asinh/asinh.c, revision cd653cf7
|
||||
sysdeps/ieee754/dbl-64/s_asinh.c
|
||||
# src/binary64/erf/erf.c, revision 384ed01d
|
||||
|
||||
@@ -9574,6 +9574,45 @@ sinh -0x6.9bbb6df7c5d08p-4
|
||||
sinh 0x1.250bfep-11
|
||||
# the next value generates larger error bounds on x86_64 (ldbl-96)
|
||||
sinh 0x2.c5d376167f4052f4p+12
|
||||
sinh -0x1.4a3f2b7dff682p+2
|
||||
sinh -0x1.41267a2ed2593p+2
|
||||
sinh 0x1.41c456749a084p+2
|
||||
sinh 0x1.41feabdae88c6p+2
|
||||
sinh -0x1.41634f0c2ae4dp+2
|
||||
sinh -0x1.4104fac728038p+2
|
||||
sinh -0x1.419af27adcfe8p+2
|
||||
sinh -0x1.41265ff261702p+2
|
||||
sinh 0x1.424c64d09a2d4p+2
|
||||
sinh 0x1.409364da45b54p+2
|
||||
sinh -0x1.409359e03b8fap+2
|
||||
sinh 0x1.4000685a224bcp+2
|
||||
sinh -0x1.49d5d039e0535p+2
|
||||
sinh -0x1.41e005cc32758p+2
|
||||
sinh -0x1.497208f2ae258p+2
|
||||
sinh -0x1.4037cd9cb7a82p+2
|
||||
sinh -0x1.41b923f7769e3p+2
|
||||
sinh -0x1.425754461bdfdp+2
|
||||
sinh -0x1.425761708c475p+2
|
||||
sinh 0x1.40e3ea0aeab5ap+2
|
||||
sinh 0x1.4144bd89e0eb8p+2
|
||||
sinh 0x1.403d4a3c61f18p+2
|
||||
sinh 0x1.41d7bcbdf9fccp+2
|
||||
sinh 0x1.42575956957ap+2
|
||||
sinh -0x1.09e94dfa477cp-2
|
||||
sinh -0x1.40a6d6ee4ba13p+2
|
||||
sinh -0x1.42575380f620bp+2
|
||||
sinh -0x1.41342809b65ffp+2
|
||||
sinh 0x1.420149be82e16p+2
|
||||
sinh -0x1.417c56d39aac4p+2
|
||||
sinh 0x1.41bc1dc9e7ebp+2
|
||||
sinh 0x1.410a9f2b0fd78p+2
|
||||
sinh 0x1.41d7ae912e7e2p+2
|
||||
sinh 0x1.403d745ee03bp+2
|
||||
sinh 0x1.42360842659f6p+2
|
||||
sinh 0x1.424c2ba022b28p+2
|
||||
sinh 0x1.78b7f7218cbbap-8
|
||||
sinh 0x1.71a7baa9570bdp-26
|
||||
sinh 0x1.0042fffffffffp-3
|
||||
sinh max
|
||||
sinh -max
|
||||
sinh min
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -11,6 +11,7 @@ ifeq ($(subdir),math)
|
||||
CFLAGS-e_cosh.c += -fexcess-precision=standard
|
||||
CFLAGS-e_lgamma_r.c += -fexcess-precision=standard
|
||||
CFLAGS-e_gamma_r.c += -fexcess-precision=standard
|
||||
CFLAGS-e_sinh.c += -fexcess-precision=standard
|
||||
CFLAGS-s_erf.c += -fexcess-precision=standard
|
||||
CFLAGS-s_erfc.c += -fexcess-precision=standard
|
||||
CFLAGS-s_erf_common.c += -fexcess-precision=standard
|
||||
|
||||
+475
-79
@@ -1,99 +1,495 @@
|
||||
/* @(#)e_sinh.c 5.1 93/09/24 */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
/* Correctly rounded hyperbolic sine for binary64 values.
|
||||
|
||||
#if defined(LIBM_SCCS) && !defined(lint)
|
||||
static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $";
|
||||
#endif
|
||||
Copyright (c) 2023-2026 Alexei Sibidanov.
|
||||
|
||||
/* __ieee754_sinh(x)
|
||||
* Method :
|
||||
* mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
|
||||
* 1. Replace x by |x| (sinh(-x) = -sinh(x)).
|
||||
* 2.
|
||||
* E + E/(E+1)
|
||||
* 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x)
|
||||
* 2
|
||||
*
|
||||
* 22 <= x <= lnovft : sinh(x) := exp(x)/2
|
||||
* lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2)
|
||||
* ln2ovft < x : sinh(x) := x*shuge (overflow)
|
||||
*
|
||||
* Special cases:
|
||||
* sinh(x) is |x| if x is +INF, -INF, or NaN.
|
||||
* only sinh(0)=0 is exact for finite x.
|
||||
*/
|
||||
The original version of this file was copied from the CORE-MATH
|
||||
project (file src/binary64/sinh/sinh.c, revision 8c2c3473).
|
||||
|
||||
#include <float.h>
|
||||
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.
|
||||
|
||||
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.
|
||||
*/
|
||||
|
||||
#include <array_length.h>
|
||||
#include <math.h>
|
||||
#include <math-narrow-eval.h>
|
||||
#include <math_private.h>
|
||||
#include <math-underflow.h>
|
||||
#include <stdint.h>
|
||||
#include <errno.h>
|
||||
#include <libm-alias-finite.h>
|
||||
|
||||
static const double one = 1.0, shuge = 1.0e307;
|
||||
#include <ddcoremath.h>
|
||||
#include "math_config.h"
|
||||
|
||||
#ifndef SECTION
|
||||
# define SECTION
|
||||
#endif
|
||||
|
||||
static inline double
|
||||
polydd_sinh (double xh, double xl, int n, const double c[][2], double *l)
|
||||
{
|
||||
int i = n - 1;
|
||||
double ch, cl, e;
|
||||
ch = fasttwosum (c[i][0], *l, &cl);
|
||||
cl += c[i][1];
|
||||
while (--i >= 0)
|
||||
{
|
||||
ch = muldd2 (xh, xl, ch, cl, &cl);
|
||||
ch = fasttwosum (c[i][0], ch, &e);
|
||||
cl = (cl + c[i][1]) + e;
|
||||
}
|
||||
*l = cl;
|
||||
return ch;
|
||||
}
|
||||
|
||||
static double __attribute__ ((noinline))
|
||||
as_exp_accurate (double x, double t, double th, double tl, double *l)
|
||||
{
|
||||
static const double ch[][2] =
|
||||
{
|
||||
{ 0x1p+0, 0x1.6c16bd194535dp-94 },
|
||||
{ 0x1p-1, -0x1.8259d904fd34fp-93 },
|
||||
{ 0x1.5555555555555p-3, 0x1.53e93e9f26e62p-57 }
|
||||
};
|
||||
const double l2h = 0x1.62e42ffp-13, l2l = 0x1.718432a1b0e26p-47,
|
||||
l2ll = 0x1.9ff0342542fc3p-102;
|
||||
double dx = x - l2h * t, dxl = l2l * t,
|
||||
dxll = l2ll * t + fma (l2l, t, -dxl);
|
||||
double dxh = dx + dxl;
|
||||
dxl = ((dx - dxh) + dxl) + dxll;
|
||||
double fl = dxh
|
||||
* (0x1.5555555555555p-5
|
||||
+ dxh * (0x1.11111113e93e9p-7 + dxh * 0x1.6c16c169400a7p-10));
|
||||
double fh = polydd_sinh (dxh, dxl, 3, ch, &fl);
|
||||
fh = muldd2 (dxh, dxl, fh, fl, &fl);
|
||||
fh = muldd2 (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;
|
||||
*l = vl;
|
||||
return vh;
|
||||
}
|
||||
|
||||
static double __attribute__ ((noinline))
|
||||
as_sinh_zero (double x)
|
||||
{
|
||||
static const double ch[][2] =
|
||||
{
|
||||
{ 0x1.5555555555555p-3, 0x1.555555555552fp-57 },
|
||||
{ 0x1.1111111111111p-7, 0x1.11111115cf00dp-63 },
|
||||
{ 0x1.a01a01a01a01ap-13, 0x1.a0011c925b85cp-73 },
|
||||
{ 0x1.71de3a556c734p-19, -0x1.b4e2835532bcdp-73 },
|
||||
{ 0x1.ae64567f54482p-26, -0x1.defcf17a6ab79p-81 }
|
||||
};
|
||||
double x2 = x * x, x2l = fma (x, x, -x2);
|
||||
double y2 = x2
|
||||
* (0x1.6124613aef206p-33
|
||||
+ x2 * (0x1.ae7f36beea815p-41 + x2 * 0x1.95785063cd974p-49));
|
||||
double y1 = polydd_sinh (x2, x2l, 5, ch, &y2);
|
||||
y1 = mulddd3 (y1, y2, x, &y2);
|
||||
y1 = muldd2 (y1, y2, x2, x2l, &y2);
|
||||
double y0 = fasttwosum (x, y1, &y1);
|
||||
y1 = fasttwosum (y1, y2, &y2);
|
||||
uint64_t t = asuint64 (y1);
|
||||
if (__glibc_unlikely (!(t & MANTISSA_MASK)))
|
||||
{
|
||||
uint64_t w = asuint64 (y2);
|
||||
if ((w ^ t) >> 63)
|
||||
t--;
|
||||
else
|
||||
t++;
|
||||
y1 = asdouble (t);
|
||||
}
|
||||
return y0 + y1;
|
||||
}
|
||||
|
||||
static __attribute__ ((noinline)) double
|
||||
as_sinh_database (double x, double f)
|
||||
{
|
||||
static const double db[][3] =
|
||||
{
|
||||
{ 0x1.364303e1ad8f6p-2, 0x1.3b07e0c779ddap-2, -0x1.bcp-106 },
|
||||
{ 0x1.4169f234f23b9p-2, 0x1.46b7b3b358f99p-2, -0x1p-56 },
|
||||
{ 0x1.616cc75d49226p-2, 0x1.687bd068c1c1ep-2, 0x1.ap-111 },
|
||||
{ 0x1.ae3773250e7d2p-2, 0x1.bafc3479fc9ccp-2, -0x1p-105 },
|
||||
{ 0x1.b7efa91915c95p-2, 0x1.c59869f17b483p-2, -0x1p-104 },
|
||||
{ 0x1.d68039861ab53p-2, 0x1.e73b46abb01e1p-2, -0x1.2p-109 },
|
||||
{ 0x1.e90f16eb88c09p-2, 0x1.fbdd4a37760b7p-2, -0x1.f8p-108 },
|
||||
{ 0x1.a3fc7e4dd47d1p-1, 0x1.d4b21ebf542fp-1, 0x1.ep-107 },
|
||||
{ 0x1.aa3b649a96091p-1, 0x1.dd32c5ed1e93p-1, 0x1.8ap-106 },
|
||||
{ 0x1.c13876341b62ep-1, 0x1.fd1d7f1c8170cp-1, 0x1.72p-105 },
|
||||
{ 0x1.2f5d3b178914ap+0, 0x1.7b8516ffd2406p+0, -0x1.28p-104 },
|
||||
{ 0x1.3ffc12b81cbc2p+0, 0x1.9a0ff413a1af3p+0, 0x1.cp-107 },
|
||||
{ 0x1.44f65dff00782p+0, 0x1.a38a3c3227609p+0, -0x1p-103 },
|
||||
{ 0x1.7346e3c591a14p+0, 0x1.01e9cfa77b855p+1, 0x1.p-102 },
|
||||
{ 0x1.b6e2c73f41415p+0, 0x1.57e377b3f0b4bp+1, -0x1p-102 },
|
||||
{ 0x1.dc5059d4e507dp+0, 0x1.9168c60ed5256p+1, 0x1.c6p-104 },
|
||||
{ 0x1.f737f1e8378c7p+0, 0x1.bffd3f94f40fbp+1, 0x1.a4p-104 },
|
||||
{ 0x1.3359640329982p+1, 0x1.5e40df3f985bep+2, 0x1.97p-102 },
|
||||
{ 0x1.58a4ff5adac35p+1, 0x1.d671928665bddp+2, 0x1p-102 },
|
||||
{ 0x1.8c0a26d055288p+1, 0x1.6056b06a21918p+3, -0x1.bep-102 },
|
||||
{ 0x1.bc3c2d0c95f52p+1, 0x1.00fef7383a978p+4, 0x1.61p-100 },
|
||||
{ 0x1.0a19aebb51e9p+3, 0x1.fee8f69c4cd25p+10, 0x1.48p-95 },
|
||||
{ 0x1.3eb8f61734227p+3, 0x1.4ab1cf45e4e26p+13, 0x1p-90 },
|
||||
{ 0x1.43a81752eabe7p+3, 0x1.81d364845ecfap+13, -0x1p-90 },
|
||||
{ 0x1.16369cd53bb69p+4, 0x1.0fbc6c02b1c9p+24, -0x1.9p-81 },
|
||||
{ 0x1.20e29ea8b51e2p+4, 0x1.08b8abba28abcp+25, 0x1.9bp-79 },
|
||||
{ 0x1.92a5c27afbe82p+4, 0x1.3c81f9a247253p+35, 0x1p-67 },
|
||||
{ 0x1.a1e4f11b513d7p+4, 0x1.9a65b6c2e2185p+36, -0x1.bcp-70 },
|
||||
{ 0x1.c089fcf166171p+4, 0x1.5c452e0e37569p+39, 0x1.4p-69 },
|
||||
{ 0x1.e42a98b3a0be5p+4, 0x1.938768ca4f8aap+42, 0x1.6dp-62 },
|
||||
{ 0x1.04db52248cbb8p+5, 0x1.0794072349523p+46, 0x1.0e8p-57 },
|
||||
{ 0x1.21bc021eeb97ep+5, 0x1.3065064a170fbp+51, 0x1.088p-52 },
|
||||
{ 0x1.39fc4d3bb711p+5, 0x1.8a4e90733b95ep+55, 0x1.6ep-50 },
|
||||
{ 0x1.3c895d86e96c9p+5, 0x1.0f33837882a6p+56, -0x1.28p-49 },
|
||||
{ 0x1.e07e71bfcf06fp+5, 0x1.91ec4412c344fp+85, 0x1p-24 },
|
||||
{ 0x1.f7216c4b435c9p+5, 0x1.a97e7be23e65ap+89, -0x1p-15 },
|
||||
{ 0x1.6474c604cc0d7p+6, 0x1.7a8f65ad009bdp+127, -0x1.08p+20 },
|
||||
{ 0x1.7a60ee15e3e9dp+6, 0x1.62e4dc3bbf53fp+135, 0x1.bp+29 },
|
||||
{ 0x1.1f0da93354198p+7, 0x1.0bd73b73fc74cp+206, 0x1.59p+102 },
|
||||
{ 0x1.54cd1fea7663ap+7, 0x1.c90810d354618p+244, 0x1.2p+135 },
|
||||
{ 0x1.556c678d5e976p+7, 0x1.37e7ac4e7f9b3p+245, 0x1.02p+141 },
|
||||
{ 0x1.7945e34b18a9ap+7, 0x1.1b0e4936a8c9bp+271, -0x1.fap+166 },
|
||||
{ 0x1.2da9e5e6af0bp+8, 0x1.27d6fe867d6f6p+434, 0x1.0ap+329 },
|
||||
{ 0x1.54ceba01331d5p+8, 0x1.9a86785b5ef3ep+490, -0x1.22p+386 },
|
||||
{ 0x1.9e7b643238a14p+8, 0x1.f5da7fe652978p+596, 0x1p+493 },
|
||||
{ 0x1.c7206c1b753e4p+8, 0x1.8670de0b68cadp+655, -0x1.78p+548 },
|
||||
{ 0x1.d6479eba7c971p+8, 0x1.62a88613629b6p+677, -0x1.4p+568 },
|
||||
{ 0x1.eb9914d4ac1c8p+8, 0x1.2b67eff65dce8p+708, -0x1.02p+603 },
|
||||
{ 0x1.0bc04af1b09f5p+9, 0x1.7b1d97c902985p+771, 0x1.56p+666 },
|
||||
{ 0x1.26ee1a46d8c8bp+9, 0x1.fbe20477df4a7p+849, -0x1.55p+745 },
|
||||
{ 0x1.4a869881f72acp+9, 0x1.9ea7540a3d1f9p+952, -0x1.2dp+848 },
|
||||
};
|
||||
int a = 0, b = array_length (db) - 1, m = (a + b) / 2;
|
||||
double ax = fabs (x);
|
||||
while (a <= b)
|
||||
{
|
||||
if (db[m][0] < ax)
|
||||
a = m + 1;
|
||||
else if (db[m][0] == ax)
|
||||
{
|
||||
f = copysign (1, x) * db[m][1] + copysign (1, x) * db[m][2];
|
||||
break;
|
||||
}
|
||||
else
|
||||
b = m - 1;
|
||||
m = (a + b) / 2;
|
||||
}
|
||||
return f;
|
||||
}
|
||||
|
||||
SECTION
|
||||
double
|
||||
__ieee754_sinh (double x)
|
||||
{
|
||||
double t, w, h;
|
||||
int32_t ix, jx;
|
||||
uint32_t lx;
|
||||
|
||||
/* High word of |x|. */
|
||||
GET_HIGH_WORD (jx, x);
|
||||
ix = jx & 0x7fffffff;
|
||||
|
||||
/* x is INF or NaN */
|
||||
if (__glibc_unlikely (ix >= 0x7ff00000))
|
||||
return x + x;
|
||||
|
||||
h = 0.5;
|
||||
if (jx < 0)
|
||||
h = -h;
|
||||
/* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
|
||||
if (ix < 0x40360000) /* |x|<22 */
|
||||
/*
|
||||
The function sinh(x) is approximated by a minimax polynomial for
|
||||
|x|<0.25. For other arguments the identity
|
||||
sinh(x)=(exp(|x|)-exp(-|x|))/2*copysign(1,x) is used. For |x|<5
|
||||
both exponents are calculated with slightly higher precision than
|
||||
double. For 5<|x|<36.736801 the exp(-|x|) is small and is
|
||||
calculated with double precision but exp(|x|) is calculated with
|
||||
higher than double precision. For 36.736801<|x|<710.47586
|
||||
exp(-|x|) becomes too small and only exp(|x|) is calculated.
|
||||
*/
|
||||
static const double t0[][2] =
|
||||
{
|
||||
if (__glibc_unlikely (ix < 0x3e300000)) { /* |x|<2**-28 */
|
||||
math_check_force_underflow (x);
|
||||
if (shuge + x > one)
|
||||
return x;
|
||||
/* sinh(tiny) = tiny with inexact */
|
||||
}
|
||||
t = __expm1 (fabs (x));
|
||||
if (ix < 0x3ff00000)
|
||||
return h * (2.0 * t - t * t / (t + one));
|
||||
return h * (t + t / (t + one));
|
||||
}
|
||||
|
||||
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
|
||||
if (ix < 0x40862e42)
|
||||
return h * __ieee754_exp (fabs (x));
|
||||
|
||||
/* |x| in [log(maxdouble), overflowthresold] */
|
||||
GET_LOW_WORD (lx, x);
|
||||
if (ix < 0x408633ce || ((ix == 0x408633ce) && (lx <= (uint32_t) 0x8fb9f87d)))
|
||||
{ 0x0p+0, 0x1p+0 },
|
||||
{ -0x1.19083535b085ep-56, 0x1.02c9a3e778061p+0 },
|
||||
{ 0x1.d73e2a475b466p-55, 0x1.059b0d3158574p+0 },
|
||||
{ 0x1.186be4bb285p-57, 0x1.0874518759bc8p+0 },
|
||||
{ 0x1.8a62e4adc610ap-54, 0x1.0b5586cf9890fp+0 },
|
||||
{ 0x1.03a1727c57b52p-59, 0x1.0e3ec32d3d1a2p+0 },
|
||||
{ -0x1.6c51039449b3ap-54, 0x1.11301d0125b51p+0 },
|
||||
{ -0x1.32fbf9af1369ep-54, 0x1.1429aaea92dep+0 },
|
||||
{ -0x1.19041b9d78a76p-55, 0x1.172b83c7d517bp+0 },
|
||||
{ 0x1.e5b4c7b4968e4p-55, 0x1.1a35beb6fcb75p+0 },
|
||||
{ 0x1.e016e00a2643cp-54, 0x1.1d4873168b9aap+0 },
|
||||
{ 0x1.dc775814a8494p-55, 0x1.2063b88628cd6p+0 },
|
||||
{ 0x1.9b07eb6c70572p-54, 0x1.2387a6e756238p+0 },
|
||||
{ 0x1.2bd339940e9dap-55, 0x1.26b4565e27cddp+0 },
|
||||
{ 0x1.612e8afad1256p-55, 0x1.29e9df51fdee1p+0 },
|
||||
{ 0x1.0024754db41d4p-54, 0x1.2d285a6e4030bp+0 },
|
||||
{ 0x1.6f46ad23182e4p-55, 0x1.306fe0a31b715p+0 },
|
||||
{ 0x1.32721843659a6p-54, 0x1.33c08b26416ffp+0 },
|
||||
{ -0x1.63aeabf42eae2p-54, 0x1.371a7373aa9cbp+0 },
|
||||
{ -0x1.5e436d661f5e2p-56, 0x1.3a7db34e59ff7p+0 },
|
||||
{ 0x1.ada0911f09ebcp-55, 0x1.3dea64c123422p+0 },
|
||||
{ -0x1.ef3691c309278p-58, 0x1.4160a21f72e2ap+0 },
|
||||
{ 0x1.89b7a04ef80dp-59, 0x1.44e086061892dp+0 },
|
||||
{ 0x1.3c1a3b69062fp-56, 0x1.486a2b5c13cdp+0 },
|
||||
{ 0x1.d4397afec42e2p-56, 0x1.4bfdad5362a27p+0 },
|
||||
{ -0x1.4b309d25957e4p-54, 0x1.4f9b2769d2ca7p+0 },
|
||||
{ -0x1.07abe1db13cacp-55, 0x1.5342b569d4f82p+0 },
|
||||
{ 0x1.9bb2c011d93acp-54, 0x1.56f4736b527dap+0 },
|
||||
{ 0x1.6324c054647acp-54, 0x1.5ab07dd485429p+0 },
|
||||
{ 0x1.ba6f93080e65ep-54, 0x1.5e76f15ad2148p+0 },
|
||||
{ -0x1.383c17e40b496p-54, 0x1.6247eb03a5585p+0 },
|
||||
{ -0x1.bb60987591c34p-54, 0x1.6623882552225p+0 },
|
||||
{ -0x1.bdd3413b26456p-54, 0x1.6a09e667f3bcdp+0 },
|
||||
{ -0x1.bbe3a683c88aap-57, 0x1.6dfb23c651a2fp+0 },
|
||||
{ -0x1.16e4786887a9ap-55, 0x1.71f75e8ec5f74p+0 },
|
||||
{ -0x1.0245957316dd4p-54, 0x1.75feb564267c9p+0 },
|
||||
{ -0x1.41577ee04993p-55, 0x1.7a11473eb0187p+0 },
|
||||
{ 0x1.05d02ba15797ep-56, 0x1.7e2f336cf4e62p+0 },
|
||||
{ -0x1.d4c1dd41532d8p-54, 0x1.82589994cce13p+0 },
|
||||
{ -0x1.fc6f89bd4f6bap-54, 0x1.868d99b4492edp+0 },
|
||||
{ 0x1.6e9f156864b26p-54, 0x1.8ace5422aa0dbp+0 },
|
||||
{ 0x1.5cc13a2e3976cp-55, 0x1.8f1ae99157736p+0 },
|
||||
{ -0x1.75fc781b57ebcp-57, 0x1.93737b0cdc5e5p+0 },
|
||||
{ -0x1.d185b7c1b85dp-54, 0x1.97d829fde4e5p+0 },
|
||||
{ 0x1.c7c46b071f2bep-56, 0x1.9c49182a3f09p+0 },
|
||||
{ -0x1.359495d1cd532p-54, 0x1.a0c667b5de565p+0 },
|
||||
{ -0x1.d2f6edb8d41e2p-54, 0x1.a5503b23e255dp+0 },
|
||||
{ 0x1.0fac90ef7fd32p-54, 0x1.a9e6b5579fdbfp+0 },
|
||||
{ 0x1.7a1cd345dcc82p-54, 0x1.ae89f995ad3adp+0 },
|
||||
{ -0x1.2805e3084d708p-57, 0x1.b33a2b84f15fbp+0 },
|
||||
{ -0x1.5584f7e54ac3ap-56, 0x1.b7f76f2fb5e47p+0 },
|
||||
{ 0x1.23dd07a2d9e84p-55, 0x1.bcc1e904bc1d2p+0 },
|
||||
{ 0x1.11065895048dep-55, 0x1.c199bdd85529cp+0 },
|
||||
{ 0x1.2884dff483cacp-54, 0x1.c67f12e57d14bp+0 },
|
||||
{ 0x1.503cbd1e949dcp-56, 0x1.cb720dcef9069p+0 },
|
||||
{ -0x1.cbc3743797a9cp-54, 0x1.d072d4a07897cp+0 },
|
||||
{ 0x1.2ed02d75b3706p-55, 0x1.d5818dcfba487p+0 },
|
||||
{ 0x1.c2300696db532p-54, 0x1.da9e603db3285p+0 },
|
||||
{ -0x1.1a5cd4f184b5cp-54, 0x1.dfc97337b9b5fp+0 },
|
||||
{ 0x1.39e8980a9cc9p-55, 0x1.e502ee78b3ff6p+0 },
|
||||
{ -0x1.e9c23179c2894p-54, 0x1.ea4afa2a490dap+0 },
|
||||
{ 0x1.dc7f486a4b6bp-54, 0x1.efa1bee615a27p+0 },
|
||||
{ 0x1.9d3e12dd8a18ap-54, 0x1.f50765b6e454p+0 },
|
||||
{ 0x1.74853f3a5931ep-55, 0x1.fa7c1819e90d8p+0 }
|
||||
};
|
||||
static const double t1[][2] =
|
||||
{
|
||||
w = __ieee754_exp (0.5 * fabs (x));
|
||||
t = h * w;
|
||||
return t * w;
|
||||
}
|
||||
{ 0x0p+0, 0x1p+0 },
|
||||
{ 0x1.ae8e38c59c72ap-54, 0x1.000b175effdc7p+0 },
|
||||
{ -0x1.7b5d0d58ea8f4p-58, 0x1.00162f3904052p+0 },
|
||||
{ 0x1.4115cb6b16a8ep-54, 0x1.0021478e11ce6p+0 },
|
||||
{ -0x1.d7c96f201bb2ep-55, 0x1.002c605e2e8cfp+0 },
|
||||
{ 0x1.84711d4c35eap-54, 0x1.003779a95f959p+0 },
|
||||
{ -0x1.0484245243778p-55, 0x1.0042936faa3d8p+0 },
|
||||
{ -0x1.4b237da2025fap-54, 0x1.004dadb113dap+0 },
|
||||
{ -0x1.5e00e62d6b30ep-56, 0x1.0058c86da1c0ap+0 },
|
||||
{ 0x1.a1d6cedbb948p-54, 0x1.0063e3a559473p+0 },
|
||||
{ -0x1.4acf197a00142p-54, 0x1.006eff583fc3dp+0 },
|
||||
{ -0x1.eaf2ea42391a6p-57, 0x1.007a1b865a8cap+0 },
|
||||
{ 0x1.da93f90835f76p-56, 0x1.0085382faef83p+0 },
|
||||
{ -0x1.6a79084ab093cp-55, 0x1.00905554425d4p+0 },
|
||||
{ 0x1.86364f8fbe8f8p-54, 0x1.009b72f41a12bp+0 },
|
||||
{ -0x1.82e8e14e3110ep-55, 0x1.00a6910f3b6fdp+0 },
|
||||
{ -0x1.4f6b2a7609f72p-55, 0x1.00b1afa5abcbfp+0 },
|
||||
{ -0x1.e1a258ea8f71ap-56, 0x1.00bcceb7707ecp+0 },
|
||||
{ 0x1.4362ca5bc26f2p-56, 0x1.00c7ee448ee02p+0 },
|
||||
{ 0x1.095a56c919d02p-54, 0x1.00d30e4d0c483p+0 },
|
||||
{ -0x1.406ac4e81a646p-57, 0x1.00de2ed0ee0f5p+0 },
|
||||
{ 0x1.b5a6902767e08p-54, 0x1.00e94fd0398ep+0 },
|
||||
{ -0x1.91b206085932p-54, 0x1.00f4714af41d3p+0 },
|
||||
{ 0x1.427068ab22306p-55, 0x1.00ff93412315cp+0 },
|
||||
{ 0x1.c1d0660524e08p-54, 0x1.010ab5b2cbd11p+0 },
|
||||
{ -0x1.e7bdfb3204be8p-54, 0x1.0115d89ff3a8bp+0 },
|
||||
{ 0x1.843aa8b9cbbc6p-55, 0x1.0120fc089ff63p+0 },
|
||||
{ -0x1.34104ee7edae8p-56, 0x1.012c1fecd613bp+0 },
|
||||
{ -0x1.2b6aeb6176892p-56, 0x1.0137444c9b5b5p+0 },
|
||||
{ 0x1.a8cd33b8a1bb2p-56, 0x1.01426927f5278p+0 },
|
||||
{ 0x1.2edc08e5da99ap-56, 0x1.014d8e7ee8d2fp+0 },
|
||||
{ 0x1.57ba2dc7e0c72p-55, 0x1.0158b4517bb88p+0 },
|
||||
{ 0x1.b61299ab8cdb8p-54, 0x1.0163da9fb3335p+0 },
|
||||
{ -0x1.90565902c5f44p-54, 0x1.016f0169949edp+0 },
|
||||
{ 0x1.70fc41c5c2d54p-55, 0x1.017a28af25567p+0 },
|
||||
{ 0x1.4b9a6e145d76cp-54, 0x1.018550706ab62p+0 },
|
||||
{ -0x1.008eff5142bfap-56, 0x1.019078ad6a19fp+0 },
|
||||
{ -0x1.77669f033c7dep-54, 0x1.019ba16628de2p+0 },
|
||||
{ -0x1.09bb78eeead0ap-54, 0x1.01a6ca9aac5f3p+0 },
|
||||
{ 0x1.371231477ece6p-54, 0x1.01b1f44af9f9ep+0 },
|
||||
{ 0x1.5e7626621eb5ap-56, 0x1.01bd1e77170b4p+0 },
|
||||
{ -0x1.bc72b100828a4p-54, 0x1.01c8491f08f08p+0 },
|
||||
{ -0x1.ce39cbbab8bbep-57, 0x1.01d37442d507p+0 },
|
||||
{ 0x1.16996709da2e2p-55, 0x1.01de9fe280ac8p+0 },
|
||||
{ -0x1.c11f5239bf536p-55, 0x1.01e9cbfe113efp+0 },
|
||||
{ 0x1.e1d4eb5edc6b4p-55, 0x1.01f4f8958c1c6p+0 },
|
||||
{ -0x1.afb99946ee3fp-54, 0x1.020025a8f6a35p+0 },
|
||||
{ -0x1.8f06d8a148a32p-54, 0x1.020b533856324p+0 },
|
||||
{ -0x1.2bf310fc54eb6p-55, 0x1.02168143b0281p+0 },
|
||||
{ -0x1.c95a035eb4176p-54, 0x1.0221afcb09e3ep+0 },
|
||||
{ -0x1.491793e46834cp-54, 0x1.022cdece68c4fp+0 },
|
||||
{ -0x1.3e8d0d9c4909p-56, 0x1.02380e4dd22adp+0 },
|
||||
{ -0x1.314aa16278aa4p-54, 0x1.02433e494b755p+0 },
|
||||
{ 0x1.48daf888e965p-55, 0x1.024e6ec0da046p+0 },
|
||||
{ 0x1.56dc8046821f4p-55, 0x1.02599fb483385p+0 },
|
||||
{ 0x1.45b42356b9d46p-54, 0x1.0264d1244c719p+0 },
|
||||
{ -0x1.082ef51b61d7ep-56, 0x1.027003103b10ep+0 },
|
||||
{ 0x1.2106ed0920a34p-56, 0x1.027b357854772p+0 },
|
||||
{ -0x1.fd4cf26ea5d0ep-54, 0x1.0286685c9e059p+0 },
|
||||
{ -0x1.09f8775e78084p-54, 0x1.02919bbd1d1d8p+0 },
|
||||
{ 0x1.64cbba902ca28p-58, 0x1.029ccf99d720ap+0 },
|
||||
{ 0x1.4383ef231d206p-54, 0x1.02a803f2d170dp+0 },
|
||||
{ 0x1.4a47a505b3a46p-54, 0x1.02b338c811703p+0 },
|
||||
{ 0x1.e47120223468p-54, 0x1.02be6e199c811p+0 },
|
||||
};
|
||||
|
||||
/* |x| > overflowthresold, sinh(x) overflow */
|
||||
return math_narrow_eval (x * shuge);
|
||||
const double s = 0x1.71547652b82fep+12;
|
||||
double ax = fabs (x), v0 = fma (ax, s, 0x1.8000002p+26);
|
||||
uint64_t jt = asuint64 (v0);
|
||||
uint64_t v = jt;
|
||||
uint64_t tt = ~((1 << 26) - 1l);
|
||||
v &= tt;
|
||||
double t = asdouble (v) - 0x1.8p26;
|
||||
uint64_t aix = asuint64 (ax);
|
||||
if (__glibc_unlikely (aix < UINT64_C(0x3fd0000000000000)))
|
||||
{ // |x| < 0x1p-2
|
||||
if (__glibc_unlikely (aix < UINT64_C(0x3e57137449123ef7)))
|
||||
{
|
||||
// |x| < x0 = 0x1.7137449123ef7p-26
|
||||
return fma (x, 0x1p-55, x);
|
||||
}
|
||||
/* With p = c[0]*x^3 + c[1]*x^5 + c[2]*x^7 + c[3]*x^9 + c[4]*x^11,
|
||||
q = x + p is a minimax approximation of sinh(x) on [x0,1/4] such that
|
||||
|q - sinh(x)|/x^3 < 2^-56.5839 */
|
||||
static const double c[] = { 0x1.5555555555555p-3, 0x1.111111111151ep-7,
|
||||
0x1.a01a019d0c767p-13, 0x1.71de444a96e11p-19,
|
||||
0x1.ae8465375242p-26 };
|
||||
double x2 = x * x, x3 = x2 * x, x4 = x2 * x2,
|
||||
p = x3
|
||||
* ((c[0] + x2 * c[1]) + x4 * ((c[2] + x2 * c[3]) + x4 * c[4]));
|
||||
// fails with e = x3*0x1.5p-53 and x=0x1.71c5b3515d069p-8 (rndz, no fma)
|
||||
double e = x3 * 0x2.1p-53, lb = x + (p - e), ub = x + (p + e);
|
||||
if (lb == ub)
|
||||
return lb;
|
||||
return as_sinh_zero (x);
|
||||
}
|
||||
if (__glibc_unlikely (aix > UINT64_C(0x408633ce8fb9f87d)))
|
||||
{ // |x| >~ 710.47586
|
||||
if (aix >= EXPONENT_MASK)
|
||||
return x + x; // nan Inf
|
||||
return copysign (0x1p1023, x) * 2.0;
|
||||
}
|
||||
// now 0.25 <= |x| < 710.47586
|
||||
// this branch was checked exhaustively with/without FMA
|
||||
int64_t il = ((uint64_t) jt << 14) >> 40, jl = -il;
|
||||
int64_t i1 = il & 0x3f, i0 = (il >> 6) & 0x3f, ie = il >> 12;
|
||||
int64_t j1 = jl & 0x3f, j0 = (jl >> 6) & 0x3f, je = jl >> 12;
|
||||
double sp = asdouble ((uint64_t) (1022 + ie) << 52),
|
||||
sm = asdouble ((uint64_t) (1022 + je) << 52);
|
||||
double t0h = t0[i0][1], t0l = t0[i0][0];
|
||||
double t1h = t1[i1][1], t1l = t1[i1][0];
|
||||
double th = t0h * t1h, tl = t0h * t1l + t1h * t0l + fma (t0h, t1h, -th);
|
||||
const double l2h = 0x1.62e42ffp-13, l2l = 0x1.718432a1b0e26p-47;
|
||||
double dx = (ax - l2h * t) + l2l * t, dx2 = dx * dx, mx = -dx;
|
||||
static const double ch[]
|
||||
= { 0x1p+0, 0x1p-1, 0x1.5555555aaaaaep-3, 0x1.55555551c98cp-5 };
|
||||
double pp = dx * ((ch[0] + dx * ch[1]) + dx2 * (ch[2] + dx * ch[3]));
|
||||
double rh, rl;
|
||||
if (__glibc_unlikely (aix > UINT64_C(0x4014000000000000)))
|
||||
{ // |x| > 5
|
||||
if (__glibc_unlikely (aix > UINT64_C(0x40425e4f7b2737fa)))
|
||||
{ // |x| >~ 36.736801
|
||||
sp = asdouble ((1021 + ie) << 52);
|
||||
rh = th;
|
||||
rl = tl + th * pp;
|
||||
rh *= copysign (1, x);
|
||||
rl *= copysign (1, x);
|
||||
double e = 0x1.1b6p-63 * th, lb = rh + (rl - e), ub = rh + (rl + e);
|
||||
if (lb == ub)
|
||||
return (lb * sp) * 2;
|
||||
|
||||
th = as_exp_accurate (ax, t, th, tl, &tl);
|
||||
th = fasttwosum (th, tl, &tl);
|
||||
th *= copysign (1, x);
|
||||
tl *= copysign (1, x);
|
||||
uint64_t uh = asuint64 (th), ul = asuint64 (tl);
|
||||
int64_t eh = (uh >> 52) & 0x7ff, el = (ul >> 52) & 0x7ff,
|
||||
ml = (ul + 8) & MANTISSA_MASK;
|
||||
th += tl;
|
||||
th *= 2;
|
||||
th *= asuint64 (sp);
|
||||
if (ml <= 16 || eh - el > 103)
|
||||
return as_sinh_database (x, th);
|
||||
return th;
|
||||
}
|
||||
// now 5 < |x| < 36.736801
|
||||
double q0h = t0[j0][1], q1h = t1[j1][1], qh = q0h * q1h;
|
||||
th *= sp;
|
||||
tl *= sp;
|
||||
qh *= sm;
|
||||
double pm = mx * ((ch[0] + mx * ch[1]) + dx2 * (ch[2] + mx * ch[3]));
|
||||
double em = qh + qh * pm;
|
||||
rh = th;
|
||||
rl = (tl - em) + th * pp;
|
||||
|
||||
rh *= copysign (1, x);
|
||||
rl *= copysign (1, x);
|
||||
// fails with e = 0x1.1dbp-63*rh and x=0x1.4971fd7b64137p+2 (rndz, no
|
||||
// fma)
|
||||
double e = 0x1.202p-63 * rh, lb = rh + (rl - e), ub = rh + (rl + e);
|
||||
if (lb == ub)
|
||||
return lb;
|
||||
|
||||
th = as_exp_accurate (ax, t, th, tl, &tl);
|
||||
if (__glibc_unlikely (aix > UINT64_C(0x403f666666666666)))
|
||||
{ // |x| > 31.4
|
||||
rh = th - qh;
|
||||
rl = ((th - rh) - qh) + tl;
|
||||
}
|
||||
else
|
||||
{ // 5 < |x| <= 31.4
|
||||
qh = q0h * q1h;
|
||||
double q0l = t0[j0][0], q1l = t1[j1][0];
|
||||
double ql = q0h * q1l + q1h * q0l + fma (q0h, q1h, -qh);
|
||||
qh *= sm;
|
||||
ql *= sm;
|
||||
qh = as_exp_accurate (-ax, -t, qh, ql, &ql);
|
||||
rh = th - qh;
|
||||
rl = (((th - rh) - qh) - ql) + tl;
|
||||
}
|
||||
}
|
||||
else
|
||||
{ // 0.25 <= |x| <= 5
|
||||
double q0h = t0[j0][1], q0l = t0[j0][0];
|
||||
double q1h = t1[j1][1], q1l = t1[j1][0];
|
||||
double qh = q0h * q1h, ql = q0h * q1l + q1h * q0l + fma (q0h, q1h, -qh);
|
||||
th *= sp;
|
||||
tl *= sp;
|
||||
qh *= sm;
|
||||
ql *= sm;
|
||||
double pm = mx * ((ch[0] + mx * ch[1]) + dx2 * (ch[2] + mx * ch[3]));
|
||||
double fph = th, fpl = tl + th * pp;
|
||||
double fmh = qh, fml = ql + qh * pm;
|
||||
|
||||
rh = fph - fmh;
|
||||
rl = ((fph - rh) - fmh) - fml + fpl;
|
||||
rh *= copysign (1, x);
|
||||
rl *= copysign (1, x);
|
||||
double e = 0x1.c0ap-62 * rh, lb = rh + (rl - e), ub = rh + (rl + e);
|
||||
if (lb == ub)
|
||||
return lb;
|
||||
th = as_exp_accurate (ax, t, th, tl, &tl);
|
||||
qh = as_exp_accurate (-ax, -t, qh, ql, &ql);
|
||||
rh = th - qh;
|
||||
rl = ((th - rh) - qh) - ql + tl;
|
||||
}
|
||||
rh = fasttwosum (rh, rl, &rl);
|
||||
uint64_t uh = asuint64 (rh), ul = asuint64 (rl);
|
||||
int64_t eh = (uh >> MANTISSA_WIDTH) & 0x7ff,
|
||||
el = (ul >> MANTISSA_WIDTH) & 0x7ff,
|
||||
ml = (ul + 8) & MANTISSA_MASK;
|
||||
rh *= copysign (1, x);
|
||||
rl *= copysign (1, x);
|
||||
rh += rl;
|
||||
// fails with ml<=14 and ul.u + 7 above with x=0x1.c13876341b62ep-1 and rndz
|
||||
if (__glibc_unlikely (ml <= 16 || eh - el > 103))
|
||||
return as_sinh_database (x, rh);
|
||||
return rh;
|
||||
}
|
||||
|
||||
#ifndef __ieee754_sinh
|
||||
|
||||
@@ -83,6 +83,18 @@ double: 0
|
||||
Function: "lgamma_upward":
|
||||
double: 0
|
||||
|
||||
Function: "sinh":
|
||||
double: 0
|
||||
|
||||
Function: "sinh_downward":
|
||||
double: 0
|
||||
|
||||
Function: "sinh_towardzero":
|
||||
double: 0
|
||||
|
||||
Function: "sinh_upward":
|
||||
double: 0
|
||||
|
||||
Function: "tgamma":
|
||||
double: 0
|
||||
|
||||
|
||||
Reference in New Issue
Block a user