611 lines · cpp
1//===-- Utils used by both MPCWrapper and MPFRWrapper----------------------===//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-exception6//7//===----------------------------------------------------------------------===//8 9#include "MPCommon.h"10 11#include "src/__support/CPP/string_view.h"12#include "src/__support/FPUtil/bfloat16.h"13#include "src/__support/FPUtil/cast.h"14#include "src/__support/macros/config.h"15#include "src/__support/macros/properties/types.h"16 17namespace LIBC_NAMESPACE_DECL {18namespace testing {19namespace mpfr {20 21MPFRNumber::MPFRNumber() : mpfr_precision(256), mpfr_rounding(MPFR_RNDN) {22 mpfr_init2(value, mpfr_precision);23}24 25MPFRNumber::MPFRNumber(const MPFRNumber &other)26 : mpfr_precision(other.mpfr_precision), mpfr_rounding(other.mpfr_rounding) {27 mpfr_init2(value, mpfr_precision);28 mpfr_set(value, other.value, mpfr_rounding);29}30 31MPFRNumber::MPFRNumber(const MPFRNumber &other, unsigned int precision)32 : mpfr_precision(precision), mpfr_rounding(other.mpfr_rounding) {33 mpfr_init2(value, mpfr_precision);34 mpfr_set(value, other.value, mpfr_rounding);35}36 37MPFRNumber::MPFRNumber(const mpfr_t x, unsigned int precision,38 RoundingMode rounding)39 : mpfr_precision(precision),40 mpfr_rounding(get_mpfr_rounding_mode(rounding)) {41 mpfr_init2(value, mpfr_precision);42 mpfr_set(value, x, mpfr_rounding);43}44 45MPFRNumber::~MPFRNumber() { mpfr_clear(value); }46 47MPFRNumber &MPFRNumber::operator=(const MPFRNumber &rhs) {48 mpfr_precision = rhs.mpfr_precision;49 mpfr_rounding = rhs.mpfr_rounding;50 mpfr_set(value, rhs.value, mpfr_rounding);51 return *this;52}53 54bool MPFRNumber::is_nan() const { return mpfr_nan_p(value); }55 56MPFRNumber MPFRNumber::abs() const {57 MPFRNumber result(*this);58 mpfr_abs(result.value, value, mpfr_rounding);59 return result;60}61 62MPFRNumber MPFRNumber::acos() const {63 MPFRNumber result(*this);64 mpfr_acos(result.value, value, mpfr_rounding);65 return result;66}67 68MPFRNumber MPFRNumber::acosh() const {69 MPFRNumber result(*this);70 mpfr_acosh(result.value, value, mpfr_rounding);71 return result;72}73 74MPFRNumber MPFRNumber::acospi() const {75 MPFRNumber result(*this);76 77#if MPFR_VERSION >= MPFR_VERSION_NUM(4, 2, 0)78 mpfr_acospi(result.value, value, mpfr_rounding);79 return result;80#else81 MPFRNumber value_acos(0.0, 1280);82 mpfr_acos(value_acos.value, value, MPFR_RNDN);83 MPFRNumber value_pi(0.0, 1280);84 mpfr_const_pi(value_pi.value, MPFR_RNDN);85 mpfr_div(result.value, value_acos.value, value_pi.value, mpfr_rounding);86 return result;87#endif88}89 90MPFRNumber MPFRNumber::add(const MPFRNumber &b) const {91 MPFRNumber result(*this);92 mpfr_add(result.value, value, b.value, mpfr_rounding);93 return result;94}95 96MPFRNumber MPFRNumber::asin() const {97 MPFRNumber result(*this);98 mpfr_asin(result.value, value, mpfr_rounding);99 return result;100}101 102MPFRNumber MPFRNumber::asinh() const {103 MPFRNumber result(*this);104 mpfr_asinh(result.value, value, mpfr_rounding);105 return result;106}107 108MPFRNumber MPFRNumber::asinpi() const {109 MPFRNumber result(*this);110#if MPFR_VERSION >= MPFR_VERSION_NUM(4, 2, 0)111 mpfr_asinpi(result.value, value, mpfr_rounding);112 return result;113#else114 MPFRNumber value_asin(0.0, 1280);115 mpfr_asin(value_asin.value, value, MPFR_RNDN);116 MPFRNumber value_pi(0.0, 1280);117 mpfr_const_pi(value_pi.value, MPFR_RNDN);118 mpfr_div(result.value, value_asin.value, value_pi.value, mpfr_rounding);119 return result;120#endif121}122 123MPFRNumber MPFRNumber::atan() const {124 MPFRNumber result(*this);125 mpfr_atan(result.value, value, mpfr_rounding);126 return result;127}128 129MPFRNumber MPFRNumber::atan2(const MPFRNumber &b) {130 MPFRNumber result(*this);131 mpfr_atan2(result.value, value, b.value, mpfr_rounding);132 return result;133}134 135MPFRNumber MPFRNumber::atanh() const {136 MPFRNumber result(*this);137 mpfr_atanh(result.value, value, mpfr_rounding);138 return result;139}140 141MPFRNumber MPFRNumber::atanpi() const {142 MPFRNumber result(*this);143#if MPFR_VERSION >= MPFR_VERSION_NUM(4, 2, 0)144 mpfr_atanpi(result.value, value, mpfr_rounding);145 return result;146#else147 MPFRNumber value_atan(0.0, mpfr_precision * 3);148 mpfr_atan(value_atan.value, value, MPFR_RNDN);149 MPFRNumber value_pi(0.0, mpfr_precision * 3);150 mpfr_const_pi(value_pi.value, MPFR_RNDN);151 mpfr_div(result.value, value_atan.value, value_pi.value, mpfr_rounding);152 return result;153#endif154}155 156MPFRNumber MPFRNumber::cbrt() const {157 MPFRNumber result(*this);158 mpfr_cbrt(result.value, value, mpfr_rounding);159 return result;160}161 162MPFRNumber MPFRNumber::ceil() const {163 MPFRNumber result(*this);164 mpfr_ceil(result.value, value);165 return result;166}167 168MPFRNumber MPFRNumber::cos() const {169 MPFRNumber result(*this);170 mpfr_cos(result.value, value, mpfr_rounding);171 return result;172}173 174MPFRNumber MPFRNumber::cosh() const {175 MPFRNumber result(*this);176 mpfr_cosh(result.value, value, mpfr_rounding);177 return result;178}179 180MPFRNumber MPFRNumber::cospi() const {181 MPFRNumber result(*this);182 183#if MPFR_VERSION >= MPFR_VERSION_NUM(4, 2, 0)184 mpfr_cospi(result.value, value, mpfr_rounding);185 return result;186#else187 if (mpfr_integer_p(value)) {188 mpz_t integer;189 mpz_init(integer);190 mpfr_get_z(integer, value, mpfr_rounding);191 192 int d = mpz_tstbit(integer, 0);193 mpfr_set_si(result.value, d ? -1 : 1, mpfr_rounding);194 mpz_clear(integer);195 return result;196 }197 198 MPFRNumber value_pi(0.0, 1280);199 mpfr_const_pi(value_pi.value, MPFR_RNDN);200 mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN);201 mpfr_cos(result.value, value_pi.value, mpfr_rounding);202 203 return result;204#endif205}206 207MPFRNumber MPFRNumber::erf() const {208 MPFRNumber result(*this);209 mpfr_erf(result.value, value, mpfr_rounding);210 return result;211}212 213MPFRNumber MPFRNumber::exp() const {214 MPFRNumber result(*this);215 mpfr_exp(result.value, value, mpfr_rounding);216 return result;217}218 219MPFRNumber MPFRNumber::exp2() const {220 MPFRNumber result(*this);221 mpfr_exp2(result.value, value, mpfr_rounding);222 return result;223}224 225MPFRNumber MPFRNumber::exp2m1() const {226 // TODO: Only use mpfr_exp2m1 once CI and buildbots get MPFR >= 4.2.0.227#if MPFR_VERSION >= MPFR_VERSION_NUM(4, 2, 0)228 MPFRNumber result(*this);229 mpfr_exp2m1(result.value, value, mpfr_rounding);230 return result;231#else232 unsigned int prec = mpfr_precision * 3;233 MPFRNumber result(*this, prec);234 235 float f = mpfr_get_flt(abs().value, mpfr_rounding);236 if (f > 0.5f && f < 0x1.0p30f) {237 mpfr_exp2(result.value, value, mpfr_rounding);238 mpfr_sub_ui(result.value, result.value, 1, mpfr_rounding);239 return result;240 }241 242 MPFRNumber ln2(2.0f, prec);243 // log(2)244 mpfr_log(ln2.value, ln2.value, mpfr_rounding);245 // x * log(2)246 mpfr_mul(result.value, value, ln2.value, mpfr_rounding);247 // e^(x * log(2)) - 1248 int ex = mpfr_expm1(result.value, result.value, mpfr_rounding);249 mpfr_subnormalize(result.value, ex, mpfr_rounding);250 return result;251#endif252}253 254MPFRNumber MPFRNumber::exp10() const {255 MPFRNumber result(*this);256 mpfr_exp10(result.value, value, mpfr_rounding);257 return result;258}259 260MPFRNumber MPFRNumber::exp10m1() const {261 // TODO: Only use mpfr_exp10m1 once CI and buildbots get MPFR >= 4.2.0.262#if MPFR_VERSION >= MPFR_VERSION_NUM(4, 2, 0)263 MPFRNumber result(*this);264 mpfr_exp10m1(result.value, value, mpfr_rounding);265 return result;266#else267 unsigned int prec = mpfr_precision * 3;268 MPFRNumber result(*this, prec);269 270 MPFRNumber ln10(10.0f, prec);271 // log(10)272 mpfr_log(ln10.value, ln10.value, mpfr_rounding);273 // x * log(10)274 mpfr_mul(result.value, value, ln10.value, mpfr_rounding);275 // e^(x * log(10)) - 1276 int ex = mpfr_expm1(result.value, result.value, mpfr_rounding);277 mpfr_subnormalize(result.value, ex, mpfr_rounding);278 return result;279#endif280}281 282MPFRNumber MPFRNumber::expm1() const {283 MPFRNumber result(*this);284 mpfr_expm1(result.value, value, mpfr_rounding);285 return result;286}287 288MPFRNumber MPFRNumber::div(const MPFRNumber &b) const {289 MPFRNumber result(*this);290 mpfr_div(result.value, value, b.value, mpfr_rounding);291 return result;292}293 294MPFRNumber MPFRNumber::floor() const {295 MPFRNumber result(*this);296 mpfr_floor(result.value, value);297 return result;298}299 300MPFRNumber MPFRNumber::fmod(const MPFRNumber &b) {301 MPFRNumber result(*this);302 mpfr_fmod(result.value, value, b.value, mpfr_rounding);303 return result;304}305 306MPFRNumber MPFRNumber::frexp(int &exp) {307 MPFRNumber result(*this);308 mpfr_exp_t resultExp;309 mpfr_frexp(&resultExp, result.value, value, mpfr_rounding);310 exp = static_cast<int>(resultExp);311 return result;312}313 314MPFRNumber MPFRNumber::hypot(const MPFRNumber &b) {315 MPFRNumber result(*this);316 mpfr_hypot(result.value, value, b.value, mpfr_rounding);317 return result;318}319 320MPFRNumber MPFRNumber::log() const {321 MPFRNumber result(*this);322 mpfr_log(result.value, value, mpfr_rounding);323 return result;324}325 326MPFRNumber MPFRNumber::log2() const {327 MPFRNumber result(*this);328 mpfr_log2(result.value, value, mpfr_rounding);329 return result;330}331 332MPFRNumber MPFRNumber::log10() const {333 MPFRNumber result(*this);334 mpfr_log10(result.value, value, mpfr_rounding);335 return result;336}337 338MPFRNumber MPFRNumber::log1p() const {339 MPFRNumber result(*this);340 mpfr_log1p(result.value, value, mpfr_rounding);341 return result;342}343 344MPFRNumber MPFRNumber::pow(const MPFRNumber &b) {345 MPFRNumber result(*this);346 mpfr_pow(result.value, value, b.value, mpfr_rounding);347 return result;348}349 350MPFRNumber MPFRNumber::remquo(const MPFRNumber &divisor, int "ient) {351 MPFRNumber remainder(*this);352 long q;353 mpfr_remquo(remainder.value, &q, value, divisor.value, mpfr_rounding);354 quotient = static_cast<int>(q);355 return remainder;356}357 358MPFRNumber MPFRNumber::round() const {359 MPFRNumber result(*this);360 mpfr_round(result.value, value);361 return result;362}363 364MPFRNumber MPFRNumber::roundeven() const {365 MPFRNumber result(*this);366#if MPFR_VERSION_MAJOR >= 4367 mpfr_roundeven(result.value, value);368#else369 mpfr_rint(result.value, value, MPFR_RNDN);370#endif371 return result;372}373 374bool MPFRNumber::round_to_long(long &result) const {375 // We first calculate the rounded value. This way, when converting376 // to long using mpfr_get_si, the rounding direction of MPFR_RNDN377 // (or any other rounding mode), does not have an influence.378 MPFRNumber roundedValue = round();379 mpfr_clear_erangeflag();380 result = mpfr_get_si(roundedValue.value, MPFR_RNDN);381 return mpfr_erangeflag_p();382}383 384bool MPFRNumber::round_to_long(mpfr_rnd_t rnd, long &result) const {385 MPFRNumber rint_result(*this);386 mpfr_rint(rint_result.value, value, rnd);387 return rint_result.round_to_long(result);388}389 390MPFRNumber MPFRNumber::rint(mpfr_rnd_t rnd) const {391 MPFRNumber result(*this);392 mpfr_rint(result.value, value, rnd);393 return result;394}395 396MPFRNumber MPFRNumber::rsqrt() const {397 MPFRNumber result(*this);398 mpfr_rec_sqrt(result.value, value, mpfr_rounding);399 return result;400}401 402MPFRNumber MPFRNumber::mod_2pi() const {403 MPFRNumber result(0.0, 1280);404 MPFRNumber _2pi(0.0, 1280);405 mpfr_const_pi(_2pi.value, MPFR_RNDN);406 mpfr_mul_si(_2pi.value, _2pi.value, 2, MPFR_RNDN);407 mpfr_fmod(result.value, value, _2pi.value, MPFR_RNDN);408 return result;409}410 411MPFRNumber MPFRNumber::mod_pi_over_2() const {412 MPFRNumber result(0.0, 1280);413 MPFRNumber pi_over_2(0.0, 1280);414 mpfr_const_pi(pi_over_2.value, MPFR_RNDN);415 mpfr_mul_d(pi_over_2.value, pi_over_2.value, 0.5, MPFR_RNDN);416 mpfr_fmod(result.value, value, pi_over_2.value, MPFR_RNDN);417 return result;418}419 420MPFRNumber MPFRNumber::mod_pi_over_4() const {421 MPFRNumber result(0.0, 1280);422 MPFRNumber pi_over_4(0.0, 1280);423 mpfr_const_pi(pi_over_4.value, MPFR_RNDN);424 mpfr_mul_d(pi_over_4.value, pi_over_4.value, 0.25, MPFR_RNDN);425 mpfr_fmod(result.value, value, pi_over_4.value, MPFR_RNDN);426 return result;427}428 429MPFRNumber MPFRNumber::sin() const {430 MPFRNumber result(*this);431 mpfr_sin(result.value, value, mpfr_rounding);432 return result;433}434 435MPFRNumber MPFRNumber::sinpi() const {436 MPFRNumber result(*this);437 438#if MPFR_VERSION >= MPFR_VERSION_NUM(4, 2, 0)439 mpfr_sinpi(result.value, value, mpfr_rounding);440 return result;441#else442 if (mpfr_integer_p(value)) {443 mpfr_set_si(result.value, 0, mpfr_rounding);444 return result;445 }446 447 MPFRNumber value_mul_two(*this);448 mpfr_mul_si(value_mul_two.value, value, 2, MPFR_RNDN);449 450 if (mpfr_integer_p(value_mul_two.value)) {451 auto d = mpfr_get_si(value, MPFR_RNDD);452 mpfr_set_si(result.value, (d & 1) ? -1 : 1, mpfr_rounding);453 return result;454 }455 456 MPFRNumber value_pi(0.0, 1280);457 mpfr_const_pi(value_pi.value, MPFR_RNDN);458 mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN);459 mpfr_sin(result.value, value_pi.value, mpfr_rounding);460 return result;461#endif462}463 464MPFRNumber MPFRNumber::sinh() const {465 MPFRNumber result(*this);466 mpfr_sinh(result.value, value, mpfr_rounding);467 return result;468}469 470MPFRNumber MPFRNumber::sqrt() const {471 MPFRNumber result(*this);472 mpfr_sqrt(result.value, value, mpfr_rounding);473 return result;474}475 476MPFRNumber MPFRNumber::sub(const MPFRNumber &b) const {477 MPFRNumber result(*this);478 mpfr_sub(result.value, value, b.value, mpfr_rounding);479 return result;480}481 482MPFRNumber MPFRNumber::tan() const {483 MPFRNumber result(*this);484 mpfr_tan(result.value, value, mpfr_rounding);485 return result;486}487 488MPFRNumber MPFRNumber::tanh() const {489 MPFRNumber result(*this);490 mpfr_tanh(result.value, value, mpfr_rounding);491 return result;492}493 494MPFRNumber MPFRNumber::tanpi() const {495 MPFRNumber result(*this);496 497#if MPFR_VERSION >= MPFR_VERSION_NUM(4, 2, 0)498 mpfr_tanpi(result.value, value, mpfr_rounding);499 return result;500#else501 MPFRNumber value_ret_exact(*this);502 MPFRNumber value_one(*this);503 mpfr_set_si(value_one.value, 1, MPFR_RNDN);504 mpfr_fmod(value_ret_exact.value, value, value_one.value, mpfr_rounding);505 mpfr_mul_si(value_ret_exact.value, value_ret_exact.value, 4, MPFR_RNDN);506 507 if (mpfr_integer_p(value_ret_exact.value)) {508 int mod = mpfr_get_si(value_ret_exact.value, MPFR_RNDN);509 mod = (mod < 0 ? -1 * mod : mod);510 511 switch (mod) {512 case 0:513 mpfr_set_si(result.value, 0, mpfr_rounding);514 break;515 case 1:516 mpfr_set_si(result.value, (mpfr_signbit(value) ? -1 : 1), mpfr_rounding);517 break;518 case 2: {519 auto d = mpfr_get_si(value, MPFR_RNDZ);520 d += mpfr_sgn(value) > 0 ? 0 : 1;521 mpfr_set_inf(result.value, (d & 1) ? -1 : 1);522 break;523 }524 case 3:525 mpfr_set_si(result.value, (mpfr_signbit(value) ? 1 : -1), mpfr_rounding);526 break;527 }528 529 return result;530 }531 532 MPFRNumber value_pi(0.0, 1280);533 mpfr_const_pi(value_pi.value, MPFR_RNDN);534 mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN);535 mpfr_tan(result.value, value_pi.value, mpfr_rounding);536 return result;537#endif538}539 540MPFRNumber MPFRNumber::trunc() const {541 MPFRNumber result(*this);542 mpfr_trunc(result.value, value);543 return result;544}545 546MPFRNumber MPFRNumber::fma(const MPFRNumber &b, const MPFRNumber &c) {547 MPFRNumber result(*this);548 mpfr_fma(result.value, value, b.value, c.value, mpfr_rounding);549 return result;550}551 552MPFRNumber MPFRNumber::mul(const MPFRNumber &b) {553 MPFRNumber result(*this);554 mpfr_mul(result.value, value, b.value, mpfr_rounding);555 return result;556}557 558cpp::string MPFRNumber::str() const {559 // 200 bytes should be more than sufficient to hold a 100-digit number560 // plus additional bytes for the decimal point, '-' sign etc.561 constexpr size_t printBufSize = 200;562 char buffer[printBufSize];563 mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);564 cpp::string_view view(buffer);565 // Trim whitespaces566 const char whitespace = ' ';567 while (!view.empty() && view.front() == whitespace)568 view.remove_prefix(1);569 while (!view.empty() && view.back() == whitespace)570 view.remove_suffix(1);571 return cpp::string(view.data());572}573 574void MPFRNumber::dump(const char *msg) const {575 mpfr_printf("%s%.128g\n", msg, value);576}577 578template <> float MPFRNumber::as<float>() const {579 return mpfr_get_flt(value, mpfr_rounding);580}581 582template <> double MPFRNumber::as<double>() const {583 return mpfr_get_d(value, mpfr_rounding);584}585 586template <> long double MPFRNumber::as<long double>() const {587 return mpfr_get_ld(value, mpfr_rounding);588}589 590#ifdef LIBC_TYPES_HAS_FLOAT16591template <> float16 MPFRNumber::as<float16>() const {592 // TODO: Either prove that this cast won't cause double-rounding errors, or593 // find a better way to get a float16.594 return fputil::cast<float16>(mpfr_get_d(value, mpfr_rounding));595}596#endif597 598#ifdef LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE599template <> float128 MPFRNumber::as<float128>() const {600 return mpfr_get_float128(value, mpfr_rounding);601}602#endif // LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE603 604template <> bfloat16 MPFRNumber::as<bfloat16>() const {605 return fputil::cast<bfloat16>(mpfr_get_flt(value, mpfr_rounding));606}607 608} // namespace mpfr609} // namespace testing610} // namespace LIBC_NAMESPACE_DECL611