1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
//===-- 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/__support/macros/optimization.h" // LIBC_UNLIKELY
#include "src/math/generic/explogxf.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 (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
return static_cast<float>(
LIBC_UNLIKELY(x_abs == 0) ? x : (x + 0.25 * x * x * x));
}
// When |x| >= 90, or x is inf or nan
if (LIBC_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 (LIBC_UNLIKELY(rounding == FE_UPWARD || rounding == FE_TOWARDZERO))
return FPBits(FPBits::MAX_NORMAL | FPBits::FloatProp::SIGN_MASK)
.get_val();
} else {
if (LIBC_UNLIKELY(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO))
return FPBits(FPBits::MAX_NORMAL).get_val();
}
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
return x + FPBits::inf(sign).get_val();
}
// |x| <= 0.078125
if (LIBC_UNLIKELY(x_abs <= 0x3da0'0000U)) {
// |x| = 0.0005589424981735646724700927734375
if (LIBC_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 static_cast<float>(fputil::multiply_add(xdbl, pe, xdbl));
}
// sinh(x) = (e^x - e^(-x)) / 2.
return static_cast<float>(exp_pm_eval</*is_sinh*/ true>(x));
}
} // namespace __llvm_libc
|
