343 lines · cpp
1//===-- Utils which wrap MPC ----------------------------------------------===//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 "MPCUtils.h"10 11#include "hdr/stdint_proxy.h"12#include "src/__support/CPP/array.h"13#include "src/__support/CPP/stringstream.h"14#include "utils/MPCWrapper/mpc_inc.h"15#include "utils/MPFRWrapper/MPCommon.h"16 17template <typename T> using FPBits = LIBC_NAMESPACE::fputil::FPBits<T>;18 19namespace LIBC_NAMESPACE_DECL {20namespace testing {21namespace mpc {22 23static inline cpp::string str(RoundingMode mode) {24 switch (mode) {25 case RoundingMode::Upward:26 return "MPFR_RNDU";27 case RoundingMode::Downward:28 return "MPFR_RNDD";29 case RoundingMode::TowardZero:30 return "MPFR_RNDZ";31 case RoundingMode::Nearest:32 return "MPFR_RNDN";33 }34}35 36class MPCNumber {37private:38 unsigned int precision;39 mpc_t value;40 mpc_rnd_t mpc_rounding;41 42public:43 explicit MPCNumber(unsigned int p) : precision(p), mpc_rounding(MPC_RNDNN) {44 mpc_init2(value, precision);45 }46 47 MPCNumber() : precision(256), mpc_rounding(MPC_RNDNN) {48 mpc_init2(value, 256);49 }50 51 MPCNumber(unsigned int p, mpc_rnd_t rnd) : precision(p), mpc_rounding(rnd) {52 mpc_init2(value, precision);53 }54 55 template <typename XType,56 cpp::enable_if_t<cpp::is_same_v<_Complex float, XType>, bool> = 0>57 MPCNumber(XType x,58 unsigned int precision = mpfr::ExtraPrecision<float>::VALUE,59 RoundingMode rnd = RoundingMode::Nearest)60 : precision(precision),61 mpc_rounding(MPC_RND(mpfr::get_mpfr_rounding_mode(rnd),62 mpfr::get_mpfr_rounding_mode(rnd))) {63 mpc_init2(value, precision);64 Complex<float> x_c = cpp::bit_cast<Complex<float>>(x);65 mpfr_t real, imag;66 mpfr_init2(real, precision);67 mpfr_init2(imag, precision);68 mpfr_set_flt(real, x_c.real, mpfr::get_mpfr_rounding_mode(rnd));69 mpfr_set_flt(imag, x_c.imag, mpfr::get_mpfr_rounding_mode(rnd));70 mpc_set_fr_fr(value, real, imag, mpc_rounding);71 mpfr_clear(real);72 mpfr_clear(imag);73 }74 75 template <typename XType,76 cpp::enable_if_t<cpp::is_same_v<_Complex double, XType>, bool> = 0>77 MPCNumber(XType x,78 unsigned int precision = mpfr::ExtraPrecision<double>::VALUE,79 RoundingMode rnd = RoundingMode::Nearest)80 : precision(precision),81 mpc_rounding(MPC_RND(mpfr::get_mpfr_rounding_mode(rnd),82 mpfr::get_mpfr_rounding_mode(rnd))) {83 mpc_init2(value, precision);84 Complex<double> x_c = cpp::bit_cast<Complex<double>>(x);85 mpc_set_d_d(value, x_c.real, x_c.imag, mpc_rounding);86 }87 88 MPCNumber(const MPCNumber &other)89 : precision(other.precision), mpc_rounding(other.mpc_rounding) {90 mpc_init2(value, precision);91 mpc_set(value, other.value, mpc_rounding);92 }93 94 ~MPCNumber() { mpc_clear(value); }95 96 MPCNumber &operator=(const MPCNumber &rhs) {97 precision = rhs.precision;98 mpc_rounding = rhs.mpc_rounding;99 mpc_init2(value, precision);100 mpc_set(value, rhs.value, mpc_rounding);101 return *this;102 }103 104 void setValue(mpc_t val) const { mpc_set(val, value, mpc_rounding); }105 106 mpc_t &getValue() { return value; }107 108 MPCNumber carg() const {109 mpfr_t res;110 MPCNumber result(precision, mpc_rounding);111 112 mpfr_init2(res, precision);113 114 mpc_arg(res, value, MPC_RND_RE(mpc_rounding));115 mpc_set_fr(result.value, res, mpc_rounding);116 117 mpfr_clear(res);118 119 return result;120 }121 122 MPCNumber cproj() const {123 MPCNumber result(precision, mpc_rounding);124 mpc_proj(result.value, value, mpc_rounding);125 return result;126 }127};128 129namespace internal {130 131template <typename InputType>132cpp::enable_if_t<cpp::is_complex_v<InputType>, MPCNumber>133unary_operation(Operation op, InputType input, unsigned int precision,134 RoundingMode rounding) {135 MPCNumber mpcInput(input, precision, rounding);136 switch (op) {137 case Operation::Carg:138 return mpcInput.carg();139 case Operation::Cproj:140 return mpcInput.cproj();141 default:142 __builtin_unreachable();143 }144}145 146template <typename InputType, typename OutputType>147bool compare_unary_operation_single_output_same_type(Operation op,148 InputType input,149 OutputType libc_result,150 double ulp_tolerance,151 RoundingMode rounding) {152 153 unsigned int precision =154 mpfr::get_precision<make_real_t<InputType>>(ulp_tolerance);155 156 MPCNumber mpc_result;157 mpc_result = unary_operation(op, input, precision, rounding);158 159 mpc_t mpc_result_val;160 mpc_init2(mpc_result_val, precision);161 mpc_result.setValue(mpc_result_val);162 163 mpfr_t real, imag;164 mpfr_init2(real, precision);165 mpfr_init2(imag, precision);166 mpc_real(real, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding));167 mpc_imag(imag, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding));168 169 mpfr::MPFRNumber mpfr_real(real, precision, rounding);170 mpfr::MPFRNumber mpfr_imag(imag, precision, rounding);171 172 double ulp_real = mpfr_real.ulp(173 (cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result)).real);174 double ulp_imag = mpfr_imag.ulp(175 (cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result)).imag);176 mpc_clear(mpc_result_val);177 mpfr_clear(real);178 mpfr_clear(imag);179 return (ulp_real <= ulp_tolerance) && (ulp_imag <= ulp_tolerance);180}181 182template bool compare_unary_operation_single_output_same_type(183 Operation, _Complex float, _Complex float, double, RoundingMode);184template bool compare_unary_operation_single_output_same_type(185 Operation, _Complex double, _Complex double, double, RoundingMode);186 187template <typename InputType, typename OutputType>188bool compare_unary_operation_single_output_different_type(189 Operation op, InputType input, OutputType libc_result, double ulp_tolerance,190 RoundingMode rounding) {191 192 unsigned int precision =193 mpfr::get_precision<make_real_t<InputType>>(ulp_tolerance);194 195 MPCNumber mpc_result;196 mpc_result = unary_operation(op, input, precision, rounding);197 198 mpc_t mpc_result_val;199 mpc_init2(mpc_result_val, precision);200 mpc_result.setValue(mpc_result_val);201 202 mpfr_t real;203 mpfr_init2(real, precision);204 mpc_real(real, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding));205 206 mpfr::MPFRNumber mpfr_real(real, precision, rounding);207 208 double ulp_real = mpfr_real.ulp(libc_result);209 mpc_clear(mpc_result_val);210 mpfr_clear(real);211 return (ulp_real <= ulp_tolerance);212}213 214template bool compare_unary_operation_single_output_different_type(215 Operation, _Complex float, float, double, RoundingMode);216template bool compare_unary_operation_single_output_different_type(217 Operation, _Complex double, double, double, RoundingMode);218 219template <typename InputType, typename OutputType>220void explain_unary_operation_single_output_different_type_error(221 Operation op, InputType input, OutputType libc_result, double ulp_tolerance,222 RoundingMode rounding) {223 224 unsigned int precision =225 mpfr::get_precision<make_real_t<InputType>>(ulp_tolerance);226 227 MPCNumber mpc_result;228 mpc_result = unary_operation(op, input, precision, rounding);229 230 mpc_t mpc_result_val;231 mpc_init2(mpc_result_val, precision);232 mpc_result.setValue(mpc_result_val);233 234 mpfr_t real;235 mpfr_init2(real, precision);236 mpc_real(real, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding));237 238 mpfr::MPFRNumber mpfr_result(real, precision, rounding);239 mpfr::MPFRNumber mpfrLibcResult(libc_result, precision, rounding);240 mpfr::MPFRNumber mpfrInputReal(241 cpp::bit_cast<Complex<make_real_t<InputType>>>(input).real, precision,242 rounding);243 mpfr::MPFRNumber mpfrInputImag(244 cpp::bit_cast<Complex<make_real_t<InputType>>>(input).imag, precision,245 rounding);246 247 cpp::array<char, 2048> msg_buf;248 cpp::StringStream msg(msg_buf);249 msg << "Match value not within tolerance value of MPFR result:\n"250 << " Input: " << mpfrInputReal.str() << " + " << mpfrInputImag.str()251 << "i\n"252 << " Rounding mode: " << str(rounding) << '\n'253 << " Libc: " << mpfrLibcResult.str() << '\n'254 << " MPC: " << mpfr_result.str() << '\n'255 << '\n'256 << " ULP error: " << mpfr_result.ulp_as_mpfr_number(libc_result).str()257 << '\n';258 tlog << msg.str();259 mpc_clear(mpc_result_val);260 mpfr_clear(real);261}262 263template void explain_unary_operation_single_output_different_type_error(264 Operation, _Complex float, float, double, RoundingMode);265template void explain_unary_operation_single_output_different_type_error(266 Operation, _Complex double, double, double, RoundingMode);267 268template <typename InputType, typename OutputType>269void explain_unary_operation_single_output_same_type_error(270 Operation op, InputType input, OutputType libc_result, double ulp_tolerance,271 RoundingMode rounding) {272 273 unsigned int precision =274 mpfr::get_precision<make_real_t<InputType>>(ulp_tolerance);275 276 MPCNumber mpc_result;277 mpc_result = unary_operation(op, input, precision, rounding);278 279 mpc_t mpc_result_val;280 mpc_init2(mpc_result_val, precision);281 mpc_result.setValue(mpc_result_val);282 283 mpfr_t real, imag;284 mpfr_init2(real, precision);285 mpfr_init2(imag, precision);286 mpc_real(real, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding));287 mpc_imag(imag, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding));288 289 mpfr::MPFRNumber mpfr_real(real, precision, rounding);290 mpfr::MPFRNumber mpfr_imag(imag, precision, rounding);291 mpfr::MPFRNumber mpfrLibcResultReal(292 cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result).real,293 precision, rounding);294 mpfr::MPFRNumber mpfrLibcResultImag(295 cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result).imag,296 precision, rounding);297 mpfr::MPFRNumber mpfrInputReal(298 cpp::bit_cast<Complex<make_real_t<InputType>>>(input).real, precision,299 rounding);300 mpfr::MPFRNumber mpfrInputImag(301 cpp::bit_cast<Complex<make_real_t<InputType>>>(input).imag, precision,302 rounding);303 304 cpp::array<char, 2048> msg_buf;305 cpp::StringStream msg(msg_buf);306 msg << "Match value not within tolerance value of MPFR result:\n"307 << " Input: " << mpfrInputReal.str() << " + " << mpfrInputImag.str()308 << "i\n"309 << " Rounding mode: " << str(rounding) << " , " << str(rounding) << '\n'310 << " Libc: " << mpfrLibcResultReal.str() << " + "311 << mpfrLibcResultImag.str() << "i\n"312 << " MPC: " << mpfr_real.str() << " + " << mpfr_imag.str() << "i\n"313 << '\n'314 << " ULP error: "315 << mpfr_real316 .ulp_as_mpfr_number(317 cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result)318 .real)319 .str()320 << " , "321 << mpfr_imag322 .ulp_as_mpfr_number(323 cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result)324 .imag)325 .str()326 << '\n';327 tlog << msg.str();328 mpc_clear(mpc_result_val);329 mpfr_clear(real);330 mpfr_clear(imag);331}332 333template void explain_unary_operation_single_output_same_type_error(334 Operation, _Complex float, _Complex float, double, RoundingMode);335template void explain_unary_operation_single_output_same_type_error(336 Operation, _Complex double, _Complex double, double, RoundingMode);337 338} // namespace internal339 340} // namespace mpc341} // namespace testing342} // namespace LIBC_NAMESPACE_DECL343