| //===-- Utils which wrap MPC ----------------------------------------------===// | |
| // | |
| // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. | |
| // See https://llvm.org/LICENSE.txt for license information. | |
| // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception | |
| // | |
| //===----------------------------------------------------------------------===// | |
| #include "MPCUtils.h" | |
| #include "hdr/stdint_proxy.h" | |
| #include "src/__support/CPP/array.h" | |
| #include "src/__support/CPP/stringstream.h" | |
| #include "utils/MPCWrapper/mpc_inc.h" | |
| #include "utils/MPFRWrapper/MPCommon.h" | |
| template <typename T> using FPBits = LIBC_NAMESPACE::fputil::FPBits<T>; | |
| namespace LIBC_NAMESPACE_DECL { | |
| namespace testing { | |
| namespace mpc { | |
| static inline cpp::string str(RoundingMode mode) { | |
| switch (mode) { | |
| case RoundingMode::Upward: | |
| return "MPFR_RNDU"; | |
| case RoundingMode::Downward: | |
| return "MPFR_RNDD"; | |
| case RoundingMode::TowardZero: | |
| return "MPFR_RNDZ"; | |
| case RoundingMode::Nearest: | |
| return "MPFR_RNDN"; | |
| } | |
| } | |
| class MPCNumber { | |
| private: | |
| unsigned int precision; | |
| mpc_t value; | |
| mpc_rnd_t mpc_rounding; | |
| public: | |
| explicit MPCNumber(unsigned int p) : precision(p), mpc_rounding(MPC_RNDNN) { | |
| mpc_init2(value, precision); | |
| } | |
| MPCNumber() : precision(256), mpc_rounding(MPC_RNDNN) { | |
| mpc_init2(value, 256); | |
| } | |
| MPCNumber(unsigned int p, mpc_rnd_t rnd) : precision(p), mpc_rounding(rnd) { | |
| mpc_init2(value, precision); | |
| } | |
| template <typename XType, | |
| cpp::enable_if_t<cpp::is_same_v<_Complex float, XType>, bool> = 0> | |
| MPCNumber(XType x, | |
| unsigned int precision = mpfr::ExtraPrecision<float>::VALUE, | |
| RoundingMode rnd = RoundingMode::Nearest) | |
| : precision(precision), | |
| mpc_rounding(MPC_RND(mpfr::get_mpfr_rounding_mode(rnd), | |
| mpfr::get_mpfr_rounding_mode(rnd))) { | |
| mpc_init2(value, precision); | |
| Complex<float> x_c = cpp::bit_cast<Complex<float>>(x); | |
| mpfr_t real, imag; | |
| mpfr_init2(real, precision); | |
| mpfr_init2(imag, precision); | |
| mpfr_set_flt(real, x_c.real, mpfr::get_mpfr_rounding_mode(rnd)); | |
| mpfr_set_flt(imag, x_c.imag, mpfr::get_mpfr_rounding_mode(rnd)); | |
| mpc_set_fr_fr(value, real, imag, mpc_rounding); | |
| mpfr_clear(real); | |
| mpfr_clear(imag); | |
| } | |
| template <typename XType, | |
| cpp::enable_if_t<cpp::is_same_v<_Complex double, XType>, bool> = 0> | |
| MPCNumber(XType x, | |
| unsigned int precision = mpfr::ExtraPrecision<double>::VALUE, | |
| RoundingMode rnd = RoundingMode::Nearest) | |
| : precision(precision), | |
| mpc_rounding(MPC_RND(mpfr::get_mpfr_rounding_mode(rnd), | |
| mpfr::get_mpfr_rounding_mode(rnd))) { | |
| mpc_init2(value, precision); | |
| Complex<double> x_c = cpp::bit_cast<Complex<double>>(x); | |
| mpc_set_d_d(value, x_c.real, x_c.imag, mpc_rounding); | |
| } | |
| MPCNumber(const MPCNumber &other) | |
| : precision(other.precision), mpc_rounding(other.mpc_rounding) { | |
| mpc_init2(value, precision); | |
| mpc_set(value, other.value, mpc_rounding); | |
| } | |
| ~MPCNumber() { mpc_clear(value); } | |
| MPCNumber &operator=(const MPCNumber &rhs) { | |
| precision = rhs.precision; | |
| mpc_rounding = rhs.mpc_rounding; | |
| mpc_init2(value, precision); | |
| mpc_set(value, rhs.value, mpc_rounding); | |
| return *this; | |
| } | |
| void setValue(mpc_t val) const { mpc_set(val, value, mpc_rounding); } | |
| mpc_t &getValue() { return value; } | |
| MPCNumber cabs() const { | |
| mpfr_t res; | |
| MPCNumber result(precision, mpc_rounding); | |
| mpfr_init2(res, precision); | |
| mpc_abs(res, value, MPC_RND_RE(mpc_rounding)); | |
| mpc_set_fr(result.value, res, mpc_rounding); | |
| mpfr_clear(res); | |
| return result; | |
| } | |
| MPCNumber carg() const { | |
| mpfr_t res; | |
| MPCNumber result(precision, mpc_rounding); | |
| mpfr_init2(res, precision); | |
| mpc_arg(res, value, MPC_RND_RE(mpc_rounding)); | |
| mpc_set_fr(result.value, res, mpc_rounding); | |
| mpfr_clear(res); | |
| return result; | |
| } | |
| MPCNumber cproj() const { | |
| MPCNumber result(precision, mpc_rounding); | |
| mpc_proj(result.value, value, mpc_rounding); | |
| return result; | |
| } | |
| }; | |
| namespace internal { | |
| template <typename InputType> | |
| cpp::enable_if_t<cpp::is_complex_v<InputType>, MPCNumber> | |
| unary_operation(Operation op, InputType input, unsigned int precision, | |
| RoundingMode rounding) { | |
| MPCNumber mpcInput(input, precision, rounding); | |
| switch (op) { | |
| case Operation::Cabs: | |
| return mpcInput.cabs(); | |
| case Operation::Carg: | |
| return mpcInput.carg(); | |
| case Operation::Cproj: | |
| return mpcInput.cproj(); | |
| default: | |
| __builtin_unreachable(); | |
| } | |
| } | |
| template <typename InputType, typename OutputType> | |
| bool compare_unary_operation_single_output_same_type(Operation op, | |
| InputType input, | |
| OutputType libc_result, | |
| double ulp_tolerance, | |
| RoundingMode rounding) { | |
| unsigned int precision = | |
| mpfr::get_precision<make_real_t<InputType>>(ulp_tolerance); | |
| MPCNumber mpc_result; | |
| mpc_result = unary_operation(op, input, precision, rounding); | |
| mpc_t mpc_result_val; | |
| mpc_init2(mpc_result_val, precision); | |
| mpc_result.setValue(mpc_result_val); | |
| mpfr_t real, imag; | |
| mpfr_init2(real, precision); | |
| mpfr_init2(imag, precision); | |
| mpc_real(real, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding)); | |
| mpc_imag(imag, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding)); | |
| mpfr::MPFRNumber mpfr_real(real, precision, rounding); | |
| mpfr::MPFRNumber mpfr_imag(imag, precision, rounding); | |
| double ulp_real = mpfr_real.ulp( | |
| (cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result)).real); | |
| double ulp_imag = mpfr_imag.ulp( | |
| (cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result)).imag); | |
| mpc_clear(mpc_result_val); | |
| mpfr_clear(real); | |
| mpfr_clear(imag); | |
| return (ulp_real <= ulp_tolerance) && (ulp_imag <= ulp_tolerance); | |
| } | |
| template bool compare_unary_operation_single_output_same_type( | |
| Operation, _Complex float, _Complex float, double, RoundingMode); | |
| template bool compare_unary_operation_single_output_same_type( | |
| Operation, _Complex double, _Complex double, double, RoundingMode); | |
| template <typename InputType, typename OutputType> | |
| bool compare_unary_operation_single_output_different_type( | |
| Operation op, InputType input, OutputType libc_result, double ulp_tolerance, | |
| RoundingMode rounding) { | |
| unsigned int precision = | |
| mpfr::get_precision<make_real_t<InputType>>(ulp_tolerance); | |
| MPCNumber mpc_result; | |
| mpc_result = unary_operation(op, input, precision, rounding); | |
| mpc_t mpc_result_val; | |
| mpc_init2(mpc_result_val, precision); | |
| mpc_result.setValue(mpc_result_val); | |
| mpfr_t real; | |
| mpfr_init2(real, precision); | |
| mpc_real(real, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding)); | |
| mpfr::MPFRNumber mpfr_real(real, precision, rounding); | |
| double ulp_real = mpfr_real.ulp(libc_result); | |
| mpc_clear(mpc_result_val); | |
| mpfr_clear(real); | |
| return (ulp_real <= ulp_tolerance); | |
| } | |
| template bool compare_unary_operation_single_output_different_type( | |
| Operation, _Complex float, float, double, RoundingMode); | |
| template bool compare_unary_operation_single_output_different_type( | |
| Operation, _Complex double, double, double, RoundingMode); | |
| template <typename InputType, typename OutputType> | |
| void explain_unary_operation_single_output_different_type_error( | |
| Operation op, InputType input, OutputType libc_result, double ulp_tolerance, | |
| RoundingMode rounding) { | |
| unsigned int precision = | |
| mpfr::get_precision<make_real_t<InputType>>(ulp_tolerance); | |
| MPCNumber mpc_result; | |
| mpc_result = unary_operation(op, input, precision, rounding); | |
| mpc_t mpc_result_val; | |
| mpc_init2(mpc_result_val, precision); | |
| mpc_result.setValue(mpc_result_val); | |
| mpfr_t real; | |
| mpfr_init2(real, precision); | |
| mpc_real(real, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding)); | |
| mpfr::MPFRNumber mpfr_result(real, precision, rounding); | |
| mpfr::MPFRNumber mpfrLibcResult(libc_result, precision, rounding); | |
| mpfr::MPFRNumber mpfrInputReal( | |
| cpp::bit_cast<Complex<make_real_t<InputType>>>(input).real, precision, | |
| rounding); | |
| mpfr::MPFRNumber mpfrInputImag( | |
| cpp::bit_cast<Complex<make_real_t<InputType>>>(input).imag, precision, | |
| rounding); | |
| cpp::array<char, 2048> msg_buf; | |
| cpp::StringStream msg(msg_buf); | |
| msg << "Match value not within tolerance value of MPFR result:\n" | |
| << " Input: " << mpfrInputReal.str() << " + " << mpfrInputImag.str() | |
| << "i\n" | |
| << " Rounding mode: " << str(rounding) << '\n' | |
| << " Libc: " << mpfrLibcResult.str() << '\n' | |
| << " MPC: " << mpfr_result.str() << '\n' | |
| << '\n' | |
| << " ULP error: " << mpfr_result.ulp_as_mpfr_number(libc_result).str() | |
| << '\n'; | |
| tlog << msg.str(); | |
| mpc_clear(mpc_result_val); | |
| mpfr_clear(real); | |
| } | |
| template void explain_unary_operation_single_output_different_type_error( | |
| Operation, _Complex float, float, double, RoundingMode); | |
| template void explain_unary_operation_single_output_different_type_error( | |
| Operation, _Complex double, double, double, RoundingMode); | |
| template <typename InputType, typename OutputType> | |
| void explain_unary_operation_single_output_same_type_error( | |
| Operation op, InputType input, OutputType libc_result, double ulp_tolerance, | |
| RoundingMode rounding) { | |
| unsigned int precision = | |
| mpfr::get_precision<make_real_t<InputType>>(ulp_tolerance); | |
| MPCNumber mpc_result; | |
| mpc_result = unary_operation(op, input, precision, rounding); | |
| mpc_t mpc_result_val; | |
| mpc_init2(mpc_result_val, precision); | |
| mpc_result.setValue(mpc_result_val); | |
| mpfr_t real, imag; | |
| mpfr_init2(real, precision); | |
| mpfr_init2(imag, precision); | |
| mpc_real(real, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding)); | |
| mpc_imag(imag, mpc_result_val, mpfr::get_mpfr_rounding_mode(rounding)); | |
| mpfr::MPFRNumber mpfr_real(real, precision, rounding); | |
| mpfr::MPFRNumber mpfr_imag(imag, precision, rounding); | |
| mpfr::MPFRNumber mpfrLibcResultReal( | |
| cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result).real, | |
| precision, rounding); | |
| mpfr::MPFRNumber mpfrLibcResultImag( | |
| cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result).imag, | |
| precision, rounding); | |
| mpfr::MPFRNumber mpfrInputReal( | |
| cpp::bit_cast<Complex<make_real_t<InputType>>>(input).real, precision, | |
| rounding); | |
| mpfr::MPFRNumber mpfrInputImag( | |
| cpp::bit_cast<Complex<make_real_t<InputType>>>(input).imag, precision, | |
| rounding); | |
| cpp::array<char, 2048> msg_buf; | |
| cpp::StringStream msg(msg_buf); | |
| msg << "Match value not within tolerance value of MPFR result:\n" | |
| << " Input: " << mpfrInputReal.str() << " + " << mpfrInputImag.str() | |
| << "i\n" | |
| << " Rounding mode: " << str(rounding) << " , " << str(rounding) << '\n' | |
| << " Libc: " << mpfrLibcResultReal.str() << " + " | |
| << mpfrLibcResultImag.str() << "i\n" | |
| << " MPC: " << mpfr_real.str() << " + " << mpfr_imag.str() << "i\n" | |
| << '\n' | |
| << " ULP error: " | |
| << mpfr_real | |
| .ulp_as_mpfr_number( | |
| cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result) | |
| .real) | |
| .str() | |
| << " , " | |
| << mpfr_imag | |
| .ulp_as_mpfr_number( | |
| cpp::bit_cast<Complex<make_real_t<InputType>>>(libc_result) | |
| .imag) | |
| .str() | |
| << '\n'; | |
| tlog << msg.str(); | |
| mpc_clear(mpc_result_val); | |
| mpfr_clear(real); | |
| mpfr_clear(imag); | |
| } | |
| template void explain_unary_operation_single_output_same_type_error( | |
| Operation, _Complex float, _Complex float, double, RoundingMode); | |
| template void explain_unary_operation_single_output_same_type_error( | |
| Operation, _Complex double, _Complex double, double, RoundingMode); | |
| } // namespace internal | |
| } // namespace mpc | |
| } // namespace testing | |
| } // namespace LIBC_NAMESPACE_DECL |