mirror of
https://sourceware.org/git/glibc.git
synced 2025-10-12 19:04:54 +03:00
b81252c4b9
The libm size improvement built with "--enable-stack-protector=strong --enable-bind-now=yes --enable-fortify-source=2": Before: text data bss dec hex filename 585192 860 12 586064 8f150 aarch64-linux-gnu/math/libm.so 960775 1068 12 961855 ead3f x86_64-linux-gnu/math/libm.so 1189174 5544 368 1195086 123c4e powerpc64le-linux-gnu/math/libm.so After: text data bss dec hex filename 584952 860 12 585824 8f060 aarch64-linux-gnu/math/libm.so 960615 1068 12 961695 eac9f x86_64-linux-gnu/math/libm.so 1189078 5544 368 1194990 123bee powerpc64le-linux-gnu/math/libm.so The are small code changes for x86_64 and powerpc64le, which do not affect performance; but on aarch64 with gcc-14 I see a slight better code generation due the usage of ldq for floating point constant loading. Reviewed-by: Andreas K. Huettel <dilfridge@gentoo.org>
109 lines
3.5 KiB
C
109 lines
3.5 KiB
C
/* Correctly-rounded hyperbolic sine function for binary32 value.
|
|
|
|
Copyright (c) 2022-2024 Alexei Sibidanov.
|
|
|
|
The original version of this file was copied from the CORE-MATH
|
|
project (file src/binary32/sinh/sinhf.c, revision 572ecec).
|
|
|
|
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 <math.h>
|
|
#include <stdint.h>
|
|
#include <libm-alias-finite.h>
|
|
#include <math-narrow-eval.h>
|
|
#include "math_config.h"
|
|
#include "e_sincoshf_data.h"
|
|
|
|
float
|
|
__ieee754_sinhf (float x)
|
|
{
|
|
static const struct
|
|
{
|
|
uint32_t uarg;
|
|
float rh, rl;
|
|
} st = { 0x74250bfeu, 0x1.250bfep-11, 0x1p-36 };
|
|
const double iln2 = 0x1.71547652b82fep+5;
|
|
double z = x;
|
|
uint32_t ux = asuint (x) << 1;
|
|
if (__glibc_unlikely (ux > 0x8565a9f8u))
|
|
{ /* |x| >~ 89.4 */
|
|
float sgn = copysignf (2.0f, x);
|
|
if (ux >= 0xff000000u)
|
|
{
|
|
if (ux << 8)
|
|
return x + x; /* nan */
|
|
return copysignf (INFINITY, x); /* +-inf */
|
|
}
|
|
float r = math_narrow_eval (sgn * 0x1.fffffep127f);
|
|
return r;
|
|
}
|
|
if (__glibc_unlikely (ux < 0x7c000000u))
|
|
{ /* |x| < 0.125 */
|
|
if (__glibc_unlikely (ux <= 0x74250bfeu))
|
|
{ /* |x| <= 0x1.250bfep-11 */
|
|
if (__glibc_unlikely (ux < 0x66000000u)) /* |x| < 0x1p-24 */
|
|
return fmaf (x, fabsf (x), x);
|
|
if (__glibc_unlikely (st.uarg == ux))
|
|
{
|
|
float sgn = copysignf (1.0f, x);
|
|
return sgn * st.rh + sgn * st.rl;
|
|
}
|
|
return (x * 0x1.555556p-3f) * (x * x) + x;
|
|
}
|
|
static const double cp[] =
|
|
{
|
|
0x1.5555555555555p-3, 0x1.11111111146e1p-7,
|
|
0x1.a01a00930dda6p-13, 0x1.71f92198aa6e9p-19
|
|
};
|
|
double z2 = z * z;
|
|
double z4 = z2 * z2;
|
|
return z + (z2 * z)
|
|
* ((cp[0] + z2 * cp[1]) + z4 * (cp[2] + z2 * (cp[3])));
|
|
}
|
|
double a = iln2 * z;
|
|
double ia = roundeven_finite (a);
|
|
double h = a - ia;
|
|
double h2 = h * h;
|
|
int64_t jp = asuint64 (ia + 0x1.8p52);
|
|
int64_t jm = -jp;
|
|
double sp = asdouble (TB[jp & 31] + ((jp >> 5) << 52));
|
|
double sm = asdouble (TB[jm & 31] + ((jm >> 5) << 52));
|
|
double te = C[0] + h2 * C[2];
|
|
double to = (C[1] + h2 * C[3]);
|
|
double rp = sp * (te + h * to);
|
|
double rm = sm * (te - h * to);
|
|
double r = rp - rm;
|
|
float ub = r;
|
|
double lb = r - 1.52e-10 * r;
|
|
if (__glibc_unlikely (ub != lb))
|
|
{
|
|
const double iln2h = 0x1.7154765p+5;
|
|
const double iln2l = 0x1.5c17f0bbbe88p-26;
|
|
h = (iln2h * z - ia) + iln2l * z;
|
|
h2 = h * h;
|
|
te = CH[0] + h2 * CH[2] + (h2 * h2) * (CH[4] + h2 * CH[6]);
|
|
to = CH[1] + h2 * (CH[3] + h2 * CH[5]);
|
|
r = sp * (te + h * to) - sm * (te - h * to);
|
|
ub = r;
|
|
}
|
|
return ub;
|
|
}
|
|
libm_alias_finite (__ieee754_sinhf, __sinhf)
|