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/StringView.h" 12 #include "src/__support/FPUtil/FPBits.h" 13 #include "src/__support/architectures.h" 14 #include "utils/UnitTest/FPMatcher.h" 15 16 #include <cmath> 17 #include <memory> 18 #include <stdint.h> 19 #include <string> 20 21 #ifdef CUSTOM_MPFR_INCLUDER 22 // Some downstream repos are monoliths carrying MPFR sources in their third 23 // party directory. In such repos, including the MPFR header as 24 // `#include <mpfr.h>` is either disallowed or not possible. If that is the 25 // case, a file named `CustomMPFRIncluder.h` should be added through which the 26 // MPFR header can be included in manner allowed in that repo. 27 #include "CustomMPFRIncluder.h" 28 #else 29 #include <mpfr.h> 30 #endif 31 32 template <typename T> using FPBits = __llvm_libc::fputil::FPBits<T>; 33 34 namespace __llvm_libc { 35 namespace testing { 36 namespace mpfr { 37 38 template <typename T> struct Precision; 39 40 template <> struct Precision<float> { 41 static constexpr unsigned int VALUE = 24; 42 }; 43 44 template <> struct Precision<double> { 45 static constexpr unsigned int VALUE = 53; 46 }; 47 48 #if !(defined(LLVM_LIBC_ARCH_X86)) 49 template <> struct Precision<long double> { 50 static constexpr unsigned int VALUE = 64; 51 }; 52 #else 53 template <> struct Precision<long double> { 54 static constexpr unsigned int VALUE = 113; 55 }; 56 #endif 57 58 class MPFRNumber { 59 // A precision value which allows sufficiently large additional 60 // precision even compared to quad-precision floating point values. 61 unsigned int mpfr_precision; 62 63 mpfr_t value; 64 65 public: 66 MPFRNumber() : mpfr_precision(256) { mpfr_init2(value, mpfr_precision); } 67 68 // We use explicit EnableIf specializations to disallow implicit 69 // conversions. Implicit conversions can potentially lead to loss of 70 // precision. 71 template <typename XType, 72 cpp::EnableIfType<cpp::IsSame<float, XType>::Value, int> = 0> 73 explicit MPFRNumber(XType x, int precision = 128) 74 : mpfr_precision(precision) { 75 mpfr_init2(value, mpfr_precision); 76 mpfr_set_flt(value, x, MPFR_RNDN); 77 } 78 79 template <typename XType, 80 cpp::EnableIfType<cpp::IsSame<double, XType>::Value, int> = 0> 81 explicit MPFRNumber(XType x, int precision = 128) 82 : mpfr_precision(precision) { 83 mpfr_init2(value, mpfr_precision); 84 mpfr_set_d(value, x, MPFR_RNDN); 85 } 86 87 template <typename XType, 88 cpp::EnableIfType<cpp::IsSame<long double, XType>::Value, int> = 0> 89 explicit MPFRNumber(XType x, int precision = 128) 90 : mpfr_precision(precision) { 91 mpfr_init2(value, mpfr_precision); 92 mpfr_set_ld(value, x, MPFR_RNDN); 93 } 94 95 template <typename XType, 96 cpp::EnableIfType<cpp::IsIntegral<XType>::Value, int> = 0> 97 explicit MPFRNumber(XType x, int precision = 128) 98 : mpfr_precision(precision) { 99 mpfr_init2(value, mpfr_precision); 100 mpfr_set_sj(value, x, MPFR_RNDN); 101 } 102 103 MPFRNumber(const MPFRNumber &other) : mpfr_precision(other.mpfr_precision) { 104 mpfr_init2(value, mpfr_precision); 105 mpfr_set(value, other.value, MPFR_RNDN); 106 } 107 108 ~MPFRNumber() { mpfr_clear(value); } 109 110 MPFRNumber &operator=(const MPFRNumber &rhs) { 111 mpfr_precision = rhs.mpfr_precision; 112 mpfr_set(value, rhs.value, MPFR_RNDN); 113 return *this; 114 } 115 116 MPFRNumber abs() const { 117 MPFRNumber result; 118 mpfr_abs(result.value, value, MPFR_RNDN); 119 return result; 120 } 121 122 MPFRNumber ceil() const { 123 MPFRNumber result; 124 mpfr_ceil(result.value, value); 125 return result; 126 } 127 128 MPFRNumber cos() const { 129 MPFRNumber result; 130 mpfr_cos(result.value, value, MPFR_RNDN); 131 return result; 132 } 133 134 MPFRNumber exp() const { 135 MPFRNumber result; 136 mpfr_exp(result.value, value, MPFR_RNDN); 137 return result; 138 } 139 140 MPFRNumber exp2() const { 141 MPFRNumber result; 142 mpfr_exp2(result.value, value, MPFR_RNDN); 143 return result; 144 } 145 146 MPFRNumber expm1() const { 147 MPFRNumber result; 148 mpfr_expm1(result.value, value, MPFR_RNDN); 149 return result; 150 } 151 152 MPFRNumber floor() const { 153 MPFRNumber result; 154 mpfr_floor(result.value, value); 155 return result; 156 } 157 158 MPFRNumber frexp(int &exp) { 159 MPFRNumber result; 160 mpfr_exp_t resultExp; 161 mpfr_frexp(&resultExp, result.value, value, MPFR_RNDN); 162 exp = resultExp; 163 return result; 164 } 165 166 MPFRNumber hypot(const MPFRNumber &b) { 167 MPFRNumber result; 168 mpfr_hypot(result.value, value, b.value, MPFR_RNDN); 169 return result; 170 } 171 172 MPFRNumber log() const { 173 MPFRNumber result; 174 mpfr_log(result.value, value, MPFR_RNDN); 175 return result; 176 } 177 178 MPFRNumber remquo(const MPFRNumber &divisor, int "ient) { 179 MPFRNumber remainder; 180 long q; 181 mpfr_remquo(remainder.value, &q, value, divisor.value, MPFR_RNDN); 182 quotient = q; 183 return remainder; 184 } 185 186 MPFRNumber round() const { 187 MPFRNumber result; 188 mpfr_round(result.value, value); 189 return result; 190 } 191 192 bool roung_to_long(long &result) const { 193 // We first calculate the rounded value. This way, when converting 194 // to long using mpfr_get_si, the rounding direction of MPFR_RNDN 195 // (or any other rounding mode), does not have an influence. 196 MPFRNumber roundedValue = round(); 197 mpfr_clear_erangeflag(); 198 result = mpfr_get_si(roundedValue.value, MPFR_RNDN); 199 return mpfr_erangeflag_p(); 200 } 201 202 bool roung_to_long(mpfr_rnd_t rnd, long &result) const { 203 MPFRNumber rint_result; 204 mpfr_rint(rint_result.value, value, rnd); 205 return rint_result.roung_to_long(result); 206 } 207 208 MPFRNumber rint(mpfr_rnd_t rnd) const { 209 MPFRNumber result; 210 mpfr_rint(result.value, value, rnd); 211 return result; 212 } 213 214 MPFRNumber mod_2pi() const { 215 MPFRNumber result(0.0, 1280); 216 MPFRNumber _2pi(0.0, 1280); 217 mpfr_const_pi(_2pi.value, MPFR_RNDN); 218 mpfr_mul_si(_2pi.value, _2pi.value, 2, MPFR_RNDN); 219 mpfr_fmod(result.value, value, _2pi.value, MPFR_RNDN); 220 return result; 221 } 222 223 MPFRNumber mod_pi_over_2() const { 224 MPFRNumber result(0.0, 1280); 225 MPFRNumber pi_over_2(0.0, 1280); 226 mpfr_const_pi(pi_over_2.value, MPFR_RNDN); 227 mpfr_mul_d(pi_over_2.value, pi_over_2.value, 0.5, MPFR_RNDN); 228 mpfr_fmod(result.value, value, pi_over_2.value, MPFR_RNDN); 229 return result; 230 } 231 232 MPFRNumber mod_pi_over_4() const { 233 MPFRNumber result(0.0, 1280); 234 MPFRNumber pi_over_4(0.0, 1280); 235 mpfr_const_pi(pi_over_4.value, MPFR_RNDN); 236 mpfr_mul_d(pi_over_4.value, pi_over_4.value, 0.25, MPFR_RNDN); 237 mpfr_fmod(result.value, value, pi_over_4.value, MPFR_RNDN); 238 return result; 239 } 240 241 MPFRNumber sin() const { 242 MPFRNumber result; 243 mpfr_sin(result.value, value, MPFR_RNDN); 244 return result; 245 } 246 247 MPFRNumber sqrt() const { 248 MPFRNumber result; 249 mpfr_sqrt(result.value, value, MPFR_RNDN); 250 return result; 251 } 252 253 MPFRNumber tan() const { 254 MPFRNumber result; 255 mpfr_tan(result.value, value, MPFR_RNDN); 256 return result; 257 } 258 259 MPFRNumber trunc() const { 260 MPFRNumber result; 261 mpfr_trunc(result.value, value); 262 return result; 263 } 264 265 MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c) { 266 MPFRNumber result(*this); 267 mpfr_fma(result.value, value, b.value, c.value, MPFR_RNDN); 268 return result; 269 } 270 271 std::string str() const { 272 // 200 bytes should be more than sufficient to hold a 100-digit number 273 // plus additional bytes for the decimal point, '-' sign etc. 274 constexpr size_t printBufSize = 200; 275 char buffer[printBufSize]; 276 mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value); 277 cpp::StringView view(buffer); 278 view = view.trim(' '); 279 return std::string(view.data()); 280 } 281 282 // These functions are useful for debugging. 283 template <typename T> T as() const; 284 285 template <> float as<float>() const { return mpfr_get_flt(value, MPFR_RNDN); } 286 template <> double as<double>() const { return mpfr_get_d(value, MPFR_RNDN); } 287 template <> long double as<long double>() const { 288 return mpfr_get_ld(value, MPFR_RNDN); 289 } 290 291 void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); } 292 293 // Return the ULP (units-in-the-last-place) difference between the 294 // stored MPFR and a floating point number. 295 // 296 // We define ULP difference as follows: 297 // If exponents of this value and the |input| are same, then: 298 // ULP(this_value, input) = abs(this_value - input) / eps(input) 299 // else: 300 // max = max(abs(this_value), abs(input)) 301 // min = min(abs(this_value), abs(input)) 302 // maxExponent = exponent(max) 303 // ULP(this_value, input) = (max - 2^maxExponent) / eps(max) + 304 // (2^maxExponent - min) / eps(min) 305 // 306 // Remarks: 307 // 1. A ULP of 0.0 will imply that the value is correctly rounded. 308 // 2. We expect that this value and the value to be compared (the [input] 309 // argument) are reasonable close, and we will provide an upper bound 310 // of ULP value for testing. Morever, most of the fractional parts of 311 // ULP value do not matter much, so using double as the return type 312 // should be good enough. 313 // 3. For close enough values (values which don't diff in their exponent by 314 // not more than 1), a ULP difference of N indicates a bit distance 315 // of N between this number and [input]. 316 // 4. A values of +0.0 and -0.0 are treated as equal. 317 template <typename T> 318 cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) { 319 T thisAsT = as<T>(); 320 if (thisAsT == input) 321 return T(0.0); 322 323 int thisExponent = fputil::FPBits<T>(thisAsT).get_exponent(); 324 int inputExponent = fputil::FPBits<T>(input).get_exponent(); 325 // Adjust the exponents for denormal numbers. 326 if (fputil::FPBits<T>(thisAsT).get_unbiased_exponent() == 0) 327 ++thisExponent; 328 if (fputil::FPBits<T>(input).get_unbiased_exponent() == 0) 329 ++inputExponent; 330 331 if (thisAsT * input < 0 || thisExponent == inputExponent) { 332 MPFRNumber inputMPFR(input); 333 mpfr_sub(inputMPFR.value, value, inputMPFR.value, MPFR_RNDN); 334 mpfr_abs(inputMPFR.value, inputMPFR.value, MPFR_RNDN); 335 mpfr_mul_2si(inputMPFR.value, inputMPFR.value, 336 -thisExponent + int(fputil::MantissaWidth<T>::VALUE), 337 MPFR_RNDN); 338 return inputMPFR.as<double>(); 339 } 340 341 // If the control reaches here, it means that this number and input are 342 // of the same sign but different exponent. In such a case, ULP error is 343 // calculated as sum of two parts. 344 thisAsT = std::abs(thisAsT); 345 input = std::abs(input); 346 T min = thisAsT > input ? input : thisAsT; 347 T max = thisAsT > input ? thisAsT : input; 348 int minExponent = fputil::FPBits<T>(min).get_exponent(); 349 int maxExponent = fputil::FPBits<T>(max).get_exponent(); 350 // Adjust the exponents for denormal numbers. 351 if (fputil::FPBits<T>(min).get_unbiased_exponent() == 0) 352 ++minExponent; 353 if (fputil::FPBits<T>(max).get_unbiased_exponent() == 0) 354 ++maxExponent; 355 356 MPFRNumber minMPFR(min); 357 MPFRNumber maxMPFR(max); 358 359 MPFRNumber pivot(uint32_t(1)); 360 mpfr_mul_2si(pivot.value, pivot.value, maxExponent, MPFR_RNDN); 361 362 mpfr_sub(minMPFR.value, pivot.value, minMPFR.value, MPFR_RNDN); 363 mpfr_mul_2si(minMPFR.value, minMPFR.value, 364 -minExponent + int(fputil::MantissaWidth<T>::VALUE), 365 MPFR_RNDN); 366 367 mpfr_sub(maxMPFR.value, maxMPFR.value, pivot.value, MPFR_RNDN); 368 mpfr_mul_2si(maxMPFR.value, maxMPFR.value, 369 -maxExponent + int(fputil::MantissaWidth<T>::VALUE), 370 MPFR_RNDN); 371 372 mpfr_add(minMPFR.value, minMPFR.value, maxMPFR.value, MPFR_RNDN); 373 return minMPFR.as<double>(); 374 } 375 }; 376 377 namespace internal { 378 379 template <typename InputType> 380 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber> 381 unary_operation(Operation op, InputType input) { 382 MPFRNumber mpfrInput(input); 383 switch (op) { 384 case Operation::Abs: 385 return mpfrInput.abs(); 386 case Operation::Ceil: 387 return mpfrInput.ceil(); 388 case Operation::Cos: 389 return mpfrInput.cos(); 390 case Operation::Exp: 391 return mpfrInput.exp(); 392 case Operation::Exp2: 393 return mpfrInput.exp2(); 394 case Operation::Expm1: 395 return mpfrInput.expm1(); 396 case Operation::Floor: 397 return mpfrInput.floor(); 398 case Operation::Log: 399 return mpfrInput.log(); 400 case Operation::Mod2PI: 401 return mpfrInput.mod_2pi(); 402 case Operation::ModPIOver2: 403 return mpfrInput.mod_pi_over_2(); 404 case Operation::ModPIOver4: 405 return mpfrInput.mod_pi_over_4(); 406 case Operation::Round: 407 return mpfrInput.round(); 408 case Operation::Sin: 409 return mpfrInput.sin(); 410 case Operation::Sqrt: 411 return mpfrInput.sqrt(); 412 case Operation::Tan: 413 return mpfrInput.tan(); 414 case Operation::Trunc: 415 return mpfrInput.trunc(); 416 default: 417 __builtin_unreachable(); 418 } 419 } 420 421 template <typename InputType> 422 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber> 423 unary_operation_two_outputs(Operation op, InputType input, int &output) { 424 MPFRNumber mpfrInput(input); 425 switch (op) { 426 case Operation::Frexp: 427 return mpfrInput.frexp(output); 428 default: 429 __builtin_unreachable(); 430 } 431 } 432 433 template <typename InputType> 434 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber> 435 binary_operation_one_output(Operation op, InputType x, InputType y) { 436 MPFRNumber inputX(x), inputY(y); 437 switch (op) { 438 case Operation::Hypot: 439 return inputX.hypot(inputY); 440 default: 441 __builtin_unreachable(); 442 } 443 } 444 445 template <typename InputType> 446 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber> 447 binary_operation_two_outputs(Operation op, InputType x, InputType y, 448 int &output) { 449 MPFRNumber inputX(x), inputY(y); 450 switch (op) { 451 case Operation::RemQuo: 452 return inputX.remquo(inputY, output); 453 default: 454 __builtin_unreachable(); 455 } 456 } 457 458 template <typename InputType> 459 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber> 460 ternary_operation_one_output(Operation op, InputType x, InputType y, 461 InputType z) { 462 // For FMA function, we just need to compare with the mpfr_fma with the same 463 // precision as InputType. Using higher precision as the intermediate results 464 // to compare might incorrectly fail due to double-rounding errors. 465 constexpr unsigned int prec = Precision<InputType>::VALUE; 466 MPFRNumber inputX(x, prec), inputY(y, prec), inputZ(z, prec); 467 switch (op) { 468 case Operation::Fma: 469 return inputX.fma(inputY, inputZ); 470 default: 471 __builtin_unreachable(); 472 } 473 } 474 475 template <typename T> 476 void explain_unary_operation_single_output_error(Operation op, T input, 477 T matchValue, 478 testutils::StreamWrapper &OS) { 479 MPFRNumber mpfrInput(input); 480 MPFRNumber mpfr_result = unary_operation(op, input); 481 MPFRNumber mpfrMatchValue(matchValue); 482 FPBits<T> inputBits(input); 483 FPBits<T> matchBits(matchValue); 484 FPBits<T> mpfr_resultBits(mpfr_result.as<T>()); 485 OS << "Match value not within tolerance value of MPFR result:\n" 486 << " Input decimal: " << mpfrInput.str() << '\n'; 487 __llvm_libc::fputil::testing::describeValue(" Input bits: ", input, OS); 488 OS << '\n' << " Match decimal: " << mpfrMatchValue.str() << '\n'; 489 __llvm_libc::fputil::testing::describeValue(" Match bits: ", matchValue, 490 OS); 491 OS << '\n' << " MPFR result: " << mpfr_result.str() << '\n'; 492 __llvm_libc::fputil::testing::describeValue( 493 " MPFR rounded: ", mpfr_result.as<T>(), OS); 494 OS << '\n'; 495 OS << " ULP error: " << std::to_string(mpfr_result.ulp(matchValue)) 496 << '\n'; 497 } 498 499 template void 500 explain_unary_operation_single_output_error<float>(Operation op, float, float, 501 testutils::StreamWrapper &); 502 template void explain_unary_operation_single_output_error<double>( 503 Operation op, double, double, testutils::StreamWrapper &); 504 template void explain_unary_operation_single_output_error<long double>( 505 Operation op, long double, long double, testutils::StreamWrapper &); 506 507 template <typename T> 508 void explain_unary_operation_two_outputs_error( 509 Operation op, T input, const BinaryOutput<T> &libc_result, 510 testutils::StreamWrapper &OS) { 511 MPFRNumber mpfrInput(input); 512 FPBits<T> inputBits(input); 513 int mpfrIntResult; 514 MPFRNumber mpfr_result = 515 unary_operation_two_outputs(op, input, mpfrIntResult); 516 517 if (mpfrIntResult != libc_result.i) { 518 OS << "MPFR integral result: " << mpfrIntResult << '\n' 519 << "Libc integral result: " << libc_result.i << '\n'; 520 } else { 521 OS << "Integral result from libc matches integral result from MPFR.\n"; 522 } 523 524 MPFRNumber mpfrMatchValue(libc_result.f); 525 OS << "Libc floating point result is not within tolerance value of the MPFR " 526 << "result.\n\n"; 527 528 OS << " Input decimal: " << mpfrInput.str() << "\n\n"; 529 530 OS << "Libc floating point value: " << mpfrMatchValue.str() << '\n'; 531 __llvm_libc::fputil::testing::describeValue( 532 " Libc floating point bits: ", libc_result.f, OS); 533 OS << "\n\n"; 534 535 OS << " MPFR result: " << mpfr_result.str() << '\n'; 536 __llvm_libc::fputil::testing::describeValue( 537 " MPFR rounded: ", mpfr_result.as<T>(), OS); 538 OS << '\n' 539 << " ULP error: " 540 << std::to_string(mpfr_result.ulp(libc_result.f)) << '\n'; 541 } 542 543 template void explain_unary_operation_two_outputs_error<float>( 544 Operation, float, const BinaryOutput<float> &, testutils::StreamWrapper &); 545 template void 546 explain_unary_operation_two_outputs_error<double>(Operation, double, 547 const BinaryOutput<double> &, 548 testutils::StreamWrapper &); 549 template void explain_unary_operation_two_outputs_error<long double>( 550 Operation, long double, const BinaryOutput<long double> &, 551 testutils::StreamWrapper &); 552 553 template <typename T> 554 void explain_binary_operation_two_outputs_error( 555 Operation op, const BinaryInput<T> &input, 556 const BinaryOutput<T> &libc_result, testutils::StreamWrapper &OS) { 557 MPFRNumber mpfrX(input.x); 558 MPFRNumber mpfrY(input.y); 559 FPBits<T> xbits(input.x); 560 FPBits<T> ybits(input.y); 561 int mpfrIntResult; 562 MPFRNumber mpfr_result = 563 binary_operation_two_outputs(op, input.x, input.y, mpfrIntResult); 564 MPFRNumber mpfrMatchValue(libc_result.f); 565 566 OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n' 567 << "MPFR integral result: " << mpfrIntResult << '\n' 568 << "Libc integral result: " << libc_result.i << '\n' 569 << "Libc floating point result: " << mpfrMatchValue.str() << '\n' 570 << " MPFR result: " << mpfr_result.str() << '\n'; 571 __llvm_libc::fputil::testing::describeValue( 572 "Libc floating point result bits: ", libc_result.f, OS); 573 __llvm_libc::fputil::testing::describeValue( 574 " MPFR rounded bits: ", mpfr_result.as<T>(), OS); 575 OS << "ULP error: " << std::to_string(mpfr_result.ulp(libc_result.f)) << '\n'; 576 } 577 578 template void explain_binary_operation_two_outputs_error<float>( 579 Operation, const BinaryInput<float> &, const BinaryOutput<float> &, 580 testutils::StreamWrapper &); 581 template void explain_binary_operation_two_outputs_error<double>( 582 Operation, const BinaryInput<double> &, const BinaryOutput<double> &, 583 testutils::StreamWrapper &); 584 template void explain_binary_operation_two_outputs_error<long double>( 585 Operation, const BinaryInput<long double> &, 586 const BinaryOutput<long double> &, testutils::StreamWrapper &); 587 588 template <typename T> 589 void explain_binary_operation_one_output_error(Operation op, 590 const BinaryInput<T> &input, 591 T libc_result, 592 testutils::StreamWrapper &OS) { 593 MPFRNumber mpfrX(input.x); 594 MPFRNumber mpfrY(input.y); 595 FPBits<T> xbits(input.x); 596 FPBits<T> ybits(input.y); 597 MPFRNumber mpfr_result = binary_operation_one_output(op, input.x, input.y); 598 MPFRNumber mpfrMatchValue(libc_result); 599 600 OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n'; 601 __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x, 602 OS); 603 __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y, 604 OS); 605 606 OS << "Libc result: " << mpfrMatchValue.str() << '\n' 607 << "MPFR result: " << mpfr_result.str() << '\n'; 608 __llvm_libc::fputil::testing::describeValue( 609 "Libc floating point result bits: ", libc_result, OS); 610 __llvm_libc::fputil::testing::describeValue( 611 " MPFR rounded bits: ", mpfr_result.as<T>(), OS); 612 OS << "ULP error: " << std::to_string(mpfr_result.ulp(libc_result)) << '\n'; 613 } 614 615 template void explain_binary_operation_one_output_error<float>( 616 Operation, const BinaryInput<float> &, float, testutils::StreamWrapper &); 617 template void explain_binary_operation_one_output_error<double>( 618 Operation, const BinaryInput<double> &, double, testutils::StreamWrapper &); 619 template void explain_binary_operation_one_output_error<long double>( 620 Operation, const BinaryInput<long double> &, long double, 621 testutils::StreamWrapper &); 622 623 template <typename T> 624 void explain_ternary_operation_one_output_error(Operation op, 625 const TernaryInput<T> &input, 626 T libc_result, 627 testutils::StreamWrapper &OS) { 628 MPFRNumber mpfrX(input.x, Precision<T>::VALUE); 629 MPFRNumber mpfrY(input.y, Precision<T>::VALUE); 630 MPFRNumber mpfrZ(input.z, Precision<T>::VALUE); 631 FPBits<T> xbits(input.x); 632 FPBits<T> ybits(input.y); 633 FPBits<T> zbits(input.z); 634 MPFRNumber mpfr_result = 635 ternary_operation_one_output(op, input.x, input.y, input.z); 636 MPFRNumber mpfrMatchValue(libc_result); 637 638 OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() 639 << " z: " << mpfrZ.str() << '\n'; 640 __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x, 641 OS); 642 __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y, 643 OS); 644 __llvm_libc::fputil::testing::describeValue("Third input bits: ", input.z, 645 OS); 646 647 OS << "Libc result: " << mpfrMatchValue.str() << '\n' 648 << "MPFR result: " << mpfr_result.str() << '\n'; 649 __llvm_libc::fputil::testing::describeValue( 650 "Libc floating point result bits: ", libc_result, OS); 651 __llvm_libc::fputil::testing::describeValue( 652 " MPFR rounded bits: ", mpfr_result.as<T>(), OS); 653 OS << "ULP error: " << std::to_string(mpfr_result.ulp(libc_result)) << '\n'; 654 } 655 656 template void explain_ternary_operation_one_output_error<float>( 657 Operation, const TernaryInput<float> &, float, testutils::StreamWrapper &); 658 template void explain_ternary_operation_one_output_error<double>( 659 Operation, const TernaryInput<double> &, double, 660 testutils::StreamWrapper &); 661 template void explain_ternary_operation_one_output_error<long double>( 662 Operation, const TernaryInput<long double> &, long double, 663 testutils::StreamWrapper &); 664 665 template <typename T> 666 bool compare_unary_operation_single_output(Operation op, T input, T libc_result, 667 double ulp_error) { 668 // If the ulp error is exactly 0.5 (i.e a tie), we would check that the result 669 // is rounded to the nearest even. 670 MPFRNumber mpfr_result = unary_operation(op, input); 671 double ulp = mpfr_result.ulp(libc_result); 672 bool bits_are_even = ((FPBits<T>(libc_result).uintval() & 1) == 0); 673 return (ulp < ulp_error) || 674 ((ulp == ulp_error) && ((ulp != 0.5) || bits_are_even)); 675 } 676 677 template bool compare_unary_operation_single_output<float>(Operation, float, 678 float, double); 679 template bool compare_unary_operation_single_output<double>(Operation, double, 680 double, double); 681 template bool compare_unary_operation_single_output<long double>(Operation, 682 long double, 683 long double, 684 double); 685 686 template <typename T> 687 bool compare_unary_operation_two_outputs(Operation op, T input, 688 const BinaryOutput<T> &libc_result, 689 double ulp_error) { 690 int mpfrIntResult; 691 MPFRNumber mpfr_result = 692 unary_operation_two_outputs(op, input, mpfrIntResult); 693 double ulp = mpfr_result.ulp(libc_result.f); 694 695 if (mpfrIntResult != libc_result.i) 696 return false; 697 698 bool bits_are_even = ((FPBits<T>(libc_result.f).uintval() & 1) == 0); 699 return (ulp < ulp_error) || 700 ((ulp == ulp_error) && ((ulp != 0.5) || bits_are_even)); 701 } 702 703 template bool 704 compare_unary_operation_two_outputs<float>(Operation, float, 705 const BinaryOutput<float> &, double); 706 template bool compare_unary_operation_two_outputs<double>( 707 Operation, double, const BinaryOutput<double> &, double); 708 template bool compare_unary_operation_two_outputs<long double>( 709 Operation, long double, const BinaryOutput<long double> &, double); 710 711 template <typename T> 712 bool compare_binary_operation_two_outputs(Operation op, 713 const BinaryInput<T> &input, 714 const BinaryOutput<T> &libc_result, 715 double ulp_error) { 716 int mpfrIntResult; 717 MPFRNumber mpfr_result = 718 binary_operation_two_outputs(op, input.x, input.y, mpfrIntResult); 719 double ulp = mpfr_result.ulp(libc_result.f); 720 721 if (mpfrIntResult != libc_result.i) { 722 if (op == Operation::RemQuo) { 723 if ((0x7 & mpfrIntResult) != (0x7 & libc_result.i)) 724 return false; 725 } else { 726 return false; 727 } 728 } 729 730 bool bits_are_even = ((FPBits<T>(libc_result.f).uintval() & 1) == 0); 731 return (ulp < ulp_error) || 732 ((ulp == ulp_error) && ((ulp != 0.5) || bits_are_even)); 733 } 734 735 template bool compare_binary_operation_two_outputs<float>( 736 Operation, const BinaryInput<float> &, const BinaryOutput<float> &, double); 737 template bool compare_binary_operation_two_outputs<double>( 738 Operation, const BinaryInput<double> &, const BinaryOutput<double> &, 739 double); 740 template bool compare_binary_operation_two_outputs<long double>( 741 Operation, const BinaryInput<long double> &, 742 const BinaryOutput<long double> &, double); 743 744 template <typename T> 745 bool compare_binary_operation_one_output(Operation op, 746 const BinaryInput<T> &input, 747 T libc_result, double ulp_error) { 748 MPFRNumber mpfr_result = binary_operation_one_output(op, input.x, input.y); 749 double ulp = mpfr_result.ulp(libc_result); 750 751 bool bits_are_even = ((FPBits<T>(libc_result).uintval() & 1) == 0); 752 return (ulp < ulp_error) || 753 ((ulp == ulp_error) && ((ulp != 0.5) || bits_are_even)); 754 } 755 756 template bool compare_binary_operation_one_output<float>( 757 Operation, const BinaryInput<float> &, float, double); 758 template bool compare_binary_operation_one_output<double>( 759 Operation, const BinaryInput<double> &, double, double); 760 template bool compare_binary_operation_one_output<long double>( 761 Operation, const BinaryInput<long double> &, long double, double); 762 763 template <typename T> 764 bool compare_ternary_operation_one_output(Operation op, 765 const TernaryInput<T> &input, 766 T libc_result, double ulp_error) { 767 MPFRNumber mpfr_result = 768 ternary_operation_one_output(op, input.x, input.y, input.z); 769 double ulp = mpfr_result.ulp(libc_result); 770 771 bool bits_are_even = ((FPBits<T>(libc_result).uintval() & 1) == 0); 772 return (ulp < ulp_error) || 773 ((ulp == ulp_error) && ((ulp != 0.5) || bits_are_even)); 774 } 775 776 template bool compare_ternary_operation_one_output<float>( 777 Operation, const TernaryInput<float> &, float, double); 778 template bool compare_ternary_operation_one_output<double>( 779 Operation, const TernaryInput<double> &, double, double); 780 template bool compare_ternary_operation_one_output<long double>( 781 Operation, const TernaryInput<long double> &, long double, double); 782 783 static mpfr_rnd_t get_mpfr_rounding_mode(RoundingMode mode) { 784 switch (mode) { 785 case RoundingMode::Upward: 786 return MPFR_RNDU; 787 break; 788 case RoundingMode::Downward: 789 return MPFR_RNDD; 790 break; 791 case RoundingMode::TowardZero: 792 return MPFR_RNDZ; 793 break; 794 case RoundingMode::Nearest: 795 return MPFR_RNDN; 796 break; 797 } 798 } 799 800 } // namespace internal 801 802 template <typename T> bool round_to_long(T x, long &result) { 803 MPFRNumber mpfr(x); 804 return mpfr.roung_to_long(result); 805 } 806 807 template bool round_to_long<float>(float, long &); 808 template bool round_to_long<double>(double, long &); 809 template bool round_to_long<long double>(long double, long &); 810 811 template <typename T> bool round_to_long(T x, RoundingMode mode, long &result) { 812 MPFRNumber mpfr(x); 813 return mpfr.roung_to_long(internal::get_mpfr_rounding_mode(mode), result); 814 } 815 816 template bool round_to_long<float>(float, RoundingMode, long &); 817 template bool round_to_long<double>(double, RoundingMode, long &); 818 template bool round_to_long<long double>(long double, RoundingMode, long &); 819 820 template <typename T> T round(T x, RoundingMode mode) { 821 MPFRNumber mpfr(x); 822 MPFRNumber result = mpfr.rint(internal::get_mpfr_rounding_mode(mode)); 823 return result.as<T>(); 824 } 825 826 template float round<float>(float, RoundingMode); 827 template double round<double>(double, RoundingMode); 828 template long double round<long double>(long double, RoundingMode); 829 830 } // namespace mpfr 831 } // namespace testing 832 } // namespace __llvm_libc 833