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_NEAREST_INTEGER_OPERATIONS_H
10 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_NEAREST_INTEGER_OPERATIONS_H
11 
12 #include "FEnvImpl.h"
13 #include "FPBits.h"
14 
15 #include "src/__support/CPP/TypeTraits.h"
16 
17 #include <errno.h>
18 #include <math.h>
19 
20 namespace __llvm_libc {
21 namespace fputil {
22 
23 template <typename T,
24           cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, int> = 0>
25 static inline T trunc(T x) {
26   FPBits<T> bits(x);
27 
28   // If x is infinity or NaN, return it.
29   // If it is zero also we should return it as is, but the logic
30   // later in this function takes care of it. But not doing a zero
31   // check, we improve the run time of non-zero values.
32   if (bits.is_inf_or_nan())
33     return x;
34 
35   int exponent = bits.get_exponent();
36 
37   // If the exponent is greater than the most negative mantissa
38   // exponent, then x is already an integer.
39   if (exponent >= static_cast<int>(MantissaWidth<T>::VALUE))
40     return x;
41 
42   // If the exponent is such that abs(x) is less than 1, then return 0.
43   if (exponent <= -1) {
44     if (bits.get_sign())
45       return T(-0.0);
46     else
47       return T(0.0);
48   }
49 
50   int trim_size = MantissaWidth<T>::VALUE - exponent;
51   bits.set_mantissa((bits.get_mantissa() >> trim_size) << trim_size);
52   return T(bits);
53 }
54 
55 template <typename T,
56           cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, int> = 0>
57 static inline T ceil(T x) {
58   FPBits<T> bits(x);
59 
60   // If x is infinity NaN or zero, return it.
61   if (bits.is_inf_or_nan() || bits.is_zero())
62     return x;
63 
64   bool is_neg = bits.get_sign();
65   int exponent = bits.get_exponent();
66 
67   // If the exponent is greater than the most negative mantissa
68   // exponent, then x is already an integer.
69   if (exponent >= static_cast<int>(MantissaWidth<T>::VALUE))
70     return x;
71 
72   if (exponent <= -1) {
73     if (is_neg)
74       return T(-0.0);
75     else
76       return T(1.0);
77   }
78 
79   uint32_t trim_size = MantissaWidth<T>::VALUE - exponent;
80   bits.set_mantissa((bits.get_mantissa() >> trim_size) << trim_size);
81   T trunc_value = T(bits);
82 
83   // If x is already an integer, return it.
84   if (trunc_value == x)
85     return x;
86 
87   // If x is negative, the ceil operation is equivalent to the trunc operation.
88   if (is_neg)
89     return trunc_value;
90 
91   return trunc_value + T(1.0);
92 }
93 
94 template <typename T,
95           cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, int> = 0>
96 static inline T floor(T x) {
97   FPBits<T> bits(x);
98   if (bits.get_sign()) {
99     return -ceil(-x);
100   } else {
101     return trunc(x);
102   }
103 }
104 
105 template <typename T,
106           cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, int> = 0>
107 static inline T round(T x) {
108   using UIntType = typename FPBits<T>::UIntType;
109   FPBits<T> bits(x);
110 
111   // If x is infinity NaN or zero, return it.
112   if (bits.is_inf_or_nan() || bits.is_zero())
113     return x;
114 
115   bool is_neg = bits.get_sign();
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>(MantissaWidth<T>::VALUE))
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     if (is_neg)
126       return T(-1.0);
127     else
128       return T(1.0);
129   }
130 
131   if (exponent <= -2) {
132     // Absolute value of x is less than 0.5.
133     if (is_neg)
134       return T(-0.0);
135     else
136       return T(0.0);
137   }
138 
139   uint32_t trim_size = MantissaWidth<T>::VALUE - exponent;
140   bool half_bit_set = bits.get_mantissa() & (UIntType(1) << (trim_size - 1));
141   bits.set_mantissa((bits.get_mantissa() >> trim_size) << trim_size);
142   T trunc_value = T(bits);
143 
144   // If x is already an integer, return it.
145   if (trunc_value == x)
146     return x;
147 
148   if (!half_bit_set) {
149     // Franctional part is less than 0.5 so round value is the
150     // same as the trunc value.
151     return trunc_value;
152   } else {
153     return is_neg ? trunc_value - T(1.0) : trunc_value + T(1.0);
154   }
155 }
156 
157 template <typename T,
158           cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, int> = 0>
159 static inline T round_using_current_rounding_mode(T x) {
160   using UIntType = typename FPBits<T>::UIntType;
161   FPBits<T> bits(x);
162 
163   // If x is infinity NaN or zero, return it.
164   if (bits.is_inf_or_nan() || bits.is_zero())
165     return x;
166 
167   bool is_neg = bits.get_sign();
168   int exponent = bits.get_exponent();
169   int rounding_mode = get_round();
170 
171   // If the exponent is greater than the most negative mantissa
172   // exponent, then x is already an integer.
173   if (exponent >= static_cast<int>(MantissaWidth<T>::VALUE))
174     return x;
175 
176   if (exponent <= -1) {
177     switch (rounding_mode) {
178     case FE_DOWNWARD:
179       return is_neg ? T(-1.0) : T(0.0);
180     case FE_UPWARD:
181       return is_neg ? T(-0.0) : T(1.0);
182     case FE_TOWARDZERO:
183       return is_neg ? T(-0.0) : T(0.0);
184     case FE_TONEAREST:
185       if (exponent <= -2 || bits.get_mantissa() == 0)
186         return is_neg ? T(-0.0) : T(0.0); // abs(x) <= 0.5
187       else
188         return is_neg ? T(-1.0) : T(1.0); // abs(x) > 0.5
189     default:
190       __builtin_unreachable();
191     }
192   }
193 
194   uint32_t trim_size = MantissaWidth<T>::VALUE - exponent;
195   FPBits<T> new_bits = bits;
196   new_bits.set_mantissa((bits.get_mantissa() >> trim_size) << trim_size);
197   T trunc_value = T(new_bits);
198 
199   // If x is already an integer, return it.
200   if (trunc_value == x)
201     return x;
202 
203   UIntType trim_value = bits.get_mantissa() & ((UIntType(1) << trim_size) - 1);
204   UIntType half_value = (UIntType(1) << (trim_size - 1));
205   // If exponent is 0, trimSize will be equal to the mantissa width, and
206   // truncIsOdd` will not be correct. So, we handle it as a special case
207   // below.
208   UIntType trunc_is_odd = new_bits.get_mantissa() & (UIntType(1) << trim_size);
209 
210   switch (rounding_mode) {
211   case FE_DOWNWARD:
212     return is_neg ? trunc_value - T(1.0) : trunc_value;
213   case FE_UPWARD:
214     return is_neg ? trunc_value : trunc_value + T(1.0);
215   case FE_TOWARDZERO:
216     return trunc_value;
217   case FE_TONEAREST:
218     if (trim_value > half_value) {
219       return is_neg ? trunc_value - T(1.0) : trunc_value + T(1.0);
220     } else if (trim_value == half_value) {
221       if (exponent == 0)
222         return is_neg ? T(-2.0) : T(2.0);
223       if (trunc_is_odd)
224         return is_neg ? trunc_value - T(1.0) : trunc_value + T(1.0);
225       else
226         return trunc_value;
227     } else {
228       return trunc_value;
229     }
230   default:
231     __builtin_unreachable();
232   }
233 }
234 
235 namespace internal {
236 
237 template <typename F, typename I,
238           cpp::EnableIfType<cpp::IsFloatingPointType<F>::Value &&
239                                 cpp::IsIntegral<I>::Value,
240                             int> = 0>
241 static inline I rounded_float_to_signed_integer(F x) {
242   constexpr I INTEGER_MIN = (I(1) << (sizeof(I) * 8 - 1));
243   constexpr I INTEGER_MAX = -(INTEGER_MIN + 1);
244   FPBits<F> bits(x);
245   auto set_domain_error_and_raise_invalid = []() {
246     if (math_errhandling & MATH_ERRNO)
247       errno = EDOM;
248     if (math_errhandling & MATH_ERREXCEPT)
249       raise_except(FE_INVALID);
250   };
251 
252   if (bits.is_inf_or_nan()) {
253     set_domain_error_and_raise_invalid();
254     return bits.get_sign() ? INTEGER_MIN : INTEGER_MAX;
255   }
256 
257   int exponent = bits.get_exponent();
258   constexpr int EXPONENT_LIMIT = sizeof(I) * 8 - 1;
259   if (exponent > EXPONENT_LIMIT) {
260     set_domain_error_and_raise_invalid();
261     return bits.get_sign() ? INTEGER_MIN : INTEGER_MAX;
262   } else if (exponent == EXPONENT_LIMIT) {
263     if (bits.get_sign() == 0 || bits.get_mantissa() != 0) {
264       set_domain_error_and_raise_invalid();
265       return bits.get_sign() ? INTEGER_MIN : INTEGER_MAX;
266     }
267     // If the control reaches here, then it means that the rounded
268     // value is the most negative number for the signed integer type I.
269   }
270 
271   // For all other cases, if `x` can fit in the integer type `I`,
272   // we just return `x`. Implicit conversion will convert the
273   // floating point value to the exact integer value.
274   return x;
275 }
276 
277 } // namespace internal
278 
279 template <typename F, typename I,
280           cpp::EnableIfType<cpp::IsFloatingPointType<F>::Value &&
281                                 cpp::IsIntegral<I>::Value,
282                             int> = 0>
283 static inline I round_to_signed_integer(F x) {
284   return internal::rounded_float_to_signed_integer<F, I>(round(x));
285 }
286 
287 template <typename F, typename I,
288           cpp::EnableIfType<cpp::IsFloatingPointType<F>::Value &&
289                                 cpp::IsIntegral<I>::Value,
290                             int> = 0>
291 static inline I round_to_signed_integer_using_current_rounding_mode(F x) {
292   return internal::rounded_float_to_signed_integer<F, I>(
293       round_using_current_rounding_mode(x));
294 }
295 
296 } // namespace fputil
297 } // namespace __llvm_libc
298 
299 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_NEAREST_INTEGER_OPERATIONS_H
300