1//===-- Utils which wrap MPFR ---------------------------------------------===//
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#include "MPFRUtils.h"
10
11#include "src/__support/CPP/string.h"
12#include "src/__support/CPP/string_view.h"
13#include "src/__support/FPUtil/FPBits.h"
14#include "src/__support/FPUtil/fpbits_str.h"
15#include "test/UnitTest/FPMatcher.h"
16
17#include "hdr/math_macros.h"
18#include <memory>
19#include <stdint.h>
20
21#include "mpfr_inc.h"
22
23template <typename T> using FPBits = LIBC_NAMESPACE::fputil::FPBits<T>;
24
25namespace LIBC_NAMESPACE {
26namespace testing {
27namespace mpfr {
28
29// A precision value which allows sufficiently large additional
30// precision compared to the floating point precision.
31template <typename T> struct ExtraPrecision;
32
33template <> struct ExtraPrecision<float> {
34 static constexpr unsigned int VALUE = 128;
35};
36
37template <> struct ExtraPrecision<double> {
38 static constexpr unsigned int VALUE = 256;
39};
40
41template <> struct ExtraPrecision<long double> {
42 static constexpr unsigned int VALUE = 256;
43};
44
45// If the ulp tolerance is less than or equal to 0.5, we would check that the
46// result is rounded correctly with respect to the rounding mode by using the
47// same precision as the inputs.
48template <typename T>
49static inline unsigned int get_precision(double ulp_tolerance) {
50 if (ulp_tolerance <= 0.5) {
51 return LIBC_NAMESPACE::fputil::FPBits<T>::FRACTION_LEN + 1;
52 } else {
53 return ExtraPrecision<T>::VALUE;
54 }
55}
56
57static inline mpfr_rnd_t get_mpfr_rounding_mode(RoundingMode mode) {
58 switch (mode) {
59 case RoundingMode::Upward:
60 return MPFR_RNDU;
61 break;
62 case RoundingMode::Downward:
63 return MPFR_RNDD;
64 break;
65 case RoundingMode::TowardZero:
66 return MPFR_RNDZ;
67 break;
68 case RoundingMode::Nearest:
69 return MPFR_RNDN;
70 break;
71 }
72 __builtin_unreachable();
73}
74
75class MPFRNumber {
76 unsigned int mpfr_precision;
77 mpfr_rnd_t mpfr_rounding;
78
79 mpfr_t value;
80
81public:
82 MPFRNumber() : mpfr_precision(256), mpfr_rounding(MPFR_RNDN) {
83 mpfr_init2(value, mpfr_precision);
84 }
85
86 // We use explicit EnableIf specializations to disallow implicit
87 // conversions. Implicit conversions can potentially lead to loss of
88 // precision.
89 template <typename XType,
90 cpp::enable_if_t<cpp::is_same_v<float, XType>, int> = 0>
91 explicit MPFRNumber(XType x,
92 unsigned int precision = ExtraPrecision<XType>::VALUE,
93 RoundingMode rounding = RoundingMode::Nearest)
94 : mpfr_precision(precision),
95 mpfr_rounding(get_mpfr_rounding_mode(mode: rounding)) {
96 mpfr_init2(value, mpfr_precision);
97 mpfr_set_flt(value, x, mpfr_rounding);
98 }
99
100 template <typename XType,
101 cpp::enable_if_t<cpp::is_same_v<double, XType>, int> = 0>
102 explicit MPFRNumber(XType x,
103 unsigned int precision = ExtraPrecision<XType>::VALUE,
104 RoundingMode rounding = RoundingMode::Nearest)
105 : mpfr_precision(precision),
106 mpfr_rounding(get_mpfr_rounding_mode(mode: rounding)) {
107 mpfr_init2(value, mpfr_precision);
108 mpfr_set_d(value, x, mpfr_rounding);
109 }
110
111 template <typename XType,
112 cpp::enable_if_t<cpp::is_same_v<long double, XType>, int> = 0>
113 explicit MPFRNumber(XType x,
114 unsigned int precision = ExtraPrecision<XType>::VALUE,
115 RoundingMode rounding = RoundingMode::Nearest)
116 : mpfr_precision(precision),
117 mpfr_rounding(get_mpfr_rounding_mode(mode: rounding)) {
118 mpfr_init2(value, mpfr_precision);
119 mpfr_set_ld(value, x, mpfr_rounding);
120 }
121
122 template <typename XType,
123 cpp::enable_if_t<cpp::is_integral_v<XType>, int> = 0>
124 explicit MPFRNumber(XType x,
125 unsigned int precision = ExtraPrecision<float>::VALUE,
126 RoundingMode rounding = RoundingMode::Nearest)
127 : mpfr_precision(precision),
128 mpfr_rounding(get_mpfr_rounding_mode(mode: rounding)) {
129 mpfr_init2(value, mpfr_precision);
130 mpfr_set_sj(value, x, mpfr_rounding);
131 }
132
133 MPFRNumber(const MPFRNumber &other)
134 : mpfr_precision(other.mpfr_precision),
135 mpfr_rounding(other.mpfr_rounding) {
136 mpfr_init2(value, mpfr_precision);
137 mpfr_set(value, other.value, mpfr_rounding);
138 }
139
140 MPFRNumber(const MPFRNumber &other, unsigned int precision)
141 : mpfr_precision(precision), mpfr_rounding(other.mpfr_rounding) {
142 mpfr_init2(value, mpfr_precision);
143 mpfr_set(value, other.value, mpfr_rounding);
144 }
145
146 ~MPFRNumber() { mpfr_clear(value); }
147
148 MPFRNumber &operator=(const MPFRNumber &rhs) {
149 mpfr_precision = rhs.mpfr_precision;
150 mpfr_rounding = rhs.mpfr_rounding;
151 mpfr_set(value, rhs.value, mpfr_rounding);
152 return *this;
153 }
154
155 bool is_nan() const { return mpfr_nan_p(value); }
156
157 MPFRNumber abs() const {
158 MPFRNumber result(*this);
159 mpfr_abs(result.value, value, mpfr_rounding);
160 return result;
161 }
162
163 MPFRNumber acos() const {
164 MPFRNumber result(*this);
165 mpfr_acos(result.value, value, mpfr_rounding);
166 return result;
167 }
168
169 MPFRNumber acosh() const {
170 MPFRNumber result(*this);
171 mpfr_acosh(result.value, value, mpfr_rounding);
172 return result;
173 }
174
175 MPFRNumber asin() const {
176 MPFRNumber result(*this);
177 mpfr_asin(result.value, value, mpfr_rounding);
178 return result;
179 }
180
181 MPFRNumber asinh() const {
182 MPFRNumber result(*this);
183 mpfr_asinh(result.value, value, mpfr_rounding);
184 return result;
185 }
186
187 MPFRNumber atan() const {
188 MPFRNumber result(*this);
189 mpfr_atan(result.value, value, mpfr_rounding);
190 return result;
191 }
192
193 MPFRNumber atan2(const MPFRNumber &b) {
194 MPFRNumber result(*this);
195 mpfr_atan2(result.value, value, b.value, mpfr_rounding);
196 return result;
197 }
198
199 MPFRNumber atanh() const {
200 MPFRNumber result(*this);
201 mpfr_atanh(result.value, value, mpfr_rounding);
202 return result;
203 }
204
205 MPFRNumber ceil() const {
206 MPFRNumber result(*this);
207 mpfr_ceil(result.value, value);
208 return result;
209 }
210
211 MPFRNumber cos() const {
212 MPFRNumber result(*this);
213 mpfr_cos(result.value, value, mpfr_rounding);
214 return result;
215 }
216
217 MPFRNumber cosh() const {
218 MPFRNumber result(*this);
219 mpfr_cosh(result.value, value, mpfr_rounding);
220 return result;
221 }
222
223 MPFRNumber erf() const {
224 MPFRNumber result(*this);
225 mpfr_erf(result.value, value, mpfr_rounding);
226 return result;
227 }
228
229 MPFRNumber exp() const {
230 MPFRNumber result(*this);
231 mpfr_exp(result.value, value, mpfr_rounding);
232 return result;
233 }
234
235 MPFRNumber exp2() const {
236 MPFRNumber result(*this);
237 mpfr_exp2(result.value, value, mpfr_rounding);
238 return result;
239 }
240
241 MPFRNumber exp2m1() const {
242 // TODO: Only use mpfr_exp2m1 once CI and buildbots get MPFR >= 4.2.0.
243#if MPFR_VERSION_MAJOR > 4 || \
244 (MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2)
245 MPFRNumber result(*this);
246 mpfr_exp2m1(result.value, value, mpfr_rounding);
247 return result;
248#else
249 unsigned int prec = mpfr_precision * 3;
250 MPFRNumber result(*this, prec);
251
252 float f = mpfr_get_flt(abs().value, mpfr_rounding);
253 if (f > 0.5f && f < 0x1.0p30f) {
254 mpfr_exp2(result.value, value, mpfr_rounding);
255 mpfr_sub_ui(result.value, result.value, 1, mpfr_rounding);
256 return result;
257 }
258
259 MPFRNumber ln2(2.0f, prec);
260 // log(2)
261 mpfr_log(ln2.value, ln2.value, mpfr_rounding);
262 // x * log(2)
263 mpfr_mul(result.value, value, ln2.value, mpfr_rounding);
264 // e^(x * log(2)) - 1
265 int ex = mpfr_expm1(result.value, result.value, mpfr_rounding);
266 mpfr_subnormalize(result.value, ex, mpfr_rounding);
267 return result;
268#endif
269 }
270
271 MPFRNumber exp10() const {
272 MPFRNumber result(*this);
273 mpfr_exp10(result.value, value, mpfr_rounding);
274 return result;
275 }
276
277 MPFRNumber expm1() const {
278 MPFRNumber result(*this);
279 mpfr_expm1(result.value, value, mpfr_rounding);
280 return result;
281 }
282
283 MPFRNumber floor() const {
284 MPFRNumber result(*this);
285 mpfr_floor(result.value, value);
286 return result;
287 }
288
289 MPFRNumber fmod(const MPFRNumber &b) {
290 MPFRNumber result(*this);
291 mpfr_fmod(result.value, value, b.value, mpfr_rounding);
292 return result;
293 }
294
295 MPFRNumber frexp(int &exp) {
296 MPFRNumber result(*this);
297 mpfr_exp_t resultExp;
298 mpfr_frexp(&resultExp, result.value, value, mpfr_rounding);
299 exp = resultExp;
300 return result;
301 }
302
303 MPFRNumber hypot(const MPFRNumber &b) {
304 MPFRNumber result(*this);
305 mpfr_hypot(result.value, value, b.value, mpfr_rounding);
306 return result;
307 }
308
309 MPFRNumber log() const {
310 MPFRNumber result(*this);
311 mpfr_log(result.value, value, mpfr_rounding);
312 return result;
313 }
314
315 MPFRNumber log2() const {
316 MPFRNumber result(*this);
317 mpfr_log2(result.value, value, mpfr_rounding);
318 return result;
319 }
320
321 MPFRNumber log10() const {
322 MPFRNumber result(*this);
323 mpfr_log10(result.value, value, mpfr_rounding);
324 return result;
325 }
326
327 MPFRNumber log1p() const {
328 MPFRNumber result(*this);
329 mpfr_log1p(result.value, value, mpfr_rounding);
330 return result;
331 }
332
333 MPFRNumber pow(const MPFRNumber &b) {
334 MPFRNumber result(*this);
335 mpfr_pow(result.value, value, b.value, mpfr_rounding);
336 return result;
337 }
338
339 MPFRNumber remquo(const MPFRNumber &divisor, int &quotient) {
340 MPFRNumber remainder(*this);
341 long q;
342 mpfr_remquo(remainder.value, &q, value, divisor.value, mpfr_rounding);
343 quotient = q;
344 return remainder;
345 }
346
347 MPFRNumber round() const {
348 MPFRNumber result(*this);
349 mpfr_round(result.value, value);
350 return result;
351 }
352
353 MPFRNumber roundeven() const {
354 MPFRNumber result(*this);
355#if MPFR_VERSION_MAJOR >= 4
356 mpfr_roundeven(result.value, value);
357#else
358 mpfr_rint(result.value, value, MPFR_RNDN);
359#endif
360 return result;
361 }
362
363 bool round_to_long(long &result) const {
364 // We first calculate the rounded value. This way, when converting
365 // to long using mpfr_get_si, the rounding direction of MPFR_RNDN
366 // (or any other rounding mode), does not have an influence.
367 MPFRNumber roundedValue = round();
368 mpfr_clear_erangeflag();
369 result = mpfr_get_si(roundedValue.value, MPFR_RNDN);
370 return mpfr_erangeflag_p();
371 }
372
373 bool round_to_long(mpfr_rnd_t rnd, long &result) const {
374 MPFRNumber rint_result(*this);
375 mpfr_rint(rint_result.value, value, rnd);
376 return rint_result.round_to_long(result);
377 }
378
379 MPFRNumber rint(mpfr_rnd_t rnd) const {
380 MPFRNumber result(*this);
381 mpfr_rint(result.value, value, rnd);
382 return result;
383 }
384
385 MPFRNumber mod_2pi() const {
386 MPFRNumber result(0.0, 1280);
387 MPFRNumber _2pi(0.0, 1280);
388 mpfr_const_pi(_2pi.value, MPFR_RNDN);
389 mpfr_mul_si(_2pi.value, _2pi.value, 2, MPFR_RNDN);
390 mpfr_fmod(result.value, value, _2pi.value, MPFR_RNDN);
391 return result;
392 }
393
394 MPFRNumber mod_pi_over_2() const {
395 MPFRNumber result(0.0, 1280);
396 MPFRNumber pi_over_2(0.0, 1280);
397 mpfr_const_pi(pi_over_2.value, MPFR_RNDN);
398 mpfr_mul_d(pi_over_2.value, pi_over_2.value, 0.5, MPFR_RNDN);
399 mpfr_fmod(result.value, value, pi_over_2.value, MPFR_RNDN);
400 return result;
401 }
402
403 MPFRNumber mod_pi_over_4() const {
404 MPFRNumber result(0.0, 1280);
405 MPFRNumber pi_over_4(0.0, 1280);
406 mpfr_const_pi(pi_over_4.value, MPFR_RNDN);
407 mpfr_mul_d(pi_over_4.value, pi_over_4.value, 0.25, MPFR_RNDN);
408 mpfr_fmod(result.value, value, pi_over_4.value, MPFR_RNDN);
409 return result;
410 }
411
412 MPFRNumber sin() const {
413 MPFRNumber result(*this);
414 mpfr_sin(result.value, value, mpfr_rounding);
415 return result;
416 }
417
418 MPFRNumber sinh() const {
419 MPFRNumber result(*this);
420 mpfr_sinh(result.value, value, mpfr_rounding);
421 return result;
422 }
423
424 MPFRNumber sqrt() const {
425 MPFRNumber result(*this);
426 mpfr_sqrt(result.value, value, mpfr_rounding);
427 return result;
428 }
429
430 MPFRNumber tan() const {
431 MPFRNumber result(*this);
432 mpfr_tan(result.value, value, mpfr_rounding);
433 return result;
434 }
435
436 MPFRNumber tanh() const {
437 MPFRNumber result(*this);
438 mpfr_tanh(result.value, value, mpfr_rounding);
439 return result;
440 }
441
442 MPFRNumber trunc() const {
443 MPFRNumber result(*this);
444 mpfr_trunc(result.value, value);
445 return result;
446 }
447
448 MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c) {
449 MPFRNumber result(*this);
450 mpfr_fma(result.value, value, b.value, c.value, mpfr_rounding);
451 return result;
452 }
453
454 cpp::string str() const {
455 // 200 bytes should be more than sufficient to hold a 100-digit number
456 // plus additional bytes for the decimal point, '-' sign etc.
457 constexpr size_t printBufSize = 200;
458 char buffer[printBufSize];
459 mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
460 cpp::string_view view(buffer);
461 // Trim whitespaces
462 const char whitespace = ' ';
463 while (!view.empty() && view.front() == whitespace)
464 view.remove_prefix(N: 1);
465 while (!view.empty() && view.back() == whitespace)
466 view.remove_suffix(N: 1);
467 return cpp::string(view.data());
468 }
469
470 // These functions are useful for debugging.
471 template <typename T> T as() const;
472
473 void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
474
475 // Return the ULP (units-in-the-last-place) difference between the
476 // stored MPFR and a floating point number.
477 //
478 // We define ULP difference as follows:
479 // If exponents of this value and the |input| are same, then:
480 // ULP(this_value, input) = abs(this_value - input) / eps(input)
481 // else:
482 // max = max(abs(this_value), abs(input))
483 // min = min(abs(this_value), abs(input))
484 // maxExponent = exponent(max)
485 // ULP(this_value, input) = (max - 2^maxExponent) / eps(max) +
486 // (2^maxExponent - min) / eps(min)
487 //
488 // Remarks:
489 // 1. A ULP of 0.0 will imply that the value is correctly rounded.
490 // 2. We expect that this value and the value to be compared (the [input]
491 // argument) are reasonable close, and we will provide an upper bound
492 // of ULP value for testing. Morever, most of the fractional parts of
493 // ULP value do not matter much, so using double as the return type
494 // should be good enough.
495 // 3. For close enough values (values which don't diff in their exponent by
496 // not more than 1), a ULP difference of N indicates a bit distance
497 // of N between this number and [input].
498 // 4. A values of +0.0 and -0.0 are treated as equal.
499 template <typename T>
500 cpp::enable_if_t<cpp::is_floating_point_v<T>, MPFRNumber>
501 ulp_as_mpfr_number(T input) {
502 T thisAsT = as<T>();
503 if (thisAsT == input)
504 return MPFRNumber(0.0);
505
506 if (is_nan()) {
507 if (FPBits<T>(input).is_nan())
508 return MPFRNumber(0.0);
509 return MPFRNumber(FPBits<T>::inf().get_val());
510 }
511
512 int thisExponent = FPBits<T>(thisAsT).get_exponent();
513 int inputExponent = FPBits<T>(input).get_exponent();
514 // Adjust the exponents for denormal numbers.
515 if (FPBits<T>(thisAsT).is_subnormal())
516 ++thisExponent;
517 if (FPBits<T>(input).is_subnormal())
518 ++inputExponent;
519
520 if (thisAsT * input < 0 || thisExponent == inputExponent) {
521 MPFRNumber inputMPFR(input);
522 mpfr_sub(inputMPFR.value, value, inputMPFR.value, MPFR_RNDN);
523 mpfr_abs(inputMPFR.value, inputMPFR.value, MPFR_RNDN);
524 mpfr_mul_2si(inputMPFR.value, inputMPFR.value,
525 -thisExponent + FPBits<T>::FRACTION_LEN, MPFR_RNDN);
526 return inputMPFR;
527 }
528
529 // If the control reaches here, it means that this number and input are
530 // of the same sign but different exponent. In such a case, ULP error is
531 // calculated as sum of two parts.
532 thisAsT = std::abs(thisAsT);
533 input = std::abs(input);
534 T min = thisAsT > input ? input : thisAsT;
535 T max = thisAsT > input ? thisAsT : input;
536 int minExponent = FPBits<T>(min).get_exponent();
537 int maxExponent = FPBits<T>(max).get_exponent();
538 // Adjust the exponents for denormal numbers.
539 if (FPBits<T>(min).is_subnormal())
540 ++minExponent;
541 if (FPBits<T>(max).is_subnormal())
542 ++maxExponent;
543
544 MPFRNumber minMPFR(min);
545 MPFRNumber maxMPFR(max);
546
547 MPFRNumber pivot(uint32_t(1));
548 mpfr_mul_2si(pivot.value, pivot.value, maxExponent, MPFR_RNDN);
549
550 mpfr_sub(minMPFR.value, pivot.value, minMPFR.value, MPFR_RNDN);
551 mpfr_mul_2si(minMPFR.value, minMPFR.value,
552 -minExponent + FPBits<T>::FRACTION_LEN, MPFR_RNDN);
553
554 mpfr_sub(maxMPFR.value, maxMPFR.value, pivot.value, MPFR_RNDN);
555 mpfr_mul_2si(maxMPFR.value, maxMPFR.value,
556 -maxExponent + FPBits<T>::FRACTION_LEN, MPFR_RNDN);
557
558 mpfr_add(minMPFR.value, minMPFR.value, maxMPFR.value, MPFR_RNDN);
559 return minMPFR;
560 }
561
562 template <typename T>
563 cpp::enable_if_t<cpp::is_floating_point_v<T>, cpp::string>
564 ulp_as_string(T input) {
565 MPFRNumber num = ulp_as_mpfr_number(input);
566 return num.str();
567 }
568
569 template <typename T>
570 cpp::enable_if_t<cpp::is_floating_point_v<T>, double> ulp(T input) {
571 MPFRNumber num = ulp_as_mpfr_number(input);
572 return num.as<double>();
573 }
574};
575
576template <> float MPFRNumber::as<float>() const {
577 return mpfr_get_flt(value, mpfr_rounding);
578}
579
580template <> double MPFRNumber::as<double>() const {
581 return mpfr_get_d(value, mpfr_rounding);
582}
583
584template <> long double MPFRNumber::as<long double>() const {
585 return mpfr_get_ld(value, mpfr_rounding);
586}
587
588namespace internal {
589
590template <typename InputType>
591cpp::enable_if_t<cpp::is_floating_point_v<InputType>, MPFRNumber>
592unary_operation(Operation op, InputType input, unsigned int precision,
593 RoundingMode rounding) {
594 MPFRNumber mpfrInput(input, precision, rounding);
595 switch (op) {
596 case Operation::Abs:
597 return mpfrInput.abs();
598 case Operation::Acos:
599 return mpfrInput.acos();
600 case Operation::Acosh:
601 return mpfrInput.acosh();
602 case Operation::Asin:
603 return mpfrInput.asin();
604 case Operation::Asinh:
605 return mpfrInput.asinh();
606 case Operation::Atan:
607 return mpfrInput.atan();
608 case Operation::Atanh:
609 return mpfrInput.atanh();
610 case Operation::Ceil:
611 return mpfrInput.ceil();
612 case Operation::Cos:
613 return mpfrInput.cos();
614 case Operation::Cosh:
615 return mpfrInput.cosh();
616 case Operation::Erf:
617 return mpfrInput.erf();
618 case Operation::Exp:
619 return mpfrInput.exp();
620 case Operation::Exp2:
621 return mpfrInput.exp2();
622 case Operation::Exp2m1:
623 return mpfrInput.exp2m1();
624 case Operation::Exp10:
625 return mpfrInput.exp10();
626 case Operation::Expm1:
627 return mpfrInput.expm1();
628 case Operation::Floor:
629 return mpfrInput.floor();
630 case Operation::Log:
631 return mpfrInput.log();
632 case Operation::Log2:
633 return mpfrInput.log2();
634 case Operation::Log10:
635 return mpfrInput.log10();
636 case Operation::Log1p:
637 return mpfrInput.log1p();
638 case Operation::Mod2PI:
639 return mpfrInput.mod_2pi();
640 case Operation::ModPIOver2:
641 return mpfrInput.mod_pi_over_2();
642 case Operation::ModPIOver4:
643 return mpfrInput.mod_pi_over_4();
644 case Operation::Round:
645 return mpfrInput.round();
646 case Operation::RoundEven:
647 return mpfrInput.roundeven();
648 case Operation::Sin:
649 return mpfrInput.sin();
650 case Operation::Sinh:
651 return mpfrInput.sinh();
652 case Operation::Sqrt:
653 return mpfrInput.sqrt();
654 case Operation::Tan:
655 return mpfrInput.tan();
656 case Operation::Tanh:
657 return mpfrInput.tanh();
658 case Operation::Trunc:
659 return mpfrInput.trunc();
660 default:
661 __builtin_unreachable();
662 }
663}
664
665template <typename InputType>
666cpp::enable_if_t<cpp::is_floating_point_v<InputType>, MPFRNumber>
667unary_operation_two_outputs(Operation op, InputType input, int &output,
668 unsigned int precision, RoundingMode rounding) {
669 MPFRNumber mpfrInput(input, precision, rounding);
670 switch (op) {
671 case Operation::Frexp:
672 return mpfrInput.frexp(exp&: output);
673 default:
674 __builtin_unreachable();
675 }
676}
677
678template <typename InputType>
679cpp::enable_if_t<cpp::is_floating_point_v<InputType>, MPFRNumber>
680binary_operation_one_output(Operation op, InputType x, InputType y,
681 unsigned int precision, RoundingMode rounding) {
682 MPFRNumber inputX(x, precision, rounding);
683 MPFRNumber inputY(y, precision, rounding);
684 switch (op) {
685 case Operation::Atan2:
686 return inputX.atan2(b: inputY);
687 case Operation::Fmod:
688 return inputX.fmod(b: inputY);
689 case Operation::Hypot:
690 return inputX.hypot(b: inputY);
691 case Operation::Pow:
692 return inputX.pow(b: inputY);
693 default:
694 __builtin_unreachable();
695 }
696}
697
698template <typename InputType>
699cpp::enable_if_t<cpp::is_floating_point_v<InputType>, MPFRNumber>
700binary_operation_two_outputs(Operation op, InputType x, InputType y,
701 int &output, unsigned int precision,
702 RoundingMode rounding) {
703 MPFRNumber inputX(x, precision, rounding);
704 MPFRNumber inputY(y, precision, rounding);
705 switch (op) {
706 case Operation::RemQuo:
707 return inputX.remquo(divisor: inputY, quotient&: output);
708 default:
709 __builtin_unreachable();
710 }
711}
712
713template <typename InputType>
714cpp::enable_if_t<cpp::is_floating_point_v<InputType>, MPFRNumber>
715ternary_operation_one_output(Operation op, InputType x, InputType y,
716 InputType z, unsigned int precision,
717 RoundingMode rounding) {
718 // For FMA function, we just need to compare with the mpfr_fma with the same
719 // precision as InputType. Using higher precision as the intermediate results
720 // to compare might incorrectly fail due to double-rounding errors.
721 MPFRNumber inputX(x, precision, rounding);
722 MPFRNumber inputY(y, precision, rounding);
723 MPFRNumber inputZ(z, precision, rounding);
724 switch (op) {
725 case Operation::Fma:
726 return inputX.fma(b: inputY, c: inputZ);
727 default:
728 __builtin_unreachable();
729 }
730}
731
732// Remark: For all the explain_*_error functions, we will use std::stringstream
733// to build the complete error messages before sending it to the outstream `OS`
734// once at the end. This will stop the error messages from interleaving when
735// the tests are running concurrently.
736template <typename T>
737void explain_unary_operation_single_output_error(Operation op, T input,
738 T matchValue,
739 double ulp_tolerance,
740 RoundingMode rounding) {
741 unsigned int precision = get_precision<T>(ulp_tolerance);
742 MPFRNumber mpfrInput(input, precision);
743 MPFRNumber mpfr_result;
744 mpfr_result = unary_operation(op, input, precision, rounding);
745 MPFRNumber mpfrMatchValue(matchValue);
746 tlog << "Match value not within tolerance value of MPFR result:\n"
747 << " Input decimal: " << mpfrInput.str() << '\n';
748 tlog << " Input bits: " << str(FPBits<T>(input)) << '\n';
749 tlog << '\n' << " Match decimal: " << mpfrMatchValue.str() << '\n';
750 tlog << " Match bits: " << str(FPBits<T>(matchValue)) << '\n';
751 tlog << '\n' << " MPFR result: " << mpfr_result.str() << '\n';
752 tlog << " MPFR rounded: " << str(FPBits<T>(mpfr_result.as<T>())) << '\n';
753 tlog << '\n';
754 tlog << " ULP error: "
755 << mpfr_result.ulp_as_mpfr_number(matchValue).str() << '\n';
756}
757
758template void explain_unary_operation_single_output_error<float>(Operation op,
759 float, float,
760 double,
761 RoundingMode);
762template void explain_unary_operation_single_output_error<double>(
763 Operation op, double, double, double, RoundingMode);
764template void explain_unary_operation_single_output_error<long double>(
765 Operation op, long double, long double, double, RoundingMode);
766
767template <typename T>
768void explain_unary_operation_two_outputs_error(
769 Operation op, T input, const BinaryOutput<T> &libc_result,
770 double ulp_tolerance, RoundingMode rounding) {
771 unsigned int precision = get_precision<T>(ulp_tolerance);
772 MPFRNumber mpfrInput(input, precision);
773 int mpfrIntResult;
774 MPFRNumber mpfr_result = unary_operation_two_outputs(op, input, mpfrIntResult,
775 precision, rounding);
776
777 if (mpfrIntResult != libc_result.i) {
778 tlog << "MPFR integral result: " << mpfrIntResult << '\n'
779 << "Libc integral result: " << libc_result.i << '\n';
780 } else {
781 tlog << "Integral result from libc matches integral result from MPFR.\n";
782 }
783
784 MPFRNumber mpfrMatchValue(libc_result.f);
785 tlog
786 << "Libc floating point result is not within tolerance value of the MPFR "
787 << "result.\n\n";
788
789 tlog << " Input decimal: " << mpfrInput.str() << "\n\n";
790
791 tlog << "Libc floating point value: " << mpfrMatchValue.str() << '\n';
792 tlog << " Libc floating point bits: " << str(FPBits<T>(libc_result.f))
793 << '\n';
794 tlog << "\n\n";
795
796 tlog << " MPFR result: " << mpfr_result.str() << '\n';
797 tlog << " MPFR rounded: " << str(FPBits<T>(mpfr_result.as<T>()))
798 << '\n';
799 tlog << '\n'
800 << " ULP error: "
801 << mpfr_result.ulp_as_mpfr_number(libc_result.f).str() << '\n';
802}
803
804template void explain_unary_operation_two_outputs_error<float>(
805 Operation, float, const BinaryOutput<float> &, double, RoundingMode);
806template void explain_unary_operation_two_outputs_error<double>(
807 Operation, double, const BinaryOutput<double> &, double, RoundingMode);
808template void explain_unary_operation_two_outputs_error<long double>(
809 Operation, long double, const BinaryOutput<long double> &, double,
810 RoundingMode);
811
812template <typename T>
813void explain_binary_operation_two_outputs_error(
814 Operation op, const BinaryInput<T> &input,
815 const BinaryOutput<T> &libc_result, double ulp_tolerance,
816 RoundingMode rounding) {
817 unsigned int precision = get_precision<T>(ulp_tolerance);
818 MPFRNumber mpfrX(input.x, precision);
819 MPFRNumber mpfrY(input.y, precision);
820 int mpfrIntResult;
821 MPFRNumber mpfr_result = binary_operation_two_outputs(
822 op, input.x, input.y, mpfrIntResult, precision, rounding);
823 MPFRNumber mpfrMatchValue(libc_result.f);
824
825 tlog << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n'
826 << "MPFR integral result: " << mpfrIntResult << '\n'
827 << "Libc integral result: " << libc_result.i << '\n'
828 << "Libc floating point result: " << mpfrMatchValue.str() << '\n'
829 << " MPFR result: " << mpfr_result.str() << '\n';
830 tlog << "Libc floating point result bits: " << str(FPBits<T>(libc_result.f))
831 << '\n';
832 tlog << " MPFR rounded bits: "
833 << str(FPBits<T>(mpfr_result.as<T>())) << '\n';
834 tlog << "ULP error: " << mpfr_result.ulp_as_mpfr_number(libc_result.f).str()
835 << '\n';
836}
837
838template void explain_binary_operation_two_outputs_error<float>(
839 Operation, const BinaryInput<float> &, const BinaryOutput<float> &, double,
840 RoundingMode);
841template void explain_binary_operation_two_outputs_error<double>(
842 Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
843 double, RoundingMode);
844template void explain_binary_operation_two_outputs_error<long double>(
845 Operation, const BinaryInput<long double> &,
846 const BinaryOutput<long double> &, double, RoundingMode);
847
848template <typename T>
849void explain_binary_operation_one_output_error(Operation op,
850 const BinaryInput<T> &input,
851 T libc_result,
852 double ulp_tolerance,
853 RoundingMode rounding) {
854 unsigned int precision = get_precision<T>(ulp_tolerance);
855 MPFRNumber mpfrX(input.x, precision);
856 MPFRNumber mpfrY(input.y, precision);
857 FPBits<T> xbits(input.x);
858 FPBits<T> ybits(input.y);
859 MPFRNumber mpfr_result =
860 binary_operation_one_output(op, input.x, input.y, precision, rounding);
861 MPFRNumber mpfrMatchValue(libc_result);
862
863 tlog << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n';
864 tlog << "First input bits: " << str(FPBits<T>(input.x)) << '\n';
865 tlog << "Second input bits: " << str(FPBits<T>(input.y)) << '\n';
866
867 tlog << "Libc result: " << mpfrMatchValue.str() << '\n'
868 << "MPFR result: " << mpfr_result.str() << '\n';
869 tlog << "Libc floating point result bits: " << str(FPBits<T>(libc_result))
870 << '\n';
871 tlog << " MPFR rounded bits: "
872 << str(FPBits<T>(mpfr_result.as<T>())) << '\n';
873 tlog << "ULP error: " << mpfr_result.ulp_as_mpfr_number(libc_result).str()
874 << '\n';
875}
876
877template void explain_binary_operation_one_output_error<float>(
878 Operation, const BinaryInput<float> &, float, double, RoundingMode);
879template void explain_binary_operation_one_output_error<double>(
880 Operation, const BinaryInput<double> &, double, double, RoundingMode);
881template void explain_binary_operation_one_output_error<long double>(
882 Operation, const BinaryInput<long double> &, long double, double,
883 RoundingMode);
884
885template <typename T>
886void explain_ternary_operation_one_output_error(Operation op,
887 const TernaryInput<T> &input,
888 T libc_result,
889 double ulp_tolerance,
890 RoundingMode rounding) {
891 unsigned int precision = get_precision<T>(ulp_tolerance);
892 MPFRNumber mpfrX(input.x, precision);
893 MPFRNumber mpfrY(input.y, precision);
894 MPFRNumber mpfrZ(input.z, precision);
895 FPBits<T> xbits(input.x);
896 FPBits<T> ybits(input.y);
897 FPBits<T> zbits(input.z);
898 MPFRNumber mpfr_result = ternary_operation_one_output(
899 op, input.x, input.y, input.z, precision, rounding);
900 MPFRNumber mpfrMatchValue(libc_result);
901
902 tlog << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str()
903 << " z: " << mpfrZ.str() << '\n';
904 tlog << " First input bits: " << str(FPBits<T>(input.x)) << '\n';
905 tlog << "Second input bits: " << str(FPBits<T>(input.y)) << '\n';
906 tlog << " Third input bits: " << str(FPBits<T>(input.z)) << '\n';
907
908 tlog << "Libc result: " << mpfrMatchValue.str() << '\n'
909 << "MPFR result: " << mpfr_result.str() << '\n';
910 tlog << "Libc floating point result bits: " << str(FPBits<T>(libc_result))
911 << '\n';
912 tlog << " MPFR rounded bits: "
913 << str(FPBits<T>(mpfr_result.as<T>())) << '\n';
914 tlog << "ULP error: " << mpfr_result.ulp_as_mpfr_number(libc_result).str()
915 << '\n';
916}
917
918template void explain_ternary_operation_one_output_error<float>(
919 Operation, const TernaryInput<float> &, float, double, RoundingMode);
920template void explain_ternary_operation_one_output_error<double>(
921 Operation, const TernaryInput<double> &, double, double, RoundingMode);
922template void explain_ternary_operation_one_output_error<long double>(
923 Operation, const TernaryInput<long double> &, long double, double,
924 RoundingMode);
925
926template <typename T>
927bool compare_unary_operation_single_output(Operation op, T input, T libc_result,
928 double ulp_tolerance,
929 RoundingMode rounding) {
930 unsigned int precision = get_precision<T>(ulp_tolerance);
931 MPFRNumber mpfr_result;
932 mpfr_result = unary_operation(op, input, precision, rounding);
933 double ulp = mpfr_result.ulp(libc_result);
934 return (ulp <= ulp_tolerance);
935}
936
937template bool compare_unary_operation_single_output<float>(Operation, float,
938 float, double,
939 RoundingMode);
940template bool compare_unary_operation_single_output<double>(Operation, double,
941 double, double,
942 RoundingMode);
943template bool compare_unary_operation_single_output<long double>(
944 Operation, long double, long double, double, RoundingMode);
945
946template <typename T>
947bool compare_unary_operation_two_outputs(Operation op, T input,
948 const BinaryOutput<T> &libc_result,
949 double ulp_tolerance,
950 RoundingMode rounding) {
951 int mpfrIntResult;
952 unsigned int precision = get_precision<T>(ulp_tolerance);
953 MPFRNumber mpfr_result = unary_operation_two_outputs(op, input, mpfrIntResult,
954 precision, rounding);
955 double ulp = mpfr_result.ulp(libc_result.f);
956
957 if (mpfrIntResult != libc_result.i)
958 return false;
959
960 return (ulp <= ulp_tolerance);
961}
962
963template bool compare_unary_operation_two_outputs<float>(
964 Operation, float, const BinaryOutput<float> &, double, RoundingMode);
965template bool compare_unary_operation_two_outputs<double>(
966 Operation, double, const BinaryOutput<double> &, double, RoundingMode);
967template bool compare_unary_operation_two_outputs<long double>(
968 Operation, long double, const BinaryOutput<long double> &, double,
969 RoundingMode);
970
971template <typename T>
972bool compare_binary_operation_two_outputs(Operation op,
973 const BinaryInput<T> &input,
974 const BinaryOutput<T> &libc_result,
975 double ulp_tolerance,
976 RoundingMode rounding) {
977 int mpfrIntResult;
978 unsigned int precision = get_precision<T>(ulp_tolerance);
979 MPFRNumber mpfr_result = binary_operation_two_outputs(
980 op, input.x, input.y, mpfrIntResult, precision, rounding);
981 double ulp = mpfr_result.ulp(libc_result.f);
982
983 if (mpfrIntResult != libc_result.i) {
984 if (op == Operation::RemQuo) {
985 if ((0x7 & mpfrIntResult) != (0x7 & libc_result.i))
986 return false;
987 } else {
988 return false;
989 }
990 }
991
992 return (ulp <= ulp_tolerance);
993}
994
995template bool compare_binary_operation_two_outputs<float>(
996 Operation, const BinaryInput<float> &, const BinaryOutput<float> &, double,
997 RoundingMode);
998template bool compare_binary_operation_two_outputs<double>(
999 Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
1000 double, RoundingMode);
1001template bool compare_binary_operation_two_outputs<long double>(
1002 Operation, const BinaryInput<long double> &,
1003 const BinaryOutput<long double> &, double, RoundingMode);
1004
1005template <typename T>
1006bool compare_binary_operation_one_output(Operation op,
1007 const BinaryInput<T> &input,
1008 T libc_result, double ulp_tolerance,
1009 RoundingMode rounding) {
1010 unsigned int precision = get_precision<T>(ulp_tolerance);
1011 MPFRNumber mpfr_result =
1012 binary_operation_one_output(op, input.x, input.y, precision, rounding);
1013 double ulp = mpfr_result.ulp(libc_result);
1014
1015 return (ulp <= ulp_tolerance);
1016}
1017
1018template bool compare_binary_operation_one_output<float>(
1019 Operation, const BinaryInput<float> &, float, double, RoundingMode);
1020template bool compare_binary_operation_one_output<double>(
1021 Operation, const BinaryInput<double> &, double, double, RoundingMode);
1022template bool compare_binary_operation_one_output<long double>(
1023 Operation, const BinaryInput<long double> &, long double, double,
1024 RoundingMode);
1025
1026template <typename T>
1027bool compare_ternary_operation_one_output(Operation op,
1028 const TernaryInput<T> &input,
1029 T libc_result, double ulp_tolerance,
1030 RoundingMode rounding) {
1031 unsigned int precision = get_precision<T>(ulp_tolerance);
1032 MPFRNumber mpfr_result = ternary_operation_one_output(
1033 op, input.x, input.y, input.z, precision, rounding);
1034 double ulp = mpfr_result.ulp(libc_result);
1035
1036 return (ulp <= ulp_tolerance);
1037}
1038
1039template bool compare_ternary_operation_one_output<float>(
1040 Operation, const TernaryInput<float> &, float, double, RoundingMode);
1041template bool compare_ternary_operation_one_output<double>(
1042 Operation, const TernaryInput<double> &, double, double, RoundingMode);
1043template bool compare_ternary_operation_one_output<long double>(
1044 Operation, const TernaryInput<long double> &, long double, double,
1045 RoundingMode);
1046
1047} // namespace internal
1048
1049template <typename T> bool round_to_long(T x, long &result) {
1050 MPFRNumber mpfr(x);
1051 return mpfr.round_to_long(result);
1052}
1053
1054template bool round_to_long<float>(float, long &);
1055template bool round_to_long<double>(double, long &);
1056template bool round_to_long<long double>(long double, long &);
1057
1058template <typename T> bool round_to_long(T x, RoundingMode mode, long &result) {
1059 MPFRNumber mpfr(x);
1060 return mpfr.round_to_long(rnd: get_mpfr_rounding_mode(mode), result);
1061}
1062
1063template bool round_to_long<float>(float, RoundingMode, long &);
1064template bool round_to_long<double>(double, RoundingMode, long &);
1065template bool round_to_long<long double>(long double, RoundingMode, long &);
1066
1067template <typename T> T round(T x, RoundingMode mode) {
1068 MPFRNumber mpfr(x);
1069 MPFRNumber result = mpfr.rint(rnd: get_mpfr_rounding_mode(mode));
1070 return result.as<T>();
1071}
1072
1073template float round<float>(float, RoundingMode);
1074template double round<double>(double, RoundingMode);
1075template long double round<long double>(long double, RoundingMode);
1076
1077} // namespace mpfr
1078} // namespace testing
1079} // namespace LIBC_NAMESPACE
1080

source code of libc/utils/MPFRWrapper/MPFRUtils.cpp