Skip to content

Commit 558975c

Browse files
bassiounixbassiounix
committed
[libc][math] Refactor acoshf16 implementation to header-only in src/__support/math folder.
1 parent 3c92206 commit 558975c

File tree

7 files changed

+195
-108
lines changed

7 files changed

+195
-108
lines changed

‎libc/shared/math.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@
1515
#include "math/acosf.h"
1616
#include "math/acosf16.h"
1717
#include "math/acoshf.h"
18+
#include "math/acoshf16.h"
1819
#include "math/erff.h"
1920
#include "math/exp.h"
2021
#include "math/exp10.h"

‎libc/shared/math/acoshf16.h

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
//===-- Shared acoshf16 function --------------------------------*- C++ -*-===//
2+
//
3+
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4+
// See https://llvm.org/LICENSE.txt for license information.
5+
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6+
//
7+
//===----------------------------------------------------------------------===//
8+
9+
#ifndef LLVM_LIBC_SHARED_MATH_ACOSHF16_H
10+
#define LLVM_LIBC_SHARED_MATH_ACOSHF16_H
11+
12+
#include "include/llvm-libc-macros/float16-macros.h"
13+
14+
#ifdef LIBC_TYPES_HAS_FLOAT16
15+
16+
#include "shared/libc_common.h"
17+
#include "src/__support/math/acoshf16.h"
18+
19+
namespace LIBC_NAMESPACE_DECL {
20+
namespace shared {
21+
22+
using math::acoshf16;
23+
24+
} // namespace shared
25+
} // namespace LIBC_NAMESPACE_DECL
26+
27+
#endif // LIBC_TYPES_HAS_FLOAT16
28+
29+
#endif // LLVM_LIBC_SHARED_MATH_ACOSHF16_H

‎libc/src/__support/math/CMakeLists.txt

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -79,6 +79,22 @@ add_header_library(
7979
libc.src.__support.macros.optimization
8080
)
8181

