1//===-- Implementation header for log2f16 -----------------------*- 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_LOG2F16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_LOG2F16_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
16#include "expxf16_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/multiply_add.h"
25#include "src/__support/common.h"
26#include "src/__support/macros/config.h"
27#include "src/__support/macros/optimization.h"
28#include "src/__support/macros/properties/cpu_features.h"
29
30namespace LIBC_NAMESPACE_DECL {
31
32namespace math {
33
34namespace log2f16_internal {
35
36#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
37#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
38LIBC_INLINE_VAR constexpr size_t N_LOG2F16_EXCEPTS = 2;
39#else
40LIBC_INLINE_VAR constexpr size_t N_LOG2F16_EXCEPTS = 9;
41#endif
42
43LIBC_INLINE_VAR constexpr fputil::ExceptValues<float16, N_LOG2F16_EXCEPTS>
44 LOG2F16_EXCEPTS = {.values: {
45// (input, RZ output, RU offset, RD offset, RN offset)
46#ifndef LIBC_TARGET_CPU_HAS_FMA_FLOAT
47 // x = 0x1.224p-1, log2f16(x) = -0x1.a34p-1 (RZ)
48 {.input: 0x3889U, .rnd_towardzero_result: 0xba8dU, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 0U},
49 // x = 0x1.e34p-1, log2f16(x) = -0x1.558p-4 (RZ)
50 {.input: 0x3b8dU, .rnd_towardzero_result: 0xad56U, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 0U},
51#endif
52 // x = 0x1.e8cp-1, log2f16(x) = -0x1.128p-4 (RZ)
53 {.input: 0x3ba3U, .rnd_towardzero_result: 0xac4aU, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 0U},
54#ifndef LIBC_TARGET_CPU_HAS_FMA_FLOAT
55 // x = 0x1.f98p-1, log2f16(x) = -0x1.2ep-6 (RZ)
56 {.input: 0x3be6U, .rnd_towardzero_result: 0xa4b8U, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 0U},
57 // x = 0x1.facp-1, log2f16(x) = -0x1.e7p-7 (RZ)
58 {.input: 0x3bebU, .rnd_towardzero_result: 0xa39cU, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 1U},
59#endif
60 // x = 0x1.fb4p-1, log2f16(x) = -0x1.b88p-7 (RZ)
61 {.input: 0x3bedU, .rnd_towardzero_result: 0xa2e2U, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 1U},
62#ifndef LIBC_TARGET_CPU_HAS_FMA_FLOAT
63 // x = 0x1.fecp-1, log2f16(x) = -0x1.cep-9 (RZ)
64 {.input: 0x3bfbU, .rnd_towardzero_result: 0x9b38U, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 1U},
65 // x = 0x1.ffcp-1, log2f16(x) = -0x1.714p-11 (RZ)
66 {.input: 0x3bffU, .rnd_towardzero_result: 0x91c5U, .rnd_upward_offset: 0U, .rnd_downward_offset: 1U, .rnd_tonearest_offset: 1U},
67 // x = 0x1.224p+0, log2f16(x) = 0x1.72cp-3 (RZ)
68 {.input: 0x3c89U, .rnd_towardzero_result: 0x31cbU, .rnd_upward_offset: 1U, .rnd_downward_offset: 0U, .rnd_tonearest_offset: 1U},
69#endif
70 }};
71#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
72
73} // namespace log2f16_internal
74
75LIBC_INLINE float16 log2f16(float16 x) {
76 using namespace math::expxf16_internal;
77 using FPBits = fputil::FPBits<float16>;
78 FPBits x_bits(x);
79
80 uint16_t x_u = x_bits.uintval();
81
82 // If x <= 0, or x is 1, or x is +inf, or x is NaN.
83 if (LIBC_UNLIKELY(x_u == 0U || x_u == 0x3c00U || x_u >= 0x7c00U)) {
84 // log2(NaN) = NaN
85 if (x_bits.is_nan()) {
86 if (x_bits.is_signaling_nan()) {
87 fputil::raise_except_if_required(FE_INVALID);
88 return FPBits::quiet_nan().get_val();
89 }
90
91 return x;
92 }
93
94 // log2(+/-0) = −inf
95 if ((x_u & 0x7fffU) == 0U) {
96 fputil::raise_except_if_required(FE_DIVBYZERO);
97 return FPBits::inf(sign: Sign::NEG).get_val();
98 }
99
100 if (x_u == 0x3c00U)
101 return FPBits::zero().get_val();
102
103 // When x < 0.
104 if (x_u > 0x8000U) {
105 fputil::set_errno_if_required(EDOM);
106 fputil::raise_except_if_required(FE_INVALID);
107 return FPBits::quiet_nan().get_val();
108 }
109
110 // log2(+inf) = +inf
111 return FPBits::inf().get_val();
112 }
113
114#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
115 if (auto r = log2f16_internal::LOG2F16_EXCEPTS.lookup(x_bits: x_u);
116 LIBC_UNLIKELY(r.has_value()))
117 return r.value();
118#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
119
120 // To compute log2(x), we perform the following range reduction:
121 // x = 2^m * 1.mant,
122 // log2(x) = m + log2(1.mant).
123 // To compute log2(1.mant), let f be the highest 6 bits including the hidden
124 // bit, and d be the difference (1.mant - f), i.e., the remaining 5 bits of
125 // the mantissa, then:
126 // log2(1.mant) = log2(f) + log2(1.mant / f)
127 // = log2(f) + log2(1 + d/f)
128 // since d/f is sufficiently small.
129 // We store log2(f) and 1/f in the lookup tables LOG2F_F and ONE_OVER_F_F
130 // respectively.
131
132 int m = -FPBits::EXP_BIAS;
133
134 // When x is subnormal, normalize it.
135 if ((x_u & FPBits::EXP_MASK) == 0U) {
136 // Can't pass an integer to fputil::cast directly.
137 constexpr float NORMALIZE_EXP = 1U << FPBits::FRACTION_LEN;
138 x_bits = FPBits(x_bits.get_val() * fputil::cast<float16>(x: NORMALIZE_EXP));
139 x_u = x_bits.uintval();
140 m -= FPBits::FRACTION_LEN;
141 }
142
143 uint16_t mant = x_bits.get_mantissa();
144 // Leading 10 - 5 = 5 bits of the mantissa.
145 int f = mant >> 5;
146 // Unbiased exponent.
147 m += x_u >> FPBits::FRACTION_LEN;
148
149 // Set bits to 1.mant instead of 2^m * 1.mant.
150 x_bits.set_biased_exponent(FPBits::EXP_BIAS);
151 float mant_f = x_bits.get_val();
152 // v = 1.mant * 1/f - 1 = d/f
153 float v = fputil::multiply_add(x: mant_f, y: ONE_OVER_F_F[f], z: -1.0f);
154
155 // Degree-3 minimax polynomial generated by Sollya with the following
156 // commands:
157 // > display = hexadecimal;
158 // > P = fpminimax(log2(1 + x)/x, 2, [|SG...|], [-2^-5, 2^-5]);
159 // > x * P;
160 float log2p1_d_over_f =
161 v * fputil::polyeval(x: v, a0: 0x1.715476p+0f, a: -0x1.71771ap-1f, a: 0x1.ecb38ep-2f);
162 // log2(1.mant) = log2(f) + log2(1 + d/f)
163 float log2_1_mant = LOG2F_F[f] + log2p1_d_over_f;
164 return fputil::cast<float16>(x: static_cast<float>(m) + log2_1_mant);
165}
166
167} // namespace math
168} // namespace LIBC_NAMESPACE_DECL
169
170#endif // LIBC_TYPES_HAS_FLOAT16
171
172#endif // LLVM_LIBC_SRC___SUPPORT_MATH_LOG2F16_H
173