1//===-- Implementation header for exp2m1f16 ----------------------*- 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_EXP2M1F16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2M1F16_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
16#include "src/__support/FPUtil/FEnvImpl.h"
17#include "src/__support/FPUtil/FPBits.h"
18#include "src/__support/FPUtil/PolyEval.h"
19#include "src/__support/FPUtil/cast.h"
20#include "src/__support/FPUtil/except_value_utils.h"
21#include "src/__support/FPUtil/multiply_add.h"
22#include "src/__support/FPUtil/rounding_mode.h"
23#include "src/__support/macros/config.h"
24#include "src/__support/macros/optimization.h"
25#include "src/__support/macros/properties/cpu_features.h"
26#include "src/__support/math/expxf16_utils.h"
27
28namespace LIBC_NAMESPACE_DECL {
29
30namespace math {
31
32LIBC_INLINE constexpr float16 exp2m1f16(float16 x) {
33#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
34 constexpr fputil::ExceptValues<float16, 6> EXP2M1F16_EXCEPTS_LO = {.values: {
35 // (input, RZ output, RU offset, RD offset, RN offset)
36 // x = 0x1.cf4p-13, exp2m1f16(x) = 0x1.41p-13 (RZ)
37 {.input: 0x0b3dU, .rnd_towardzero_result: 0x0904U, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 1U},
38 // x = 0x1.4fcp-12, exp2m1f16(x) = 0x1.d14p-13 (RZ)
39 {.input: 0x0d3fU, .rnd_towardzero_result: 0x0b45U, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 1U},
40 // x = 0x1.63p-11, exp2m1f16(x) = 0x1.ec4p-12 (RZ)
41 {.input: 0x118cU, .rnd_towardzero_result: 0x0fb1U, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 0U},
42 // x = 0x1.6fp-7, exp2m1f16(x) = 0x1.fe8p-8 (RZ)
43 {.input: 0x21bcU, .rnd_towardzero_result: 0x1ffaU, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 1U},
44 // x = -0x1.c6p-10, exp2m1f16(x) = -0x1.3a8p-10 (RZ)
45 {.input: 0x9718U, .rnd_towardzero_result: 0x94eaU, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 0U},
46 // x = -0x1.cfcp-10, exp2m1f16(x) = -0x1.414p-10 (RZ)
47 {.input: 0x973fU, .rnd_towardzero_result: 0x9505U, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 0U},
48 }};
49
50#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
51 constexpr size_t N_EXP2M1F16_EXCEPTS_HI = 6;
52#else
53 constexpr size_t N_EXP2M1F16_EXCEPTS_HI = 7;
54#endif
55
56 constexpr fputil::ExceptValues<float16, N_EXP2M1F16_EXCEPTS_HI>
57 EXP2M1F16_EXCEPTS_HI = {.values: {
58 // (input, RZ output, RU offset, RD offset, RN offset)
59 // x = 0x1.e58p-3, exp2m1f16(x) = 0x1.6dcp-3 (RZ)
60 {.input: 0x3396U, .rnd_towardzero_result: 0x31b7U, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 0U},
61#ifndef LIBC_TARGET_CPU_HAS_FMA_FLOAT
62 // x = 0x1.2e8p-2, exp2m1f16(x) = 0x1.d14p-3 (RZ)
63 {.input: 0x34baU, .rnd_towardzero_result: 0x3345U, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 0U},
64#endif
65 // x = 0x1.ad8p-2, exp2m1f16(x) = 0x1.598p-2 (RZ)
66 {.input: 0x36b6U, .rnd_towardzero_result: 0x3566U, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 0U},
67#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
68 // x = 0x1.edcp-2, exp2m1f16(x) = 0x1.964p-2 (RZ)
69 {0x37b7U, 0x3659U, 1U, 0U, 1U},
70#endif
71 // x = -0x1.804p-3, exp2m1f16(x) = -0x1.f34p-4 (RZ)
72 {.input: 0xb201U, .rnd_towardzero_result: 0xafcdU, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 1U},
73 // x = -0x1.f3p-3, exp2m1f16(x) = -0x1.3e4p-3 (RZ)
74 {.input: 0xb3ccU, .rnd_towardzero_result: 0xb0f9U, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 0U},
75 // x = -0x1.294p-1, exp2m1f16(x) = -0x1.53p-2 (RZ)
76 {.input: 0xb8a5U, .rnd_towardzero_result: 0xb54cU, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 1U},
77#ifndef LIBC_TARGET_CPU_HAS_FMA_FLOAT
78 // x = -0x1.a34p-1, exp2m1f16(x) = -0x1.bb4p-2 (RZ)
79 {.input: 0xba8dU, .rnd_towardzero_result: 0xb6edU, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 1U},
80#endif
81 }};
82#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
83
84 using namespace math::expxf16_internal;
85 using FPBits = fputil::FPBits<float16>;
86 FPBits x_bits(x);
87
88 uint16_t x_u = x_bits.uintval();
89 uint16_t x_abs = x_u & 0x7fffU;
90
91 // When |x| <= 2^(-3), or |x| >= 11, or x is NaN.
92 if (LIBC_UNLIKELY(x_abs <= 0x3000U || x_abs >= 0x4980U)) {
93 // exp2m1(NaN) = NaN
94 if (x_bits.is_nan()) {
95 if (x_bits.is_signaling_nan()) {
96 fputil::raise_except_if_required(FE_INVALID);
97 return FPBits::quiet_nan().get_val();
98 }
99
100 return x;
101 }
102
103 // When x >= 16.
104 if (x_u >= 0x4c00 && x_bits.is_pos()) {
105 // exp2m1(+inf) = +inf
106 if (x_bits.is_inf())
107 return FPBits::inf().get_val();
108
109#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
110 fputil::set_errno_if_required(ERANGE);
111 fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
112 return FPBits::inf().get_val();
113#else
114 switch (fputil::quick_get_round()) {
115 case FE_TONEAREST:
116 case FE_UPWARD:
117 fputil::set_errno_if_required(ERANGE);
118 fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
119 return FPBits::inf().get_val();
120 default:
121 return FPBits::max_normal().get_val();
122 }
123#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
124 }
125
126 // When x < -11.
127 if (x_u > 0xc980U) {
128 // exp2m1(-inf) = -1
129 if (x_bits.is_inf())
130 return FPBits::one(sign: Sign::NEG).get_val();
131
132 // When -12 < x < -11, round(2^x - 1, HP, RN) = -0x1.ffcp-1.
133 if (x_u < 0xca00U)
134 return fputil::round_result_slightly_down(
135 value_rn: fputil::cast<float16>(x: -0x1.ffcp-1));
136
137 // When x <= -12, round(2^x - 1, HP, RN) = -1.
138#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
139 fputil::raise_except_if_required(FE_INEXACT);
140 return FPBits::one(Sign::NEG).get_val();
141#else
142 switch (fputil::quick_get_round()) {
143 case FE_TONEAREST:
144 case FE_DOWNWARD:
145 fputil::raise_except_if_required(FE_INEXACT);
146 return FPBits::one(sign: Sign::NEG).get_val();
147 default:
148 fputil::raise_except_if_required(FE_INEXACT);
149 return fputil::cast<float16>(x: -0x1.ffcp-1);
150 }
151#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
152 }
153
154 // When |x| <= 2^(-3).
155 if (x_abs <= 0x3000U) {
156#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
157 if (auto r = EXP2M1F16_EXCEPTS_LO.lookup(x_bits: x_u);
158 LIBC_UNLIKELY(r.has_value()))
159 return r.value();
160#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
161
162 float xf = x;
163 // Degree-5 minimax polynomial generated by Sollya with the following
164 // commands:
165 // > display = hexadecimal;
166 // > P = fpminimax((2^x - 1)/x, 4, [|SG...|], [-2^-3, 2^-3]);
167 // > x * P;
168 return fputil::cast<float16>(
169 x: xf * fputil::polyeval(x: xf, a0: 0x1.62e43p-1f, a: 0x1.ebfbdep-3f,
170 a: 0x1.c6af88p-5f, a: 0x1.3b45d6p-7f,
171 a: 0x1.641e7cp-10f));
172 }
173 }
174
175#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
176 if (auto r = EXP2M1F16_EXCEPTS_HI.lookup(x_bits: x_u); LIBC_UNLIKELY(r.has_value()))
177 return r.value();
178#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
179
180 // exp2(x) = exp2(hi + mid) * exp2(lo)
181 auto [exp2_hi_mid, exp2_lo] = exp2_range_reduction(x);
182 // exp2m1(x) = exp2(hi + mid) * exp2(lo) - 1
183 return fputil::cast<float16>(
184 x: fputil::multiply_add(x: exp2_hi_mid, y: exp2_lo, z: -1.0f));
185}
186
187} // namespace math
188
189} // namespace LIBC_NAMESPACE_DECL
190
191#endif // LIBC_TYPES_HAS_FLOAT16
192
193#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP2M1F16_H
194