Skip to content

[libc][math][c23] Add atanhf16 C23 math function. #132612

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 17 commits into from
Apr 25, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions libc/config/linux/x86_64/entrypoints.txt
Original file line number Diff line number Diff line change
Expand Up @@ -664,6 +664,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.acoshf16
libc.src.math.asinf16
libc.src.math.asinhf16
libc.src.math.atanhf16
libc.src.math.canonicalizef16
libc.src.math.ceilf16
libc.src.math.copysignf16
Expand Down
2 changes: 1 addition & 1 deletion libc/docs/headers/math/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -267,7 +267,7 @@ Higher Math Functions
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| atan2pi | | | | | | 7.12.4.11 | F.10.1.11 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| atanh | |check| | | | | | 7.12.5.3 | F.10.2.3 |
| atanh | |check| | | | |check| | | 7.12.5.3 | F.10.2.3 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| atanpi | | | | | | 7.12.4.10 | F.10.1.10 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
Expand Down
7 changes: 7 additions & 0 deletions libc/include/math.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,13 @@ functions:
return_type: float
arguments:
- type: float
- name: atanhf16
standards:
- stdc
return_type: _Float16
arguments:
- type: _Float16
guard: LIBC_TYPES_HAS_FLOAT16
- name: canonicalize
standards:
- stdc
Expand Down
1 change: 1 addition & 0 deletions libc/src/math/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,7 @@ add_math_entrypoint_object(atan2f128)

add_math_entrypoint_object(atanh)
add_math_entrypoint_object(atanhf)
add_math_entrypoint_object(atanhf16)

add_math_entrypoint_object(canonicalize)
add_math_entrypoint_object(canonicalizef)
Expand Down
21 changes: 21 additions & 0 deletions libc/src/math/atanhf16.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
//===-- Implementation header for atanhf16 ----------------------*- 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_ATANHF16_H
#define LLVM_LIBC_SRC_MATH_ATANHF16_H

#include "src/__support/macros/config.h"
#include "src/__support/macros/properties/types.h"

namespace LIBC_NAMESPACE_DECL {

float16 atanhf16(float16 x);

} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC_MATH_ATANHF16_H
20 changes: 20 additions & 0 deletions libc/src/math/generic/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -4029,6 +4029,26 @@ add_entrypoint_object(
libc.src.__support.macros.optimization
)

add_entrypoint_object(
atanhf16
SRCS
atanhf16.cpp
HDRS
../atanhf16.h
DEPENDS
.explogxf
libc.hdr.errno_macros
libc.hdr.fenv_macros
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.FPUtil.polyeval
libc.src.__support.macros.optimization
libc.src.__support.macros.properties.types
)

