[libc][math] Added sinhf function.

Differential Revision: https://reviews.llvm.org/D129278
This commit is contained in:
Kirill Okhotnikov
2022-07-07 13:48:11 +02:00
parent fcb9d7e2cf
commit a7f55f0805
15 changed files with 332 additions and 0 deletions

View File

@@ -187,6 +187,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.roundf
libc.src.math.roundl
libc.src.math.sincosf
libc.src.math.sinhf
libc.src.math.sinf
libc.src.math.sqrt
libc.src.math.sqrtf

View File

@@ -206,6 +206,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.roundf
libc.src.math.roundl
libc.src.math.sincosf
libc.src.math.sinhf
libc.src.math.sinf
libc.src.math.sqrt
libc.src.math.sqrtf

View File

@@ -213,6 +213,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.roundl
libc.src.math.sin
libc.src.math.sincosf
libc.src.math.sinhf
libc.src.math.sinf
libc.src.math.sqrt
libc.src.math.sqrtf

View File

@@ -191,6 +191,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sin
libc.src.math.sincosf
libc.src.math.sinf
libc.src.math.sinhf
libc.src.math.sqrt
libc.src.math.sqrtf
libc.src.math.sqrtl

View File

@@ -472,6 +472,7 @@ def StdC : StandardSpec<"stdc"> {
FunctionSpec<"nextafterl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoubleType>]>,
FunctionSpec<"coshf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
FunctionSpec<"sinhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
]
>;

View File

@@ -182,6 +182,7 @@ add_math_entrypoint_object(sincosf)
add_math_entrypoint_object(sin)
add_math_entrypoint_object(sinf)
add_math_entrypoint_object(sinhf)
add_math_entrypoint_object(sqrt)
add_math_entrypoint_object(sqrtf)

View File

@@ -1147,3 +1147,21 @@ add_entrypoint_object(
-O3
)
add_entrypoint_object(
sinhf
SRCS
sinhf.cpp
HDRS
../sinhf.h
expxf.h
DEPENDS
.common_constants
libc.src.__support.FPUtil.fputil
libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.nearest_integer
libc.src.__support.FPUtil.polyeval
libc.include.math
COMPILE_OPTIONS
-O3
)

View File

