1 //===-- Floating-point manipulation functions -------------------*- 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_FPUTIL_MANIPULATIONFUNCTIONS_H
10 #define LLVM_LIBC_SRC___SUPPORT_FPUTIL_MANIPULATIONFUNCTIONS_H
11
12 #include "FPBits.h"
13 #include "NearestIntegerOperations.h"
14 #include "NormalFloat.h"
15 #include "dyadic_float.h"
16 #include "rounding_mode.h"
17
18 #include "hdr/math_macros.h"
19 #include "src/__support/CPP/bit.h"
20 #include "src/__support/CPP/limits.h" // INT_MAX, INT_MIN
21 #include "src/__support/CPP/type_traits.h"
22 #include "src/__support/FPUtil/FEnvImpl.h"
23 #include "src/__support/macros/attributes.h"
24 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
25
26 namespace LIBC_NAMESPACE {
27 namespace fputil {
28
29 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
frexp(T x,int & exp)30 LIBC_INLINE T frexp(T x, int &exp) {
31 FPBits<T> bits(x);
32 if (bits.is_inf_or_nan())
33 return x;
34 if (bits.is_zero()) {
35 exp = 0;
36 return x;
37 }
38
39 NormalFloat<T> normal(bits);
40 exp = normal.exponent + 1;
41 normal.exponent = -1;
42 return normal;
43 }
44
45 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
modf(T x,T & iptr)46 LIBC_INLINE T modf(T x, T &iptr) {
47 FPBits<T> bits(x);
48 if (bits.is_zero() || bits.is_nan()) {
49 iptr = x;
50 return x;
51 } else if (bits.is_inf()) {
52 iptr = x;
53 return FPBits<T>::zero(bits.sign()).get_val();
54 } else {
55 iptr = trunc(x);
56 if (x == iptr) {
57 // If x is already an integer value, then return zero with the right
58 // sign.
59 return FPBits<T>::zero(bits.sign()).get_val();
60 } else {
61 return x - iptr;
62 }
63 }
64 }
65
66 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
copysign(T x,T y)67 LIBC_INLINE T copysign(T x, T y) {
68 FPBits<T> xbits(x);
69 xbits.set_sign(FPBits<T>(y).sign());
70 return xbits.get_val();
71 }
72
73 template <typename T> struct IntLogbConstants;
74
75 template <> struct IntLogbConstants<int> {
76 LIBC_INLINE_VAR static constexpr int FP_LOGB0 = FP_ILOGB0;
77 LIBC_INLINE_VAR static constexpr int FP_LOGBNAN = FP_ILOGBNAN;
78 LIBC_INLINE_VAR static constexpr int T_MAX = INT_MAX;
79 LIBC_INLINE_VAR static constexpr int T_MIN = INT_MIN;
80 };
81
82 template <> struct IntLogbConstants<long> {
83 LIBC_INLINE_VAR static constexpr long FP_LOGB0 = FP_ILOGB0;
84 LIBC_INLINE_VAR static constexpr long FP_LOGBNAN = FP_ILOGBNAN;
85 LIBC_INLINE_VAR static constexpr long T_MAX = LONG_MAX;
86 LIBC_INLINE_VAR static constexpr long T_MIN = LONG_MIN;
87 };
88
89 template <typename T, typename U>
90 LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_floating_point_v<U>, T>
91 intlogb(U x) {
92 FPBits<U> bits(x);
93 if (LIBC_UNLIKELY(bits.is_zero() || bits.is_inf_or_nan())) {
94 set_errno_if_required(EDOM);
95 raise_except_if_required(FE_INVALID);
96
97 if (bits.is_zero())
98 return IntLogbConstants<T>::FP_LOGB0;
99 if (bits.is_nan())
100 return IntLogbConstants<T>::FP_LOGBNAN;
101 // bits is inf.
102 return IntLogbConstants<T>::T_MAX;
103 }
104
105 DyadicFloat<FPBits<U>::STORAGE_LEN> normal(bits.get_val());
106 int exponent = normal.get_unbiased_exponent();
107 // The C standard does not specify the return value when an exponent is
108 // out of int range. However, XSI conformance required that INT_MAX or
109 // INT_MIN are returned.
110 // NOTE: It is highly unlikely that exponent will be out of int range as
111 // the exponent is only 15 bits wide even for the 128-bit floating point
112 // format.
113 if (LIBC_UNLIKELY(exponent > IntLogbConstants<T>::T_MAX ||
114 exponent < IntLogbConstants<T>::T_MIN)) {
115 set_errno_if_required(ERANGE);
116 raise_except_if_required(FE_INVALID);
117 return exponent > 0 ? IntLogbConstants<T>::T_MAX
118 : IntLogbConstants<T>::T_MIN;
119 }
120
121 return static_cast<T>(exponent);
122 }
123
124 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
125 LIBC_INLINE constexpr T logb(T x) {
126 FPBits<T> bits(x);
127 if (LIBC_UNLIKELY(bits.is_zero() || bits.is_inf_or_nan())) {
128 if (bits.is_nan())
129 return x;
130
131 raise_except_if_required(FE_DIVBYZERO);
132
133 if (bits.is_zero()) {
134 set_errno_if_required(ERANGE);
135 return FPBits<T>::inf(Sign::NEG).get_val();
136 }
137 // bits is inf.
138 return FPBits<T>::inf().get_val();
139 }
140
141 DyadicFloat<FPBits<T>::STORAGE_LEN> normal(bits.get_val());
142 return static_cast<T>(normal.get_unbiased_exponent());
143 }
144
145 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
146 LIBC_INLINE constexpr T ldexp(T x, int exp) {
147 FPBits<T> bits(x);
148 if (LIBC_UNLIKELY((exp == 0) || bits.is_zero() || bits.is_inf_or_nan()))
149 return x;
150
151 // NormalFloat uses int32_t to store the true exponent value. We should ensure
152 // that adding |exp| to it does not lead to integer rollover. But, if |exp|
153 // value is larger the exponent range for type T, then we can return infinity
154 // early. Because the result of the ldexp operation can be a subnormal number,
155 // we need to accommodate the (mantissaWidth + 1) worth of shift in
156 // calculating the limit.
157 constexpr int EXP_LIMIT =
158 FPBits<T>::MAX_BIASED_EXPONENT + FPBits<T>::FRACTION_LEN + 1;
159 if (LIBC_UNLIKELY(exp > EXP_LIMIT)) {
160 int rounding_mode = quick_get_round();
161 Sign sign = bits.sign();
162
163 if ((sign == Sign::POS && rounding_mode == FE_DOWNWARD) ||
164 (sign == Sign::NEG && rounding_mode == FE_UPWARD) ||
165 (rounding_mode == FE_TOWARDZERO))
166 return FPBits<T>::max_normal(sign).get_val();
167
168 set_errno_if_required(ERANGE);
169 raise_except_if_required(FE_OVERFLOW);
170 return FPBits<T>::inf(sign).get_val();
171 }
172
173 // Similarly on the negative side we return zero early if |exp| is too small.
174 if (LIBC_UNLIKELY(exp < -EXP_LIMIT)) {
175 int rounding_mode = quick_get_round();
176 Sign sign = bits.sign();
177
178 if ((sign == Sign::POS && rounding_mode == FE_UPWARD) ||
179 (sign == Sign::NEG && rounding_mode == FE_DOWNWARD))
180 return FPBits<T>::min_subnormal(sign).get_val();
181
182 set_errno_if_required(ERANGE);
183 raise_except_if_required(FE_UNDERFLOW);
184 return FPBits<T>::zero(sign).get_val();
185 }
186
187 // For all other values, NormalFloat to T conversion handles it the right way.
188 DyadicFloat<FPBits<T>::STORAGE_LEN> normal(bits.get_val());
189 normal.exponent += exp;
190 return static_cast<T>(normal);
191 }
192
193 template <typename T, typename U,
194 cpp::enable_if_t<cpp::is_floating_point_v<T> &&
195 cpp::is_floating_point_v<U> &&
196 (sizeof(T) <= sizeof(U)),
197 int> = 0>
198 LIBC_INLINE T nextafter(T from, U to) {
199 FPBits<T> from_bits(from);
200 if (from_bits.is_nan())
201 return from;
202
203 FPBits<U> to_bits(to);
204 if (to_bits.is_nan())
205 return static_cast<T>(to);
206
207 // NOTE: This would work only if `U` has a greater or equal precision than
208 // `T`. Otherwise `from` could loose its precision and the following statement
209 // could incorrectly evaluate to `true`.
210 if (static_cast<U>(from) == to)
211 return static_cast<T>(to);
212
213 using StorageType = typename FPBits<T>::StorageType;
214 if (from != T(0)) {
215 if ((static_cast<U>(from) < to) == (from > T(0))) {
216 from_bits = FPBits<T>(StorageType(from_bits.uintval() + 1));
217 } else {
218 from_bits = FPBits<T>(StorageType(from_bits.uintval() - 1));
219 }
220 } else {
221 from_bits = FPBits<T>::min_subnormal(to_bits.sign());
222 }
223
224 if (from_bits.is_subnormal())
225 raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
226 else if (from_bits.is_inf())
227 raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
228
229 return from_bits.get_val();
230 }
231
232 template <bool IsDown, typename T,
233 cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
234 LIBC_INLINE constexpr T nextupdown(T x) {
235 constexpr Sign sign = IsDown ? Sign::NEG : Sign::POS;
236
237 FPBits<T> xbits(x);
238 if (xbits.is_nan() || xbits == FPBits<T>::max_normal(sign) ||
239 xbits == FPBits<T>::inf(sign))
240 return x;
241
242 using StorageType = typename FPBits<T>::StorageType;
243 if (x != T(0)) {
244 if (xbits.sign() == sign) {
245 xbits = FPBits<T>(StorageType(xbits.uintval() + 1));
246 } else {
247 xbits = FPBits<T>(StorageType(xbits.uintval() - 1));
248 }
249 } else {
250 xbits = FPBits<T>::min_subnormal(sign);
251 }
252
253 return xbits.get_val();
254 }
255
256 } // namespace fputil
257 } // namespace LIBC_NAMESPACE
258
259 #ifdef LIBC_TYPES_LONG_DOUBLE_IS_X86_FLOAT80
260 #include "x86_64/NextAfterLongDouble.h"
261 #include "x86_64/NextUpDownLongDouble.h"
262 #endif // LIBC_TYPES_LONG_DOUBLE_IS_X86_FLOAT80
263
264 #endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_MANIPULATIONFUNCTIONS_H
265