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.isInfOrNaN())
35     return x;
36 
37   int exponent = bits.getExponent();
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.getSign())
47       return T(-0.0);
48     else
49       return T(0.0);
50   }
51 
52   int trimSize = MantissaWidth<T>::value - exponent;
53   bits.setMantissa((bits.getMantissa() >> trimSize) << trimSize);
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.isInfOrNaN() || bits.isZero())
64     return x;
65 
66   bool isNeg = bits.getSign();
67   int exponent = bits.getExponent();
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 (isNeg)
76       return T(-0.0);
77     else
78       return T(1.0);
79   }
80 
81   uint32_t trimSize = MantissaWidth<T>::value - exponent;
82   bits.setMantissa((bits.getMantissa() >> trimSize) << trimSize);
83   T truncValue = T(bits);
84 
85   // If x is already an integer, return it.
86   if (truncValue == x)
87     return x;
88 
89   // If x is negative, the ceil operation is equivalent to the trunc operation.
90   if (isNeg)
91     return truncValue;
92 
93   return truncValue + 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.getSign()) {
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.isInfOrNaN() || bits.isZero())
115     return x;
116 
117   bool isNeg = bits.getSign();
118   int exponent = bits.getExponent();
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 (isNeg)
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 (isNeg)
136       return T(-0.0);
137     else
138       return T(0.0);
139   }
140 
141   uint32_t trimSize = MantissaWidth<T>::value - exponent;
142   bool halfBitSet = bits.getMantissa() & (UIntType(1) << (trimSize - 1));
143   bits.setMantissa((bits.getMantissa() >> trimSize) << trimSize);
144   T truncValue = T(bits);
145 
146   // If x is already an integer, return it.
147   if (truncValue == x)
148     return x;
149 
150   if (!halfBitSet) {
151     // Franctional part is less than 0.5 so round value is the
152     // same as the trunc value.
153     return truncValue;
154   } else {
155     return isNeg ? truncValue - T(1.0) : truncValue + T(1.0);
156   }
157 }
158 
159 template <typename T,
160           cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, int> = 0>
161 static inline T roundUsingCurrentRoundingMode(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.isInfOrNaN() || bits.isZero())
167     return x;
168 
169   bool isNeg = bits.getSign();
170   int exponent = bits.getExponent();
171   int roundingMode = getRound();
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 (roundingMode) {
180     case FE_DOWNWARD:
181       return isNeg ? T(-1.0) : T(0.0);
182     case FE_UPWARD:
183       return isNeg ? T(-0.0) : T(1.0);
184     case FE_TOWARDZERO:
185       return isNeg ? T(-0.0) : T(0.0);
186     case FE_TONEAREST:
187       if (exponent <= -2 || bits.getMantissa() == 0)
188         return isNeg ? T(-0.0) : T(0.0); // abs(x) <= 0.5
189       else
190         return isNeg ? T(-1.0) : T(1.0); // abs(x) > 0.5
191     default:
192       __builtin_unreachable();
193     }
194   }
195 
196   uint32_t trimSize = MantissaWidth<T>::value - exponent;
197   FPBits<T> newBits = bits;
198   newBits.setMantissa((bits.getMantissa() >> trimSize) << trimSize);
199   T truncValue = T(newBits);
200 
201   // If x is already an integer, return it.
202   if (truncValue == x)
203     return x;
204 
205   UIntType trimValue = bits.getMantissa() & ((UIntType(1) << trimSize) - 1);
206   UIntType halfValue = (UIntType(1) << (trimSize - 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 truncIsOdd = newBits.getMantissa() & (UIntType(1) << trimSize);
211 
212   switch (roundingMode) {
213   case FE_DOWNWARD:
214     return isNeg ? truncValue - T(1.0) : truncValue;
215   case FE_UPWARD:
216     return isNeg ? truncValue : truncValue + T(1.0);
217   case FE_TOWARDZERO:
218     return truncValue;
219   case FE_TONEAREST:
220     if (trimValue > halfValue) {
221       return isNeg ? truncValue - T(1.0) : truncValue + T(1.0);
222     } else if (trimValue == halfValue) {
223       if (exponent == 0)
224         return isNeg ? T(-2.0) : T(2.0);
225       if (truncIsOdd)
226         return isNeg ? truncValue - T(1.0) : truncValue + T(1.0);
227       else
228         return truncValue;
229     } else {
230       return truncValue;
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 roundedFloatToSignedInteger(F x) {
244   constexpr I IntegerMin = (I(1) << (sizeof(I) * 8 - 1));
245   constexpr I IntegerMax = -(IntegerMin + 1);
246   FPBits<F> bits(x);
247   auto setDomainErrorAndRaiseInvalid = []() {
248 #if math_errhandling & MATH_ERRNO
249     errno = EDOM;
250 #endif
251 #if math_errhandling & MATH_ERREXCEPT
252     raiseExcept(FE_INVALID);
253 #endif
254   };
255 
256   if (bits.isInfOrNaN()) {
257     setDomainErrorAndRaiseInvalid();
258     return bits.getSign() ? IntegerMin : IntegerMax;
259   }
260 
261   int exponent = bits.getExponent();
262   constexpr int exponentLimit = sizeof(I) * 8 - 1;
263   if (exponent > exponentLimit) {
264     setDomainErrorAndRaiseInvalid();
265     return bits.getSign() ? IntegerMin : IntegerMax;
266   } else if (exponent == exponentLimit) {
267     if (bits.getSign() == 0 || bits.getMantissa() != 0) {
268       setDomainErrorAndRaiseInvalid();
269       return bits.getSign() ? IntegerMin : IntegerMax;
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 roundToSignedInteger(F x) {
288   return internal::roundedFloatToSignedInteger<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 roundToSignedIntegerUsingCurrentRoundingMode(F x) {
296   return internal::roundedFloatToSignedInteger<F, I>(
297       roundUsingCurrentRoundingMode(x));
298 }
299 
300 } // namespace fputil
301 } // namespace __llvm_libc
302 
303 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_NEAREST_INTEGER_OPERATIONS_H
304