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