@@ -0,0 +1,84 @@
//===-- Single-precision sinh function ------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
#include "src/math/sinhf.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/math/generic/expxf.h"
namespace __llvm_libc {
LLVM_LIBC_FUNCTION(float, sinhf, (float x)) {
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
bool sign = xbits.get_sign();
uint32_t x_abs = xbits.uintval() & FPBits::FloatProp::EXP_MANT_MASK;
// |x| <= 2^-26
if (unlikely(x_abs <= 0x3280'0000U)) {
return unlikely(x_abs == 0) ? x : (x + 0.25 * x * x * x);
}
// When |x| >= 90, or x is inf or nan
if (unlikely(x_abs >= 0x42b4'0000U)) {
if (xbits.is_nan())
return x + 1.0f; // sNaN to qNaN + signal
if (xbits.is_inf())
return x;
int rounding = fputil::get_round();
if (sign) {
if (unlikely(rounding == FE_UPWARD || rounding == FE_TOWARDZERO))
return FPBits(FPBits::MAX_NORMAL | FPBits::FloatProp::SIGN_MASK)
.get_val();
} else {
if (unlikely(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO))
return FPBits(FPBits::MAX_NORMAL).get_val();
}
errno = ERANGE;
return x + FPBits::inf(sign).get_val();
}
// |x| <= 0.078125
if (unlikely(x_abs <= 0x3da0'0000U)) {
// |x| = 0.0005589424981735646724700927734375
if (unlikely(x_abs == 0x3a12'85ffU)) {
if (fputil::get_round() == FE_TONEAREST)
return x;
}
double xdbl = x;
double x2 = xdbl * xdbl;
// Sollya: fpminimax(sinh(x),[|3,5,7|],[|D...|],[-1/16-1/64;1/16+1/64],x);
// Sollya output: x * (0x1p0 + x^0x1p1 * (0x1.5555555556583p-3 + x^0x1p1
// * (0x1.111110d239f1fp-7
// + x^0x1p1 * 0x1.a02b5a284013cp-13)))
// Therefore, output of Sollya = x * pe;
double pe = fputil::polyeval(x2, 0.0, 0x1.5555555556583p-3,
0x1.111110d239f1fp-7, 0x1.a02b5a284013cp-13);
return fputil::multiply_add(xdbl, pe, xdbl);
}
// MULT_POWER2 = -1
auto ep_p = exp_eval<-1>(x); // 0.5 * exp(x)
auto ep_m = exp_eval<-1>(-x); // 0.5 * exp(-x)
// 0.5 * expm1(x) = ep_p.mult_exp * (ep_p.r + 1) - 0.5
// = ep_p.mult_exp * ep_p.r + ep_p.mult_exp - 0.5
// 0.5 * expm1(-x) = ep_m.mult_exp * (ep_m.r + 1) - 0.5
// = ep_m.mult_exp * ep_m.r + ep_m.mult_exp - 0.5
// sinh(x) = 0.5 * expm1(x) - 0.5 * expm1(-x)
// Using expm1 instead of exp improved precision around zero.
double ep = fputil::multiply_add(ep_p.mult_exp, ep_p.r, ep_p.mult_exp - 0.5) -
fputil::multiply_add(ep_m.mult_exp, ep_m.r, ep_m.mult_exp - 0.5);
return ep;
}
} // namespace __llvm_libc

18
libc/src/math/sinhf.h Normal file
View File

@@ -0,0 +1,18 @@
//===-- Implementation header for sinhf -------------------------*- C++ -*-===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC_MATH_SINHF_H
#define LLVM_LIBC_SRC_MATH_SINHF_H
namespace __llvm_libc {
float sinhf(float x);
} // namespace __llvm_libc
#endif // LLVM_LIBC_SRC_MATH_SINHF_H

View File

@@ -1330,6 +1330,22 @@ add_fp_unittest(
libc.src.__support.FPUtil.fputil
)
add_fp_unittest(
sinhf_test
NEED_MPFR
SUITE
libc_math_unittests
SRCS
sinhf_test.cpp
HDRS
sdcomp26094.h
DEPENDS
libc.include.errno
libc.src.math.sinhf
libc.src.__support.CPP.array
libc.src.__support.FPUtil.fputil
)
add_subdirectory(generic)
add_subdirectory(exhaustive)
add_subdirectory(differential_testing)

View File

@@ -219,3 +219,20 @@ add_fp_unittest(
-lpthread
)
add_fp_unittest(
sinhf_test
NO_RUN_POSTBUILD
NEED_MPFR
SUITE
libc_math_exhaustive_tests
SRCS
sinhf_test.cpp
DEPENDS
.exhaustive_test
libc.include.math
libc.src.math.sinhf
libc.src.__support.FPUtil.fputil
LINK_LIBRARIES
-lpthread
)

View File

@@ -0,0 +1,75 @@
//===-- Exhaustive test for sinhf -----------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
#include "exhaustive_test.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/math/sinhf.h"
#include "utils/MPFRWrapper/MPFRUtils.h"
#include "utils/UnitTest/FPMatcher.h"
#include <thread>
using FPBits = __llvm_libc::fputil::FPBits<float>;
namespace mpfr = __llvm_libc::testing::mpfr;
struct LlvmLibcSinhfExhaustiveTest : public LlvmLibcExhaustiveTest<uint32_t> {
bool check(uint32_t start, uint32_t stop,
mpfr::RoundingMode rounding) override {
mpfr::ForceRoundingMode r(rounding);
uint32_t bits = start;
bool result = true;
do {
FPBits xbits(bits);
float x = float(xbits);
result &= EXPECT_MPFR_MATCH(mpfr::Operation::Sinh, x,
__llvm_libc::sinhf(x), 0.5, rounding);
} while (bits++ < stop);
return result;
}
};
// Range: [0, 90];
static constexpr uint32_t POS_START = 0x0000'0000U;
static constexpr uint32_t POS_STOP = 0x42b4'0000U;
TEST_F(LlvmLibcSinhfExhaustiveTest, PostiveRangeRoundNearestTieToEven) {
test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Nearest);
}
TEST_F(LlvmLibcSinhfExhaustiveTest, PostiveRangeRoundUp) {
test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Upward);
}
TEST_F(LlvmLibcSinhfExhaustiveTest, PostiveRangeRoundDown) {
test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::Downward);
}
TEST_F(LlvmLibcSinhfExhaustiveTest, PostiveRangeRoundTowardZero) {
test_full_range(POS_START, POS_STOP, mpfr::RoundingMode::TowardZero);
}
// Range: [-90, 0];
static constexpr uint32_t NEG_START = 0x8000'0000U;
static constexpr uint32_t NEG_STOP = 0xc2b4'0000U;
TEST_F(LlvmLibcSinhfExhaustiveTest, NegativeRangeRoundNearestTieToEven) {
test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Nearest);
}
TEST_F(LlvmLibcSinhfExhaustiveTest, NegativeRangeRoundUp) {
test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Upward);
}
TEST_F(LlvmLibcSinhfExhaustiveTest, NegativeRangeRoundDown) {
test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::Downward);
}
TEST_F(LlvmLibcSinhfExhaustiveTest, NegativeRangeRoundTowardZero) {
test_full_range(NEG_START, NEG_STOP, mpfr::RoundingMode::TowardZero);
}

