Skip to content

Commit 9102fa2

Browse files
overmightyovermighty
committed
[libc][math][c23] Add sinhf16 and coshf16 C23 math functions
Part of #95250.
1 parent 3f11a3c commit 9102fa2

File tree

17 files changed

+756
-2
lines changed

17 files changed

+756
-2
lines changed

‎libc/config/gpu/entrypoints.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -485,6 +485,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
485485
libc.src.math.canonicalizef16
486486
libc.src.math.ceilf16
487487
libc.src.math.copysignf16
488+
libc.src.math.coshf16
488489
libc.src.math.exp10f16
489490
libc.src.math.exp10m1f16
490491
libc.src.math.exp2f16
@@ -549,6 +550,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
549550
libc.src.math.scalbnf16
550551
libc.src.math.setpayloadf16
551552
libc.src.math.setpayloadsigf16
553+
libc.src.math.sinhf16
552554
libc.src.math.totalorderf16
553555
libc.src.math.totalordermagf16
554556
libc.src.math.truncf16

‎libc/config/linux/x86_64/entrypoints.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -594,6 +594,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
594594
libc.src.math.canonicalizef16
595595
libc.src.math.ceilf16
596596
libc.src.math.copysignf16
597+
libc.src.math.coshf16
597598
libc.src.math.exp10f16
598599
libc.src.math.exp10m1f16
599600
libc.src.math.exp2f16
@@ -660,6 +661,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
660661
libc.src.math.scalbnf16
661662
libc.src.math.setpayloadf16
662663
libc.src.math.setpayloadsigf16
664+
libc.src.math.sinhf16
663665
libc.src.math.totalorderf16
664666
libc.src.math.totalordermagf16
665667
libc.src.math.truncf16

‎libc/docs/math/index.rst‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -274,7 +274,7 @@ Higher Math Functions
274274
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
275275
| cos | |check| | |check| | | | | 7.12.4.5 | F.10.1.5 |
276276
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
277-
| cosh | |check| | | | | | 7.12.5.4 | F.10.2.4 |
277+
| cosh | |check| | | | |check| | | 7.12.5.4 | F.10.2.4 |
278278
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
279279
| cospi | |check| | | | | | 7.12.4.12 | F.10.1.12 |
280280
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
@@ -336,7 +336,7 @@ Higher Math Functions
336336
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
337337
| sincos | |check| | |check| | | | | | |
338338
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
339-
| sinh | |check| | | | | | 7.12.5.5 | F.10.2.5 |
339+
| sinh | |check| | | | |check| | | 7.12.5.5 | F.10.2.5 |
340340
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
341341
| sinpi | |check| | | | | | 7.12.4.13 | F.10.1.13 |
342342
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+

