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