xref: /llvm-project/libc/src/__support/FPUtil/NearestIntegerOperations.h (revision 111b062f63ba52552f41e425449ba1db048dc51a) 
1 //===-- Nearest integer floating-point operations ---------------*- 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_NEARESTINTEGEROPERATIONS_H
10 #define LLVM_LIBC_SRC___SUPPORT_FPUTIL_NEARESTINTEGEROPERATIONS_H
11 
12 #include "FEnvImpl.h"
13 #include "FPBits.h"
14 #include "rounding_mode.h"
15 
16 #include "hdr/math_macros.h"
17 #include "src/__support/CPP/type_traits.h"
18 #include "src/__support/common.h"
19 #include "src/__support/macros/config.h"
20 
21 namespace LIBC_NAMESPACE_DECL {
22 namespace fputil {
23 
24 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
25 LIBC_INLINE T trunc(T x) {
26   using StorageType = typename FPBits<T>::StorageType;
27   FPBits<T> bits(x);
28 
29   // If x is infinity or NaN, return it.
30   // If it is zero also we should return it as is, but the logic
31   // later in this function takes care of it. But not doing a zero
32   // check, we improve the run time of non-zero values.
33   if (bits.is_inf_or_nan())
34     return x;
35 
36   int exponent = bits.get_exponent();
37 
38   // If the exponent is greater than the most negative mantissa
39   // exponent, then x is already an integer.
40   if (exponent >= static_cast<int>(FPBits<T>::FRACTION_LEN))
41     return x;
42 
43   // If the exponent is such that abs(x) is less than 1, then return 0.
44   if (exponent <= -1)
45     return FPBits<T>::zero(bits.sign()).get_val();
46 
47   int trim_size = FPBits<T>::FRACTION_LEN - exponent;
48   StorageType trunc_mantissa =
49       static_cast<StorageType>((bits.get_mantissa() >> trim_size) << trim_size);
50   bits.set_mantissa(trunc_mantissa);
51   return bits.get_val();
52 }
53 
54 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
55 LIBC_INLINE T ceil(T x) {
56   using StorageType = typename FPBits<T>::StorageType;
57   FPBits<T> bits(x);
58 
59   // If x is infinity NaN or zero, return it.
60   if (bits.is_inf_or_nan() || bits.is_zero())
61     return x;
62 
63   bool is_neg = bits.is_neg();
64   int exponent = bits.get_exponent();
65 
66   // If the exponent is greater than the most negative mantissa
67   // exponent, then x is already an integer.
68   if (exponent >= static_cast<int>(FPBits<T>::FRACTION_LEN))
69     return x;
70 
71   if (exponent <= -1) {
72     if (is_neg)
73       return T(-0.0);
74     else
75       return T(1.0);
76   }
77 
78   uint32_t trim_size = FPBits<T>::FRACTION_LEN - exponent;
79   StorageType x_u = bits.uintval();
80   StorageType trunc_u =
81       static_cast<StorageType>((x_u >> trim_size) << trim_size);
82 
83   // If x is already an integer, return it.
84   if (trunc_u == x_u)
85     return x;
86 
87   bits.set_uintval(trunc_u);
88   T trunc_value = bits.get_val();
89 
90   // If x is negative, the ceil operation is equivalent to the trunc operation.
91   if (is_neg)
92     return trunc_value;
93 
94   return trunc_value + T(1.0);
95 }
96 
97 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
98 LIBC_INLINE T floor(T x) {
99   FPBits<T> bits(x);
100   if (bits.is_neg()) {
101     return -ceil(-x);
102   } else {
103     return trunc(x);
104   }
105 }
106 
107 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
108 LIBC_INLINE T round(T x) {
109   using StorageType = typename FPBits<T>::StorageType;
110   FPBits<T> bits(x);
111 
112   // If x is infinity NaN or zero, return it.
113   if (bits.is_inf_or_nan() || bits.is_zero())
114     return x;
115 
116   int exponent = bits.get_exponent();
117 
118   // If the exponent is greater than the most negative mantissa
119   // exponent, then x is already an integer.
120   if (exponent >= static_cast<int>(FPBits<T>::FRACTION_LEN))
121     return x;
122 
123   if (exponent == -1) {
124     // Absolute value of x is greater than equal to 0.5 but less than 1.
125     return FPBits<T>::one(bits.sign()).get_val();
126   }
127 
128   if (exponent <= -2) {
129     // Absolute value of x is less than 0.5.
130     return FPBits<T>::zero(bits.sign()).get_val();
131   }
132 
133   uint32_t trim_size = FPBits<T>::FRACTION_LEN - exponent;
134   bool half_bit_set =
135       bool(bits.get_mantissa() & (StorageType(1) << (trim_size - 1)));
136   StorageType x_u = bits.uintval();
137   StorageType trunc_u =
138       static_cast<StorageType>((x_u >> trim_size) << trim_size);
139 
140   // If x is already an integer, return it.
141   if (trunc_u == x_u)
142     return x;
143 
144   bits.set_uintval(trunc_u);
145   T trunc_value = bits.get_val();
146 
147   if (!half_bit_set) {
148     // Franctional part is less than 0.5 so round value is the
149     // same as the trunc value.
150     return trunc_value;
151   } else {
152     return bits.is_neg() ? trunc_value - T(1.0) : trunc_value + T(1.0);
153   }
154 }
155 
156 template <typename T>
157 LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_floating_point_v<T>, T>
158 round_using_specific_rounding_mode(T x, int rnd) {
159   using StorageType = typename FPBits<T>::StorageType;
160   FPBits<T> bits(x);
161 
162   // If x is infinity NaN or zero, return it.
163   if (bits.is_inf_or_nan() || bits.is_zero())
164     return x;
165 
166   bool is_neg = bits.is_neg();
167   int exponent = bits.get_exponent();
168 
169   // If the exponent is greater than the most negative mantissa
170   // exponent, then x is already an integer.
171   if (exponent >= static_cast<int>(FPBits<T>::FRACTION_LEN))
172     return x;
173 
174   if (exponent <= -1) {
175     switch (rnd) {
176     case FP_INT_DOWNWARD:
177       return is_neg ? T(-1.0) : T(0.0);
178     case FP_INT_UPWARD:
179       return is_neg ? T(-0.0) : T(1.0);
180     case FP_INT_TOWARDZERO:
181       return is_neg ? T(-0.0) : T(0.0);
182     case FP_INT_TONEARESTFROMZERO:
183       if (exponent < -1)
184         return is_neg ? T(-0.0) : T(0.0); // abs(x) < 0.5
185       return is_neg ? T(-1.0) : T(1.0);   // abs(x) >= 0.5
186     case FP_INT_TONEAREST:
187     default:
188       if (exponent <= -2 || bits.get_mantissa() == 0)
189         return is_neg ? T(-0.0) : T(0.0); // abs(x) <= 0.5
190       else
191         return is_neg ? T(-1.0) : T(1.0); // abs(x) > 0.5
192     }
193   }
194 
195   uint32_t trim_size = FPBits<T>::FRACTION_LEN - exponent;
196   StorageType x_u = bits.uintval();
197   StorageType trunc_u =
198       static_cast<StorageType>((x_u >> trim_size) << trim_size);
199 
200   // If x is already an integer, return it.
201   if (trunc_u == x_u)
202     return x;
203 
204   FPBits<T> new_bits(trunc_u);
205   T trunc_value = new_bits.get_val();
206 
207   StorageType trim_value =
208       bits.get_mantissa() &
209       static_cast<StorageType>(((StorageType(1) << trim_size) - 1));
210   StorageType half_value =
211       static_cast<StorageType>((StorageType(1) << (trim_size - 1)));
212   // If exponent is 0, trimSize will be equal to the mantissa width, and
213   // truncIsOdd` will not be correct. So, we handle it as a special case
214   // below.
215   StorageType trunc_is_odd =
216       new_bits.get_mantissa() & (StorageType(1) << trim_size);
217 
218   switch (rnd) {
219   case FP_INT_DOWNWARD:
220     return is_neg ? trunc_value - T(1.0) : trunc_value;
221   case FP_INT_UPWARD:
222     return is_neg ? trunc_value : trunc_value + T(1.0);
223   case FP_INT_TOWARDZERO:
224     return trunc_value;
225   case FP_INT_TONEARESTFROMZERO:
226     if (trim_value >= half_value)
227       return is_neg ? trunc_value - T(1.0) : trunc_value + T(1.0);
228     return trunc_value;
229   case FP_INT_TONEAREST:
230   default:
231     if (trim_value > half_value) {
232       return is_neg ? trunc_value - T(1.0) : trunc_value + T(1.0);
233     } else if (trim_value == half_value) {
234       if (exponent == 0)
235         return is_neg ? T(-2.0) : T(2.0);
236       if (trunc_is_odd)
237         return is_neg ? trunc_value - T(1.0) : trunc_value + T(1.0);
238       else
239         return trunc_value;
240     } else {
241       return trunc_value;
242     }
243   }
244 }
245 
246 template <typename T>
247 LIBC_INLINE cpp::enable_if_t<cpp::is_floating_point_v<T>, T>
248 round_using_current_rounding_mode(T x) {
249   int rounding_mode = quick_get_round();
250 
251   switch (rounding_mode) {
252   case FE_DOWNWARD:
253     return round_using_specific_rounding_mode(x, FP_INT_DOWNWARD);
254   case FE_UPWARD:
255     return round_using_specific_rounding_mode(x, FP_INT_UPWARD);
256   case FE_TOWARDZERO:
257     return round_using_specific_rounding_mode(x, FP_INT_TOWARDZERO);
258   case FE_TONEAREST:
259     return round_using_specific_rounding_mode(x, FP_INT_TONEAREST);
260   default:
261     __builtin_unreachable();
262   }
263 }
264 
265 template <bool IsSigned, typename T>
266 LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_floating_point_v<T>, T>
267 fromfp(T x, int rnd, unsigned int width) {
268   using StorageType = typename FPBits<T>::StorageType;
269 
270   constexpr StorageType EXPLICIT_BIT =
271       FPBits<T>::SIG_MASK - FPBits<T>::FRACTION_MASK;
272 
273   if (width == 0U) {
274     raise_except_if_required(FE_INVALID);
275     return FPBits<T>::quiet_nan().get_val();
276   }
277 
278   FPBits<T> bits(x);
279 
280   if (bits.is_inf_or_nan()) {
281     raise_except_if_required(FE_INVALID);
282     return FPBits<T>::quiet_nan().get_val();
283   }
284 
285   T rounded_value = round_using_specific_rounding_mode(x, rnd);
286 
287   if constexpr (IsSigned) {
288     // T can't hold a finite number >= 2.0 * 2^EXP_BIAS.
289     if (width - 1 > FPBits<T>::EXP_BIAS)
290       return rounded_value;
291 
292     StorageType range_exp =
293         static_cast<StorageType>(width - 1 + FPBits<T>::EXP_BIAS);
294     // rounded_value < -2^(width - 1)
295     T range_min =
296         FPBits<T>::create_value(Sign::NEG, range_exp, EXPLICIT_BIT).get_val();
297     if (rounded_value < range_min) {
298       raise_except_if_required(FE_INVALID);
299       return FPBits<T>::quiet_nan().get_val();
300     }
301     // rounded_value > 2^(width - 1) - 1
302     T range_max =
303         FPBits<T>::create_value(Sign::POS, range_exp, EXPLICIT_BIT).get_val() -
304         T(1.0);
305     if (rounded_value > range_max) {
306       raise_except_if_required(FE_INVALID);
307       return FPBits<T>::quiet_nan().get_val();
308     }
309 
310     return rounded_value;
311   }
312 
313   if (rounded_value < T(0.0)) {
314     raise_except_if_required(FE_INVALID);
315     return FPBits<T>::quiet_nan().get_val();
316   }
317 
318   // T can't hold a finite number >= 2.0 * 2^EXP_BIAS.
319   if (width > FPBits<T>::EXP_BIAS)
320     return rounded_value;
321 
322   StorageType range_exp = static_cast<StorageType>(width + FPBits<T>::EXP_BIAS);
323   // rounded_value > 2^width - 1
324   T range_max =
325       FPBits<T>::create_value(Sign::POS, range_exp, EXPLICIT_BIT).get_val() -
326       T(1.0);
327   if (rounded_value > range_max) {
328     raise_except_if_required(FE_INVALID);
329     return FPBits<T>::quiet_nan().get_val();
330   }
331 
332   return rounded_value;
333 }
334 
335 template <bool IsSigned, typename T>
336 LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_floating_point_v<T>, T>
337 fromfpx(T x, int rnd, unsigned int width) {
338   T rounded_value = fromfp<IsSigned>(x, rnd, width);
339   FPBits<T> bits(rounded_value);
340 
341   if (!bits.is_nan() && rounded_value != x)
342     raise_except_if_required(FE_INEXACT);
343 
344   return rounded_value;
345 }
346 
347 namespace internal {
348 
349 template <typename FloatType, typename IntType,
350           cpp::enable_if_t<cpp::is_floating_point_v<FloatType> &&
351                                cpp::is_integral_v<IntType>,
352                            int> = 0>
353 LIBC_INLINE IntType rounded_float_to_signed_integer(FloatType x) {
354   constexpr IntType INTEGER_MIN = (IntType(1) << (sizeof(IntType) * 8 - 1));
355   constexpr IntType INTEGER_MAX = -(INTEGER_MIN + 1);
356   FPBits<FloatType> bits(x);
357   auto set_domain_error_and_raise_invalid = []() {
358     set_errno_if_required(EDOM);
359     raise_except_if_required(FE_INVALID);
360   };
361 
362   if (bits.is_inf_or_nan()) {
363     set_domain_error_and_raise_invalid();
364     return bits.is_neg() ? INTEGER_MIN : INTEGER_MAX;
365   }
366 
367   int exponent = bits.get_exponent();
368   constexpr int EXPONENT_LIMIT = sizeof(IntType) * 8 - 1;
369   if (exponent > EXPONENT_LIMIT) {
370     set_domain_error_and_raise_invalid();
371     return bits.is_neg() ? INTEGER_MIN : INTEGER_MAX;
372   } else if (exponent == EXPONENT_LIMIT) {
373     if (bits.is_pos() || bits.get_mantissa() != 0) {
374       set_domain_error_and_raise_invalid();
375       return bits.is_neg() ? INTEGER_MIN : INTEGER_MAX;
376     }
377     // If the control reaches here, then it means that the rounded
378     // value is the most negative number for the signed integer type IntType.
379   }
380 
381   // For all other cases, if `x` can fit in the integer type `IntType`,
382   // we just return `x`. static_cast will convert the floating
383   // point value to the exact integer value.
384   return static_cast<IntType>(x);
385 }
386 
387 } // namespace internal
388 
389 template <typename FloatType, typename IntType,
390           cpp::enable_if_t<cpp::is_floating_point_v<FloatType> &&
391                                cpp::is_integral_v<IntType>,
392                            int> = 0>
393 LIBC_INLINE IntType round_to_signed_integer(FloatType x) {
394   return internal::rounded_float_to_signed_integer<FloatType, IntType>(
395       round(x));
396 }
397 
398 template <typename FloatType, typename IntType,
399           cpp::enable_if_t<cpp::is_floating_point_v<FloatType> &&
400                                cpp::is_integral_v<IntType>,
401                            int> = 0>
402 LIBC_INLINE IntType
403 round_to_signed_integer_using_current_rounding_mode(FloatType x) {
404   return internal::rounded_float_to_signed_integer<FloatType, IntType>(
405       round_using_current_rounding_mode(x));
406 }
407 
408 } // namespace fputil
409 } // namespace LIBC_NAMESPACE_DECL
410 
411 #endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_NEARESTINTEGEROPERATIONS_H
412