‎libc/spec/stdc.td‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -706,7 +706,11 @@ def StdC : StandardSpec<"stdc"> {
706706
FunctionSpec<"pow", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
707707

708708
FunctionSpec<"coshf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
709+
GuardedFunctionSpec<"coshf16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
710+
709711
FunctionSpec<"sinhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
712+
GuardedFunctionSpec<"sinhf16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
713+
710714
FunctionSpec<"tanhf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
711715

712716
FunctionSpec<"acosf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,

‎libc/src/math/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -83,8 +83,11 @@ add_math_entrypoint_object(copysignf128)
8383

8484
add_math_entrypoint_object(cos)
8585
add_math_entrypoint_object(cosf)
86+
8687
add_math_entrypoint_object(cosh)
8788
add_math_entrypoint_object(coshf)
89+
add_math_entrypoint_object(coshf16)
90+
8891
add_math_entrypoint_object(cospif)
8992

9093
add_math_entrypoint_object(daddl)
@@ -468,6 +471,7 @@ add_math_entrypoint_object(sinpif)
468471

469472
add_math_entrypoint_object(sinh)
470473
add_math_entrypoint_object(sinhf)
474+
add_math_entrypoint_object(sinhf16)
471475

472476
add_math_entrypoint_object(sqrt)
473477
add_math_entrypoint_object(sqrtf)

‎libc/src/math/coshf16.h‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
//===-- Implementation header for coshf16 -----------------------*- 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_MATH_COSHF16_H
10+
#define LLVM_LIBC_SRC_MATH_COSHF16_H
11+
12+
#include "src/__support/macros/config.h"
13+
#include "src/__support/macros/properties/types.h"
14+
15+
namespace LIBC_NAMESPACE_DECL {
16+
17+
float16 coshf16(float16 x);
18+
19+
} // namespace LIBC_NAMESPACE_DECL
20+
21+
#endif // LLVM_LIBC_SRC_MATH_COSHF16_H

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

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4084,6 +4084,25 @@ add_entrypoint_object(
40844084
-O3
40854085
)
40864086

4087+
add_entrypoint_object(
4088+
coshf16
4089+
SRCS
4090+
coshf16.cpp
4091+
HDRS
4092+
../coshf16.h
4093+
DEPENDS
4094+
.expxf16
4095+
libc.hdr.errno_macros
4096+
libc.hdr.fenv_macros
4097+
libc.src.__support.FPUtil.except_value_utils
4098+
libc.src.__support.FPUtil.fenv_impl
4099+
libc.src.__support.FPUtil.fp_bits
4100+
libc.src.__support.FPUtil.rounding_mode
4101+
libc.src.__support.macros.optimization
4102+
COMPILE_OPTIONS
4103+
-O3
4104+
)
4105+
40874106
add_entrypoint_object(
40884107
sinhf
40894108
SRCS
@@ -4099,6 +4118,25 @@ add_entrypoint_object(
40994118
-O3
41004119
)
41014120

4121+
add_entrypoint_object(
4122+
sinhf16
4123+
SRCS
4124+
sinhf16.cpp
4125+
HDRS
4126+
../sinhf16.h
4127+
DEPENDS
4128+
.expxf16
4129+
libc.hdr.errno_macros
4130+
libc.hdr.fenv_macros
4131+
libc.src.__support.FPUtil.except_value_utils
4132+
libc.src.__support.FPUtil.fenv_impl
4133+
libc.src.__support.FPUtil.fp_bits
4134+
libc.src.__support.FPUtil.rounding_mode
4135+
libc.src.__support.macros.optimization
4136+
COMPILE_OPTIONS
4137+
-O3
4138+
)
4139+
41024140
add_entrypoint_object(
41034141
tanhf
41044142
SRCS

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

Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,103 @@
1+
//===-- Half-precision cosh(x) function -----------------------------------===//
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+
#include "src/math/coshf16.h"
10+
#include "expxf16.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/except_value_utils.h"
16+
#include "src/__support/FPUtil/rounding_mode.h"
17+
#include "src/__support/common.h"
18+
#include "src/__support/macros/config.h"
19+
#include "src/__support/macros/optimization.h"
20+
21+
namespace LIBC_NAMESPACE_DECL {
22+
23+
static constexpr fputil::ExceptValues<float16, 9> COSHF16_EXCEPTS_POS = {{
24+
// x = 0x1.6ap-5, coshf16(x) = 0x1p+0 (RZ)
25+
{0x29a8U, 0x3c00U, 1U, 0U, 1U},
26+
// x = 0x1.8c4p+0, coshf16(x) = 0x1.3a8p+1 (RZ)
27+
{0x3e31U, 0x40eaU, 1U, 0U, 0U},
28+
// x = 0x1.994p+0, coshf16(x) = 0x1.498p+1 (RZ)
29+
{0x3e65U, 0x4126U, 1U, 0U, 0U},
30+
// x = 0x1.b6p+0, coshf16(x) = 0x1.6d8p+1 (RZ)
31+
{0x3ed8U, 0x41b6U, 1U, 0U, 1U},
32+
// x = 0x1.aap+1, coshf16(x) = 0x1.be8p+3 (RZ)
33+
{0x42a8U, 0x4afaU, 1U, 0U, 1U},
34+
// x = 0x1.cc4p+1, coshf16(x) = 0x1.23cp+4 (RZ)
35+
{0x4331U, 0x4c8fU, 1U, 0U, 0U},
36+
// x = 0x1.288p+2, coshf16(x) = 0x1.9b4p+5 (RZ)
37+
{0x44a2U, 0x526dU, 1U, 0U, 0U},
38+
// x = 0x1.958p+2, coshf16(x) = 0x1.1a4p+8 (RZ)
39+
{0x4656U, 0x5c69U, 1U, 0U, 0U},
40+
// x = 0x1.5fp+3, coshf16(x) = 0x1.c54p+14 (RZ)
41+
{0x497cU, 0x7715U, 1U, 0U, 1U},
42+
}};
43+
44+
static constexpr fputil::ExceptValues<float16, 4> COSHF16_EXCEPTS_NEG = {{
45+
// x = -0x1.6ap-5, coshf16(x) = 0x1p+0 (RZ)
46+
{0xa9a8U, 0x3c00U, 1U, 0U, 1U},
47+
// x = -0x1.b6p+0, coshf16(x) = 0x1.6d8p+1 (RZ)
48+
{0xbed8U, 0x41b6U, 1U, 0U, 1U},
49+
// x = -0x1.288p+2, coshf16(x) = 0x1.9b4p+5 (RZ)
50+
{0xc4a2U, 0x526dU, 1U, 0U, 0U},
51+
// x = -0x1.5fp+3, coshf16(x) = 0x1.c54p+14 (RZ)
52+
{0xc97cU, 0x7715U, 1U, 0U, 1U},
53+
}};
54+
55+
LLVM_LIBC_FUNCTION(float16, coshf16, (float16 x)) {
56+
using FPBits = fputil::FPBits<float16>;
57+
FPBits x_bits(x);
58+
59+
uint16_t x_u = x_bits.uintval();
60+
uint16_t x_abs = x_u & 0x7fffU;
61+
62+
// When |x| >= acosh(2^16), or x is NaN.
63+
if (LIBC_UNLIKELY(x_abs >= 0x49e5U)) {
64+
// cosh(NaN) = NaN
65+
if (x_bits.is_nan()) {
66+
if (x_bits.is_signaling_nan()) {
67+
fputil::raise_except_if_required(FE_INVALID);
68+
return FPBits::quiet_nan().get_val();
69+
}
70+
71+
return x;
72+
}
73+
74+
// When |x| >= acosh(2^16).
75+
if (x_abs >= 0x49e5U) {
76+
// cosh(+/-inf) = +inf
77+
if (x_bits.is_inf())
78+
return FPBits::inf().get_val();
79+
80+
switch (fputil::quick_get_round()) {
81+
case FE_TONEAREST:
82+
case FE_UPWARD:
83+
fputil::set_errno_if_required(ERANGE);
84+
fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
85+
return FPBits::inf().get_val();
86+
default:
87+
return FPBits::max_normal().get_val();
88+
}
89+
}
90+
}
91+
92+
if (x_bits.is_pos()) {
93+
if (auto r = COSHF16_EXCEPTS_POS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
94+
return r.value();
95+
} else {
96+
if (auto r = COSHF16_EXCEPTS_NEG.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
97+
return r.value();
98+
}
99+
100+
return eval_sinh_or_cosh</*IsSinh=*/false>(x);
101+
}
102+
103+
} // namespace LIBC_NAMESPACE_DECL

‎libc/src/math/generic/expxf16.h‎

Lines changed: 113 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -174,6 +174,119 @@ LIBC_INLINE ExpRangeReduction exp10_range_reduction(float16 x) {
174174
return {exp2_hi_mid, exp10_lo};
175175
}
176176

177+
// Generated by Sollya with the following commands:
178+
// > display = hexadecimal;
179+
// > round(log2(exp(1)), SG, RN);
180+
static constexpr float LOG2F_E = 0x1.715476p+0f;
181+
182+
// Generated by Sollya with the following commands:
183+
// > display = hexadecimal;
184+
// > round(log(2), SG, RN);
185+
static constexpr float LOGF_2 = 0x1.62e43p-1f;
186+
187+
// Generated by Sollya with the following commands:
188+
// > display = hexadecimal;
189+
// > for i from 0 to 31 do printsingle(round(2^(i * 2^-5), SG, RN));
190+
static constexpr cpp::array<uint32_t, 32> EXP2_MID_5_BITS = {
191+
0x3f80'0000U, 0x3f82'cd87U, 0x3f85'aac3U, 0x3f88'980fU, 0x3f8b'95c2U,
192+
0x3f8e'a43aU, 0x3f91'c3d3U, 0x3f94'f4f0U, 0x3f98'37f0U, 0x3f9b'8d3aU,
193+
0x3f9e'f532U, 0x3fa2'7043U, 0x3fa5'fed7U, 0x3fa9'a15bU, 0x3fad'583fU,
194+
0x3fb1'23f6U, 0x3fb5'04f3U, 0x3fb8'fbafU, 0x3fbd'08a4U, 0x3fc1'2c4dU,
195+
0x3fc5'672aU, 0x3fc9'b9beU, 0x3fce'248cU, 0x3fd2'a81eU, 0x3fd7'44fdU,
196+
0x3fdb'fbb8U, 0x3fe0'ccdfU, 0x3fe5'b907U, 0x3fea'c0c7U, 0x3fef'e4baU,
197+
0x3ff5'257dU, 0x3ffa'83b3U,
198+
};
199+
200+
// This function correctly calculates sinh(x) and cosh(x) by calculating exp(x)
201+
// and exp(-x) simultaneously.
202+
// To compute e^x, we perform the following range reduction:
203+
// find hi, mid, lo such that:
204+
// x = (hi + mid) * log(2) + lo, in which
205+
// hi is an integer,
206+
// 0 <= mid * 2^5 < 32 is an integer
207+
// -2^(-5) <= lo * log2(e) <= 2^-5.
208+
// In particular,
209+
// hi + mid = round(x * log2(e) * 2^5) * 2^(-5).
210+
// Then,
211+
// e^x = 2^(hi + mid) * e^lo = 2^hi * 2^mid * e^lo.
212+
// We store 2^mid in the lookup table EXP2_MID_5_BITS, and compute 2^hi * 2^mid
213+
// by adding hi to the exponent field of 2^mid.
214+
// e^lo is computed using a degree-3 minimax polynomial generated by Sollya:
215+
// e^lo ~ P(lo)
216+
// = 1 + lo + c2 * lo^2 + ... + c5 * lo^5
217+
// = (1 + c2*lo^2 + c4*lo^4) + lo * (1 + c3*lo^2 + c5*lo^4)
218+
// = P_even + lo * P_odd
219+
// To compute e^(-x), notice that:
220+
// e^(-x) = 2^(-(hi + mid)) * e^(-lo)
221+
// ~ 2^(-(hi + mid)) * P(-lo)
222+
// = 2^(-(hi + mid)) * (P_even - lo * P_odd)
223+
// So:
224+
// sinh(x) = (e^x - e^(-x)) / 2
225+
// ~ 0.5 * (2^(hi + mid) * (P_even + lo * P_odd) -
226+
// 2^(-(hi + mid)) * (P_even - lo * P_odd))
227+
// = 0.5 * (P_even * (2^(hi + mid) - 2^(-(hi + mid))) +
228+
// lo * P_odd * (2^(hi + mid) + 2^(-(hi + mid))))
229+
// And similarly:
230+
// cosh(x) = (e^x + e^(-x)) / 2
231+
// ~ 0.5 * (P_even * (2^(hi + mid) + 2^(-(hi + mid))) +
232+
// lo * P_odd * (2^(hi + mid) - 2^(-(hi + mid))))
233+
// The main point of these formulas is that the expensive part of calculating
234+
// the polynomials approximating lower parts of e^x and e^(-x) is shared and
235+
// only done once.
236+
template <bool IsSinh> LIBC_INLINE float16 eval_sinh_or_cosh(float16 x) {
237+
float xf = x;
238+
float kf = fputil::nearest_integer(xf * (LOG2F_E * 0x1.0p+5f));
239+
int x_hi_mid_p = static_cast<int>(kf);
240+
int x_hi_mid_m = -x_hi_mid_p;
241+
242+
int x_hi_p = x_hi_mid_p >> 5;
243+
int x_hi_m = x_hi_mid_m >> 5;
244+
int x_mid_p = x_hi_mid_p & 0x1f;
245+
int x_mid_m = x_hi_mid_m & 0x1f;
246+
247+
uint32_t exp2_hi_mid_bits_p =
248+
EXP2_MID_5_BITS[x_mid_p] +
249+
static_cast<uint32_t>(x_hi_p << fputil::FPBits<float>::FRACTION_LEN);
250+
uint32_t exp2_hi_mid_bits_m =
251+
EXP2_MID_5_BITS[x_mid_m] +
252+
static_cast<uint32_t>(x_hi_m << fputil::FPBits<float>::FRACTION_LEN);
253+
// exp2_hi_mid_p = 2^(hi + mid)
254+
float exp2_hi_mid_p = fputil::FPBits<float>(exp2_hi_mid_bits_p).get_val();
255+
// exp2_hi_mid_m = 2^(-(hi + mid))
256+
float exp2_hi_mid_m = fputil::FPBits<float>(exp2_hi_mid_bits_m).get_val();
257+
258+
// exp2_hi_mid_sum = 2^(hi + mid) + 2^(-(hi + mid))
259+
float exp2_hi_mid_sum = exp2_hi_mid_p + exp2_hi_mid_m;
260+
// exp2_hi_mid_diff = 2^(hi + mid) - 2^(-(hi + mid))
261+
float exp2_hi_mid_diff = exp2_hi_mid_p - exp2_hi_mid_m;
262+
263+
// lo = x - (hi + mid) = round(x * log2(e) * 2^5) * log(2) * (-2^(-5)) + x
264+
float lo = fputil::multiply_add(kf, LOGF_2 * -0x1.0p-5f, xf);
265+
float lo_sq = lo * lo;
266+
267+
// Degree-3 minimax polynomial generated by Sollya with the following
268+
// commands:
269+
// > display = hexadecimal;
270+
// > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-5, 2^-5]);
271+
// > 1 + x * P;
272+
constexpr cpp::array<float, 4> COEFFS = {0x1p+0f, 0x1p+0f, 0x1.0004p-1f,
273+
0x1.555778p-3f};
274+
float half_p_odd =
275+
fputil::polyeval(lo_sq, COEFFS[1] * 0.5f, COEFFS[3] * 0.5f);
276+
float half_p_even =
277+
fputil::polyeval(lo_sq, COEFFS[0] * 0.5f, COEFFS[2] * 0.5f);
278+
279+
// sinh(x) = lo * (0.5 * P_odd * (2^(hi + mid) + 2^(-(hi + mid)))) +
280+
// (0.5 * P_even * (2^(hi + mid) - 2^(-(hi + mid))))
281+
if constexpr (IsSinh)
282+
return static_cast<float16>(fputil::multiply_add(
283+
lo, half_p_odd * exp2_hi_mid_sum, half_p_even * exp2_hi_mid_diff));
284+
// cosh(x) = lo * (0.5 * P_odd * (2^(hi + mid) - 2^(-(hi + mid)))) +
285+
// (0.5 * P_even * (2^(hi + mid) + 2^(-(hi + mid))))
286+
return static_cast<float16>(fputil::multiply_add(
287+
lo, half_p_odd * exp2_hi_mid_diff, half_p_even * exp2_hi_mid_sum));
288+
}
289+
177290
} // namespace LIBC_NAMESPACE_DECL
178291

179292
#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPXF16_H

0 commit comments

Comments
 (0)