libc/src/math/generic/asinhf.cpp

   1 //===-- Single-precision asinh 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/asinhf.h"
  10 #include "src/__support/FPUtil/FPBits.h"
  11 #include "src/__support/FPUtil/PolyEval.h"
  12 #include "src/__support/FPUtil/multiply_add.h"
  13 #include "src/__support/FPUtil/sqrt.h"
  14 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
  15 #include "src/math/generic/common_constants.h"
  16 #include "src/math/generic/explogxf.h"
  17
  18 namespace __llvm_libc {
  19
  20 LLVM_LIBC_FUNCTION(float, asinhf, (float x)) {
  21   using FPBits_t = typename fputil::FPBits<float>;
  22   FPBits_t xbits(x);
  23   uint32_t x_u = xbits.uintval();
  24   uint32_t x_abs = x_u & FPBits_t::FloatProp::EXP_MANT_MASK;
  25
  26   // |x| <= 2^-3
  27   if (LIBC_UNLIKELY(x_abs <= 0x3e80'0000U)) {
  28     // |x| <= 2^-26
  29     if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
  30       return static_cast<float>(LIBC_UNLIKELY(x_abs == 0)
  31                                     ? x
  32                                     : (x - 0x1.5555555555555p-3 * x * x * x));
  33     }
  34
  35     double x_d = x;
  36     double x_sq = x_d * x_d;
  37     // Generated by Sollya with:
  38     // > P = fpminimax(asinh(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16], [|D...|],
  39     //                 [0, 2^-2]);
  40     double p = fputil::polyeval(
  41         x_sq, 0.0, -0x1.555555555551ep-3, 0x1.3333333325495p-4,
  42         -0x1.6db6db5a7622bp-5, 0x1.f1c70f82928c6p-6, -0x1.6e893934266b7p-6,
  43         0x1.1c0b41d3fbe78p-6, -0x1.c0f47810b3c4fp-7, 0x1.2c8602690143dp-7);
  44     return static_cast<float>(fputil::multiply_add(x_d, p, x_d));
  45   }
  46
  47   const double SIGN[2] = {1.0, -1.0};
  48   double x_sign = SIGN[x_u >> 31];
  49   double x_d = x;
  50
  51   // Helper functions to set results for exceptional cases.
  52   auto round_result_slightly_down = [x_sign](float r) -> float {
  53     return fputil::multiply_add(static_cast<float>(x_sign), r,
  54                                 static_cast<float>(x_sign) * (-0x1.0p-24f));
  55   };
  56   auto round_result_slightly_up = [x_sign](float r) -> float {
  57     return fputil::multiply_add(static_cast<float>(x_sign), r,
  58                                 static_cast<float>(x_sign) * 0x1.0p-24f);
  59   };
  60
  61   if (LIBC_UNLIKELY(x_abs >= 0x4bdd'65a5U)) {
  62     if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) {
  63       // x is +-inf or nan
  64       return x;
  65     }
  66
  67     // Exceptional cases when x > 2^24.
  68     switch (x_abs) {
  69     case 0x4bdd65a5: // |x| = 0x1.bacb4ap24f
  70       return round_result_slightly_down(0x1.1e0696p4f);
  71     case 0x4c803f2c: // |x| = 0x1.007e58p26f
  72       return round_result_slightly_down(0x1.2b786cp4f);
  73     case 0x4f8ffb03: // |x| = 0x1.1ff606p32f
  74       return round_result_slightly_up(0x1.6fdd34p4f);
  75     case 0x5c569e88: // |x| = 0x1.ad3d1p57f
  76       return round_result_slightly_up(0x1.45c146p5f);
  77     case 0x5e68984e: // |x| = 0x1.d1309cp61f
  78       return round_result_slightly_up(0x1.5c9442p5f);
  79     case 0x655890d3: // |x| = 0x1.b121a6p75f
  80       return round_result_slightly_down(0x1.a9a3f2p5f);
  81     case 0x65de7ca6: // |x| = 0x1.bcf94cp76f
  82       return round_result_slightly_up(0x1.af66cp5f);
  83     case 0x6eb1a8ec: // |x| = 0x1.6351d8p94f
  84       return round_result_slightly_down(0x1.08b512p6f);
  85     case 0x7997f30a: // |x| = 0x1.2fe614p116f
  86       return round_result_slightly_up(0x1.451436p6f);
  87     }
  88   } else {
  89     // Exceptional cases when x < 2^24.
  90     if (LIBC_UNLIKELY(x_abs == 0x45abaf26)) {
  91       // |x| = 0x1.575e4cp12f
  92       return round_result_slightly_down(0x1.29becap3f);
  93     }
  94     if (LIBC_UNLIKELY(x_abs == 0x49d29048)) {
  95       // |x| = 0x1.a5209p20f
  96       return round_result_slightly_down(0x1.e1b92p3f);
  97     }
  98   }
  99
 100   // asinh(x) = log(x + sqrt(x^2 + 1))
 101   return static_cast<float>(
 102       x_sign *
 103       log_eval(fputil::multiply_add(
 104           x_d, x_sign, fputil::sqrt(fputil::multiply_add(x_d, x_d, 1.0)))));
 105 }
 106
 107 } // namespace __llvm_libc