1//===-- Single-precision sinh function ------------------------------------===//
2//
3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4// See https://llvm.org/LICENSE.txt for license information.
5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6//
7//===----------------------------------------------------------------------===//
8
9#include "src/math/sinhf.h"
10#include "src/__support/FPUtil/FPBits.h"
11#include "src/__support/FPUtil/rounding_mode.h"
12#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
13#include "src/math/generic/explogxf.h"
14
15namespace LIBC_NAMESPACE {
16
17LLVM_LIBC_FUNCTION(float, sinhf, (float x)) {
18 using FPBits = typename fputil::FPBits<float>;
19 FPBits xbits(x);
20 uint32_t x_abs = xbits.abs().uintval();
21
22 // When |x| >= 90, or x is inf or nan
23 if (LIBC_UNLIKELY(x_abs >= 0x42b4'0000U || x_abs <= 0x3da0'0000U)) {
24 // |x| <= 0.078125
25 if (x_abs <= 0x3da0'0000U) {
26 // |x| = 0.0005589424981735646724700927734375
27 if (LIBC_UNLIKELY(x_abs == 0x3a12'85ffU)) {
28 if (fputil::fenv_is_round_to_nearest())
29 return x;
30 }
31
32 // |x| <= 2^-26
33 if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
34 return static_cast<float>(
35 LIBC_UNLIKELY(x_abs == 0) ? x : (x + 0.25 * x * x * x));
36 }
37
38 double xdbl = x;
39 double x2 = xdbl * xdbl;
40 // Sollya: fpminimax(sinh(x),[|3,5,7|],[|D...|],[-1/16-1/64;1/16+1/64],x);
41 // Sollya output: x * (0x1p0 + x^0x1p1 * (0x1.5555555556583p-3 + x^0x1p1
42 // * (0x1.111110d239f1fp-7
43 // + x^0x1p1 * 0x1.a02b5a284013cp-13)))
44 // Therefore, output of Sollya = x * pe;
45 double pe = fputil::polyeval(x: x2, a0: 0.0, a: 0x1.5555555556583p-3,
46 a: 0x1.111110d239f1fp-7, a: 0x1.a02b5a284013cp-13);
47 return static_cast<float>(fputil::multiply_add(x: xdbl, y: pe, z: xdbl));
48 }
49
50 if (xbits.is_nan())
51 return x + 1.0f; // sNaN to qNaN + signal
52
53 if (xbits.is_inf())
54 return x;
55
56 int rounding = fputil::quick_get_round();
57 if (xbits.is_neg()) {
58 if (LIBC_UNLIKELY(rounding == FE_UPWARD || rounding == FE_TOWARDZERO))
59 return -FPBits::max_normal().get_val();
60 } else {
61 if (LIBC_UNLIKELY(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO))
62 return FPBits::max_normal().get_val();
63 }
64
65 fputil::set_errno_if_required(ERANGE);
66 fputil::raise_except_if_required(FE_OVERFLOW);
67
68 return x + FPBits::inf(sign: xbits.sign()).get_val();
69 }
70
71 // sinh(x) = (e^x - e^(-x)) / 2.
72 return static_cast<float>(exp_pm_eval</*is_sinh*/ true>(x));
73}
74
75} // namespace LIBC_NAMESPACE
76

source code of libc/src/math/generic/sinhf.cpp