libc/src/math/generic/sinhf.cpp

   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/macros/optimization.h" // LIBC_UNLIKELY
  12 #include "src/math/generic/explogxf.h"
  13
  14 namespace __llvm_libc {
  15
  16 LLVM_LIBC_FUNCTION(float, sinhf, (float x)) {
  17   using FPBits = typename fputil::FPBits<float>;
  18   FPBits xbits(x);
  19   bool sign = xbits.get_sign();
  20   uint32_t x_abs = xbits.uintval() & FPBits::FloatProp::EXP_MANT_MASK;
  21
  22   // |x| <= 2^-26
  23   if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
  24     return static_cast<float>(
  25         LIBC_UNLIKELY(x_abs == 0) ? x : (x + 0.25 * x * x * x));
  26   }
  27
  28   // When |x| >= 90, or x is inf or nan
  29   if (LIBC_UNLIKELY(x_abs >= 0x42b4'0000U)) {
  30     if (xbits.is_nan())
  31       return x + 1.0f; // sNaN to qNaN + signal
  32
  33     if (xbits.is_inf())
  34       return x;
  35
  36     int rounding = fputil::get_round();
  37     if (sign) {
  38       if (LIBC_UNLIKELY(rounding == FE_UPWARD || rounding == FE_TOWARDZERO))
  39         return FPBits(FPBits::MAX_NORMAL | FPBits::FloatProp::SIGN_MASK)
  40             .get_val();
  41     } else {
  42       if (LIBC_UNLIKELY(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO))
  43         return FPBits(FPBits::MAX_NORMAL).get_val();
  44     }
  45
  46     fputil::set_errno_if_required(ERANGE);
  47     fputil::raise_except_if_required(FE_OVERFLOW);
  48
  49     return x + FPBits::inf(sign).get_val();
  50   }
  51
  52   // |x| <= 0.078125
  53   if (LIBC_UNLIKELY(x_abs <= 0x3da0'0000U)) {
  54     // |x| = 0.0005589424981735646724700927734375
  55     if (LIBC_UNLIKELY(x_abs == 0x3a12'85ffU)) {
  56       if (fputil::get_round() == FE_TONEAREST)
  57         return x;
  58     }
  59
  60     double xdbl = x;
  61     double x2 = xdbl * xdbl;
  62     // Sollya: fpminimax(sinh(x),[|3,5,7|],[|D...|],[-1/16-1/64;1/16+1/64],x);
  63     // Sollya output: x * (0x1p0 + x^0x1p1 * (0x1.5555555556583p-3 + x^0x1p1
  64     //                  * (0x1.111110d239f1fp-7
  65     //                  + x^0x1p1 * 0x1.a02b5a284013cp-13)))
  66     // Therefore, output of Sollya = x * pe;
  67     double pe = fputil::polyeval(x2, 0.0, 0x1.5555555556583p-3,
  68                                  0x1.111110d239f1fp-7, 0x1.a02b5a284013cp-13);
  69     return static_cast<float>(fputil::multiply_add(xdbl, pe, xdbl));
  70   }
  71
  72   // sinh(x) = (e^x - e^(-x)) / 2.
  73   return static_cast<float>(exp_pm_eval</*is_sinh*/ true>(x));
  74 }
  75
  76 } // namespace __llvm_libc