1//===----------------------------------------------------------------------===//
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/// \file
10/// Implementation of the exp2bf16 function, which computes 2^x for bfloat16.
11///
12//===----------------------------------------------------------------------===//
13
14#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP2BF16_H
15#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2BF16_H
16
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/bfloat16.h"
23#include "src/__support/FPUtil/cast.h"
24#include "src/__support/FPUtil/multiply_add.h"
25#include "src/__support/FPUtil/nearest_integer.h"
26#include "src/__support/FPUtil/rounding_mode.h"
27#include "src/__support/macros/attributes.h"
28#include "src/__support/macros/config.h"
29#include "src/__support/macros/optimization.h"
30#include "src/__support/macros/properties/types.h"
31
32namespace LIBC_NAMESPACE_DECL {
33namespace math {
34
35namespace exp2bf16_internal {
36// Generated by Sollya with the following commands:
37// > display = hexadecimal;
38// > for i in [|-126.5, -127, -127.5, -128, -129, -130, -131, -132, -133|]
39// do print(round(2^i, SG, RN) @ "f,");
40LIBC_INLINE_VAR constexpr float EXP2_DENORMAL[9] = {
41 0x1.6a09e8p-127f, 0x1p-127f, 0x1.6a09e8p-128f, 0x1p-128f, 0x1p-129f,
42 0x1p-130f, 0x1p-131f, 0x1p-132f, 0x1p-133f,
43};
44
45// Generated by Sollya with the following commands:
46// > display = hexadecimal;
47// > for i from 0 to 15 do print(round(2^(i/16), SG, RN) @ "f,");
48LIBC_INLINE_VAR constexpr float EXP2_MID[16] = {
49 0x1p0f, 0x1.0b5586p0f, 0x1.172b84p0f, 0x1.2387a6p0f,
50 0x1.306fep0f, 0x1.3dea64p0f, 0x1.4bfdaep0f, 0x1.5ab07ep0f,
51 0x1.6a09e6p0f, 0x1.7a1148p0f, 0x1.8ace54p0f, 0x1.9c4918p0f,
52 0x1.ae89fap0f, 0x1.c199bep0f, 0x1.d5818ep0f, 0x1.ea4afap0f,
53};
54} // namespace exp2bf16_internal
55
56LIBC_INLINE bfloat16 exp2bf16(bfloat16 x) {
57 using FPBits = fputil::FPBits<bfloat16>;
58 FPBits x_bits(x);
59 uint16_t x_uint = x_bits.uintval();
60
61 // We first check for NaN
62 if (LIBC_UNLIKELY(x_bits.is_nan())) {
63 // If it is a signaling NaN, then we raise an FE_INVALID exception
64 // and return a quiet NaN
65 if (x_bits.is_signaling_nan()) {
66 fputil::raise_except_if_required(FE_INVALID);
67 return FPBits::quiet_nan().get_val();
68 }
69 // If not, just return the NaN
70 return x;
71 }
72
73 // exp2bf16 may only give numeric output for input x in (-134.0, 128.0).
74 // All the remaining exceptional cases (overflow, underflow and denormal
75 // outputs) have |x| >= 126.5, so we filter them with a single check.
76 constexpr uint16_t EXCEPTIONAL_ABS_BASE = 0x42FD; // 0x1.fap6 = 126.5
77 if (LIBC_UNLIKELY(x_bits.abs().uintval() >= EXCEPTIONAL_ABS_BASE)) {
78 if (x_bits.is_pos()) {
79 // Overflow: x >= 128.0 (0x4300 = 0x1p7)
80 if (x_uint >= 0x4300) {
81 // exp2(Inf) = Inf
82 if (x_bits.is_inf())
83 return x;
84
85#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
86 // If we can only RN, then it'll go towards INF
87 fputil::set_errno_if_required(ERANGE);
88 fputil::raise_overflow_except_if_required<float>();
89 return FPBits::inf().get_val();
90#else
91 // If not, get what other rounding modes are available
92 switch (fputil::quick_get_round()) {
93 case FE_TONEAREST:
94 case FE_UPWARD:
95 // Both of these go to Inf
96 fputil::set_errno_if_required(ERANGE);
97 fputil::raise_overflow_except_if_required<float>();
98 return FPBits::inf().get_val();
99 default:
100 // RD and RZ go to the maximum possible value in the
101 // range of the function, being the maximum normal.
102 return FPBits::max_normal().get_val();
103 }
104#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
105 }
106 } else {
107 // Underflow: x <= -134.0 (0xC306 = -0x1.0cp7)
108 if (x_uint >= 0xC306) {
109 // exp2(-Inf) = 0
110 if (x_bits.is_inf())
111 return FPBits::zero().get_val();
112
113 // Need these for every output now.
114 fputil::set_errno_if_required(ERANGE);
115 fputil::raise_underflow_except_if_required<float>();
116#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
117 return FPBits::zero().get_val();
118#else
119 switch (fputil::quick_get_round()) {
120 case FE_UPWARD:
121 // RU rounds up to the minimum subnormal
122 return FPBits::min_subnormal().get_val();
123 default:
124 // All other rounding modes go down to zero
125 return FPBits::zero().get_val();
126 }
127#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
128 }
129
130 // Handle denormal output separately, by using a LUT
131 // We are handling these separately because the trick of adding hi
132 // directly to the exponent of the result will not work if the
133 // output is denormal. There are only 9 such inputs and they are
134 // all back-to-back when expressed as an uint16_t, with
135 // x_uint >= 0xC2FD (-0x1.fap6 = -126.5) and
136 // x_uint <= 0xC305 (-0x1.0ap7 = -133.0).
137 constexpr uint16_t DENORMAL_LUT_BASE = 0xC2FD;
138 return fputil::cast<bfloat16>(
139 x: exp2bf16_internal::EXP2_DENORMAL[x_uint - DENORMAL_LUT_BASE]);
140 }
141 }
142
143 // For inputs not handled above in -134.0 < x < 128.0,
144 // we perform the following range reduction:
145 // x = hi + mid + lo
146 // where,
147 // hi is an integer
148 // mid is a grid point of the form j / 2^n in [0, 1)
149 // lo is what's leftover ie x - k / 2^n
150 // and |lo| <= 2^(-n-1)
151 // with,
152 // n = 4
153 // k = round(2^n * x)
154 // hi = k >> n
155 // j = k & (2^n - 1)
156 // then,
157 // exp2(x) = exp2(hi + mid + lo)
158 // = exp2(hi) * exp2(mid) * exp2(lo)
159 //
160 // To get exp2(hi), we'll just add hi to the exponent bits of the
161 // result of exp2(mid) * exp2(lo), which is exact.
162 // exp2(mid) we shall obtain by consulting our LUT,
163 // exp2(lo) by the polynomial.
164
165 float x_float = fputil::cast<float>(x);
166 float k_float = fputil::nearest_integer(x: x_float * 16.0f);
167
168 int k_int = static_cast<int>(k_float);
169 int hi = k_int >> 4;
170 int j = k_int & (16 - 1);
171
172 // lo = x - (hi + mid) = x - k / 2^n
173 // = k * (-2^(-n)) + x
174 float lo = fputil::multiply_add(x: k_float, y: -0.0625f, z: x_float);
175
176 float exp2_mid = exp2bf16_internal::EXP2_MID[j];
177
178 // Degree-3 minimax polynomial generated by Sollya with the following
179 // commands:
180 // > display = hexadecimal;
181 // > P = fpminimax(2^x, [|1, 2, 3|], [|SG...|], [-2^-5, 2^-5], relative,
182 // floating, 1);
183 // 0x1p0 + x * (0x1.62e43p-1 + x * (0x1.ebfff4p-3 + x * 0x1.c6aceep-5))
184 float exp2_lo = fputil::polyeval(x: lo, a0: 0x1p0f, a: 0x1.62e43p-1f, a: 0x1.ebfff4p-3f,
185 a: 0x1.c6aceep-5f);
186
187 // We tried using Estrin's scheme as follows,
188 // float lo2 = lo * lo;
189 // float p0 = fputil::multiply_add(lo, 0x1.62e43p-1f, 0x1p0f);
190 // float p1 = fputil::multiply_add(lo, 0x1.c6aceep-5f, 0x1.ebfff4p-3f);
191 // float exp2_lo = fputil::multiply_add(lo2, p1, p0);
192 // however that led to a 1-ULP error in all RU, negative, small input cases.
193 // The latency increase wasn't significant (~1%) so we reverted back
194 // to Horner's scheme.
195
196 float res = exp2_mid * exp2_lo;
197
198 fputil::FPBits<float> res_bits(res);
199 uint32_t res_int = res_bits.uintval();
200
201 uint32_t exp_hi = hi << fputil::FPBits<float>::FRACTION_LEN;
202 res_int += exp_hi;
203
204 res = fputil::FPBits<float>(res_int).get_val();
205 return fputil::cast<bfloat16>(x: res);
206}
207
208} // namespace math
209} // namespace LIBC_NAMESPACE_DECL
210
211#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP2BF16_H
212