82+
add_header_library(
83+
acoshf16
84+
HDRS
85+
acoshf16.h
86+
DEPENDS
87+
.acoshf_utils
88+
libc.src.__support.FPUtil.cast
89+
libc.src.__support.FPUtil.except_value_utils
90+
libc.src.__support.FPUtil.fenv_impl
91+
libc.src.__support.FPUtil.fp_bits
92+
libc.src.__support.FPUtil.multiply_add
93+
libc.src.__support.FPUtil.polyeval
94+
libc.src.__support.FPUtil.sqrt
95+
libc.src.__support.macros.optimization
96+
)
97+
8298
add_header_library(
8399
asin_utils
84100
HDRS

‎libc/src/__support/math/acoshf16.h

Lines changed: 121 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,121 @@
1+
//===-- Implementation header for acoshf16 ----------------------*- C++ -*-===//
2+
//
3+
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4+
// See https://llvm.org/LICENSE.txt for license information.
5+
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6+
//
7+
//===----------------------------------------------------------------------===//
8+
9+
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_ACOSHF_H
10+
#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSHF_H
11+
12+
#include "include/llvm-libc-macros/float16-macros.h"
13+
14+
#ifdef LIBC_TYPES_HAS_FLOAT16
15+
16+
#include "acoshf_utils.h"
17+
#include "src/__support/FPUtil/FEnvImpl.h"
18+
#include "src/__support/FPUtil/FPBits.h"
19+
#include "src/__support/FPUtil/PolyEval.h"
20+
#include "src/__support/FPUtil/cast.h"
21+
#include "src/__support/FPUtil/except_value_utils.h"
22+
#include "src/__support/FPUtil/multiply_add.h"
23+
#include "src/__support/FPUtil/sqrt.h"
24+
#include "src/__support/macros/config.h"
25+
#include "src/__support/macros/optimization.h"
26+
27+
namespace LIBC_NAMESPACE_DECL {
28+
29+
namespace math {
30+
31+
static constexpr size_t N_EXCEPTS = 2;
32+
static constexpr fputil::ExceptValues<float16, N_EXCEPTS> ACOSHF16_EXCEPTS{{
33+
// (input, RZ output, RU offset, RD offset, RN offset)
34+
// x = 0x1.6dcp+1, acoshf16(x) = 0x1.b6p+0 (RZ)
35+
{0x41B7, 0x3ED8, 1, 0, 0},
36+
// x = 0x1.39p+0, acoshf16(x) = 0x1.4f8p-1 (RZ)
37+
{0x3CE4, 0x393E, 1, 0, 1},
38+
}};
39+
40+
static constexpr float16 acoshf16(float16 x) {
41+
using FPBits = fputil::FPBits<float16>;
42+
FPBits xbits(x);
43+
uint16_t x_u = xbits.uintval();
44+
45+
// Check for NaN input first.
46+
if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) {
47+
if (xbits.is_signaling_nan()) {
48+
fputil::raise_except_if_required(FE_INVALID);
49+
return FPBits::quiet_nan().get_val();
50+
}
51+
if (xbits.is_neg()) {
52+
fputil::set_errno_if_required(EDOM);
53+
fputil::raise_except_if_required(FE_INVALID);
54+
return FPBits::quiet_nan().get_val();
55+
}
56+
return x;
57+
}
58+
59+
// Domain error for inputs less than 1.0.
60+
if (LIBC_UNLIKELY(x <= 1.0f)) {
61+
if (x == 1.0f)
62+
return FPBits::zero().get_val();
63+
fputil::set_errno_if_required(EDOM);
64+
fputil::raise_except_if_required(FE_INVALID);
65+
return FPBits::quiet_nan().get_val();
66+
}
67+
68+
if (auto r = ACOSHF16_EXCEPTS.lookup(xbits.uintval());
69+
LIBC_UNLIKELY(r.has_value()))
70+
return r.value();
71+
72+
float xf = x;
73+
// High-precision polynomial approximation for inputs close to 1.0
74+
// ([1, 1.25)).
75+
//
76+
// Brief derivation:
77+
// 1. Expand acosh(1 + delta) using Taylor series around delta=0:
78+
// acosh(1 + delta) ≈ sqrt(2 * delta) * [1 - delta/12 + 3*delta^2/160
79+
// - 5*delta^3/896 + 35*delta^4/18432 + ...]
80+
// 2. Truncate the series to fit accurately for delta in [0, 0.25].
81+
// 3. Polynomial coefficients (from sollya) used here are:
82+
// P(delta) ≈ 1 - 0x1.555556p-4 * delta + 0x1.333334p-6 * delta^2
83+
// - 0x1.6db6dcp-8 * delta^3 + 0x1.f1c71cp-10 * delta^4
84+
// 4. The Sollya commands used to generate these coefficients were:
85+
// > display = hexadecimal;
86+
// > round(1/12, SG, RN);
87+
// > round(3/160, SG, RN);
88+
// > round(5/896, SG, RN);
89+
// > round(35/18432, SG, RN);
90+
// With hexadecimal display mode enabled, the outputs were:
91+
// 0x1.555556p-4
92+
// 0x1.333334p-6
93+
// 0x1.6db6dcp-8
94+
// 0x1.f1c71cp-10
95+
// 5. The maximum absolute error, estimated using:
96+
// dirtyinfnorm(acosh(1 + x) - sqrt(2*x) * P(x), [0, 0.25])
97+
// is:
98+
// 0x1.d84281p-22
99+
if (LIBC_UNLIKELY(x_u < 0x3D00U)) {
100+
float delta = xf - 1.0f;
101+
float sqrt_2_delta = fputil::sqrt<float>(2.0 * delta);
102+
float pe = fputil::polyeval(delta, 0x1p+0f, -0x1.555556p-4f, 0x1.333334p-6f,
103+
-0x1.6db6dcp-8f, 0x1.f1c71cp-10f);
104+
float approx = sqrt_2_delta * pe;
105+
return fputil::cast<float16>(approx);
106+
}
107+
108+
// acosh(x) = log(x + sqrt(x^2 - 1))
109+
float sqrt_term = fputil::sqrt<float>(fputil::multiply_add(xf, xf, -1.0f));
110+
float result = static_cast<float>(log_eval(xf + sqrt_term));
111+
112+
return fputil::cast<float16>(result);
113+
}
114+
115+
} // namespace math
116+
117+
} // namespace LIBC_NAMESPACE_DECL
118+
119+
#endif // LIBC_TYPES_HAS_FLOAT16
120+
121+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOSHF_H

