libc/src/__support/math/erfcf16.h

//===-- Implementation header for erfcf16 -----------------------*- 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_ERFCF16_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_ERFCF16_H

#include "include/llvm-libc-macros/float16-macros.h"

#ifdef LIBC_TYPES_HAS_FLOAT16

#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/cast.h"
#include "src/__support/FPUtil/except_value_utils.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY

namespace LIBC_NAMESPACE_DECL {

namespace math {

LIBC_INLINE float16 erfcf16(float16 x) {
  // Polynomials approximating erfc(|x|) on ( k/2, (k+1)/2 ) generated by
  // Sollya with: > P = fpminimax(erfc(x), [|0, 1, 2, 3, 4, 5, 6, 7|], [|D...|],
  //                             [k/2, (k+1)/2]);
  //
  // for k = 0..7.
  constexpr double COEFFS[8][8] = {
      {0x1.00000006e5898p0, -0x1.20dd7b4435481p0, 0x1.e2ffc0264c37cp-17,
       0x1.80ef5566d3135p-2, 0x1.96ef459387c13p-10, -0x1.e6f847a52d3fep-4,
       0x1.9a59af64a5dfap-7, 0x1.f652f3b14963bp-7},
      {0x1.001a9f1d3c0b4p0, -0x1.221a42637ae8p0, 0x1.9c3ddc1e5d176p-6,
       0x1.347481936316bp-2, 0x1.1db66a2746b8bp-3, -0x1.1d7a264182166p-2,
       0x1.ef858095e1609p-4, -0x1.26936198d6b07p-6},
      {0x1.f08a38741577bp-1, -0x1.e26ae90822fp-1, -0x1.f300a88478724p-2,
       0x1.115a0580ba3f8p0, -0x1.1a1c2bc3ffcdap-1, 0x1.888d0c9a082b8p-4,
       0x1.f41c9f1877fdfp-8, -0x1.a72fb878df30bp-9},
      {0x1.c1cabe622be64p-1, -0x1.ee19257d1037ap-2, -0x1.7a725229a24b8p0,
       0x1.2089bc62b1069p1, -0x1.673f1b9fbced5p0, 0x1.da7db686b0475p-2,
       -0x1.49a292662ca6ap-4, 0x1.7e2b298da504dp-8},
      {0x1.a466e958ca525p0, -0x1.8a850d50339a9p1, 0x1.29b741ccfdd3ap1,
       -0x1.b250e97982f77p-1, 0x1.ea6896c5fa419p-4, 0x1.b332978509bf1p-7,
       -0x1.9f1c5d9122108p-8, 0x1.2fb4d83883203p-11},
      {0x1.12fb16acc5644p1, -0x1.25003c37e1b24p2, 0x1.0e0dc863b2f12p2,
       -0x1.16f179d8797a6p1, 0x1.5c8b140bf43f8p-1, -0x1.07473d994c2edp-3,
       0x1.bd0c1c2d2a9e3p-7, -0x1.448767255aabbp-11},
      {0x1.13a691333e22ap0, -0x1.0e2a642d8d318p1, 0x1.c7d88193df94bp0,
       -0x1.acf7c7b79498fp-1, 0x1.e626f6202a127p-3, -0x1.4babcc8609859p-5,
       0x1.f860060a3658p-9, -0x1.49a4190580d4p-13},
      {0x1.cf3a84bf655afp-3, -0x1.968e6988b5cep-2, 0x1.326f1f9499739p-2,
       -0x1.011785d112c9bp-3, 0x1.0343a0c05e285p-5, -0x1.3a3ae5e48130ep-8,
       0x1.a7c5af0484892p-12, -0x1.ea82a062010c3p-17},
  };

  static constexpr size_t N_ERFCF16_EXCEPTS = 3;
  static constexpr fputil::ExceptValues<float16, N_ERFCF16_EXCEPTS>
      ERFCF16_EXCEPTS = {{
          // (input, RZ output, RU offset, RD offset, RN offset)
          // x = 0x0.000216, erfc(x) = 0x1.000
          {0x0B17, 0x3BFF, 1, 0, 0},
          // x = 0x0.00346, erfc(x) = 0x0.ff8
          {0x1B17, 0x3BF7, 1, 0, 1},
          // x = -0x0.00346, erfc(x) = 0x1.00c
          {0x9B17, 0x3C04, 1, 0, 0},
      }};
  using FPBits = typename fputil::FPBits<float16>;
  FPBits xbits(x);
  uint16_t x_u = xbits.uintval();

  if (auto r = ERFCF16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
    return r.value();

  uint16_t x_abs = xbits.abs().uintval();

  // Special cases: NaN and Inf
  if (LIBC_UNLIKELY(x_abs >= 0x7c00U)) {
    if (x_abs > 0x7c00U) {
      if (xbits.is_signaling_nan()) {
        fputil::raise_except_if_required(FE_INVALID);
        return FPBits::quiet_nan().get_val();
      }
      return x;
    }
    // erfc(+Inf) = 0, erfc(-Inf) = 2
    return xbits.is_neg() ? fputil::cast<float16>(2.0)
                          : fputil::cast<float16>(0.0);
  }

  if (LIBC_UNLIKELY(x_abs == 0))
    return 1.0f16;

  // Asymptotic behavior: erfc(x) rounds to 0 or 2 for |x| >= 4.0.
  if (LIBC_UNLIKELY(x_abs >= 0x4400U)) { // |x| >= 4.0
    if (xbits.is_neg()) {
      // 0x1.0p-12 is ~0.25 ULP of 2.0 in float16, small enough to round
      // to 2.0 in RN, but large enough to round down in RD/RZ.
      float xf = fputil::cast<float>(x);
      return fputil::cast<float16>(2.0 - 0x1.0p-12 * (-4.0 / xf));
    }
    fputil::set_errno_if_required(ERANGE);
    fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
    if (fputil::fenv_is_round_up())
      return FPBits::min_subnormal().get_val();
#endif
    return 0.0f16;
  }

  // Polynomial approximation:
  //   erfc(x) ~ erfc(|x|)      if x >= 0
  //   erfc(x) ~ 2 - erfc(|x|)  if x < 0
  // erfc(|x|) is evaluated using a degree-7 polynomial on each sub-interval.

  int idx = static_cast<int>(xbits.abs().get_val() * 2.0f);
  double xd = fputil::cast<double>(xbits.abs().get_val());
  double xsq = xd * xd;
  double x4 = xsq * xsq;

  double p01 = fputil::multiply_add(xd, COEFFS[idx][1], COEFFS[idx][0]);
  double p23 = fputil::multiply_add(xd, COEFFS[idx][3], COEFFS[idx][2]);
  double p45 = fputil::multiply_add(xd, COEFFS[idx][5], COEFFS[idx][4]);
  double p67 = fputil::multiply_add(xd, COEFFS[idx][7], COEFFS[idx][6]);

  double p03 = fputil::multiply_add(xsq, p23, p01);
  double p47 = fputil::multiply_add(xsq, p67, p45);

  double erfc_abs = fputil::multiply_add(x4, p47, p03);

  if (xbits.is_neg())
    return fputil::cast<float16>(2.0 - erfc_abs);

  return fputil::cast<float16>(erfc_abs);
}

} // namespace math

} // namespace LIBC_NAMESPACE_DECL

#endif // LIBC_TYPES_HAS_FLOAT16

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ERFCF16_H