libc/utils/MPFRWrapper/MPFRUtils.cpp

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