[libc][complex] Testing infra for MPC (#121261)

This PR aims to add the groundwork to test the precision of libc complex
functions against MPC. I took `cargf` as a test to verify that the infra
works fine.
This commit is contained in:
Shourya Goel 2025-01-28 11:01:16 +05:30 committed by GitHub
parent 7109f52197
commit 7f37b34d31
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
31 changed files with 1716 additions and 908 deletions

View File

@ -262,6 +262,7 @@ if(LIBC_TARGET_OS_IS_GPU)
endif()
include(LLVMLibCCheckMPFR)
include(LLVMLibCCheckMPC)
if(LLVM_LIBC_CLANG_TIDY)
set(LLVM_LIBC_ENABLE_LINTING ON)

View File

@ -0,0 +1,22 @@
if(LIBC_TESTS_CAN_USE_MPFR)
set(LLVM_LIBC_MPC_INSTALL_PATH "" CACHE PATH "Path to where MPC is installed (e.g. C:/src/install or ~/src/install)")
if(LLVM_LIBC_MPC_INSTALL_PATH)
set(LIBC_TESTS_CAN_USE_MPC TRUE)
elseif(LIBC_TARGET_OS_IS_GPU OR LLVM_LIBC_FULL_BUILD)
# In full build mode, the MPC library should be built using our own facilities,
# which is currently not possible.
set(LIBC_TESTS_CAN_USE_MPC FALSE)
else()
try_compile(
LIBC_TESTS_CAN_USE_MPC
${CMAKE_CURRENT_BINARY_DIR}
SOURCES
${LIBC_SOURCE_DIR}/utils/MPCWrapper/check_mpc.cpp
COMPILE_DEFINITIONS
${LIBC_COMPILE_OPTIONS_DEFAULT}
LINK_LIBRARIES
-lmpc -lmpfr -lgmp -latomic
)
endif()
endif()

View File

@ -243,12 +243,20 @@ add_header_library(
HDRS
complex_type.h
DEPENDS
libc.src.__support.CPP.bit
libc.src.__support.FPUtil.fp_bits
libc.src.__support.macros.properties.types
libc.src.__support.macros.properties.complex_types
)
add_header_library(
complex_basic_ops
HDRS
complex_basic_ops.h
DEPENDS
.complex_type
libc.src.__support.CPP.bit
libc.src.__support.FPUtil.fp_bits
)
add_header_library(
integer_operations
HDRS

View File

@ -26,6 +26,7 @@
#include "src/__support/CPP/type_traits/is_array.h"
#include "src/__support/CPP/type_traits/is_base_of.h"
#include "src/__support/CPP/type_traits/is_class.h"
#include "src/__support/CPP/type_traits/is_complex.h"
#include "src/__support/CPP/type_traits/is_const.h"
#include "src/__support/CPP/type_traits/is_constant_evaluated.h"
#include "src/__support/CPP/type_traits/is_convertible.h"

View File

@ -0,0 +1,36 @@
//===-- complex basic operations --------------------------------*- C++ -*-===//
//
// 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
//
//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_COMPLEX_BASIC_OPERATIONS_H
#define LLVM_LIBC_SRC___SUPPORT_COMPLEX_BASIC_OPERATIONS_H
#include "complex_type.h"
#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FPBits.h"
namespace LIBC_NAMESPACE_DECL {
template <typename T> LIBC_INLINE constexpr T conjugate(T c) {
Complex<make_real_t<T>> c_c = cpp::bit_cast<Complex<make_real_t<T>>>(c);
c_c.imag = -c_c.imag;
return cpp::bit_cast<T>(c_c);
}
template <typename T> LIBC_INLINE constexpr T project(T c) {
using real_t = make_real_t<T>;
Complex<real_t> c_c = cpp::bit_cast<Complex<real_t>>(c);
if (fputil::FPBits<real_t>(c_c.real).is_inf() ||
fputil::FPBits<real_t>(c_c.imag).is_inf())
return cpp::bit_cast<T>(
Complex<real_t>{(fputil::FPBits<real_t>::inf(Sign::POS).get_val()),
static_cast<real_t>(c_c.imag > 0 ? 0.0 : -0.0)});
return c;
}
} // namespace LIBC_NAMESPACE_DECL
#endif // LLVM_LIBC_SRC___SUPPORT_COMPLEX_BASIC_OPERATIONS_H

View File

@ -9,8 +9,6 @@
#ifndef LLVM_LIBC_SRC___SUPPORT_COMPLEX_TYPE_H
#define LLVM_LIBC_SRC___SUPPORT_COMPLEX_TYPE_H
#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/properties/complex_types.h"
#include "src/__support/macros/properties/types.h"
@ -71,24 +69,5 @@ template <> struct make_real<cfloat128> {
template <typename T> using make_real_t = typename make_real<T>::type;
template <typename T> LIBC_INLINE constexpr T conjugate(T c) {
Complex<make_real_t<T>> c_c = cpp::bit_cast<Complex<make_real_t<T>>>(c);
c_c.imag = -c_c.imag;
return cpp::bit_cast<T>(c_c);
}
template <typename T> LIBC_INLINE constexpr T project(T c) {
using real_t = make_real_t<T>;
Complex<real_t> c_c = cpp::bit_cast<Complex<real_t>>(c);
if (fputil::FPBits<real_t>(c_c.real).is_inf() ||
fputil::FPBits<real_t>(c_c.imag).is_inf()) {
return cpp::bit_cast<T>(
Complex<real_t>{(fputil::FPBits<real_t>::inf(Sign::POS).get_val()),
static_cast<real_t>(c_c.imag > 0 ? 0.0 : -0.0)});
} else {
return c;
}
}
} // namespace LIBC_NAMESPACE_DECL
#endif // LLVM_LIBC_SRC___SUPPORT_COMPLEX_TYPE_H

View File

@ -7,7 +7,7 @@ add_entrypoint_object(
COMPILE_OPTIONS
${libc_opt_high_flag}
DEPENDS
libc.src.__support.complex_type
libc.src.__support.complex_basic_ops
)
add_entrypoint_object(
@ -19,7 +19,7 @@ add_entrypoint_object(
COMPILE_OPTIONS
${libc_opt_high_flag}
DEPENDS
libc.src.__support.complex_type
libc.src.__support.complex_basic_ops
)
add_entrypoint_object(
@ -31,7 +31,7 @@ add_entrypoint_object(
COMPILE_OPTIONS
${libc_opt_high_flag}
DEPENDS
libc.src.__support.complex_type
libc.src.__support.complex_basic_ops
)
add_entrypoint_object(
@ -43,7 +43,7 @@ add_entrypoint_object(
COMPILE_OPTIONS
${libc_opt_high_flag}
DEPENDS
libc.src.__support.complex_type
libc.src.__support.complex_basic_ops
libc.src.__support.macros.properties.types
libc.src.__support.macros.properties.complex_types
)
@ -57,7 +57,7 @@ add_entrypoint_object(
COMPILE_OPTIONS
${libc_opt_high_flag}
DEPENDS
libc.src.__support.complex_type
libc.src.__support.complex_basic_ops
libc.src.__support.macros.properties.types
libc.src.__support.macros.properties.complex_types
)
@ -71,7 +71,7 @@ add_entrypoint_object(
COMPILE_OPTIONS
${libc_opt_high_flag}
DEPENDS
libc.src.__support.complex_type
libc.src.__support.complex_basic_ops
)
add_entrypoint_object(
@ -83,7 +83,7 @@ add_entrypoint_object(
COMPILE_OPTIONS
${libc_opt_high_flag}
DEPENDS
libc.src.__support.complex_type
libc.src.__support.complex_basic_ops
)
add_entrypoint_object(
@ -95,7 +95,7 @@ add_entrypoint_object(
COMPILE_OPTIONS
${libc_opt_high_flag}
DEPENDS
libc.src.__support.complex_type
libc.src.__support.complex_basic_ops
)
add_entrypoint_object(
@ -107,7 +107,7 @@ add_entrypoint_object(
COMPILE_OPTIONS
${libc_opt_high_flag}
DEPENDS
libc.src.__support.complex_type
libc.src.__support.complex_basic_ops
libc.src.__support.macros.properties.types
libc.src.__support.macros.properties.complex_types
)
@ -121,7 +121,7 @@ add_entrypoint_object(
COMPILE_OPTIONS
${libc_opt_high_flag}
DEPENDS
libc.src.__support.complex_type
libc.src.__support.complex_basic_ops
libc.src.__support.macros.properties.types
libc.src.__support.macros.properties.complex_types
)

View File

@ -8,7 +8,7 @@
#include "src/complex/conj.h"
#include "src/__support/common.h"
#include "src/__support/complex_type.h"
#include "src/__support/complex_basic_ops.h"
namespace LIBC_NAMESPACE_DECL {

View File

@ -8,7 +8,7 @@
#include "src/complex/conjf.h"
#include "src/__support/common.h"
#include "src/__support/complex_type.h"
#include "src/__support/complex_basic_ops.h"
namespace LIBC_NAMESPACE_DECL {

View File

@ -8,7 +8,7 @@
#include "src/complex/conjf128.h"
#include "src/__support/common.h"
#include "src/__support/complex_type.h"
#include "src/__support/complex_basic_ops.h"
namespace LIBC_NAMESPACE_DECL {

View File

@ -8,7 +8,7 @@
#include "src/complex/conjf16.h"
#include "src/__support/common.h"
#include "src/__support/complex_type.h"
#include "src/__support/complex_basic_ops.h"
namespace LIBC_NAMESPACE_DECL {

View File

@ -8,7 +8,7 @@
#include "src/complex/conjl.h"
#include "src/__support/common.h"
#include "src/__support/complex_type.h"
#include "src/__support/complex_basic_ops.h"
namespace LIBC_NAMESPACE_DECL {

View File

@ -8,7 +8,7 @@
#include "src/complex/cproj.h"
#include "src/__support/common.h"
#include "src/__support/complex_type.h"
#include "src/__support/complex_basic_ops.h"
namespace LIBC_NAMESPACE_DECL {

View File

@ -8,7 +8,7 @@
#include "src/complex/cprojf.h"
#include "src/__support/common.h"
#include "src/__support/complex_type.h"
#include "src/__support/complex_basic_ops.h"
namespace LIBC_NAMESPACE_DECL {

View File

@ -8,7 +8,7 @@
#include "src/complex/cprojf128.h"
#include "src/__support/common.h"
#include "src/__support/complex_type.h"
#include "src/__support/complex_basic_ops.h"
namespace LIBC_NAMESPACE_DECL {

View File

@ -8,7 +8,7 @@
#include "src/complex/cprojf16.h"
#include "src/__support/common.h"
#include "src/__support/complex_type.h"
#include "src/__support/complex_basic_ops.h"
namespace LIBC_NAMESPACE_DECL {

View File

@ -8,7 +8,7 @@
#include "src/complex/cprojl.h"
#include "src/__support/common.h"
#include "src/__support/complex_type.h"
#include "src/__support/complex_basic_ops.h"
namespace LIBC_NAMESPACE_DECL {

View File

@ -11,7 +11,6 @@
#include "src/__support/CPP/array.h"
#include "src/__support/CPP/type_traits.h"
#include "src/__support/CPP/type_traits/is_complex.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/fpbits_str.h"

View File

@ -1,12 +1,21 @@
function(add_fp_unittest name)
cmake_parse_arguments(
"MATH_UNITTEST"
"NEED_MPFR;UNIT_TEST_ONLY;HERMETIC_TEST_ONLY" # Optional arguments
"NEED_MPFR;NEED_MPC;UNIT_TEST_ONLY;HERMETIC_TEST_ONLY" # Optional arguments
"" # Single value arguments
"LINK_LIBRARIES;DEPENDS" # Multi-value arguments
${ARGN}
)
if(MATH_UNITTEST_NEED_MPC)
set(MATH_UNITTEST_NEED_MPFR TRUE)
if(NOT LIBC_TESTS_CAN_USE_MPC)
message(VERBOSE "Complex test ${name} will be skipped as MPC library is not available.")
return()
endif()
list(APPEND MATH_UNITTEST_LINK_LIBRARIES libcMPCWrapper)
endif()
if(MATH_UNITTEST_NEED_MPFR)
if(NOT LIBC_TESTS_CAN_USE_MPFR)
message(VERBOSE "Math test ${name} will be skipped as MPFR library is not available.")

View File

@ -1,6 +1,6 @@
add_custom_target(libc-complex-unittests)
add_libc_test(
add_fp_unittest(
conj_test
SUITE
libc-complex-unittests
@ -8,11 +8,9 @@ add_libc_test(
conj_test.cpp
DEPENDS
libc.src.complex.conj
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
conjf_test
SUITE
libc-complex-unittests
@ -20,11 +18,9 @@ add_libc_test(
conjf_test.cpp
DEPENDS
libc.src.complex.conjf
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
conjl_test
SUITE
libc-complex-unittests
@ -32,11 +28,9 @@ add_libc_test(
conjl_test.cpp
DEPENDS
libc.src.complex.conjl
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
conjf16_test
SUITE
libc-complex-unittests
@ -44,11 +38,9 @@ add_libc_test(
conjf16_test.cpp
DEPENDS
libc.src.complex.conjf16
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
conjf128_test
SUITE
libc-complex-unittests
@ -56,11 +48,9 @@ add_libc_test(
conjf128_test.cpp
DEPENDS
libc.src.complex.conjf128
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
cproj_test
SUITE
libc-complex-unittests
@ -68,23 +58,20 @@ add_libc_test(
cproj_test.cpp
DEPENDS
libc.src.complex.cproj
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
cprojf_test
NEED_MPC
SUITE
libc-complex-unittests
SRCS
cprojf_test.cpp
DEPENDS
libc.src.complex.cprojf
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
cprojl_test
SUITE
libc-complex-unittests
@ -92,11 +79,9 @@ add_libc_test(
cprojl_test.cpp
DEPENDS
libc.src.complex.cprojl
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
cprojf16_test
SUITE
libc-complex-unittests
@ -104,11 +89,9 @@ add_libc_test(
cprojf16_test.cpp
DEPENDS
libc.src.complex.cprojf16
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
cprojf128_test
SUITE
libc-complex-unittests
@ -116,11 +99,9 @@ add_libc_test(
cprojf128_test.cpp
DEPENDS
libc.src.complex.cprojf128
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
creal_test
SUITE
libc-complex-unittests
@ -128,11 +109,9 @@ add_libc_test(
creal_test.cpp
DEPENDS
libc.src.complex.creal
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
crealf_test
SUITE
libc-complex-unittests
@ -140,11 +119,9 @@ add_libc_test(
crealf_test.cpp
DEPENDS
libc.src.complex.crealf
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
creall_test
SUITE
libc-complex-unittests
@ -152,11 +129,9 @@ add_libc_test(
creall_test.cpp
DEPENDS
libc.src.complex.creall
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
crealf16_test
SUITE
libc-complex-unittests
@ -164,11 +139,9 @@ add_libc_test(
crealf16_test.cpp
DEPENDS
libc.src.complex.crealf16
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
crealf128_test
SUITE
libc-complex-unittests
@ -176,11 +149,9 @@ add_libc_test(
crealf128_test.cpp
DEPENDS
libc.src.complex.crealf128
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
cimag_test
SUITE
libc-complex-unittests
@ -188,11 +159,9 @@ add_libc_test(
cimag_test.cpp
DEPENDS
libc.src.complex.cimag
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
cimagf_test
SUITE
libc-complex-unittests
@ -200,11 +169,9 @@ add_libc_test(
cimagf_test.cpp
DEPENDS
libc.src.complex.cimagf
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
cimagl_test
SUITE
libc-complex-unittests
@ -212,11 +179,9 @@ add_libc_test(
cimagl_test.cpp
DEPENDS
libc.src.complex.cimagl
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
cimagf16_test
SUITE
libc-complex-unittests
@ -224,11 +189,9 @@ add_libc_test(
cimagf16_test.cpp
DEPENDS
libc.src.complex.cimagf16
LINK_LIBRARIES
LibcFPTestHelpers
)
add_libc_test(
add_fp_unittest(
cimagf128_test
SUITE
libc-complex-unittests
@ -236,6 +199,4 @@ add_libc_test(
cimagf128_test.cpp
DEPENDS
libc.src.complex.cimagf128
LINK_LIBRARIES
LibcFPTestHelpers
)

View File

@ -10,4 +10,16 @@
#include "src/complex/cprojf.h"
#include "utils/MPCWrapper/MPCUtils.h"
using LlvmLibcCprojTestMPC = LIBC_NAMESPACE::testing::FPTest<float>;
namespace mpc = LIBC_NAMESPACE::testing::mpc;
TEST_F(LlvmLibcCprojTestMPC, MPCRND) {
_Complex float test = 5.0 + 10.0i;
EXPECT_MPC_MATCH_ALL_ROUNDING(mpc::Operation::Cproj, test,
LIBC_NAMESPACE::cprojf(test), 0.5);
}
LIST_CPROJ_TESTS(_Complex float, float, LIBC_NAMESPACE::cprojf)

View File

@ -2,4 +2,5 @@ add_subdirectory(hdrgen)
if(LLVM_INCLUDE_TESTS)
add_subdirectory(MPFRWrapper)
add_subdirectory(MPCWrapper)
endif()

View File

@ -0,0 +1,24 @@
if(LIBC_TESTS_CAN_USE_MPC)
add_library(libcMPCWrapper STATIC
MPCUtils.cpp
MPCUtils.h
)
_get_common_test_compile_options(compile_options "" "")
list(REMOVE_ITEM compile_options "-ffreestanding")
target_compile_options(libcMPCWrapper PRIVATE -O3 ${compile_options})
add_dependencies(
libcMPCWrapper
libcMPCommon
libc.src.__support.CPP.array
libc.src.__support.CPP.string
libc.src.__support.CPP.stringstream
libc.src.__support.CPP.type_traits
libc.src.__support.FPUtil.fp_bits
libc.src.__support.complex_type
LibcTest.unit
)
target_include_directories(libcMPCWrapper PUBLIC ${LIBC_SOURCE_DIR})
target_link_libraries(libcMPCWrapper PUBLIC libcMPCommon LibcFPTestHelpers.unit LibcTest.unit mpc)
elseif(NOT LIBC_TARGET_OS_IS_GPU AND NOT LLVM_LIBC_FULL_BUILD)
message(WARNING "Math tests using MPC will be skipped.")
endif()

View File

@ -0,0 +1,344 @@
//===-- 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 "src/__support/CPP/array.h"
#include "src/__support/CPP/stringstream.h"
#include "utils/MPFRWrapper/MPCommon.h"
#include <stdint.h>
#include "mpc.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 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::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

View File

@ -0,0 +1,270 @@
//===-- MPCUtils.h ----------------------------------------------*- C++ -*-===//
//
// 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
//
//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_UTILS_MPCWRAPPER_MPCUTILS_H
#define LLVM_LIBC_UTILS_MPCWRAPPER_MPCUTILS_H
#include "src/__support/CPP/type_traits.h"
#include "src/__support/complex_type.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/properties/complex_types.h"
#include "src/__support/macros/properties/types.h"
#include "test/UnitTest/RoundingModeUtils.h"
#include "test/UnitTest/Test.h"
#include <stdint.h>
namespace LIBC_NAMESPACE_DECL {
namespace testing {
namespace mpc {
enum class Operation {
// Operations which take a single complex floating point number as input
// and produce a single floating point number as output which has the same
// floating point type as the real/imaginary part of the input.
BeginUnaryOperationsSingleOutputDifferentOutputType,
Carg,
Cabs,
EndUnaryOperationsSingleOutputDifferentOutputType,
// Operations which take a single complex floating point number as input
// and produce a single complex floating point number of the same kind
// as output.
BeginUnaryOperationsSingleOutputSameOutputType,
Cproj,
Csqrt,
Clog,
Cexp,
Csinh,
Ccosh,
Ctanh,
Casinh,
Cacosh,
Catanh,
Csin,
Ccos,
Ctan,
Casin,
Cacos,
Catan,
EndUnaryOperationsSingleOutputSameOutputType,
// Operations which take two complex floating point numbers as input
// and produce a single complex floating point number of the same kind
// as output.
BeginBinaryOperationsSingleOutput,
Cpow,
EndBinaryOperationsSingleOutput,
};
using LIBC_NAMESPACE::fputil::testing::RoundingMode;
template <typename T> struct BinaryInput {
static_assert(LIBC_NAMESPACE::cpp::is_complex_v<T>,
"Template parameter of BinaryInput must be a complex floating "
"point type.");
using Type = T;
T x, y;
};
namespace internal {
template <typename InputType, typename OutputType>
bool compare_unary_operation_single_output_same_type(Operation op,
InputType input,
OutputType libc_output,
double ulp_tolerance,
RoundingMode rounding);
template <typename InputType, typename OutputType>
bool compare_unary_operation_single_output_different_type(
Operation op, InputType input, OutputType libc_output, double ulp_tolerance,
RoundingMode rounding);
template <typename InputType, typename OutputType>
bool compare_binary_operation_one_output(Operation op,
const BinaryInput<InputType> &input,
OutputType libc_output,
double ulp_tolerance,
RoundingMode rounding);
template <typename InputType, typename OutputType>
void explain_unary_operation_single_output_same_type_error(
Operation op, InputType input, OutputType match_value, double ulp_tolerance,
RoundingMode rounding);
template <typename InputType, typename OutputType>
void explain_unary_operation_single_output_different_type_error(
Operation op, InputType input, OutputType match_value, double ulp_tolerance,
RoundingMode rounding);
template <typename InputType, typename OutputType>
void explain_binary_operation_one_output_error(
Operation op, const BinaryInput<InputType> &input, OutputType match_value,
double ulp_tolerance, RoundingMode rounding);
template <Operation op, typename InputType, typename OutputType>
class MPCMatcher : public testing::Matcher<OutputType> {
private:
InputType input;
OutputType match_value;
double ulp_tolerance;
RoundingMode rounding;
public:
MPCMatcher(InputType testInput, double ulp_tolerance, RoundingMode rounding)
: input(testInput), ulp_tolerance(ulp_tolerance), rounding(rounding) {}
bool match(OutputType libcResult) {
match_value = libcResult;
return match(input, match_value);
}
void explainError() override { // NOLINT
explain_error(input, match_value);
}
private:
template <typename InType, typename OutType>
bool match(InType in, OutType out) {
if (cpp::is_same_v<InType, OutType>) {
return compare_unary_operation_single_output_same_type(
op, in, out, ulp_tolerance, rounding);
} else {
return compare_unary_operation_single_output_different_type(
op, in, out, ulp_tolerance, rounding);
}
}
template <typename T, typename U>
bool match(const BinaryInput<T> &in, U out) {
return compare_binary_operation_one_output(op, in, out, ulp_tolerance,
rounding);
}
template <typename InType, typename OutType>
void explain_error(InType in, OutType out) {
if (cpp::is_same_v<InType, OutType>) {
explain_unary_operation_single_output_same_type_error(
op, in, out, ulp_tolerance, rounding);
} else {
explain_unary_operation_single_output_different_type_error(
op, in, out, ulp_tolerance, rounding);
}
}
template <typename T, typename U>
void explain_error(const BinaryInput<T> &in, U out) {
explain_binary_operation_one_output_error(op, in, out, ulp_tolerance,
rounding);
}
};
} // namespace internal
// Return true if the input and ouput types for the operation op are valid
// types.
template <Operation op, typename InputType, typename OutputType>
constexpr bool is_valid_operation() {
return (Operation::BeginBinaryOperationsSingleOutput < op &&
op < Operation::EndBinaryOperationsSingleOutput &&
cpp::is_complex_type_same<InputType, OutputType> &&
cpp::is_complex_v<InputType>) ||
(Operation::BeginUnaryOperationsSingleOutputSameOutputType < op &&
op < Operation::EndUnaryOperationsSingleOutputSameOutputType &&
cpp::is_complex_type_same<InputType, OutputType> &&
cpp::is_complex_v<InputType>) ||
(Operation::BeginUnaryOperationsSingleOutputDifferentOutputType < op &&
op < Operation::EndUnaryOperationsSingleOutputDifferentOutputType &&
cpp::is_same_v<make_real_t<InputType>, OutputType> &&
cpp::is_complex_v<InputType>);
}
template <Operation op, typename InputType, typename OutputType>
cpp::enable_if_t<is_valid_operation<op, InputType, OutputType>(),
internal::MPCMatcher<op, InputType, OutputType>>
get_mpc_matcher(InputType input, [[maybe_unused]] OutputType output,
double ulp_tolerance, RoundingMode rounding) {
return internal::MPCMatcher<op, InputType, OutputType>(input, ulp_tolerance,
rounding);
}
} // namespace mpc
} // namespace testing
} // namespace LIBC_NAMESPACE_DECL
#define EXPECT_MPC_MATCH_DEFAULT(op, input, match_value, ulp_tolerance) \
EXPECT_THAT(match_value, \
LIBC_NAMESPACE::testing::mpc::get_mpc_matcher<op>( \
input, match_value, ulp_tolerance, \
LIBC_NAMESPACE::fputil::testing::RoundingMode::Nearest))
#define EXPECT_MPC_MATCH_ROUNDING(op, input, match_value, ulp_tolerance, \
rounding) \
EXPECT_THAT(match_value, LIBC_NAMESPACE::testing::mpc::get_mpc_matcher<op>( \
input, match_value, ulp_tolerance, rounding))
#define EXPECT_MPC_MATCH_ALL_ROUNDING_HELPER(op, input, match_value, \
ulp_tolerance, rounding) \
{ \
MPCRND::ForceRoundingMode __r(rounding); \
if (__r.success) { \
EXPECT_MPC_MATCH_ROUNDING(op, input, match_value, ulp_tolerance, \
rounding); \
} \
}
#define EXPECT_MPC_MATCH_ALL_ROUNDING(op, input, match_value, ulp_tolerance) \
{ \
namespace MPCRND = LIBC_NAMESPACE::fputil::testing; \
for (int i = 0; i < 4; i++) { \
MPCRND::RoundingMode r_mode = static_cast<MPCRND::RoundingMode>(i); \
EXPECT_MPC_MATCH_ALL_ROUNDING_HELPER(op, input, match_value, \
ulp_tolerance, r_mode); \
} \
}
#define TEST_MPC_MATCH_ROUNDING(op, input, match_value, ulp_tolerance, \
rounding) \
LIBC_NAMESPACE::testing::mpc::get_mpc_matcher<op>(input, match_value, \
ulp_tolerance, rounding) \
.match(match_value)
#define ASSERT_MPC_MATCH_DEFAULT(op, input, match_value, ulp_tolerance) \
ASSERT_THAT(match_value, \
LIBC_NAMESPACE::testing::mpc::get_mpc_matcher<op>( \
input, match_value, ulp_tolerance, \
LIBC_NAMESPACE::fputil::testing::RoundingMode::Nearest))
#define ASSERT_MPC_MATCH_ROUNDING(op, input, match_value, ulp_tolerance, \
rounding) \
ASSERT_THAT(match_value, LIBC_NAMESPACE::testing::mpc::get_mpc_matcher<op>( \
input, match_value, ulp_tolerance, rounding))
#define ASSERT_MPC_MATCH_ALL_ROUNDING_HELPER(op, input, match_value, \
ulp_tolerance, rounding) \
{ \
MPCRND::ForceRoundingMode __r(rounding); \
if (__r.success) { \
ASSERT_MPC_MATCH_ROUNDING(op, input, match_value, ulp_tolerance, \
rounding); \
} \
}
#define ASSERT_MPC_MATCH_ALL_ROUNDING(op, input, match_value, ulp_tolerance) \
{ \
namespace MPCRND = LIBC_NAMESPACE::fputil::testing; \
for (int i = 0; i < 4; i++) { \
MPCRND::RoundingMode r_mode = static_cast<MPCRND::RoundingMode>(i); \
ASSERT_MPC_MATCH_ALL_ROUNDING_HELPER(op, input, match_value, \
ulp_tolerance, r_mode); \
} \
}
#endif // LLVM_LIBC_UTILS_MPCWRAPPER_MPCUTILS_H

View File

@ -0,0 +1,8 @@
#include <mpc.h>
int main() {
mpc_t x;
mpc_init2(x, 256);
mpc_clear(x);
return 0;
}

View File

@ -1,21 +1,40 @@
if(LIBC_TESTS_CAN_USE_MPFR)
add_library(libcMPFRWrapper STATIC
MPFRUtils.cpp
MPFRUtils.h
if(LIBC_TESTS_CAN_USE_MPFR OR LIBC_TESTS_CAN_USE_MPC)
add_library(libcMPCommon STATIC
MPCommon.cpp
MPCommon.h
mpfr_inc.h
)
_get_common_test_compile_options(compile_options "" "")
# mpfr/gmp headers do not work with -ffreestanding flag.
list(REMOVE_ITEM compile_options "-ffreestanding")
target_compile_options(libcMPFRWrapper PRIVATE -O3 ${compile_options})
target_compile_options(libcMPCommon PRIVATE -O3 ${compile_options})
add_dependencies(
libcMPFRWrapper
libc.src.__support.CPP.array
libc.src.__support.CPP.stringstream
libcMPCommon
libc.src.__support.CPP.string
libc.src.__support.CPP.string_view
libc.src.__support.CPP.type_traits
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fp_bits
)
target_include_directories(libcMPCommon PUBLIC ${LIBC_SOURCE_DIR})
target_link_libraries(libcMPCommon PUBLIC LibcFPTestHelpers.unit mpfr gmp)
elseif(NOT LIBC_TARGET_OS_IS_GPU AND NOT LLVM_LIBC_FULL_BUILD)
message(WARNING "Math tests using MPFR will be skipped.")
endif()
if(LIBC_TESTS_CAN_USE_MPFR)
add_library(libcMPFRWrapper STATIC
MPFRUtils.cpp
MPFRUtils.h
)
_get_common_test_compile_options(compile_options "" "")
target_compile_options(libcMPFRWrapper PRIVATE -O3 ${compile_options})
add_dependencies(
libcMPFRWrapper
libcMPCommon
libc.src.__support.CPP.array
libc.src.__support.CPP.stringstream
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.fpbits_str
LibcTest.unit
)
@ -24,7 +43,7 @@ if(LIBC_TESTS_CAN_USE_MPFR)
target_link_directories(libcMPFRWrapper PUBLIC ${LLVM_LIBC_MPFR_INSTALL_PATH}/lib)
endif()
target_include_directories(libcMPFRWrapper PUBLIC ${LIBC_SOURCE_DIR})
target_link_libraries(libcMPFRWrapper PUBLIC LibcFPTestHelpers.unit LibcTest.unit mpfr gmp)
target_link_libraries(libcMPFRWrapper PUBLIC libcMPCommon LibcFPTestHelpers.unit LibcTest.unit)
elseif(NOT LIBC_TARGET_OS_IS_GPU AND NOT LLVM_LIBC_FULL_BUILD)
message(WARNING "Math tests using MPFR will be skipped.")
endif()

View File

@ -0,0 +1,561 @@
//===-- Utils used by both MPCWrapper and MPFRWrapper----------------------===//
//
// 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 "MPCommon.h"
#include "src/__support/CPP/string_view.h"
#include "src/__support/FPUtil/cast.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/properties/types.h"
namespace LIBC_NAMESPACE_DECL {
namespace testing {
namespace mpfr {
MPFRNumber::MPFRNumber() : mpfr_precision(256), mpfr_rounding(MPFR_RNDN) {
mpfr_init2(value, mpfr_precision);
}
MPFRNumber::MPFRNumber(const MPFRNumber &other)
: mpfr_precision(other.mpfr_precision), mpfr_rounding(other.mpfr_rounding) {
mpfr_init2(value, mpfr_precision);
mpfr_set(value, other.value, mpfr_rounding);
}
MPFRNumber::MPFRNumber(const MPFRNumber &other, unsigned int precision)
: mpfr_precision(precision), mpfr_rounding(other.mpfr_rounding) {
mpfr_init2(value, mpfr_precision);
mpfr_set(value, other.value, mpfr_rounding);
}
MPFRNumber::MPFRNumber(const mpfr_t x, unsigned int precision,
RoundingMode rounding)
: mpfr_precision(precision),
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
mpfr_init2(value, mpfr_precision);
mpfr_set(value, x, mpfr_rounding);
}
MPFRNumber::~MPFRNumber() { mpfr_clear(value); }
MPFRNumber &MPFRNumber::operator=(const MPFRNumber &rhs) {
mpfr_precision = rhs.mpfr_precision;
mpfr_rounding = rhs.mpfr_rounding;
mpfr_set(value, rhs.value, mpfr_rounding);
return *this;
}
bool MPFRNumber::is_nan() const { return mpfr_nan_p(value); }
MPFRNumber MPFRNumber::abs() const {
MPFRNumber result(*this);
mpfr_abs(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::acos() const {
MPFRNumber result(*this);
mpfr_acos(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::acosh() const {
MPFRNumber result(*this);
mpfr_acosh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::add(const MPFRNumber &b) const {
MPFRNumber result(*this);
mpfr_add(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::asin() const {
MPFRNumber result(*this);
mpfr_asin(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::asinh() const {
MPFRNumber result(*this);
mpfr_asinh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::atan() const {
MPFRNumber result(*this);
mpfr_atan(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::atan2(const MPFRNumber &b) {
MPFRNumber result(*this);
mpfr_atan2(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::atanh() const {
MPFRNumber result(*this);
mpfr_atanh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::cbrt() const {
MPFRNumber result(*this);
mpfr_cbrt(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::ceil() const {
MPFRNumber result(*this);
mpfr_ceil(result.value, value);
return result;
}
MPFRNumber MPFRNumber::cos() const {
MPFRNumber result(*this);
mpfr_cos(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::cosh() const {
MPFRNumber result(*this);
mpfr_cosh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::cospi() const {
MPFRNumber result(*this);
#if MPFR_VERSION_MAJOR > 4 || \
(MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2)
mpfr_cospi(result.value, value, mpfr_rounding);
return result;
#else
if (mpfr_integer_p(value)) {
mpz_t integer;
mpz_init(integer);
mpfr_get_z(integer, value, mpfr_rounding);
int d = mpz_tstbit(integer, 0);
mpfr_set_si(result.value, d ? -1 : 1, mpfr_rounding);
mpz_clear(integer);
return result;
}
MPFRNumber value_pi(0.0, 1280);
mpfr_const_pi(value_pi.value, MPFR_RNDN);
mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN);
mpfr_cos(result.value, value_pi.value, mpfr_rounding);
return result;
#endif
}
MPFRNumber MPFRNumber::erf() const {
MPFRNumber result(*this);
mpfr_erf(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::exp() const {
MPFRNumber result(*this);
mpfr_exp(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::exp2() const {
MPFRNumber result(*this);
mpfr_exp2(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::exp2m1() const {
// TODO: Only use mpfr_exp2m1 once CI and buildbots get MPFR >= 4.2.0.
#if MPFR_VERSION_MAJOR > 4 || \
(MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2)
MPFRNumber result(*this);
mpfr_exp2m1(result.value, value, mpfr_rounding);
return result;
#else
unsigned int prec = mpfr_precision * 3;
MPFRNumber result(*this, prec);
float f = mpfr_get_flt(abs().value, mpfr_rounding);
if (f > 0.5f && f < 0x1.0p30f) {
mpfr_exp2(result.value, value, mpfr_rounding);
mpfr_sub_ui(result.value, result.value, 1, mpfr_rounding);
return result;
}
MPFRNumber ln2(2.0f, prec);
// log(2)
mpfr_log(ln2.value, ln2.value, mpfr_rounding);
// x * log(2)
mpfr_mul(result.value, value, ln2.value, mpfr_rounding);
// e^(x * log(2)) - 1
int ex = mpfr_expm1(result.value, result.value, mpfr_rounding);
mpfr_subnormalize(result.value, ex, mpfr_rounding);
return result;
#endif
}
MPFRNumber MPFRNumber::exp10() const {
MPFRNumber result(*this);
mpfr_exp10(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::exp10m1() const {
// TODO: Only use mpfr_exp10m1 once CI and buildbots get MPFR >= 4.2.0.
#if MPFR_VERSION_MAJOR > 4 || \
(MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2)
MPFRNumber result(*this);
mpfr_exp10m1(result.value, value, mpfr_rounding);
return result;
#else
unsigned int prec = mpfr_precision * 3;
MPFRNumber result(*this, prec);
MPFRNumber ln10(10.0f, prec);
// log(10)
mpfr_log(ln10.value, ln10.value, mpfr_rounding);
// x * log(10)
mpfr_mul(result.value, value, ln10.value, mpfr_rounding);
// e^(x * log(10)) - 1
int ex = mpfr_expm1(result.value, result.value, mpfr_rounding);
mpfr_subnormalize(result.value, ex, mpfr_rounding);
return result;
#endif
}
MPFRNumber MPFRNumber::expm1() const {
MPFRNumber result(*this);
mpfr_expm1(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::div(const MPFRNumber &b) const {
MPFRNumber result(*this);
mpfr_div(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::floor() const {
MPFRNumber result(*this);
mpfr_floor(result.value, value);
return result;
}
MPFRNumber MPFRNumber::fmod(const MPFRNumber &b) {
MPFRNumber result(*this);
mpfr_fmod(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::frexp(int &exp) {
MPFRNumber result(*this);
mpfr_exp_t resultExp;
mpfr_frexp(&resultExp, result.value, value, mpfr_rounding);
exp = resultExp;
return result;
}
MPFRNumber MPFRNumber::hypot(const MPFRNumber &b) {
MPFRNumber result(*this);
mpfr_hypot(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::log() const {
MPFRNumber result(*this);
mpfr_log(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::log2() const {
MPFRNumber result(*this);
mpfr_log2(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::log10() const {
MPFRNumber result(*this);
mpfr_log10(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::log1p() const {
MPFRNumber result(*this);
mpfr_log1p(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::pow(const MPFRNumber &b) {
MPFRNumber result(*this);
mpfr_pow(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::remquo(const MPFRNumber &divisor, int &quotient) {
MPFRNumber remainder(*this);
long q;
mpfr_remquo(remainder.value, &q, value, divisor.value, mpfr_rounding);
quotient = q;
return remainder;
}
MPFRNumber MPFRNumber::round() const {
MPFRNumber result(*this);
mpfr_round(result.value, value);
return result;
}
MPFRNumber MPFRNumber::roundeven() const {
MPFRNumber result(*this);
#if MPFR_VERSION_MAJOR >= 4
mpfr_roundeven(result.value, value);
#else
mpfr_rint(result.value, value, MPFR_RNDN);
#endif
return result;
}
bool MPFRNumber::round_to_long(long &result) const {
// We first calculate the rounded value. This way, when converting
// to long using mpfr_get_si, the rounding direction of MPFR_RNDN
// (or any other rounding mode), does not have an influence.
MPFRNumber roundedValue = round();
mpfr_clear_erangeflag();
result = mpfr_get_si(roundedValue.value, MPFR_RNDN);
return mpfr_erangeflag_p();
}
bool MPFRNumber::round_to_long(mpfr_rnd_t rnd, long &result) const {
MPFRNumber rint_result(*this);
mpfr_rint(rint_result.value, value, rnd);
return rint_result.round_to_long(result);
}
MPFRNumber MPFRNumber::rint(mpfr_rnd_t rnd) const {
MPFRNumber result(*this);
mpfr_rint(result.value, value, rnd);
return result;
}
MPFRNumber MPFRNumber::mod_2pi() const {
MPFRNumber result(0.0, 1280);
MPFRNumber _2pi(0.0, 1280);
mpfr_const_pi(_2pi.value, MPFR_RNDN);
mpfr_mul_si(_2pi.value, _2pi.value, 2, MPFR_RNDN);
mpfr_fmod(result.value, value, _2pi.value, MPFR_RNDN);
return result;
}
MPFRNumber MPFRNumber::mod_pi_over_2() const {
MPFRNumber result(0.0, 1280);
MPFRNumber pi_over_2(0.0, 1280);
mpfr_const_pi(pi_over_2.value, MPFR_RNDN);
mpfr_mul_d(pi_over_2.value, pi_over_2.value, 0.5, MPFR_RNDN);
mpfr_fmod(result.value, value, pi_over_2.value, MPFR_RNDN);
return result;
}
MPFRNumber MPFRNumber::mod_pi_over_4() const {
MPFRNumber result(0.0, 1280);
MPFRNumber pi_over_4(0.0, 1280);
mpfr_const_pi(pi_over_4.value, MPFR_RNDN);
mpfr_mul_d(pi_over_4.value, pi_over_4.value, 0.25, MPFR_RNDN);
mpfr_fmod(result.value, value, pi_over_4.value, MPFR_RNDN);
return result;
}
MPFRNumber MPFRNumber::sin() const {
MPFRNumber result(*this);
mpfr_sin(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::sinpi() const {
MPFRNumber result(*this);
#if MPFR_VERSION_MAJOR > 4 || \
(MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2)
mpfr_sinpi(result.value, value, mpfr_rounding);
return result;
#else
if (mpfr_integer_p(value)) {
mpfr_set_si(result.value, 0, mpfr_rounding);
return result;
}
MPFRNumber value_mul_two(*this);
mpfr_mul_si(value_mul_two.value, value, 2, MPFR_RNDN);
if (mpfr_integer_p(value_mul_two.value)) {
auto d = mpfr_get_si(value, MPFR_RNDD);
mpfr_set_si(result.value, (d & 1) ? -1 : 1, mpfr_rounding);
return result;
}
MPFRNumber value_pi(0.0, 1280);
mpfr_const_pi(value_pi.value, MPFR_RNDN);
mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN);
mpfr_sin(result.value, value_pi.value, mpfr_rounding);
return result;
#endif
}
MPFRNumber MPFRNumber::sinh() const {
MPFRNumber result(*this);
mpfr_sinh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::sqrt() const {
MPFRNumber result(*this);
mpfr_sqrt(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::sub(const MPFRNumber &b) const {
MPFRNumber result(*this);
mpfr_sub(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::tan() const {
MPFRNumber result(*this);
mpfr_tan(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::tanh() const {
MPFRNumber result(*this);
mpfr_tanh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::tanpi() const {
MPFRNumber result(*this);
#if MPFR_VERSION_MAJOR > 4 || \
(MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2)
mpfr_tanpi(result.value, value, mpfr_rounding);
return result;
#else
MPFRNumber value_ret_exact(*this);
MPFRNumber value_one(*this);
mpfr_set_si(value_one.value, 1, MPFR_RNDN);
mpfr_fmod(value_ret_exact.value, value, value_one.value, mpfr_rounding);
mpfr_mul_si(value_ret_exact.value, value_ret_exact.value, 4, MPFR_RNDN);
if (mpfr_integer_p(value_ret_exact.value)) {
int mod = mpfr_get_si(value_ret_exact.value, MPFR_RNDN);
mod = (mod < 0 ? -1 * mod : mod);
switch (mod) {
case 0:
mpfr_set_si(result.value, 0, mpfr_rounding);
break;
case 1:
mpfr_set_si(result.value, (mpfr_signbit(value) ? -1 : 1), mpfr_rounding);
break;
case 2: {
auto d = mpfr_get_si(value, MPFR_RNDZ);
d += mpfr_sgn(value) > 0 ? 0 : 1;
mpfr_set_inf(result.value, (d & 1) ? -1 : 1);
break;
}
case 3:
mpfr_set_si(result.value, (mpfr_signbit(value) ? 1 : -1), mpfr_rounding);
break;
}
return result;
}
MPFRNumber value_pi(0.0, 1280);
mpfr_const_pi(value_pi.value, MPFR_RNDN);
mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN);
mpfr_tan(result.value, value_pi.value, mpfr_rounding);
return result;
#endif
}
MPFRNumber MPFRNumber::trunc() const {
MPFRNumber result(*this);
mpfr_trunc(result.value, value);
return result;
}
MPFRNumber MPFRNumber::fma(const MPFRNumber &b, const MPFRNumber &c) {
MPFRNumber result(*this);
mpfr_fma(result.value, value, b.value, c.value, mpfr_rounding);
return result;
}
MPFRNumber MPFRNumber::mul(const MPFRNumber &b) {
MPFRNumber result(*this);
mpfr_mul(result.value, value, b.value, mpfr_rounding);
return result;
}
cpp::string MPFRNumber::str() const {
// 200 bytes should be more than sufficient to hold a 100-digit number
// plus additional bytes for the decimal point, '-' sign etc.
constexpr size_t printBufSize = 200;
char buffer[printBufSize];
mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
cpp::string_view view(buffer);
// Trim whitespaces
const char whitespace = ' ';
while (!view.empty() && view.front() == whitespace)
view.remove_prefix(1);
while (!view.empty() && view.back() == whitespace)
view.remove_suffix(1);
return cpp::string(view.data());
}
void MPFRNumber::dump(const char *msg) const {
mpfr_printf("%s%.128g\n", msg, value);
}
template <> float MPFRNumber::as<float>() const {
return mpfr_get_flt(value, mpfr_rounding);
}
template <> double MPFRNumber::as<double>() const {
return mpfr_get_d(value, mpfr_rounding);
}
template <> long double MPFRNumber::as<long double>() const {
return mpfr_get_ld(value, mpfr_rounding);
}
#ifdef LIBC_TYPES_HAS_FLOAT16
template <> float16 MPFRNumber::as<float16>() const {
// TODO: Either prove that this cast won't cause double-rounding errors, or
// find a better way to get a float16.
return fputil::cast<float16>(mpfr_get_d(value, mpfr_rounding));
}
#endif
#ifdef LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
template <> float128 MPFRNumber::as<float128>() const {
return mpfr_get_float128(value, mpfr_rounding);
}
#endif // LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
} // namespace mpfr
} // namespace testing
} // namespace LIBC_NAMESPACE_DECL

View File

@ -0,0 +1,342 @@
//===-- MPCommon.h ----------------------------------------------*- C++ -*-===//
//
// 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
//
//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_UTILS_MPFRWRAPPER_MPCOMMON_H
#define LLVM_LIBC_UTILS_MPFRWRAPPER_MPCOMMON_H
#include "src/__support/CPP/string.h"
#include "src/__support/CPP/type_traits.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/macros/config.h"
#include "test/UnitTest/RoundingModeUtils.h"
#include <stdint.h>
#include "mpfr_inc.h"
#ifdef LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
extern "C" {
int mpfr_set_float128(mpfr_ptr, float128, mpfr_rnd_t);
float128 mpfr_get_float128(mpfr_srcptr, mpfr_rnd_t);
}
#endif
namespace LIBC_NAMESPACE_DECL {
namespace testing {
namespace mpfr {
template <typename T> using FPBits = LIBC_NAMESPACE::fputil::FPBits<T>;
using LIBC_NAMESPACE::fputil::testing::RoundingMode;
// A precision value which allows sufficiently large additional
// precision compared to the floating point precision.
template <typename T> struct ExtraPrecision;
#ifdef LIBC_TYPES_HAS_FLOAT16
template <> struct ExtraPrecision<float16> {
static constexpr unsigned int VALUE = 128;
};
#endif
template <> struct ExtraPrecision<float> {
static constexpr unsigned int VALUE = 128;
};
template <> struct ExtraPrecision<double> {
static constexpr unsigned int VALUE = 256;
};
template <> struct ExtraPrecision<long double> {
#ifdef LIBC_TYPES_LONG_DOUBLE_IS_FLOAT128
static constexpr unsigned int VALUE = 512;
#else
static constexpr unsigned int VALUE = 256;
#endif
};
#if defined(LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE)
template <> struct ExtraPrecision<float128> {
static constexpr unsigned int VALUE = 512;
};
#endif // LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
// If the ulp tolerance is less than or equal to 0.5, we would check that the
// result is rounded correctly with respect to the rounding mode by using the
// same precision as the inputs.
template <typename T>
static inline unsigned int get_precision(double ulp_tolerance) {
if (ulp_tolerance <= 0.5) {
return LIBC_NAMESPACE::fputil::FPBits<T>::FRACTION_LEN + 1;
} else {
return ExtraPrecision<T>::VALUE;
}
}
static inline mpfr_rnd_t get_mpfr_rounding_mode(RoundingMode mode) {
switch (mode) {
case RoundingMode::Upward:
return MPFR_RNDU;
break;
case RoundingMode::Downward:
return MPFR_RNDD;
break;
case RoundingMode::TowardZero:
return MPFR_RNDZ;
break;
case RoundingMode::Nearest:
return MPFR_RNDN;
break;
}
__builtin_unreachable();
}
class MPFRNumber {
unsigned int mpfr_precision;
mpfr_rnd_t mpfr_rounding;
mpfr_t value;
public:
MPFRNumber();
// We use explicit EnableIf specializations to disallow implicit
// conversions. Implicit conversions can potentially lead to loss of
// precision. We exceptionally allow implicit conversions from float16
// to float, as the MPFR API does not support float16, thus requiring
// conversion to a higher-precision format.
template <typename XType,
cpp::enable_if_t<cpp::is_same_v<float, XType>
#ifdef LIBC_TYPES_HAS_FLOAT16
|| cpp::is_same_v<float16, XType>
#endif
,
int> = 0>
explicit MPFRNumber(XType x,
unsigned int precision = ExtraPrecision<XType>::VALUE,
RoundingMode rounding = RoundingMode::Nearest)
: mpfr_precision(precision),
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
mpfr_init2(value, mpfr_precision);
mpfr_set_flt(value, x, mpfr_rounding);
}
template <typename XType,
cpp::enable_if_t<cpp::is_same_v<double, XType>, int> = 0>
explicit MPFRNumber(XType x,
unsigned int precision = ExtraPrecision<XType>::VALUE,
RoundingMode rounding = RoundingMode::Nearest)
: mpfr_precision(precision),
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
mpfr_init2(value, mpfr_precision);
mpfr_set_d(value, x, mpfr_rounding);
}
template <typename XType,
cpp::enable_if_t<cpp::is_same_v<long double, XType>, int> = 0>
explicit MPFRNumber(XType x,
unsigned int precision = ExtraPrecision<XType>::VALUE,
RoundingMode rounding = RoundingMode::Nearest)
: mpfr_precision(precision),
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
mpfr_init2(value, mpfr_precision);
mpfr_set_ld(value, x, mpfr_rounding);
}
#ifdef LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
template <typename XType,
cpp::enable_if_t<cpp::is_same_v<float128, XType>, int> = 0>
explicit MPFRNumber(XType x,
unsigned int precision = ExtraPrecision<XType>::VALUE,
RoundingMode rounding = RoundingMode::Nearest)
: mpfr_precision(precision),
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
mpfr_init2(value, mpfr_precision);
mpfr_set_float128(value, x, mpfr_rounding);
}
#endif // LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
template <typename XType,
cpp::enable_if_t<cpp::is_integral_v<XType>, int> = 0>
explicit MPFRNumber(XType x,
unsigned int precision = ExtraPrecision<float>::VALUE,
RoundingMode rounding = RoundingMode::Nearest)
: mpfr_precision(precision),
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
mpfr_init2(value, mpfr_precision);
mpfr_set_sj(value, x, mpfr_rounding);
}
MPFRNumber(const MPFRNumber &other);
MPFRNumber(const MPFRNumber &other, unsigned int precision);
MPFRNumber(const mpfr_t x, unsigned int precision, RoundingMode rounding);
~MPFRNumber();
MPFRNumber &operator=(const MPFRNumber &rhs);
bool is_nan() const;
MPFRNumber abs() const;
MPFRNumber acos() const;
MPFRNumber acosh() const;
MPFRNumber add(const MPFRNumber &b) const;
MPFRNumber asin() const;
MPFRNumber asinh() const;
MPFRNumber atan() const;
MPFRNumber atan2(const MPFRNumber &b);
MPFRNumber atanh() const;
MPFRNumber cbrt() const;
MPFRNumber ceil() const;
MPFRNumber cos() const;
MPFRNumber cosh() const;
MPFRNumber cospi() const;
MPFRNumber erf() const;
MPFRNumber exp() const;
MPFRNumber exp2() const;
MPFRNumber exp2m1() const;
MPFRNumber exp10() const;
MPFRNumber exp10m1() const;
MPFRNumber expm1() const;
MPFRNumber div(const MPFRNumber &b) const;
MPFRNumber floor() const;
MPFRNumber fmod(const MPFRNumber &b);
MPFRNumber frexp(int &exp);
MPFRNumber hypot(const MPFRNumber &b);
MPFRNumber log() const;
MPFRNumber log2() const;
MPFRNumber log10() const;
MPFRNumber log1p() const;
MPFRNumber pow(const MPFRNumber &b);
MPFRNumber remquo(const MPFRNumber &divisor, int &quotient);
MPFRNumber round() const;
MPFRNumber roundeven() const;
bool round_to_long(long &result) const;
bool round_to_long(mpfr_rnd_t rnd, long &result) const;
MPFRNumber rint(mpfr_rnd_t rnd) const;
MPFRNumber mod_2pi() const;
MPFRNumber mod_pi_over_2() const;
MPFRNumber mod_pi_over_4() const;
MPFRNumber sin() const;
MPFRNumber sinpi() const;
MPFRNumber sinh() const;
MPFRNumber sqrt() const;
MPFRNumber sub(const MPFRNumber &b) const;
MPFRNumber tan() const;
MPFRNumber tanh() const;
MPFRNumber tanpi() const;
MPFRNumber trunc() const;
MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c);
MPFRNumber mul(const MPFRNumber &b);
cpp::string str() const;
template <typename T> T as() const;
void dump(const char *msg) const;
// Return the ULP (units-in-the-last-place) difference between the
// stored MPFR and a floating point number.
//
// We define ULP difference as follows:
// If exponents of this value and the |input| are same, then:
// ULP(this_value, input) = abs(this_value - input) / eps(input)
// else:
// max = max(abs(this_value), abs(input))
// min = min(abs(this_value), abs(input))
// maxExponent = exponent(max)
// ULP(this_value, input) = (max - 2^maxExponent) / eps(max) +
// (2^maxExponent - min) / eps(min)
//
// Remarks:
// 1. A ULP of 0.0 will imply that the value is correctly rounded.
// 2. We expect that this value and the value to be compared (the [input]
// argument) are reasonable close, and we will provide an upper bound
// of ULP value for testing. Morever, most of the fractional parts of
// ULP value do not matter much, so using double as the return type
// should be good enough.
// 3. For close enough values (values which don't diff in their exponent by
// not more than 1), a ULP difference of N indicates a bit distance
// of N between this number and [input].
// 4. A values of +0.0 and -0.0 are treated as equal.
template <typename T>
cpp::enable_if_t<cpp::is_floating_point_v<T>, MPFRNumber>
ulp_as_mpfr_number(T input) {
T thisAsT = as<T>();
if (thisAsT == input)
return MPFRNumber(0.0);
if (is_nan()) {
if (FPBits<T>(input).is_nan())
return MPFRNumber(0.0);
return MPFRNumber(FPBits<T>::inf().get_val());
}
int thisExponent = FPBits<T>(thisAsT).get_exponent();
int inputExponent = FPBits<T>(input).get_exponent();
// Adjust the exponents for denormal numbers.
if (FPBits<T>(thisAsT).is_subnormal())
++thisExponent;
if (FPBits<T>(input).is_subnormal())
++inputExponent;
if (thisAsT * input < 0 || thisExponent == inputExponent) {
MPFRNumber inputMPFR(input);
mpfr_sub(inputMPFR.value, value, inputMPFR.value, MPFR_RNDN);
mpfr_abs(inputMPFR.value, inputMPFR.value, MPFR_RNDN);
mpfr_mul_2si(inputMPFR.value, inputMPFR.value,
-thisExponent + FPBits<T>::FRACTION_LEN, MPFR_RNDN);
return inputMPFR;
}
// If the control reaches here, it means that this number and input are
// of the same sign but different exponent. In such a case, ULP error is
// calculated as sum of two parts.
thisAsT = FPBits<T>(thisAsT).abs().get_val();
input = FPBits<T>(input).abs().get_val();
T min = thisAsT > input ? input : thisAsT;
T max = thisAsT > input ? thisAsT : input;
int minExponent = FPBits<T>(min).get_exponent();
int maxExponent = FPBits<T>(max).get_exponent();
// Adjust the exponents for denormal numbers.
if (FPBits<T>(min).is_subnormal())
++minExponent;
if (FPBits<T>(max).is_subnormal())
++maxExponent;
MPFRNumber minMPFR(min);
MPFRNumber maxMPFR(max);
MPFRNumber pivot(uint32_t(1));
mpfr_mul_2si(pivot.value, pivot.value, maxExponent, MPFR_RNDN);
mpfr_sub(minMPFR.value, pivot.value, minMPFR.value, MPFR_RNDN);
mpfr_mul_2si(minMPFR.value, minMPFR.value,
-minExponent + FPBits<T>::FRACTION_LEN, MPFR_RNDN);
mpfr_sub(maxMPFR.value, maxMPFR.value, pivot.value, MPFR_RNDN);
mpfr_mul_2si(maxMPFR.value, maxMPFR.value,
-maxExponent + FPBits<T>::FRACTION_LEN, MPFR_RNDN);
mpfr_add(minMPFR.value, minMPFR.value, maxMPFR.value, MPFR_RNDN);
return minMPFR;
}
template <typename T>
cpp::enable_if_t<cpp::is_floating_point_v<T>, cpp::string>
ulp_as_string(T input) {
MPFRNumber num = ulp_as_mpfr_number(input);
return num.str();
}
template <typename T>
cpp::enable_if_t<cpp::is_floating_point_v<T>, double> ulp(T input) {
MPFRNumber num = ulp_as_mpfr_number(input);
return num.as<double>();
}
};
} // namespace mpfr
} // namespace testing
} // namespace LIBC_NAMESPACE_DECL
#endif // LLVM_LIBC_UTILS_MPFRWRAPPER_MPCOMMON_H

View File

@ -7,806 +7,17 @@
//===----------------------------------------------------------------------===//
#include "MPFRUtils.h"
#include "MPCommon.h"
#include "src/__support/CPP/array.h"
#include "src/__support/CPP/string.h"
#include "src/__support/CPP/string_view.h"
#include "src/__support/CPP/stringstream.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/cast.h"
#include "src/__support/FPUtil/fpbits_str.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/properties/types.h"
#include <stdint.h>
#include "mpfr_inc.h"
#ifdef LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
extern "C" {
int mpfr_set_float128(mpfr_ptr, float128, mpfr_rnd_t);
float128 mpfr_get_float128(mpfr_srcptr, mpfr_rnd_t);
}
#endif
template <typename T> using FPBits = LIBC_NAMESPACE::fputil::FPBits<T>;
namespace LIBC_NAMESPACE_DECL {
namespace testing {
namespace mpfr {
// A precision value which allows sufficiently large additional
// precision compared to the floating point precision.
template <typename T> struct ExtraPrecision;
#ifdef LIBC_TYPES_HAS_FLOAT16
template <> struct ExtraPrecision<float16> {
static constexpr unsigned int VALUE = 128;
};
#endif
template <> struct ExtraPrecision<float> {
static constexpr unsigned int VALUE = 128;
};
template <> struct ExtraPrecision<double> {
static constexpr unsigned int VALUE = 256;
};
template <> struct ExtraPrecision<long double> {
#ifdef LIBC_TYPES_LONG_DOUBLE_IS_FLOAT128
static constexpr unsigned int VALUE = 512;
#else
static constexpr unsigned int VALUE = 256;
#endif
};
#if defined(LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE)
template <> struct ExtraPrecision<float128> {
static constexpr unsigned int VALUE = 512;
};
#endif // LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
// If the ulp tolerance is less than or equal to 0.5, we would check that the
// result is rounded correctly with respect to the rounding mode by using the
// same precision as the inputs.
template <typename T>
static inline unsigned int get_precision(double ulp_tolerance) {
if (ulp_tolerance <= 0.5) {
return LIBC_NAMESPACE::fputil::FPBits<T>::FRACTION_LEN + 1;
} else {
return ExtraPrecision<T>::VALUE;
}
}
static inline mpfr_rnd_t get_mpfr_rounding_mode(RoundingMode mode) {
switch (mode) {
case RoundingMode::Upward:
return MPFR_RNDU;
break;
case RoundingMode::Downward:
return MPFR_RNDD;
break;
case RoundingMode::TowardZero:
return MPFR_RNDZ;
break;
case RoundingMode::Nearest:
return MPFR_RNDN;
break;
}
__builtin_unreachable();
}
class MPFRNumber {
unsigned int mpfr_precision;
mpfr_rnd_t mpfr_rounding;
mpfr_t value;
public:
MPFRNumber() : mpfr_precision(256), mpfr_rounding(MPFR_RNDN) {
mpfr_init2(value, mpfr_precision);
}
// We use explicit EnableIf specializations to disallow implicit
// conversions. Implicit conversions can potentially lead to loss of
// precision. We exceptionally allow implicit conversions from float16
// to float, as the MPFR API does not support float16, thus requiring
// conversion to a higher-precision format.
template <typename XType,
cpp::enable_if_t<cpp::is_same_v<float, XType>
#ifdef LIBC_TYPES_HAS_FLOAT16
|| cpp::is_same_v<float16, XType>
#endif
,
int> = 0>
explicit MPFRNumber(XType x,
unsigned int precision = ExtraPrecision<XType>::VALUE,
RoundingMode rounding = RoundingMode::Nearest)
: mpfr_precision(precision),
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
mpfr_init2(value, mpfr_precision);
mpfr_set_flt(value, x, mpfr_rounding);
}
template <typename XType,
cpp::enable_if_t<cpp::is_same_v<double, XType>, int> = 0>
explicit MPFRNumber(XType x,
unsigned int precision = ExtraPrecision<XType>::VALUE,
RoundingMode rounding = RoundingMode::Nearest)
: mpfr_precision(precision),
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
mpfr_init2(value, mpfr_precision);
mpfr_set_d(value, x, mpfr_rounding);
}
template <typename XType,
cpp::enable_if_t<cpp::is_same_v<long double, XType>, int> = 0>
explicit MPFRNumber(XType x,
unsigned int precision = ExtraPrecision<XType>::VALUE,
RoundingMode rounding = RoundingMode::Nearest)
: mpfr_precision(precision),
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
mpfr_init2(value, mpfr_precision);
mpfr_set_ld(value, x, mpfr_rounding);
}
#ifdef LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
template <typename XType,
cpp::enable_if_t<cpp::is_same_v<float128, XType>, int> = 0>
explicit MPFRNumber(XType x,
unsigned int precision = ExtraPrecision<XType>::VALUE,
RoundingMode rounding = RoundingMode::Nearest)
: mpfr_precision(precision),
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
mpfr_init2(value, mpfr_precision);
mpfr_set_float128(value, x, mpfr_rounding);
}
#endif // LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
template <typename XType,
cpp::enable_if_t<cpp::is_integral_v<XType>, int> = 0>
explicit MPFRNumber(XType x,
unsigned int precision = ExtraPrecision<float>::VALUE,
RoundingMode rounding = RoundingMode::Nearest)
: mpfr_precision(precision),
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
mpfr_init2(value, mpfr_precision);
mpfr_set_sj(value, x, mpfr_rounding);
}
MPFRNumber(const MPFRNumber &other)
: mpfr_precision(other.mpfr_precision),
mpfr_rounding(other.mpfr_rounding) {
mpfr_init2(value, mpfr_precision);
mpfr_set(value, other.value, mpfr_rounding);
}
MPFRNumber(const MPFRNumber &other, unsigned int precision)
: mpfr_precision(precision), mpfr_rounding(other.mpfr_rounding) {
mpfr_init2(value, mpfr_precision);
mpfr_set(value, other.value, mpfr_rounding);
}
~MPFRNumber() { mpfr_clear(value); }
MPFRNumber &operator=(const MPFRNumber &rhs) {
mpfr_precision = rhs.mpfr_precision;
mpfr_rounding = rhs.mpfr_rounding;
mpfr_set(value, rhs.value, mpfr_rounding);
return *this;
}
bool is_nan() const { return mpfr_nan_p(value); }
MPFRNumber abs() const {
MPFRNumber result(*this);
mpfr_abs(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber acos() const {
MPFRNumber result(*this);
mpfr_acos(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber acosh() const {
MPFRNumber result(*this);
mpfr_acosh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber add(const MPFRNumber &b) const {
MPFRNumber result(*this);
mpfr_add(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber asin() const {
MPFRNumber result(*this);
mpfr_asin(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber asinh() const {
MPFRNumber result(*this);
mpfr_asinh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber atan() const {
MPFRNumber result(*this);
mpfr_atan(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber atan2(const MPFRNumber &b) {
MPFRNumber result(*this);
mpfr_atan2(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber atanh() const {
MPFRNumber result(*this);
mpfr_atanh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber cbrt() const {
MPFRNumber result(*this);
mpfr_cbrt(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber ceil() const {
MPFRNumber result(*this);
mpfr_ceil(result.value, value);
return result;
}
MPFRNumber cos() const {
MPFRNumber result(*this);
mpfr_cos(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber cosh() const {
MPFRNumber result(*this);
mpfr_cosh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber cospi() const {
MPFRNumber result(*this);
#if MPFR_VERSION_MAJOR > 4 || \
(MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2)
mpfr_cospi(result.value, value, mpfr_rounding);
return result;
#else
if (mpfr_integer_p(value)) {
mpz_t integer;
mpz_init(integer);
mpfr_get_z(integer, value, mpfr_rounding);
int d = mpz_tstbit(integer, 0);
mpfr_set_si(result.value, d ? -1 : 1, mpfr_rounding);
mpz_clear(integer);
return result;
}
MPFRNumber value_pi(0.0, 1280);
mpfr_const_pi(value_pi.value, MPFR_RNDN);
mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN);
mpfr_cos(result.value, value_pi.value, mpfr_rounding);
return result;
#endif
}
MPFRNumber erf() const {
MPFRNumber result(*this);
mpfr_erf(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber exp() const {
MPFRNumber result(*this);
mpfr_exp(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber exp2() const {
MPFRNumber result(*this);
mpfr_exp2(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber exp2m1() const {
// TODO: Only use mpfr_exp2m1 once CI and buildbots get MPFR >= 4.2.0.
#if MPFR_VERSION_MAJOR > 4 || \
(MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2)
MPFRNumber result(*this);
mpfr_exp2m1(result.value, value, mpfr_rounding);
return result;
#else
unsigned int prec = mpfr_precision * 3;
MPFRNumber result(*this, prec);
float f = mpfr_get_flt(abs().value, mpfr_rounding);
if (f > 0.5f && f < 0x1.0p30f) {
mpfr_exp2(result.value, value, mpfr_rounding);
mpfr_sub_ui(result.value, result.value, 1, mpfr_rounding);
return result;
}
MPFRNumber ln2(2.0f, prec);
// log(2)
mpfr_log(ln2.value, ln2.value, mpfr_rounding);
// x * log(2)
mpfr_mul(result.value, value, ln2.value, mpfr_rounding);
// e^(x * log(2)) - 1
int ex = mpfr_expm1(result.value, result.value, mpfr_rounding);
mpfr_subnormalize(result.value, ex, mpfr_rounding);
return result;
#endif
}
MPFRNumber exp10() const {
MPFRNumber result(*this);
mpfr_exp10(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber exp10m1() const {
// TODO: Only use mpfr_exp10m1 once CI and buildbots get MPFR >= 4.2.0.
#if MPFR_VERSION_MAJOR > 4 || \
(MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2)
MPFRNumber result(*this);
mpfr_exp10m1(result.value, value, mpfr_rounding);
return result;
#else
unsigned int prec = mpfr_precision * 3;
MPFRNumber result(*this, prec);
MPFRNumber ln10(10.0f, prec);
// log(10)
mpfr_log(ln10.value, ln10.value, mpfr_rounding);
// x * log(10)
mpfr_mul(result.value, value, ln10.value, mpfr_rounding);
// e^(x * log(10)) - 1
int ex = mpfr_expm1(result.value, result.value, mpfr_rounding);
mpfr_subnormalize(result.value, ex, mpfr_rounding);
return result;
#endif
}
MPFRNumber expm1() const {
MPFRNumber result(*this);
mpfr_expm1(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber div(const MPFRNumber &b) const {
MPFRNumber result(*this);
mpfr_div(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber floor() const {
MPFRNumber result(*this);
mpfr_floor(result.value, value);
return result;
}
MPFRNumber fmod(const MPFRNumber &b) {
MPFRNumber result(*this);
mpfr_fmod(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber frexp(int &exp) {
MPFRNumber result(*this);
mpfr_exp_t resultExp;
mpfr_frexp(&resultExp, result.value, value, mpfr_rounding);
exp = resultExp;
return result;
}
MPFRNumber hypot(const MPFRNumber &b) {
MPFRNumber result(*this);
mpfr_hypot(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber log() const {
MPFRNumber result(*this);
mpfr_log(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber log2() const {
MPFRNumber result(*this);
mpfr_log2(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber log10() const {
MPFRNumber result(*this);
mpfr_log10(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber log1p() const {
MPFRNumber result(*this);
mpfr_log1p(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber pow(const MPFRNumber &b) {
MPFRNumber result(*this);
mpfr_pow(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber remquo(const MPFRNumber &divisor, int &quotient) {
MPFRNumber remainder(*this);
long q;
mpfr_remquo(remainder.value, &q, value, divisor.value, mpfr_rounding);
quotient = q;
return remainder;
}
MPFRNumber round() const {
MPFRNumber result(*this);
mpfr_round(result.value, value);
return result;
}
MPFRNumber roundeven() const {
MPFRNumber result(*this);
#if MPFR_VERSION_MAJOR >= 4
mpfr_roundeven(result.value, value);
#else
mpfr_rint(result.value, value, MPFR_RNDN);
#endif
return result;
}
bool round_to_long(long &result) const {
// We first calculate the rounded value. This way, when converting
// to long using mpfr_get_si, the rounding direction of MPFR_RNDN
// (or any other rounding mode), does not have an influence.
MPFRNumber roundedValue = round();
mpfr_clear_erangeflag();
result = mpfr_get_si(roundedValue.value, MPFR_RNDN);
return mpfr_erangeflag_p();
}
bool round_to_long(mpfr_rnd_t rnd, long &result) const {
MPFRNumber rint_result(*this);
mpfr_rint(rint_result.value, value, rnd);
return rint_result.round_to_long(result);
}
MPFRNumber rint(mpfr_rnd_t rnd) const {
MPFRNumber result(*this);
mpfr_rint(result.value, value, rnd);
return result;
}
MPFRNumber mod_2pi() const {
MPFRNumber result(0.0, 1280);
MPFRNumber _2pi(0.0, 1280);
mpfr_const_pi(_2pi.value, MPFR_RNDN);
mpfr_mul_si(_2pi.value, _2pi.value, 2, MPFR_RNDN);
mpfr_fmod(result.value, value, _2pi.value, MPFR_RNDN);
return result;
}
MPFRNumber mod_pi_over_2() const {
MPFRNumber result(0.0, 1280);
MPFRNumber pi_over_2(0.0, 1280);
mpfr_const_pi(pi_over_2.value, MPFR_RNDN);
mpfr_mul_d(pi_over_2.value, pi_over_2.value, 0.5, MPFR_RNDN);
mpfr_fmod(result.value, value, pi_over_2.value, MPFR_RNDN);
return result;
}
MPFRNumber mod_pi_over_4() const {
MPFRNumber result(0.0, 1280);
MPFRNumber pi_over_4(0.0, 1280);
mpfr_const_pi(pi_over_4.value, MPFR_RNDN);
mpfr_mul_d(pi_over_4.value, pi_over_4.value, 0.25, MPFR_RNDN);
mpfr_fmod(result.value, value, pi_over_4.value, MPFR_RNDN);
return result;
}
MPFRNumber sin() const {
MPFRNumber result(*this);
mpfr_sin(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber sinpi() const {
MPFRNumber result(*this);
#if MPFR_VERSION_MAJOR > 4 || \
(MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2)
mpfr_sinpi(result.value, value, mpfr_rounding);
return result;
#else
if (mpfr_integer_p(value)) {
mpfr_set_si(result.value, 0, mpfr_rounding);
return result;
}
MPFRNumber value_mul_two(*this);
mpfr_mul_si(value_mul_two.value, value, 2, MPFR_RNDN);
if (mpfr_integer_p(value_mul_two.value)) {
auto d = mpfr_get_si(value, MPFR_RNDD);
mpfr_set_si(result.value, (d & 1) ? -1 : 1, mpfr_rounding);
return result;
}
MPFRNumber value_pi(0.0, 1280);
mpfr_const_pi(value_pi.value, MPFR_RNDN);
mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN);
mpfr_sin(result.value, value_pi.value, mpfr_rounding);
return result;
#endif
}
MPFRNumber sinh() const {
MPFRNumber result(*this);
mpfr_sinh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber sqrt() const {
MPFRNumber result(*this);
mpfr_sqrt(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber sub(const MPFRNumber &b) const {
MPFRNumber result(*this);
mpfr_sub(result.value, value, b.value, mpfr_rounding);
return result;
}
MPFRNumber tan() const {
MPFRNumber result(*this);
mpfr_tan(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber tanh() const {
MPFRNumber result(*this);
mpfr_tanh(result.value, value, mpfr_rounding);
return result;
}
MPFRNumber tanpi() const {
MPFRNumber result(*this);
#if MPFR_VERSION_MAJOR > 4 || \
(MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2)
mpfr_tanpi(result.value, value, mpfr_rounding);
return result;
#else
MPFRNumber value_ret_exact(*this);
MPFRNumber value_one(*this);
mpfr_set_si(value_one.value, 1, MPFR_RNDN);
mpfr_fmod(value_ret_exact.value, value, value_one.value, mpfr_rounding);
mpfr_mul_si(value_ret_exact.value, value_ret_exact.value, 4, MPFR_RNDN);
if (mpfr_integer_p(value_ret_exact.value)) {
int mod = mpfr_get_si(value_ret_exact.value, MPFR_RNDN);
mod = (mod < 0 ? -1 * mod : mod);
switch (mod) {
case 0:
mpfr_set_si(result.value, 0, mpfr_rounding);
break;
case 1:
mpfr_set_si(result.value, (mpfr_signbit(value) ? -1 : 1),
mpfr_rounding);
break;
case 2: {
auto d = mpfr_get_si(value, MPFR_RNDZ);
d += mpfr_sgn(value) > 0 ? 0 : 1;
mpfr_set_inf(result.value, (d & 1) ? -1 : 1);
break;
}
case 3:
mpfr_set_si(result.value, (mpfr_signbit(value) ? 1 : -1),
mpfr_rounding);
break;
}
return result;
}
MPFRNumber value_pi(0.0, 1280);
mpfr_const_pi(value_pi.value, MPFR_RNDN);
mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN);
mpfr_tan(result.value, value_pi.value, mpfr_rounding);
return result;
#endif
}
MPFRNumber trunc() const {
MPFRNumber result(*this);
mpfr_trunc(result.value, value);
return result;
}
MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c) {
MPFRNumber result(*this);
mpfr_fma(result.value, value, b.value, c.value, mpfr_rounding);
return result;
}
MPFRNumber mul(const MPFRNumber &b) {
MPFRNumber result(*this);
mpfr_mul(result.value, value, b.value, mpfr_rounding);
return result;
}
cpp::string str() const {
// 200 bytes should be more than sufficient to hold a 100-digit number
// plus additional bytes for the decimal point, '-' sign etc.
constexpr size_t printBufSize = 200;
char buffer[printBufSize];
mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
cpp::string_view view(buffer);
// Trim whitespaces
const char whitespace = ' ';
while (!view.empty() && view.front() == whitespace)
view.remove_prefix(1);
while (!view.empty() && view.back() == whitespace)
view.remove_suffix(1);
return cpp::string(view.data());
}
// These functions are useful for debugging.
template <typename T> T as() const;
void dump(const char *msg) const { mpfr_printf("%s%.128g\n", msg, value); }
// Return the ULP (units-in-the-last-place) difference between the
// stored MPFR and a floating point number.
//
// We define ULP difference as follows:
// If exponents of this value and the |input| are same, then:
// ULP(this_value, input) = abs(this_value - input) / eps(input)
// else:
// max = max(abs(this_value), abs(input))
// min = min(abs(this_value), abs(input))
// maxExponent = exponent(max)
// ULP(this_value, input) = (max - 2^maxExponent) / eps(max) +
// (2^maxExponent - min) / eps(min)
//
// Remarks:
// 1. A ULP of 0.0 will imply that the value is correctly rounded.
// 2. We expect that this value and the value to be compared (the [input]
// argument) are reasonable close, and we will provide an upper bound
// of ULP value for testing. Morever, most of the fractional parts of
// ULP value do not matter much, so using double as the return type
// should be good enough.
// 3. For close enough values (values which don't diff in their exponent by
// not more than 1), a ULP difference of N indicates a bit distance
// of N between this number and [input].
// 4. A values of +0.0 and -0.0 are treated as equal.
template <typename T>
cpp::enable_if_t<cpp::is_floating_point_v<T>, MPFRNumber>
ulp_as_mpfr_number(T input) {
T thisAsT = as<T>();
if (thisAsT == input)
return MPFRNumber(0.0);
if (is_nan()) {
if (FPBits<T>(input).is_nan())
return MPFRNumber(0.0);
return MPFRNumber(FPBits<T>::inf().get_val());
}
int thisExponent = FPBits<T>(thisAsT).get_exponent();
int inputExponent = FPBits<T>(input).get_exponent();
// Adjust the exponents for denormal numbers.
if (FPBits<T>(thisAsT).is_subnormal())
++thisExponent;
if (FPBits<T>(input).is_subnormal())
++inputExponent;
if (thisAsT * input < 0 || thisExponent == inputExponent) {
MPFRNumber inputMPFR(input);
mpfr_sub(inputMPFR.value, value, inputMPFR.value, MPFR_RNDN);
mpfr_abs(inputMPFR.value, inputMPFR.value, MPFR_RNDN);
mpfr_mul_2si(inputMPFR.value, inputMPFR.value,
-thisExponent + FPBits<T>::FRACTION_LEN, MPFR_RNDN);
return inputMPFR;
}
// If the control reaches here, it means that this number and input are
// of the same sign but different exponent. In such a case, ULP error is
// calculated as sum of two parts.
thisAsT = FPBits<T>(thisAsT).abs().get_val();
input = FPBits<T>(input).abs().get_val();
T min = thisAsT > input ? input : thisAsT;
T max = thisAsT > input ? thisAsT : input;
int minExponent = FPBits<T>(min).get_exponent();
int maxExponent = FPBits<T>(max).get_exponent();
// Adjust the exponents for denormal numbers.
if (FPBits<T>(min).is_subnormal())
++minExponent;
if (FPBits<T>(max).is_subnormal())
++maxExponent;
MPFRNumber minMPFR(min);
MPFRNumber maxMPFR(max);
MPFRNumber pivot(uint32_t(1));
mpfr_mul_2si(pivot.value, pivot.value, maxExponent, MPFR_RNDN);
mpfr_sub(minMPFR.value, pivot.value, minMPFR.value, MPFR_RNDN);
mpfr_mul_2si(minMPFR.value, minMPFR.value,
-minExponent + FPBits<T>::FRACTION_LEN, MPFR_RNDN);
mpfr_sub(maxMPFR.value, maxMPFR.value, pivot.value, MPFR_RNDN);
mpfr_mul_2si(maxMPFR.value, maxMPFR.value,
-maxExponent + FPBits<T>::FRACTION_LEN, MPFR_RNDN);
mpfr_add(minMPFR.value, minMPFR.value, maxMPFR.value, MPFR_RNDN);
return minMPFR;
}
template <typename T>
cpp::enable_if_t<cpp::is_floating_point_v<T>, cpp::string>
ulp_as_string(T input) {
MPFRNumber num = ulp_as_mpfr_number(input);
return num.str();
}
template <typename T>
cpp::enable_if_t<cpp::is_floating_point_v<T>, double> ulp(T input) {
MPFRNumber num = ulp_as_mpfr_number(input);
return num.as<double>();
}
};
template <> float MPFRNumber::as<float>() const {
return mpfr_get_flt(value, mpfr_rounding);
}
template <> double MPFRNumber::as<double>() const {
return mpfr_get_d(value, mpfr_rounding);
}
template <> long double MPFRNumber::as<long double>() const {
return mpfr_get_ld(value, mpfr_rounding);
}
#ifdef LIBC_TYPES_HAS_FLOAT16
template <> float16 MPFRNumber::as<float16>() const {
// TODO: Either prove that this cast won't cause double-rounding errors, or
// find a better way to get a float16.
return fputil::cast<float16>(mpfr_get_d(value, mpfr_rounding));
}
#endif
#ifdef LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
template <> float128 MPFRNumber::as<float128>() const {
return mpfr_get_float128(value, mpfr_rounding);
}
#endif // LIBC_TYPES_FLOAT128_IS_NOT_LONG_DOUBLE
namespace internal {
template <typename InputType>

View File

@ -21,7 +21,7 @@ namespace testing {
namespace mpfr {
enum class Operation : int {
// Operations with take a single floating point number as input
// Operations which take a single floating point number as input
// and produce a single floating point number as output. The input
// and output floating point numbers are of the same kind.
BeginUnaryOperationsSingleOutput,
@ -87,10 +87,10 @@ enum class Operation : int {
EndBinaryOperationsSingleOutput,
// Operations which take two floating point numbers of the same type as
// input and produce two outputs. The first output is a floating nubmer of
// the same type as the inputs. The second output is af type 'int'.
// input and produce two outputs. The first output is a floating point number
// of the same type as the inputs. The second output is of type 'int'.
BeginBinaryOperationsTwoOutputs,
RemQuo, // The first output, the floating point output, is the remainder.
RemQuo, // The first output(floating point) is the remainder.
EndBinaryOperationsTwoOutputs,
// Operations which take three floating point nubmers of the same type as