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