‎libc/src/math/generic/CMakeLists.txt

Lines changed: 2 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -3880,18 +3880,8 @@ add_entrypoint_object(
38803880
HDRS
38813881
../acoshf16.h
38823882
DEPENDS
3883-
.explogxf
3884-
libc.hdr.errno_macros
3885-
libc.hdr.fenv_macros
3886-
libc.src.__support.FPUtil.cast
3887-
libc.src.__support.FPUtil.except_value_utils
3888-
libc.src.__support.FPUtil.fenv_impl
3889-
libc.src.__support.FPUtil.fp_bits
3890-
libc.src.__support.FPUtil.multiply_add
3891-
libc.src.__support.FPUtil.polyeval
3892-
libc.src.__support.FPUtil.sqrt
3893-
libc.src.__support.macros.optimization
3894-
libc.src.__support.macros.properties.types
3883+
libc.src.__support.math.acoshf16
3884+
libc.src.errno.errno
38953885
)
38963886

38973887
add_entrypoint_object(

‎libc/src/math/generic/acoshf16.cpp

Lines changed: 2 additions & 96 deletions
Original file line numberDiff line numberDiff line change
@@ -7,104 +7,10 @@
77
//===----------------------------------------------------------------------===//
88

99
#include "src/math/acoshf16.h"
10-
#include "explogxf.h"
11-
#include "hdr/errno_macros.h"
12-
#include "hdr/fenv_macros.h"
13-
#include "src/__support/FPUtil/FEnvImpl.h"
14-
#include "src/__support/FPUtil/FPBits.h"
15-
#include "src/__support/FPUtil/PolyEval.h"
16-
#include "src/__support/FPUtil/cast.h"
17-
#include "src/__support/FPUtil/except_value_utils.h"
18-
#include "src/__support/FPUtil/multiply_add.h"
19-
#include "src/__support/FPUtil/sqrt.h"
20-
#include "src/__support/common.h"
21-
#include "src/__support/macros/config.h"
22-
#include "src/__support/macros/optimization.h"
10+
#include "src/__support/math/acoshf16.h"
2311

2412
namespace LIBC_NAMESPACE_DECL {
2513

26-
static constexpr size_t N_EXCEPTS = 2;
27-
static constexpr fputil::ExceptValues<float16, N_EXCEPTS> ACOSHF16_EXCEPTS{{
28-
// (input, RZ output, RU offset, RD offset, RN offset)
29-
// x = 0x1.6dcp+1, acoshf16(x) = 0x1.b6p+0 (RZ)
30-
{0x41B7, 0x3ED8, 1, 0, 0},
31-
// x = 0x1.39p+0, acoshf16(x) = 0x1.4f8p-1 (RZ)
32-
{0x3CE4, 0x393E, 1, 0, 1},
33-
}};
34-
35-
LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) {
36-
using FPBits = fputil::FPBits<float16>;
37-
FPBits xbits(x);
38-
uint16_t x_u = xbits.uintval();
39-
40-
// Check for NaN input first.
41-
if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) {
42-
if (xbits.is_signaling_nan()) {
43-
fputil::raise_except_if_required(FE_INVALID);
44-
return FPBits::quiet_nan().get_val();
45-
}
46-
if (xbits.is_neg()) {
47-
fputil::set_errno_if_required(EDOM);
48-
fputil::raise_except_if_required(FE_INVALID);
49-
return FPBits::quiet_nan().get_val();
50-
}
51-
return x;
52-
}
53-
54-
// Domain error for inputs less than 1.0.
55-
if (LIBC_UNLIKELY(x <= 1.0f)) {
56-
if (x == 1.0f)
57-
return FPBits::zero().get_val();
58-
fputil::set_errno_if_required(EDOM);
59-
fputil::raise_except_if_required(FE_INVALID);
60-
return FPBits::quiet_nan().get_val();
61-
}
62-
63-
if (auto r = ACOSHF16_EXCEPTS.lookup(xbits.uintval());
64-
LIBC_UNLIKELY(r.has_value()))
65-
return r.value();
66-
67-
float xf = x;
68-
// High-precision polynomial approximation for inputs close to 1.0
69-
// ([1, 1.25)).
70-
//
71-
// Brief derivation:
72-
// 1. Expand acosh(1 + delta) using Taylor series around delta=0:
73-
// acosh(1 + delta) ≈ sqrt(2 * delta) * [1 - delta/12 + 3*delta^2/160
74-
// - 5*delta^3/896 + 35*delta^4/18432 + ...]
75-
// 2. Truncate the series to fit accurately for delta in [0, 0.25].
76-
// 3. Polynomial coefficients (from sollya) used here are:
77-
// P(delta) ≈ 1 - 0x1.555556p-4 * delta + 0x1.333334p-6 * delta^2
78-
// - 0x1.6db6dcp-8 * delta^3 + 0x1.f1c71cp-10 * delta^4
79-
// 4. The Sollya commands used to generate these coefficients were:
80-
// > display = hexadecimal;
81-
// > round(1/12, SG, RN);
82-
// > round(3/160, SG, RN);
83-
// > round(5/896, SG, RN);
84-
// > round(35/18432, SG, RN);
85-
// With hexadecimal display mode enabled, the outputs were:
86-
// 0x1.555556p-4
87-
// 0x1.333334p-6
88-
// 0x1.6db6dcp-8
89-
// 0x1.f1c71cp-10
90-
// 5. The maximum absolute error, estimated using:
91-
// dirtyinfnorm(acosh(1 + x) - sqrt(2*x) * P(x), [0, 0.25])
92-
// is:
93-
// 0x1.d84281p-22
94-
if (LIBC_UNLIKELY(x_u < 0x3D00U)) {
95-
float delta = xf - 1.0f;
96-
float sqrt_2_delta = fputil::sqrt<float>(2.0 * delta);
97-
float pe = fputil::polyeval(delta, 0x1p+0f, -0x1.555556p-4f, 0x1.333334p-6f,
98-
-0x1.6db6dcp-8f, 0x1.f1c71cp-10f);
99-
float approx = sqrt_2_delta * pe;
100-
return fputil::cast<float16>(approx);
101-
}
102-
103-
// acosh(x) = log(x + sqrt(x^2 - 1))
104-
float sqrt_term = fputil::sqrt<float>(fputil::multiply_add(xf, xf, -1.0f));
105-
float result = static_cast<float>(log_eval(xf + sqrt_term));
106-
107-
return fputil::cast<float16>(result);
108-
}
14+
LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { return math::acoshf16(x); }
10915