add_object_library(
inv_trigf_utils
HDRS
Expand Down
98 changes: 98 additions & 0 deletions libc/src/math/generic/atanhf16.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,98 @@
//===-- Half-precision atanh(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/atanhf16.h"
#include "explogxf.h"
#include "hdr/errno_macros.h"
#include "hdr/fenv_macros.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.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/common.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"

namespace LIBC_NAMESPACE_DECL {

static constexpr size_t N_EXCEPTS = 1;
static constexpr fputil::ExceptValues<float16, N_EXCEPTS> ATANHF16_EXCEPTS{{
// (input, RZ output, RU offset, RD offset, RN offset)
// x = 0x1.a5cp-4, atanhf16(x) = 0x1.a74p-4 (RZ)
{0x2E97, 0x2E9D, 1, 0, 0},
}};

LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
using FPBits = fputil::FPBits<float16>;

FPBits xbits(x);
Sign sign = xbits.sign();
uint16_t x_abs = xbits.abs().uintval();

// |x| >= 1
if (LIBC_UNLIKELY(x_abs >= 0x3c00U)) {
if (xbits.is_nan()) {
if (xbits.is_signaling_nan()) {
fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}
return x;
}

// |x| == 1.0
if (x_abs == 0x3c00U) {
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_DIVBYZERO);
return FPBits::inf(sign).get_val();
}
// |x| > 1.0
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}

if (auto r = ATANHF16_EXCEPTS.lookup(xbits.uintval());
LIBC_UNLIKELY(r.has_value()))
return r.value();

// For |x| less than approximately 0.24
if (LIBC_UNLIKELY(x_abs <= 0x33f3U)) {
// atanh(+/-0) = +/-0
if (LIBC_UNLIKELY(x_abs == 0U))
return x;
// The Taylor expansion of atanh(x) is:
// atanh(x) = x + x^3/3 + x^5/5 + x^7/7 + x^9/9 + x^11/11
// = x * [1 + x^2/3 + x^4/5 + x^6/7 + x^8/9 + x^10/11]
// When |x| < 2^-5 (0x0800U), this can be approximated by:
// atanh(x) ≈ x + (1/3)*x^3
if (LIBC_UNLIKELY(x_abs < 0x0800U)) {
float xf = x;
return fputil::cast<float16>(xf + 0x1.555556p-2f * xf * xf * xf);
}

// For 2^-5 <= |x| <= 0x1.fccp-3 (~0.24):
// Let t = x^2.
// Define P(t) ≈ (1/3)*t + (1/5)*t^2 + (1/7)*t^3 + (1/9)*t^4 + (1/11)*t^5.
// Coefficients (from Sollya, RN, hexadecimal):
// 1/3 = 0x1.555556p-2, 1/5 = 0x1.99999ap-3, 1/7 = 0x1.24924ap-3,
// 1/9 = 0x1.c71c72p-4, 1/11 = 0x1.745d18p-4
// Thus, atanh(x) ≈ x * (1 + P(x^2)).
float xf = x;
float x2 = xf * xf;
float pe = fputil::polyeval(x2, 0.0f, 0x1.555556p-2f, 0x1.99999ap-3f,
0x1.24924ap-3f, 0x1.c71c72p-4f, 0x1.745d18p-4f);
return fputil::cast<float16>(fputil::multiply_add(xf, pe, xf));
}

float xf = x;
return fputil::cast<float16>(0.5 * log_eval_f((xf + 1.0f) / (xf - 1.0f)));
}

} // namespace LIBC_NAMESPACE_DECL
78 changes: 78 additions & 0 deletions libc/src/math/generic/common_constants.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,45 @@

