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

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