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>
trunc(T x)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>
ceil(T x)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>
floor(T x)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>
round(T x)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>
round_using_current_rounding_mode(T x)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>
rounded_float_to_signed_integer(F x)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>
round_to_signed_integer(F x)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>
round_to_signed_integer_using_current_rounding_mode(F x)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