1//===-- Implementation header for expf16 ------------------------*- 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_EXPF16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
16#include "expf16_utils.h"
17#include "hdr/errno_macros.h"
18#include "hdr/fenv_macros.h"
19#include "src/__support/FPUtil/FEnvImpl.h"
20#include "src/__support/FPUtil/FPBits.h"
21#include "src/__support/FPUtil/PolyEval.h"
22#include "src/__support/FPUtil/cast.h"
23#include "src/__support/FPUtil/except_value_utils.h"
24#include "src/__support/FPUtil/rounding_mode.h"
25#include "src/__support/common.h"
26#include "src/__support/macros/config.h"
27#include "src/__support/macros/optimization.h"
28
29namespace LIBC_NAMESPACE_DECL {
30
31namespace math {
32
33LIBC_INLINE constexpr float16 expf16(float16 x) {
34#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
35 constexpr fputil::ExceptValues<float16, 2> EXPF16_EXCEPTS_LO = {.values: {
36 // (input, RZ output, RU offset, RD offset, RN offset)
37 // x = 0x1.de4p-8, expf16(x) = 0x1.01cp+0 (RZ)
38 {.input: 0x1f79U, .rnd_towardzero_result: 0x3c07U, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 0U},
39 // x = 0x1.73cp-6, expf16(x) = 0x1.05cp+0 (RZ)
40 {.input: 0x25cfU, .rnd_towardzero_result: 0x3c17U, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 0U},
41 }};
42
43 constexpr fputil::ExceptValues<float16, 3> EXPF16_EXCEPTS_HI = {.values: {
44 // (input, RZ output, RU offset, RD offset, RN offset)
45 // x = 0x1.c34p+0, expf16(x) = 0x1.74cp+2 (RZ)
46 {.input: 0x3f0dU, .rnd_towardzero_result: 0x45d3U, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 1U},
47 // x = -0x1.488p-5, expf16(x) = 0x1.ebcp-1 (RZ)
48 {.input: 0xa922U, .rnd_towardzero_result: 0x3bafU, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 0U},
49 // x = -0x1.55p-5, expf16(x) = 0x1.ebp-1 (RZ)
50 {.input: 0xa954U, .rnd_towardzero_result: 0x3bacU, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 0U},
51 }};
52#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
53
54 using FPBits = fputil::FPBits<float16>;
55 FPBits x_bits(x);
56
57 uint16_t x_u = x_bits.uintval();
58 uint16_t x_abs = x_u & 0x7fffU;
59
60 // When 0 < |x| <= 2^(-5), or |x| >= 12, or x is NaN.
61 if (LIBC_UNLIKELY(x_abs <= 0x2800U || x_abs >= 0x4a00U)) {
62 // exp(NaN) = NaN
63 if (x_bits.is_nan()) {
64 if (x_bits.is_signaling_nan()) {
65 fputil::raise_except_if_required(FE_INVALID);
66 return FPBits::quiet_nan().get_val();
67 }
68
69 return x;
70 }
71
72 // When x >= 12.
73 if (x_bits.is_pos() && x_abs >= 0x4a00U) {
74 // exp(+inf) = +inf
75 if (x_bits.is_inf())
76 return FPBits::inf().get_val();
77
78#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
79 fputil::set_errno_if_required(ERANGE);
80 fputil::raise_except_if_required(FE_OVERFLOW);
81 return FPBits::inf().get_val();
82#else
83 switch (fputil::quick_get_round()) {
84 case FE_TONEAREST:
85 case FE_UPWARD:
86 fputil::set_errno_if_required(ERANGE);
87 fputil::raise_except_if_required(FE_OVERFLOW);
88 return FPBits::inf().get_val();
89 default:
90 return FPBits::max_normal().get_val();
91 }
92#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
93 }
94
95 // When x <= -18.
96 if (x_u >= 0xcc80U) {
97 // exp(-inf) = +0
98 if (x_bits.is_inf())
99 return FPBits::zero().get_val();
100
101 fputil::set_errno_if_required(ERANGE);
102 fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
103
104#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
105 return FPBits::zero().get_val();
106#else
107 switch (fputil::quick_get_round()) {
108 case FE_UPWARD:
109 return FPBits::min_subnormal().get_val();
110 default:
111 return FPBits::zero().get_val();
112 }
113#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
114 }
115
116 // When 0 < |x| <= 2^(-5).
117 if (x_abs <= 0x2800U && !x_bits.is_zero()) {
118#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
119 if (auto r = EXPF16_EXCEPTS_LO.lookup(x_bits: x_u); LIBC_UNLIKELY(r.has_value()))
120 return r.value();
121#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
122
123 float xf = x;
124 // Degree-3 minimax polynomial generated by Sollya with the following
125 // commands:
126 // > display = hexadecimal;
127 // > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-5, 2^-5]);
128 // > 1 + x * P;
129 return fputil::cast<float16>(
130 x: fputil::polyeval(x: xf, a0: 0x1p+0f, a: 0x1p+0f, a: 0x1.0004p-1f, a: 0x1.555778p-3f));
131 }
132 }
133
134#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
135 if (auto r = EXPF16_EXCEPTS_HI.lookup(x_bits: x_u); LIBC_UNLIKELY(r.has_value()))
136 return r.value();
137#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
138
139 // exp(x) = exp(hi + mid) * exp(lo)
140 auto [exp_hi_mid, exp_lo] = exp_range_reduction(x);
141 return fputil::cast<float16>(x: exp_hi_mid * exp_lo);
142}
143
144} // namespace math
145
146} // namespace LIBC_NAMESPACE_DECL
147
148#endif // LIBC_TYPES_HAS_FLOAT16
149
150#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_H
151