1//===-- MPCommon.h ----------------------------------------------*- 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_MPFRWRAPPER_MPCOMMON_H
10#define LLVM_LIBC_UTILS_MPFRWRAPPER_MPCOMMON_H
11
12#include "src/__support/CPP/string.h"
13#include "src/__support/CPP/type_traits.h"
14#include "src/__support/FPUtil/FPBits.h"
15#include "src/__support/macros/config.h"
16#include "src/__support/macros/properties/types.h"
17#include "test/UnitTest/RoundingModeUtils.h"
18
19#include <stdint.h>
20
21#include "mpfr_inc.h"
22
23#ifdef LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
24extern "C" {
25int mpfr_set_float128(mpfr_ptr, float128, mpfr_rnd_t);
26float128 mpfr_get_float128(mpfr_srcptr, mpfr_rnd_t);
27}
28#endif
29
30namespace LIBC_NAMESPACE_DECL {
31namespace testing {
32namespace mpfr {
33
34template <typename T> using FPBits = LIBC_NAMESPACE::fputil::FPBits<T>;
35using LIBC_NAMESPACE::fputil::testing::RoundingMode;
36
37// A precision value which allows sufficiently large additional
38// precision compared to the floating point precision.
39template <typename T> struct ExtraPrecision;
40
41#ifdef LIBC_TYPES_HAS_FLOAT16
42template <> struct ExtraPrecision<float16> {
43 static constexpr unsigned int VALUE = 128;
44};
45#endif
46
47template <> struct ExtraPrecision<float> {
48 static constexpr unsigned int VALUE = 128;
49};
50
51template <> struct ExtraPrecision<double> {
52 static constexpr unsigned int VALUE = 256;
53};
54
55template <> struct ExtraPrecision<long double> {
56#ifdef LIBC_TYPES_LONG_DOUBLE_IS_FLOAT128
57 static constexpr unsigned int VALUE = 512;
58#else
59 static constexpr unsigned int VALUE = 256;
60#endif
61};
62
63#if defined(LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE)
64template <> struct ExtraPrecision<float128> {
65 static constexpr unsigned int VALUE = 512;
66};
67#endif // LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
68
69template <> struct ExtraPrecision<bfloat16> {
70 static constexpr unsigned int VALUE = 64;
71};
72
73// If the ulp tolerance is less than or equal to 0.5, we would check that the
74// result is rounded correctly with respect to the rounding mode by using the
75// same precision as the inputs.
76template <typename T>
77static inline unsigned int get_precision(double ulp_tolerance) {
78 if (ulp_tolerance <= 0.5) {
79 return LIBC_NAMESPACE::fputil::FPBits<T>::FRACTION_LEN + 1;
80 } else {
81 return ExtraPrecision<T>::VALUE;
82 }
83}
84
85static inline mpfr_rnd_t get_mpfr_rounding_mode(RoundingMode mode) {
86 switch (mode) {
87 case RoundingMode::Upward:
88 return MPFR_RNDU;
89 break;
90 case RoundingMode::Downward:
91 return MPFR_RNDD;
92 break;
93 case RoundingMode::TowardZero:
94 return MPFR_RNDZ;
95 break;
96 case RoundingMode::Nearest:
97 return MPFR_RNDN;
98 break;
99 }
100 __builtin_unreachable();
101}
102
103class MPFRNumber {
104 unsigned int mpfr_precision;
105 mpfr_rnd_t mpfr_rounding;
106 mpfr_t value;
107
108public:
109 MPFRNumber();
110 // We use explicit EnableIf specializations to disallow implicit
111 // conversions. Implicit conversions can potentially lead to loss of
112 // precision. We exceptionally allow implicit conversions from float16
113 // to float, as the MPFR API does not support float16, thus requiring
114 // conversion to a higher-precision format.
115 template <typename XType,
116 cpp::enable_if_t<cpp::is_same_v<float, XType>
117#ifdef LIBC_TYPES_HAS_FLOAT16
118 || cpp::is_same_v<float16, XType>
119#endif
120 || cpp::is_same_v<bfloat16, XType>,
121 int> = 0>
122 explicit MPFRNumber(XType x,
123 unsigned int precision = ExtraPrecision<XType>::VALUE,
124 RoundingMode rounding = RoundingMode::Nearest)
125 : mpfr_precision(precision),
126 mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
127 mpfr_init2(value, mpfr_precision);
128 mpfr_set_flt(value, x, mpfr_rounding);
129 }
130
131 template <typename XType,
132 cpp::enable_if_t<cpp::is_same_v<double, XType>, int> = 0>
133 explicit MPFRNumber(XType x,
134 unsigned int precision = ExtraPrecision<XType>::VALUE,
135 RoundingMode rounding = RoundingMode::Nearest)
136 : mpfr_precision(precision),
137 mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
138 mpfr_init2(value, mpfr_precision);
139 mpfr_set_d(value, x, mpfr_rounding);
140 }
141
142 template <typename XType,
143 cpp::enable_if_t<cpp::is_same_v<long double, XType>, int> = 0>
144 explicit MPFRNumber(XType x,
145 unsigned int precision = ExtraPrecision<XType>::VALUE,
146 RoundingMode rounding = RoundingMode::Nearest)
147 : mpfr_precision(precision),
148 mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
149 mpfr_init2(value, mpfr_precision);
150 mpfr_set_ld(value, x, mpfr_rounding);
151 }
152
153#ifdef LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
154 template <typename XType,
155 cpp::enable_if_t<cpp::is_same_v<float128, XType>, int> = 0>
156 explicit MPFRNumber(XType x,
157 unsigned int precision = ExtraPrecision<XType>::VALUE,
158 RoundingMode rounding = RoundingMode::Nearest)
159 : mpfr_precision(precision),
160 mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
161 mpfr_init2(value, mpfr_precision);
162 mpfr_set_float128(value, x, mpfr_rounding);
163 }
164#endif // LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
165
166 template <typename XType,
167 cpp::enable_if_t<cpp::is_integral_v<XType>, int> = 0>
168 explicit MPFRNumber(XType x,
169 unsigned int precision = ExtraPrecision<float>::VALUE,
170 RoundingMode rounding = RoundingMode::Nearest)
171 : mpfr_precision(precision),
172 mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
173 mpfr_init2(value, mpfr_precision);
174 mpfr_set_sj(value, x, mpfr_rounding);
175 }
176
177 MPFRNumber(const MPFRNumber &other);
178 MPFRNumber(const MPFRNumber &other, unsigned int precision);
179 MPFRNumber(const mpfr_t x, unsigned int precision, RoundingMode rounding);
180
181 ~MPFRNumber();
182
183 MPFRNumber &operator=(const MPFRNumber &rhs);
184
185 bool is_nan() const;
186 MPFRNumber abs() const;
187 MPFRNumber acos() const;
188 MPFRNumber acosh() const;
189 MPFRNumber acospi() const;
190 MPFRNumber add(const MPFRNumber &b) const;
191 MPFRNumber asin() const;
192 MPFRNumber asinh() const;
193 MPFRNumber atan() const;
194 MPFRNumber atan2(const MPFRNumber &b);
195 MPFRNumber atanh() const;
196 MPFRNumber cbrt() const;
197 MPFRNumber ceil() const;
198 MPFRNumber cos() const;
199 MPFRNumber cosh() const;
200 MPFRNumber cospi() const;
201 MPFRNumber erf() const;
202 MPFRNumber exp() const;
203 MPFRNumber exp2() const;
204 MPFRNumber exp2m1() const;
205 MPFRNumber exp10() const;
206 MPFRNumber exp10m1() const;
207 MPFRNumber expm1() const;
208 MPFRNumber div(const MPFRNumber &b) const;
209 MPFRNumber floor() const;
210 MPFRNumber fmod(const MPFRNumber &b);
211 MPFRNumber frexp(int &exp);
212 MPFRNumber hypot(const MPFRNumber &b);
213 MPFRNumber log() const;
214 MPFRNumber log2() const;
215 MPFRNumber log10() const;
216 MPFRNumber log1p() const;
217 MPFRNumber pow(const MPFRNumber &b);
218 MPFRNumber remquo(const MPFRNumber &divisor, int &quotient);
219 MPFRNumber round() const;
220 MPFRNumber roundeven() const;
221 bool round_to_long(long &result) const;
222 bool round_to_long(mpfr_rnd_t rnd, long &result) const;
223 MPFRNumber rint(mpfr_rnd_t rnd) const;
224 MPFRNumber mod_2pi() const;
225 MPFRNumber mod_pi_over_2() const;
226 MPFRNumber mod_pi_over_4() const;
227 MPFRNumber sin() const;
228 MPFRNumber sinpi() const;
229 MPFRNumber sinh() const;
230 MPFRNumber sqrt() const;
231 MPFRNumber sub(const MPFRNumber &b) const;
232 MPFRNumber tan() const;
233 MPFRNumber tanh() const;
234 MPFRNumber tanpi() const;
235 MPFRNumber trunc() const;
236 MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c);
237 MPFRNumber mul(const MPFRNumber &b);
238 cpp::string str() const;
239
240 template <typename T> T as() const;
241 void dump(const char *msg) const;
242
243 // Return the ULP (units-in-the-last-place) difference between the
244 // stored MPFR and a floating point number.
245 //
246 // We define ULP difference as follows:
247 // If exponents of this value and the |input| are same, then:
248 // ULP(this_value, input) = abs(this_value - input) / eps(input)
249 // else:
250 // max = max(abs(this_value), abs(input))
251 // min = min(abs(this_value), abs(input))
252 // maxExponent = exponent(max)
253 // ULP(this_value, input) = (max - 2^maxExponent) / eps(max) +
254 // (2^maxExponent - min) / eps(min)
255 //
256 // Remarks:
257 // 1. A ULP of 0.0 will imply that the value is correctly rounded.
258 // 2. We expect that this value and the value to be compared (the [input]
259 // argument) are reasonable close, and we will provide an upper bound
260 // of ULP value for testing. Morever, most of the fractional parts of
261 // ULP value do not matter much, so using double as the return type
262 // should be good enough.
263 // 3. For close enough values (values which don't diff in their exponent by
264 // not more than 1), a ULP difference of N indicates a bit distance
265 // of N between this number and [input].
266 // 4. A values of +0.0 and -0.0 are treated as equal.
267 template <typename T>
268 cpp::enable_if_t<cpp::is_floating_point_v<T>, MPFRNumber>
269 ulp_as_mpfr_number(T input) {
270 T thisAsT = as<T>();
271 if (thisAsT == input)
272 return MPFRNumber(0.0);
273
274 if (is_nan()) {
275 if (FPBits<T>(input).is_nan())
276 return MPFRNumber(0.0);
277 return MPFRNumber(FPBits<T>::inf().get_val());
278 }
279
280 int thisExponent = FPBits<T>(thisAsT).get_exponent();
281 int inputExponent = FPBits<T>(input).get_exponent();
282 // Adjust the exponents for denormal numbers.
283 if (FPBits<T>(thisAsT).is_subnormal())
284 ++thisExponent;
285 if (FPBits<T>(input).is_subnormal())
286 ++inputExponent;
287
288 if (thisAsT * input < 0 || thisExponent == inputExponent) {
289 MPFRNumber inputMPFR(input);
290 mpfr_sub(inputMPFR.value, value, inputMPFR.value, MPFR_RNDN);
291 mpfr_abs(inputMPFR.value, inputMPFR.value, MPFR_RNDN);
292 mpfr_mul_2si(inputMPFR.value, inputMPFR.value,
293 -thisExponent + FPBits<T>::FRACTION_LEN, MPFR_RNDN);
294 return inputMPFR;
295 }
296
297 // If the control reaches here, it means that this number and input are
298 // of the same sign but different exponent. In such a case, ULP error is
299 // calculated as sum of two parts.
300 thisAsT = FPBits<T>(thisAsT).abs().get_val();
301 input = FPBits<T>(input).abs().get_val();
302 T min = thisAsT > input ? input : thisAsT;
303 T max = thisAsT > input ? thisAsT : input;
304 int minExponent = FPBits<T>(min).get_exponent();
305 int maxExponent = FPBits<T>(max).get_exponent();
306 // Adjust the exponents for denormal numbers.
307 if (FPBits<T>(min).is_subnormal())
308 ++minExponent;
309 if (FPBits<T>(max).is_subnormal())
310 ++maxExponent;
311
312 MPFRNumber minMPFR(min);
313 MPFRNumber maxMPFR(max);
314
315 MPFRNumber pivot(uint32_t(1));
316 mpfr_mul_2si(pivot.value, pivot.value, maxExponent, MPFR_RNDN);
317
318 mpfr_sub(minMPFR.value, pivot.value, minMPFR.value, MPFR_RNDN);
319 mpfr_mul_2si(minMPFR.value, minMPFR.value,
320 -minExponent + FPBits<T>::FRACTION_LEN, MPFR_RNDN);
321
322 mpfr_sub(maxMPFR.value, maxMPFR.value, pivot.value, MPFR_RNDN);
323 mpfr_mul_2si(maxMPFR.value, maxMPFR.value,
324 -maxExponent + FPBits<T>::FRACTION_LEN, MPFR_RNDN);
325
326 mpfr_add(minMPFR.value, minMPFR.value, maxMPFR.value, MPFR_RNDN);
327 return minMPFR;
328 }
329
330 template <typename T>
331 cpp::enable_if_t<cpp::is_floating_point_v<T>, cpp::string>
332 ulp_as_string(T input) {
333 MPFRNumber num = ulp_as_mpfr_number(input);
334 return num.str();
335 }
336
337 template <typename T>
338 cpp::enable_if_t<cpp::is_floating_point_v<T>, double> ulp(T input) {
339 MPFRNumber num = ulp_as_mpfr_number(input);
340 return num.as<double>();
341 }
342};
343
344} // namespace mpfr
345} // namespace testing
346} // namespace LIBC_NAMESPACE_DECL
347
348#endif // LLVM_LIBC_UTILS_MPFRWRAPPER_MPCOMMON_H
349

source code of libc/utils/MPFRWrapper/MPCommon.h