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