libc/src/__support/math/acoshf_utils.h

//===-- Collection of utils for acoshf --------------------------*- 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_ACOSHF_UTILS_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSHF_UTILS_H

#include "acosh_float_constants.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/multiply_add.h"

namespace LIBC_NAMESPACE_DECL {

namespace acoshf_internal {

// x should be positive, normal finite value
LIBC_INLINE static double log_eval(double x) {
  // For x = 2^ex * (1 + mx)
  //   log(x) = ex * log(2) + log(1 + mx)
  using FPB = fputil::FPBits<double>;
  FPB bs(x);

  double ex = static_cast<double>(bs.get_exponent());

  // p1 is the leading 7 bits of mx, i.e.
  // p1 * 2^(-7) <= m_x < (p1 + 1) * 2^(-7).
  int p1 = static_cast<int>(bs.get_mantissa() >> (FPB::FRACTION_LEN - 7));

  // Set bs to (1 + (mx - p1*2^(-7))
  bs.set_uintval(bs.uintval() & (FPB::FRACTION_MASK >> 7));
  bs.set_biased_exponent(FPB::EXP_BIAS);
  // dx = (mx - p1*2^(-7)) / (1 + p1*2^(-7)).
  double dx = (bs.get_val() - 1.0) * ONE_OVER_F[p1];

  // Minimax polynomial of log(1 + dx) generated by Sollya with:
  // > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]);
  const double COEFFS[6] = {-0x1.ffffffffffffcp-2, 0x1.5555555552ddep-2,
                            -0x1.ffffffefe562dp-3, 0x1.9999817d3a50fp-3,
                            -0x1.554317b3f67a5p-3, 0x1.1dc5c45e09c18p-3};
  double dx2 = dx * dx;
  double c1 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]);
  double c2 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]);
  double c3 = fputil::multiply_add(dx, COEFFS[5], COEFFS[4]);

  double p = fputil::polyeval(dx2, dx, c1, c2, c3);
  double result =
      fputil::multiply_add(ex, /*log(2)*/ 0x1.62e42fefa39efp-1, LOG_F[p1] + p);
  return result;
}

} // namespace acoshf_internal

} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOSHF_UTILS_H