11016
} // namespace LIBC_NAMESPACE_DECL

‎utils/bazel/llvm-project-overlay/libc/BUILD.bazel

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2142,6 +2142,22 @@ libc_support_library(
21422142
],
21432143
)
21442144

2145+
libc_support_library(
2146+
name = "__support_math_acoshf16",
2147+
hdrs = ["src/__support/math/acoshf16.h"],
2148+
deps = [
2149+
":__support_math_acoshf_utils",
2150+
":__support_fputil_cast",
2151+
":__support_fputil_except_value_utils",
2152+
":__support_fputil_fenv_impl",
2153+
":__support_fputil_fp_bits",
2154+
":__support_fputil_multiply_add",
2155+
":__support_fputil_polyeval",
2156+
":__support_fputil_sqrt",
2157+
":__support_macros_optimization",
2158+
],
2159+
)
2160+
21452161
libc_support_library(
21462162
name = "__support_math_asin_utils",
21472163
hdrs = ["src/__support/math/asin_utils.h"],
@@ -2683,6 +2699,14 @@ libc_math_function(
26832699
],
26842700
)
26852701

2702+
libc_math_function(
2703+
name = "acoshf16",
2704+
additional_deps = [
2705+
":__support_math_acoshf16",
2706+
":errno",
2707+
],
2708+
)
2709+
26862710
libc_math_function(
26872711
name = "asinf",
26882712
additional_deps = [

0 commit comments

Comments
 (0)