[libc][math] Implement C23 half precision erfc function (#180930)

Add support for the half-precision complementary error function
`erfcf16``, using a Sollya generated polynomial implementation with
proper handling of special cases.

Extend the MPFR utilities with erfc support to allow tests.

closes: #180927
This commit is contained in:
Anonmiraj 2026-03-27 16:42:03 +02:00 committed by GitHub
parent c3bffc8e82
commit e06b5e53a2
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
25 changed files with 414 additions and 1 deletions

View File

@ -532,6 +532,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.coshf16
libc.src.math.cospif16
libc.src.math.erff16
libc.src.math.erfcf16
libc.src.math.exp10f16
libc.src.math.exp10m1f16
libc.src.math.exp2f16

View File

@ -534,6 +534,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.coshf16
libc.src.math.cospif16
libc.src.math.erff16
libc.src.math.erfcf16
libc.src.math.exp10f16
libc.src.math.exp10m1f16
libc.src.math.exp2f16

View File

@ -691,6 +691,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.copysignf16
libc.src.math.cospif16
libc.src.math.erff16
libc.src.math.erfcf16
libc.src.math.expf16
libc.src.math.f16add
libc.src.math.f16addf

View File

@ -705,6 +705,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.coshf16
libc.src.math.cospif16
libc.src.math.erff16
libc.src.math.erfcf16
libc.src.math.exp10f16
libc.src.math.exp10m1f16
libc.src.math.exp2f16

View File

@ -759,6 +759,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.coshf16
libc.src.math.cospif16
libc.src.math.erff16
libc.src.math.erfcf16
libc.src.math.exp10f16
libc.src.math.exp10m1f16
libc.src.math.exp2f16

View File

@ -403,6 +403,13 @@ functions:
arguments:
- type: _Float16
guard: LIBC_TYPES_HAS_FLOAT16
- name: erfcf16
standards:
- stdc
return_type: _Float16
arguments:
- type: _Float16
guard: LIBC_TYPES_HAS_FLOAT16
- name: exp
standards:
- stdc

View File

@ -81,6 +81,7 @@
#include "math/dfmaf128.h"
#include "math/dfmal.h"
#include "math/dsqrtl.h"
#include "math/erfcf16.h"
#include "math/erff.h"
#include "math/erff16.h"
#include "math/exp.h"

View File

@ -0,0 +1,27 @@
//===-- Shared erfcf16 function ---------------------------------*- 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_SHARED_MATH_ERFCF16_H
#define LLVM_LIBC_SHARED_MATH_ERFCF16_H
#ifdef LIBC_TYPES_HAS_FLOAT16
#include "shared/libc_common.h"
#include "src/__support/math/erfcf16.h"
namespace LIBC_NAMESPACE_DECL {
namespace shared {
using math::erfcf16;
} // namespace shared
} // namespace LIBC_NAMESPACE_DECL
#endif // LIBC_TYPES_HAS_FLOAT16
#endif // LLVM_LIBC_SHARED_MATH_ERFCF16_H

View File

@ -960,6 +960,20 @@ add_header_library(
libc.src.__support.macros.properties.cpu_features
)
add_header_library(
erfcf16
HDRS
erfcf16.h
DEPENDS
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.multiply_add
libc.src.__support.macros.config
libc.src.__support.macros.optimization
)
add_header_library(
erff
HDRS

View File

@ -0,0 +1,145 @@
//===-- 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);
if (fputil::fenv_is_round_up())
return FPBits::min_subnormal().get_val();
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

View File

@ -142,6 +142,8 @@ add_math_entrypoint_object(erf)
add_math_entrypoint_object(erff)
add_math_entrypoint_object(erff16)
add_math_entrypoint_object(erfcf16)
add_math_entrypoint_object(exp)
add_math_entrypoint_object(expf)
add_math_entrypoint_object(expf16)

21
libc/src/math/erfcf16.h Normal file
View File

@ -0,0 +1,21 @@
//===-- 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_MATH_ERFCF16_H
#define LLVM_LIBC_SRC_MATH_ERFCF16_H
#include "src/__support/macros/config.h"
#include "src/__support/macros/properties/types.h"
namespace LIBC_NAMESPACE_DECL {
float16 erfcf16(float16 x);
} // namespace LIBC_NAMESPACE_DECL
#endif // LLVM_LIBC_SRC_MATH_ERFCF16_H

View File

@ -1292,7 +1292,16 @@ add_entrypoint_object(
../erff16.h
DEPENDS
libc.src.__support.math.erff16
libc.src.errno.errno
)
add_entrypoint_object(
erfcf16
SRCS
erfcf16.cpp
HDRS
../erfcf16.h
DEPENDS
libc.src.__support.math.erfcf16
)
add_entrypoint_object(

View File

@ -0,0 +1,16 @@
//===-- Half-precision erfc(x) function -----------------------------------===//
//
// 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
//
//===----------------------------------------------------------------------===//
#include "src/math/erfcf16.h"
#include "src/__support/math/erfcf16.h"
namespace LIBC_NAMESPACE_DECL {
LLVM_LIBC_FUNCTION(float16, erfcf16, (float16 x)) { return math::erfcf16(x); }
} // namespace LIBC_NAMESPACE_DECL

View File

@ -80,6 +80,7 @@ add_fp_unittest(
libc.src.__support.math.dsqrtl
libc.src.__support.math.exp10m1f
libc.src.__support.math.exp10m1f16
libc.src.__support.math.erfcf16
libc.src.__support.math.erff
libc.src.__support.math.erff16
libc.src.__support.math.exp

View File

@ -32,6 +32,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat16) {
EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::coshf16(0.0f16));
EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::cospif16(0.0f16));
EXPECT_FP_EQ(0x0p+0f16, LIBC_NAMESPACE::shared::erff16(0.0f16));
EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::erfcf16(0.0f));
EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::exp10f16(0.0f16));
EXPECT_FP_EQ(0x0p+0f16, LIBC_NAMESPACE::shared::exp10m1f16(0.0f16));
EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::exp2f16(0.0f16));

View File

@ -2723,6 +2723,18 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
add_fp_unittest(
erfcf16_test
NEED_MPFR
SUITE
libc-math-unittests
SRCS
erfcf16_test.cpp
DEPENDS
libc.src.math.erfcf16
libc.src.__support.FPUtil.fp_bits
)
add_fp_unittest(
pow_test
NEED_MPFR

View File

@ -0,0 +1,43 @@
//===-- Exhaustive test for erfcf16 ---------------------------------------===//
//
// 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
//
//===----------------------------------------------------------------------===//
#include "src/__support/macros/optimization.h"
#include "src/math/erfcf16.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
#include "utils/MPFRWrapper/MPFRUtils.h"
using LlvmLibcErfcf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
// Range: [0, Inf];
// 0x0000 is +0.0, 0x7c00 is +Inf.
static constexpr uint16_t POS_START = 0x0000U;
static constexpr uint16_t POS_STOP = 0x7c00U;
// Range: [-0, -Inf];
// 0x8000 is -0.0, 0xfc00 is -Inf.
static constexpr uint16_t NEG_START = 0x8000U;
static constexpr uint16_t NEG_STOP = 0xfc00U;
TEST_F(LlvmLibcErfcf16Test, PositiveRange) {
for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
float16 x = FPBits(v).get_val();
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Erfc, x,
LIBC_NAMESPACE::erfcf16(x), 0.5);
}
}
TEST_F(LlvmLibcErfcf16Test, NegativeRange) {
for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
float16 x = FPBits(v).get_val();
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Erfc, x,
LIBC_NAMESPACE::erfcf16(x), 0.5);
}
}

View File

@ -5183,6 +5183,17 @@ add_fp_unittest(
libc.src.math.erff16
)
add_fp_unittest(
erfcf16_test
SUITE
libc-math-smoke-tests
SRCS
erfcf16_test.cpp
DEPENDS
libc.hdr.errno_macros
libc.src.math.erfcf16
)
add_fp_unittest(
pow_test
SUITE

View File

@ -0,0 +1,66 @@
//===-- Unittests for erfcf16 ---------------------------------------------===//
//
// 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
//
//===----------------------------------------------------------------------===//
#include "hdr/errno_macros.h"
#include "hdr/math_macros.h"
#include "hdr/stdint_proxy.h"
#include "src/math/erfcf16.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
using LlvmLibcErfcTest = LIBC_NAMESPACE::testing::FPTest<float16>;
TEST_F(LlvmLibcErfcTest, SpecialNumbers) {
EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::erfcf16(sNaN), FE_INVALID);
EXPECT_MATH_ERRNO(0);
EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::erfcf16(aNaN));
EXPECT_FP_EQ_ALL_ROUNDING(0.0f16, LIBC_NAMESPACE::erfcf16(inf));
EXPECT_FP_EQ_ALL_ROUNDING(2.0f16, LIBC_NAMESPACE::erfcf16(neg_inf));
EXPECT_FP_EQ_ALL_ROUNDING(1.0f16, LIBC_NAMESPACE::erfcf16(zero));
EXPECT_FP_EQ_ALL_ROUNDING(1.0f16, LIBC_NAMESPACE::erfcf16(neg_zero));
}
TEST_F(LlvmLibcErfcTest, Underflow) {
EXPECT_FP_EQ_WITH_EXCEPTION(zero, LIBC_NAMESPACE::erfcf16(4.0f16),
FE_UNDERFLOW | FE_INEXACT);
EXPECT_MATH_ERRNO(ERANGE);
EXPECT_FP_EQ_WITH_EXCEPTION(zero, LIBC_NAMESPACE::erfcf16(5.75f16),
FE_UNDERFLOW | FE_INEXACT);
EXPECT_MATH_ERRNO(ERANGE);
EXPECT_FP_EQ_WITH_EXCEPTION(zero, LIBC_NAMESPACE::erfcf16(11.25f16),
FE_UNDERFLOW | FE_INEXACT);
EXPECT_MATH_ERRNO(ERANGE);
}
#ifdef LIBC_TEST_FTZ_DAZ
using namespace LIBC_NAMESPACE::testing;
TEST_F(LlvmLibcErfcTest, FTZMode) {
ModifyMXCSR mxcsr(FTZ);
EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::erfcf16(min_denormal));
EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::erfcf16(max_denormal));
}
TEST_F(LlvmLibcErfcTest, DAZMode) {
ModifyMXCSR mxcsr(DAZ);
EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::erfcf16(min_denormal));
EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::erfcf16(max_denormal));
}
TEST_F(LlvmLibcErfcTest, FTZDAZMode) {
ModifyMXCSR mxcsr(FTZ | DAZ);
EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::erfcf16(min_denormal));
EXPECT_FP_EQ(1.0f16, LIBC_NAMESPACE::erfcf16(max_denormal));
}
#endif

View File

@ -210,6 +210,12 @@ MPFRNumber MPFRNumber::erf() const {
return result;
}
MPFRNumber MPFRNumber::erfc() const {
MPFRNumber result(*this);
mpfr_erfc(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::exp() const {
MPFRNumber result(*this);
mpfr_exp(result.value, value, mpfr_rounding);

View File

@ -200,6 +200,7 @@ public:
MPFRNumber cosh() const;
MPFRNumber cospi() const;
MPFRNumber erf() const;
MPFRNumber erfc() const;
MPFRNumber exp() const;
MPFRNumber exp2() const;
MPFRNumber exp2m1() const;

View File

@ -59,6 +59,8 @@ unary_operation(Operation op, InputType input, unsigned int precision,
return mpfrInput.cospi();
case Operation::Erf:
return mpfrInput.erf();
case Operation::Erfc:
return mpfrInput.erfc();
case Operation::Exp:
return mpfrInput.exp();
case Operation::Exp2:

View File

@ -40,6 +40,7 @@ enum class Operation : int {
Cosh,
Cospi,
Erf,
Erfc,
Exp,
Exp2,
Exp2m1,

View File

@ -3830,6 +3830,21 @@ libc_support_library(
],
)
libc_support_library(
name = "__support_math_erfcf16",
hdrs = ["src/__support/math/erfcf16.h"],
deps = [
":__support_fputil_cast",
":__support_fputil_except_value_utils",
":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
":__support_fputil_multiply_add",
":__support_macros_config",
":__support_macros_optimization",
":llvm_libc_macros_float16_macros",
],
)
libc_support_library(
name = "__support_math_exp_float_constants",
hdrs = ["src/__support/math/exp_float_constants.h"],
@ -6769,6 +6784,13 @@ libc_math_function(
],
)
libc_math_function(
name = "erfcf16",
additional_deps = [
":__support_math_erfcf16",
],
)
libc_math_function(
name = "exp",
additional_deps = [