libc/src/__support/math/asinbf16.h

//===-- Implementation header for asinbf16 ----------------------*- 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___SUPPORT_MATH_ASINBF16_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINBF16_H

#include "inv_trigf_utils.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/bfloat16.h"
#include "src/__support/FPUtil/cast.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/sqrt.h"
#include "src/__support/macros/optimization.h"

namespace LIBC_NAMESPACE_DECL {
namespace math {

LIBC_INLINE LIBC_CONSTEXPR bfloat16 asinbf16(bfloat16 x) {
  // Generated by Sollya using the following command:
  // > display = hexadecimal;
  // > round(pi/2, SG, RN);
  constexpr float PI_2 = 0x1.921fb6p0f;

  using FPBits = fputil::FPBits<bfloat16>;
  FPBits xbits(x);

  uint16_t x_u = xbits.uintval();
  uint16_t x_abs = x_u & 0x7fff;
  float x_sign = (x_u >> 15) ? -1 : 1;
  float xf = x;

  float xf_abs = (xf < 0 ? -xf : xf);
  float x_sq = xf_abs * xf_abs;

  // Case 1: |x| <= 0.5
  if (x_abs <= 0x3F00) { // x_abs <= 0.5
    // |x| = {0}
    if (LIBC_UNLIKELY(x_abs == 0))
      return x; // with sign

    if (LIBC_UNLIKELY(x_abs <= 0x3D00)) {
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
      int rounding = fputil::quick_get_round();
      if ((xbits.is_pos() && rounding == FE_UPWARD) ||
          (xbits.is_neg() && rounding == FE_DOWNWARD)) {
        return fputil::cast<bfloat16>(fputil::multiply_add(xf, 0x1.0p-9f, xf));
      }
#endif
      return x;
    }

    float xp = fputil::cast<float>(inv_trigf_utils_internal::asin_eval(x_sq));
    float result = xf * (fputil::multiply_add<float>(x_sq, xp, 1.0f));
    return fputil::cast<bfloat16>(result);
  }

  // Case 2: 0.5 <|x| <= 1
  //  using reduction: asin(x) = pi/2 - 2*asin(sqrt((1-x)/2))
  if (x_abs <= 0x3F80) { // x_abs <= 1
    // |x| = {1}
    if (LIBC_UNLIKELY(x_abs == 0x3F80)) {
      return fputil::cast<bfloat16>(x_sign * PI_2);
    }

    float t = fputil::multiply_add<float>(xf_abs, -0.5f, 0.5f);
    float t_sqrt = fputil::sqrt<float>(t);
    float tp = fputil::cast<float>(inv_trigf_utils_internal::asin_eval(t));
    float asin_sqrt_t = t_sqrt * (fputil::multiply_add<float>(t, tp, 1.0f));
    float result = fputil::multiply_add<float>(-2.0f, asin_sqrt_t, PI_2);
    return fputil::cast<bfloat16>(x_sign * result);
  }

  // Case 3: NaN and Inf
  // NaN
  if (xbits.is_nan()) {
    if (xbits.is_signaling_nan()) {
      fputil::raise_except_if_required(FE_INVALID);
      return FPBits::quiet_nan().get_val();
    }
    return x; // quiet NaN
  }
  // |x|>1 & inf
  fputil::raise_except_if_required(FE_INVALID);
  fputil::set_errno_if_required(EDOM); // Domain is bounded
  return FPBits::quiet_nan().get_val();
}

} // namespace math
} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINBF16_H