View File

@@ -0,0 +1,89 @@
//===-- Unittests for sinhf -----------------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
#include "src/__support/CPP/Array.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/math/sinhf.h"
#include "utils/MPFRWrapper/MPFRUtils.h"
#include "utils/UnitTest/FPMatcher.h"
#include "utils/UnitTest/Test.h"
#include <math.h>
#include <errno.h>
#include <stdint.h>
using FPBits = __llvm_libc::fputil::FPBits<float>;
namespace mpfr = __llvm_libc::testing::mpfr;
DECLARE_SPECIAL_CONSTANTS(float)
TEST(LlvmLibcSinhfTest, SpecialNumbers) {
errno = 0;
EXPECT_FP_EQ(aNaN, __llvm_libc::sinhf(aNaN));
EXPECT_MATH_ERRNO(0);
EXPECT_FP_EQ(0.0f, __llvm_libc::sinhf(0.0f));
EXPECT_MATH_ERRNO(0);
EXPECT_FP_EQ(-0.0f, __llvm_libc::sinhf(-0.0f));
EXPECT_MATH_ERRNO(0);
EXPECT_FP_EQ(inf, __llvm_libc::sinhf(inf));
EXPECT_MATH_ERRNO(0);
EXPECT_FP_EQ(neg_inf, __llvm_libc::sinhf(neg_inf));
EXPECT_MATH_ERRNO(0);
}
TEST(LlvmLibcSinhfTest, InFloatRange) {
constexpr uint32_t COUNT = 1000000;
constexpr uint32_t STEP = UINT32_MAX / COUNT;
for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
float x = float(FPBits(v));
if (isnan(x) || isinf(x))
continue;
ASSERT_MPFR_MATCH(mpfr::Operation::Sinh, x, __llvm_libc::sinhf(x), 0.5);
}
}
// For small values, sinh(x) is x.
TEST(LlvmLibcSinhfTest, SmallValues) {
float x = float(FPBits(uint32_t(0x17800000)));
float result = __llvm_libc::sinhf(x);
EXPECT_MPFR_MATCH(mpfr::Operation::Sinh, x, result, 0.5);
EXPECT_FP_EQ(x, result);
x = float(FPBits(uint32_t(0x00400000)));
result = __llvm_libc::sinhf(x);
EXPECT_MPFR_MATCH(mpfr::Operation::Sinh, x, result, 0.5);
EXPECT_FP_EQ(x, result);
}
TEST(LlvmLibcSinhfTest, Overflow) {
errno = 0;
EXPECT_FP_EQ(inf, __llvm_libc::sinhf(float(FPBits(0x7f7fffffU))));
EXPECT_MATH_ERRNO(ERANGE);
EXPECT_FP_EQ(inf, __llvm_libc::sinhf(float(FPBits(0x42cffff8U))));
EXPECT_MATH_ERRNO(ERANGE);
EXPECT_FP_EQ(inf, __llvm_libc::sinhf(float(FPBits(0x42d00008U))));
EXPECT_MATH_ERRNO(ERANGE);
}
TEST(LlvmLibcSinhfTest, ExceptionalValues) {
float x = float(FPBits(uint32_t(0x3a12'85ffU)));
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinh, x,
__llvm_libc::sinhf(x), 0.5);
x = -float(FPBits(uint32_t(0x3a12'85ffU)));
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinh, x,
__llvm_libc::sinhf(x), 0.5);
}

View File

@@ -337,6 +337,12 @@ public:
return result;
}
MPFRNumber sinh() const {
MPFRNumber result(*this);
mpfr_sinh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber sqrt() const {
MPFRNumber result(*this);
mpfr_sqrt(result.value, value, mpfr_rounding);
@@ -515,6 +521,8 @@ unary_operation(Operation op, InputType input, unsigned int precision,
return mpfrInput.round();
case Operation::Sin:
return mpfrInput.sin();
case Operation::Sinh:
return mpfrInput.sinh();
case Operation::Sqrt:
return mpfrInput.sqrt();
case Operation::Tan:

View File

@@ -41,6 +41,7 @@ enum class Operation : int {
ModPIOver4,
Round,
Sin,
Sinh,
Sqrt,
Tan,
Trunc,