namespace LIBC_NAMESPACE_DECL {

// Lookup table for logf(f) = logf(1 + n*2^(-7)) where n = 0..127,
// computed and stored as float precision constants.
// Generated by Sollya with the following commands:
// display = hexadecimal;
// for n from 0 to 127 do { print(single(1 / (1 + n / 128.0))); };
const float ONE_OVER_F_FLOAT[128] = {
0x1p0f, 0x1.fc07fp-1f, 0x1.f81f82p-1f, 0x1.f4465ap-1f,
0x1.f07c2p-1f, 0x1.ecc07cp-1f, 0x1.e9131ap-1f, 0x1.e573acp-1f,
0x1.e1e1e2p-1f, 0x1.de5d6ep-1f, 0x1.dae608p-1f, 0x1.d77b66p-1f,
0x1.d41d42p-1f, 0x1.d0cb58p-1f, 0x1.cd8568p-1f, 0x1.ca4b3p-1f,
0x1.c71c72p-1f, 0x1.c3f8fp-1f, 0x1.c0e07p-1f, 0x1.bdd2b8p-1f,
0x1.bacf92p-1f, 0x1.b7d6c4p-1f, 0x1.b4e81cp-1f, 0x1.b20364p-1f,
0x1.af286cp-1f, 0x1.ac5702p-1f, 0x1.a98ef6p-1f, 0x1.a6d01ap-1f,
0x1.a41a42p-1f, 0x1.a16d4p-1f, 0x1.9ec8eap-1f, 0x1.9c2d14p-1f,
0x1.99999ap-1f, 0x1.970e5p-1f, 0x1.948b1p-1f, 0x1.920fb4p-1f,
0x1.8f9c18p-1f, 0x1.8d3018p-1f, 0x1.8acb9p-1f, 0x1.886e6p-1f,
0x1.861862p-1f, 0x1.83c978p-1f, 0x1.818182p-1f, 0x1.7f406p-1f,
0x1.7d05f4p-1f, 0x1.7ad22p-1f, 0x1.78a4c8p-1f, 0x1.767dcep-1f,
0x1.745d18p-1f, 0x1.724288p-1f, 0x1.702e06p-1f, 0x1.6e1f76p-1f,
0x1.6c16c2p-1f, 0x1.6a13cep-1f, 0x1.681682p-1f, 0x1.661ec6p-1f,
0x1.642c86p-1f, 0x1.623fa8p-1f, 0x1.605816p-1f, 0x1.5e75bcp-1f,
0x1.5c9882p-1f, 0x1.5ac056p-1f, 0x1.58ed24p-1f, 0x1.571ed4p-1f,
0x1.555556p-1f, 0x1.539094p-1f, 0x1.51d07ep-1f, 0x1.501502p-1f,
0x1.4e5e0ap-1f, 0x1.4cab88p-1f, 0x1.4afd6ap-1f, 0x1.49539ep-1f,
0x1.47ae14p-1f, 0x1.460cbcp-1f, 0x1.446f86p-1f, 0x1.42d662p-1f,
0x1.414142p-1f, 0x1.3fb014p-1f, 0x1.3e22ccp-1f, 0x1.3c995ap-1f,
0x1.3b13b2p-1f, 0x1.3991c2p-1f, 0x1.381382p-1f, 0x1.3698ep-1f,
0x1.3521dp-1f, 0x1.33ae46p-1f, 0x1.323e34p-1f, 0x1.30d19p-1f,
0x1.2f684cp-1f, 0x1.2e025cp-1f, 0x1.2c9fb4p-1f, 0x1.2b404ap-1f,
0x1.29e412p-1f, 0x1.288b02p-1f, 0x1.27350cp-1f, 0x1.25e228p-1f,
0x1.24924ap-1f, 0x1.234568p-1f, 0x1.21fb78p-1f, 0x1.20b47p-1f,
0x1.1f7048p-1f, 0x1.1e2ef4p-1f, 0x1.1cf06ap-1f, 0x1.1bb4a4p-1f,
0x1.1a7b96p-1f, 0x1.194538p-1f, 0x1.181182p-1f, 0x1.16e068p-1f,
0x1.15b1e6p-1f, 0x1.1485fp-1f, 0x1.135c82p-1f, 0x1.12358ep-1f,
0x1.111112p-1f, 0x1.0fef02p-1f, 0x1.0ecf56p-1f, 0x1.0db20ap-1f,
0x1.0c9714p-1f, 0x1.0b7e6ep-1f, 0x1.0a681p-1f, 0x1.0953f4p-1f,
0x1.08421p-1f, 0x1.07326p-1f, 0x1.0624dep-1f, 0x1.05198p-1f,
0x1.041042p-1f, 0x1.03091cp-1f, 0x1.020408p-1f, 0x1.010102p-1f};

// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127.
const double ONE_OVER_F[128] = {
0x1.0000000000000p+0, 0x1.fc07f01fc07f0p-1, 0x1.f81f81f81f820p-1,
Expand Down Expand Up @@ -59,6 +98,45 @@ const double ONE_OVER_F[128] = {
0x1.05197f7d73404p-1, 0x1.0410410410410p-1, 0x1.03091b51f5e1ap-1,
0x1.0204081020408p-1, 0x1.0101010101010p-1};

// Lookup table for log(f) = log(1 + n*2^(-7)) where n = 0..127,
// computed and stored as float precision constants.
// Generated by Sollya with the following commands:
// display = hexadecimal;
// for n from 0 to 127 do { print(single(log(1 + n / 128.0))); };
const float LOG_F_FLOAT[128] = {
0.0f, 0x1.fe02a6p-8f, 0x1.fc0a8cp-7f, 0x1.7b91bp-6f,
0x1.f829bp-6f, 0x1.39e87cp-5f, 0x1.77459p-5f, 0x1.b42dd8p-5f,
0x1.f0a30cp-5f, 0x1.16536ep-4f, 0x1.341d7ap-4f, 0x1.51b074p-4f,
0x1.6f0d28p-4f, 0x1.8c345ep-4f, 0x1.a926d4p-4f, 0x1.c5e548p-4f,
0x1.e27076p-4f, 0x1.fec914p-4f, 0x1.0d77e8p-3f, 0x1.1b72aep-3f,
0x1.29553p-3f, 0x1.371fc2p-3f, 0x1.44d2b6p-3f, 0x1.526e5ep-3f,
0x1.5ff308p-3f, 0x1.6d60fep-3f, 0x1.7ab89p-3f, 0x1.87fa06p-3f,
0x1.9525aap-3f, 0x1.a23bc2p-3f, 0x1.af3c94p-3f, 0x1.bc2868p-3f,
0x1.c8ff7cp-3f, 0x1.d5c216p-3f, 0x1.e27076p-3f, 0x1.ef0adcp-3f,
0x1.fb9186p-3f, 0x1.04025ap-2f, 0x1.0a324ep-2f, 0x1.1058cp-2f,
0x1.1675cap-2f, 0x1.1c898cp-2f, 0x1.22942p-2f, 0x1.2895a2p-2f,
0x1.2e8e2cp-2f, 0x1.347ddap-2f, 0x1.3a64c6p-2f, 0x1.404308p-2f,
0x1.4618bcp-2f, 0x1.4be5fap-2f, 0x1.51aad8p-2f, 0x1.576772p-2f,
0x1.5d1bdcp-2f, 0x1.62c83p-2f, 0x1.686c82p-2f, 0x1.6e08eap-2f,
0x1.739d8p-2f, 0x1.792a56p-2f, 0x1.7eaf84p-2f, 0x1.842d1ep-2f,
0x1.89a338p-2f, 0x1.8f11e8p-2f, 0x1.947942p-2f, 0x1.99d958p-2f,
0x1.9f323ep-2f, 0x1.a4840ap-2f, 0x1.a9cecap-2f, 0x1.af1294p-2f,
0x1.b44f78p-2f, 0x1.b9858ap-2f, 0x1.beb4dap-2f, 0x1.c3dd7ap-2f,
0x1.c8ff7cp-2f, 0x1.ce1afp-2f, 0x1.d32fe8p-2f, 0x1.d83e72p-2f,
0x1.dd46ap-2f, 0x1.e24882p-2f, 0x1.e74426p-2f, 0x1.ec399ep-2f,
0x1.f128f6p-2f, 0x1.f6124p-2f, 0x1.faf588p-2f, 0x1.ffd2ep-2f,
0x1.02552ap-1f, 0x1.04bdfap-1f, 0x1.0723e6p-1f, 0x1.0986f4p-1f,
0x1.0be72ep-1f, 0x1.0e4498p-1f, 0x1.109f3ap-1f, 0x1.12f71ap-1f,
0x1.154c3ep-1f, 0x1.179eacp-1f, 0x1.19ee6cp-1f, 0x1.1c3b82p-1f,
0x1.1e85f6p-1f, 0x1.20cdcep-1f, 0x1.23130ep-1f, 0x1.2555bcp-1f,
0x1.2795e2p-1f, 0x1.29d38p-1f, 0x1.2c0e9ep-1f, 0x1.2e4744p-1f,
0x1.307d74p-1f, 0x1.32b134p-1f, 0x1.34e28ap-1f, 0x1.37117cp-1f,
0x1.393e0ep-1f, 0x1.3b6844p-1f, 0x1.3d9026p-1f, 0x1.3fb5b8p-1f,
0x1.41d8fep-1f, 0x1.43f9fep-1f, 0x1.4618bcp-1f, 0x1.48353ep-1f,
0x1.4a4f86p-1f, 0x1.4c679ap-1f, 0x1.4e7d82p-1f, 0x1.50913cp-1f,
0x1.52a2d2p-1f, 0x1.54b246p-1f, 0x1.56bf9ep-1f, 0x1.58cadcp-1f,
0x1.5ad404p-1f, 0x1.5cdb1ep-1f, 0x1.5ee02ap-1f, 0x1.60e33p-1f};

// Lookup table for log(f) = log(1 + n*2^(-7)) where n = 0..127.
const double LOG_F[128] = {
0x0.0000000000000p+0, 0x1.fe02a6b106788p-8, 0x1.fc0a8b0fc03e3p-7,
Expand Down
8 changes: 8 additions & 0 deletions libc/src/math/generic/common_constants.h
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,17 @@

namespace LIBC_NAMESPACE_DECL {

// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127,
// computed and stored as float precision constants.
extern const float ONE_OVER_F_FLOAT[128];

// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127.
extern const double ONE_OVER_F[128];

// Lookup table for log(f) = log(1 + n*2^(-7)) where n = 0..127,
// computed and stored as float precision constants.
extern const float LOG_F_FLOAT[128];

// Lookup table for log(f) = log(1 + n*2^(-7)) where n = 0..127.
extern const double LOG_F[128];

Expand Down
43 changes: 43 additions & 0 deletions libc/src/math/generic/explogxf.h
Original file line number Diff line number Diff line change
Expand Up @@ -297,6 +297,49 @@ LIBC_INLINE static double log2_eval(double x) {
return result;
}

// x should be positive, normal finite value
// TODO: Simplify range reduction and polynomial degree for float16.
// See issue #137190.
LIBC_INLINE static float log_eval_f(float x) {
// For x = 2^ex * (1 + mx), logf(x) = ex * logf(2) + logf(1 + mx).
using FPBits = fputil::FPBits<float>;
FPBits xbits(x);

float ex = static_cast<float>(xbits.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>(xbits.get_mantissa() >> (FPBits::FRACTION_LEN - 7));

// Set bits to (1 + (mx - p1*2^(-7)))
xbits.set_uintval(xbits.uintval() & (FPBits::FRACTION_MASK >> 7));
xbits.set_biased_exponent(FPBits::EXP_BIAS);
// dx = (mx - p1*2^(-7)) / (1 + p1*2^(-7)).
float dx = (xbits.get_val() - 1.0f) * ONE_OVER_F_FLOAT[p1];

// Minimax polynomial for log(1 + dx), generated using Sollya:
// > P = fpminimax(log(1 + x)/x, 6, [|SG...|], [0, 2^-7]);
// > Q = (P - 1) / x;
// > for i from 0 to degree(Q) do print(coeff(Q, i));
constexpr float COEFFS[6] = {-0x1p-1f, 0x1.555556p-2f, -0x1.00022ep-2f,
0x1.9ea056p-3f, -0x1.e50324p-2f, 0x1.c018fp3f};

float dx2 = dx * dx;

float c1 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]);
float c2 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]);
float c3 = fputil::multiply_add(dx, COEFFS[5], COEFFS[4]);

float p = fputil::polyeval(dx2, dx, c1, c2, c3);

// Generated by Sollya with the following commands:
// > display = hexadecimal;
// > round(log(2), SG, RN);
constexpr float LOGF_2 = 0x1.62e43p-1f;

float result = fputil::multiply_add(ex, LOGF_2, LOG_F_FLOAT[p1] + p);
return result;
}

// x should be positive, normal finite value
LIBC_INLINE static double log_eval(double x) {
// For x = 2^ex * (1 + mx)
Expand Down
11 changes: 11 additions & 0 deletions libc/test/src/math/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -2143,6 +2143,17 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)

add_fp_unittest(
atanhf16_test
NEED_MPFR
SUITE
libc-math-unittests
SRCS
atanhf16_test.cpp
DEPENDS
libc.src.math.atanhf16
)

add_fp_unittest(
fmul_test
NEED_MPFR
Expand Down
Loading
Loading