Files
multiprecision/test/test_pow.cpp
2025-07-19 22:30:57 +02:00

484 lines
16 KiB
C++

///////////////////////////////////////////////////////////////
// Copyright 2011 - 2025 John Maddock.
// Copyright Christopher Kormanyos 2021 - 2025.
// Distributed under the Boost
// Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
//
// This work is based on an earlier work:
// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
#ifdef _MSC_VER
#define _SCL_SECURE_NO_WARNINGS
#endif
#include <test.hpp>
#include <test_pow_data.hpp>
#include <test_pow_data_reduced.hpp>
#include <boost/detail/lightweight_test.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#if !defined(TEST_MPF_50) && !defined(TEST_MPF) && !defined(TEST_BACKEND) && !defined(TEST_CPP_DEC_FLOAT) && !defined(TEST_MPFR) && !defined(TEST_MPFR_50) && !defined(TEST_MPFI_50) && !defined(TEST_FLOAT128) && !defined(TEST_CPP_BIN_FLOAT) && !defined(TEST_CPP_DOUBLE_FLOAT)
#define TEST_MPF_50
//# define TEST_MPF
#define TEST_BACKEND
#define TEST_CPP_DEC_FLOAT
#define TEST_FLOAT128
#define TEST_MPFR_50
#define TEST_MPFI_50
#define TEST_CPP_BIN_FLOAT
#define TEST_CPP_DOUBLE_FLOAT
#ifdef _MSC_VER
#pragma message("CAUTION!!: No backend type specified so testing everything.... this will take some time!!")
#endif
#ifdef __GNUC__
#pragma warning "CAUTION!!: No backend type specified so testing everything.... this will take some time!!"
#endif
#endif
#include <test_traits.hpp> // Note: include this AFTER the test-backends are defined
#if defined(TEST_MPF_50)
#include <boost/multiprecision/gmp.hpp>
#endif
#if defined(TEST_MPFR_50)
#include <boost/multiprecision/mpfr.hpp>
#endif
#if defined(TEST_MPFI_50)
#include <boost/multiprecision/mpfi.hpp>
#endif
#ifdef TEST_BACKEND
#include <boost/multiprecision/concepts/mp_number_archetypes.hpp>
#endif
#ifdef TEST_CPP_DEC_FLOAT
#include <boost/multiprecision/cpp_dec_float.hpp>
#endif
#ifdef TEST_FLOAT128
#include <boost/multiprecision/float128.hpp>
#endif
#ifdef TEST_CPP_BIN_FLOAT
#include <boost/multiprecision/cpp_bin_float.hpp>
#endif
#ifdef TEST_CPP_DOUBLE_FLOAT
#include <boost/multiprecision/cpp_double_fp.hpp>
#endif
#if (defined(TEST_CPP_DEC_FLOAT) || defined(TEST_CPP_BIN_FLOAT) || defined(TEST_CPP_DOUBLE_FLOAT) || defined(TEST_MPF_50) || defined(TEST_FLOAT128))
// These were the only ones I checked locally at the moment that use random tests.
#include <array>
#include <ctime>
#include <random>
#endif // (defined(TEST_CPP_DEC_FLOAT) || defined(TEST_CPP_BIN_FLOAT) || defined(TEST_CPP_DOUBLE_FLOAT) || defined(TEST_MPF_50) || defined(TEST_FLOAT128))
template<typename FloatType> auto my_zero() -> FloatType&;
namespace local {
template <class DataArrayType>
struct data_maker
{
private:
using local_array_type = DataArrayType;
public:
static auto get_data() -> const local_array_type&
{
// We do not expect to ever instantiate this struct.
static const local_array_type instance { };
return instance;
}
};
template <>
struct data_maker<test_pow_data::test_pow_data_array_type_reduced> final
{
private:
using local_array_type = test_pow_data::test_pow_data_array_type_reduced;
public:
static auto get_data() -> const local_array_type&
{
static const local_array_type instance { test_pow_data::data_reduced };
return instance;
}
};
template <>
struct data_maker<test_pow_data::test_pow_data_array_type_default> final
{
private:
using local_array_type = test_pow_data::test_pow_data_array_type_default;
public:
static auto get_data() -> const local_array_type&
{
static const local_array_type instance { test_pow_data::data };
return instance;
}
};
} // namespace local
#if (defined(TEST_CPP_DEC_FLOAT) || defined(TEST_CPP_BIN_FLOAT) || defined(TEST_CPP_DOUBLE_FLOAT) || defined(TEST_MPF_50) || defined(TEST_FLOAT128))
// These were the only ones I checked locally at the moment.
// TEST_MPF_50 is included (even though it won't be tested)
// because its value of has_quiet_NaN is false. This has been
// done intentionally. Its inclusion here is intended to verify
// the constexpr logic implemented in test_issue722().
template <class T>
void test_issue722()
{
BOOST_IF_CONSTEXPR (std::numeric_limits<T>::is_specialized && std::numeric_limits<T>::has_infinity && std::numeric_limits<T>::has_quiet_NaN)
{
std::mt19937_64 gen { };
gen.seed(static_cast<typename std::mt19937_64::result_type>(std::clock()));
std::uniform_real_distribution<float>
dist
(
static_cast<float>(1.01L),
static_cast<float>(1.04L)
);
for (int index = 0; index < 8; ++index)
{
static_cast<void>(index);
const T neg_inf = T(-std::numeric_limits<T>::infinity() * dist(gen));
const T pow_neg_inf_to_the_nan = pow(neg_inf, std::numeric_limits<T>::quiet_NaN() * dist(gen));
BOOST_CHECK((boost::multiprecision::isnan)(pow_neg_inf_to_the_nan));
}
for (int index = 0; index < 8; ++index)
{
static_cast<void>(index);
const T val_x_zero { ::my_zero<T>() * dist(gen) };
const T val_a_zero { ::my_zero<T>() * dist(gen) };
const T pow_zero_to_the_zero = pow(val_x_zero, val_a_zero);
BOOST_CHECK(pow_zero_to_the_zero == 1);
}
}
}
template <class T>
void test_small_a_in_default_ops()
{
// Table[N[(789/1000)^(10^(n)), 101], {n, -12, -6, 1}]
using small_a_data_array_type = std::array<T, std::size_t { UINT8_C(7) }>;
const small_a_data_array_type small_a_data =
{{
T("0.99999999999976301104186376527231033639051940385775656688844414344774455623761946001430627233999497696"),
T("0.99999999999763011041864018009258589471870605197334212386735372254907219885974849025436678031854699987"),
T("0.99999999997630110418665453787411005197468225239383193650701563365941781771201272808876889833419359752"),
T("0.99999999976301104189181907356423443492950586157450272561606631208248254028299217451384055629333440182"),
T("0.99999999763011042144556021597917228924335227306284422554015278215963516103516330474932028841209422025"),
T("0.99999997630110446719254821691095026610571257635164394962451615449422537744479656992450152929555193362"),
T("0.99999976301106994561811131758260769485991144451216442896857591003045319114657385381607667593523204524")
}};
unsigned max_err = 0U;
const T x_val(T(789U) / 1000);
const T local_tol { std::numeric_limits<T>::epsilon() / 8 };
T small_a = pow(T(10), -12);
for (std::size_t k { UINT8_C(0) }; k < std::tuple_size<small_a_data_array_type>::value; ++k)
{
const T val = pow(x_val, small_a);
small_a *= 10;
using std::fabs;
const T delta { val / small_a_data[k] };
const T closeness { fabs(1 - delta) };
const T rel = closeness / local_tol;
unsigned err = rel.template convert_to<unsigned>();
if (err > max_err)
{
max_err = err;
}
}
std::cout << "Max error data having small a: " << max_err << std::endl;
BOOST_TEST(max_err < 64U);
}
#endif // (defined(TEST_CPP_DEC_FLOAT) || defined(TEST_CPP_BIN_FLOAT) || defined(TEST_CPP_DOUBLE_FLOAT) || defined(TEST_MPF_50) || defined(TEST_FLOAT128))
template <class T>
void test()
{
std::cout << "Testing type: " << typeid(T).name() << std::endl;
using local_test_pow_data_array_type =
typename std::conditional<::has_poor_exp_range_or_precision_support<T>::value,
test_pow_data::test_pow_data_array_type_reduced,
test_pow_data::test_pow_data_array_type_default>::type;
const auto& data = local::data_maker<local_test_pow_data_array_type>::get_data();
unsigned max_err = 0U;
for (unsigned k = 0; k < data.size(); k++)
{
T val = pow(T(data[k][0]), T(data[k][1]));
T e = relative_error(val, T(data[k][2]));
unsigned err = e.template convert_to<unsigned>();
if (err > max_err)
{
max_err = err;
}
val = pow(T(data[k][0]), -T(data[k][1]));
e = relative_error(val, T(1 / T(data[k][2])));
err = e.template convert_to<unsigned>();
if (err > max_err)
{
max_err = err;
}
}
std::cout << "Max error from data table was: " << max_err << std::endl;
#if defined(BOOST_INTEL) && defined(TEST_FLOAT128)
BOOST_TEST(max_err < 40000U);
#else
BOOST_TEST(max_err < 8000U);
#endif
//
// Special cases and edge-cases come last:
//
BOOST_CHECK_EQUAL(pow(T(0), T(0)), 1);
BOOST_CHECK_EQUAL(pow(T(0), 0), 1);
BOOST_CHECK_EQUAL(pow(T(0), T(2)), 0);
BOOST_CHECK_EQUAL(pow(T(0), 2), 0);
BOOST_CHECK_EQUAL(pow(T(1), T(0)), 1);
BOOST_CHECK_EQUAL(pow(T(1), 0), 1);
BOOST_CHECK_EQUAL(pow(T(1), T(2)), 1);
BOOST_CHECK_EQUAL(pow(T(1), 2), 1);
#if (defined(TEST_CPP_DEC_FLOAT) || defined(TEST_CPP_BIN_FLOAT) || defined(TEST_CPP_DOUBLE_FLOAT) || defined(TEST_MPF_50) || defined(TEST_FLOAT128))
// These were the only ones I checked locally at the moment.
std::mt19937_64 gen { };
gen.seed(static_cast<typename std::mt19937_64::result_type>(std::clock()));
std::uniform_real_distribution<float>
dist
(
static_cast<float>(1.01L),
static_cast<float>(1.04L)
);
#endif // (defined(TEST_CPP_DEC_FLOAT) || defined(TEST_CPP_BIN_FLOAT) || defined(TEST_CPP_DOUBLE_FLOAT) || defined(TEST_MPF_50) || defined(TEST_FLOAT128))
BOOST_IF_CONSTEXPR(std::numeric_limits<T>::is_specialized && std::numeric_limits<T>::has_infinity)
{
BOOST_CHECK_EQUAL(pow(T(1.25F), std::numeric_limits<T>::infinity()), std::numeric_limits<T>::infinity());
BOOST_CHECK_EQUAL(pow(T(1.25F), -std::numeric_limits<T>::infinity()), T(0));
#if (defined(TEST_CPP_DEC_FLOAT) || defined(TEST_CPP_BIN_FLOAT) || defined(TEST_CPP_DOUBLE_FLOAT) || defined(TEST_MPF_50) || defined(TEST_FLOAT128))
// These were the only ones I checked locally at the moment.
for (int npow = -8; npow < 8; ++npow)
{
T zero_arg = my_zero<T>() * dist(gen);
const T pow_zero_to_a = pow(zero_arg, T(static_cast<float>(npow) * dist(gen)));
if (npow < 0)
{
BOOST_CHECK((boost::multiprecision::isinf)(pow_zero_to_a));
}
else if (npow > 0)
{
BOOST_CHECK_EQUAL(pow_zero_to_a, T(0));
}
else
{
BOOST_CHECK_EQUAL(pow_zero_to_a, T(1));
}
}
for (int narg = 2; narg < 8; ++narg)
{
T finite_arg = T(static_cast<float>(narg) * dist(gen));
const T pow_finite_to_zero = pow(finite_arg, my_zero<T>() * dist(gen));
BOOST_CHECK_EQUAL(pow_finite_to_zero, T(1));
}
for (int nexp = -8; nexp < 8; ++nexp)
{
T finite_exp = T(static_cast<float>(nexp) * dist(gen));
const T pow_infinite_to_finite = pow(std::numeric_limits<T>::infinity() * dist(gen), finite_exp);
if (nexp < 0)
{
BOOST_CHECK_EQUAL(pow_infinite_to_finite, T(0));
}
else if (nexp > 0)
{
BOOST_CHECK((boost::multiprecision::isinf)(pow_infinite_to_finite));
}
else
{
BOOST_CHECK_EQUAL(pow_infinite_to_finite, T(1));
}
}
BOOST_IF_CONSTEXPR(std::numeric_limits<T>::is_specialized && std::numeric_limits<T>::has_quiet_NaN)
{
for (int index = 0; index < 8; ++index)
{
static_cast<void>(index);
T inf_or_zero_arg = T(::my_zero<T>() * dist(gen));
T pow_inf_or_zero_to_the_nan = pow(inf_or_zero_arg, std::numeric_limits<T>::quiet_NaN() * dist(gen));
BOOST_CHECK((boost::multiprecision::isnan)(pow_inf_or_zero_to_the_nan));
inf_or_zero_arg = T(std::numeric_limits<T>::infinity() * dist(gen));
pow_inf_or_zero_to_the_nan = pow(inf_or_zero_arg, std::numeric_limits<T>::quiet_NaN() * dist(gen));
BOOST_CHECK((boost::multiprecision::isnan)(pow_inf_or_zero_to_the_nan));
}
}
#endif // (defined(TEST_CPP_DEC_FLOAT) || defined(TEST_CPP_BIN_FLOAT) || defined(TEST_CPP_DOUBLE_FLOAT) || defined(TEST_MPF_50) || defined(TEST_FLOAT128))
}
BOOST_IF_CONSTEXPR((!boost::multiprecision::is_interval_number<T>::value) && std::numeric_limits<T>::is_specialized && std::numeric_limits<T>::has_quiet_NaN)
{
BOOST_CHECK((boost::math::isnan)(pow(T(1.25F), std::numeric_limits<T>::quiet_NaN())));
BOOST_CHECK_EQUAL(pow(std::numeric_limits<T>::quiet_NaN(), T(0)), T(1));
BOOST_CHECK((boost::math::isnan)(pow(std::numeric_limits<T>::quiet_NaN(), T(1.25F))));
}
if (!boost::multiprecision::is_interval_number<T>::value)
{
T bug_case = -1.05 * log((std::numeric_limits<T>::max)()) / log(T(1.01));
for (unsigned i = 0; i < 100; ++i, bug_case *= 1.05)
{
if (std::numeric_limits<T>::has_infinity)
{
BOOST_CHECK_EQUAL(pow(T(1.01), bug_case), 0);
}
else
{
BOOST_CHECK_LE(pow(T(1.01), bug_case), (std::numeric_limits<T>::min)());
}
}
}
}
int main()
{
#ifdef TEST_BACKEND
test<boost::multiprecision::number<boost::multiprecision::concepts::number_backend_float_architype> >();
#endif
#ifdef TEST_MPF_50
test<boost::multiprecision::mpf_float_50>();
test<boost::multiprecision::mpf_float_100>();
test_issue722<boost::multiprecision::mpf_float_50>();
test_small_a_in_default_ops<boost::multiprecision::mpf_float_50>();
#endif
#ifdef TEST_MPFR_50
test<boost::multiprecision::mpfr_float_50>();
test<boost::multiprecision::mpfr_float_100>();
#endif
#ifdef TEST_MPFI_50
test<boost::multiprecision::mpfi_float_50>();
test<boost::multiprecision::mpfi_float_100>();
#endif
#ifdef TEST_CPP_DEC_FLOAT
test<boost::multiprecision::cpp_dec_float_50>();
test<boost::multiprecision::cpp_dec_float_100>();
test_issue722<boost::multiprecision::cpp_dec_float_50>();
test_small_a_in_default_ops<boost::multiprecision::cpp_dec_float_50>();
#ifndef SLOW_COMPLER
// Some "peculiar" digit counts which stress our code:
test<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<65> > >();
test<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<64> > >();
test<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<63> > >();
test<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<62> > >();
test<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<61, long long> > >();
test<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<60, long long> > >();
test<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<59, long long, std::allocator<char> > > >();
test<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<58, long long, std::allocator<char> > > >();
// Check low multiprecision digit counts.
test<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<9> > >();
test<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<18> > >();
#endif
#endif
#ifdef TEST_FLOAT128
test<boost::multiprecision::float128>();
test_issue722<boost::multiprecision::float128>();
test_small_a_in_default_ops<boost::multiprecision::float128>();
#endif
#ifdef TEST_CPP_BIN_FLOAT
test<boost::multiprecision::cpp_bin_float_50>();
test<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<35, boost::multiprecision::digit_base_10, std::allocator<char>, long long> > >();
test_issue722<boost::multiprecision::cpp_bin_float_50>();
test_small_a_in_default_ops<boost::multiprecision::cpp_bin_float_50>();
#endif
#ifdef TEST_CPP_DOUBLE_FLOAT
test<boost::multiprecision::cpp_double_double>();
test<boost::multiprecision::cpp_double_long_double>();
#if defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128)
test<boost::multiprecision::cpp_double_float128>();
#endif
test_issue722<boost::multiprecision::cpp_double_double>();
test_small_a_in_default_ops<boost::multiprecision::cpp_double_double>();
test_small_a_in_default_ops<boost::multiprecision::cpp_double_long_double>();
#if defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128)
test_small_a_in_default_ops<boost::multiprecision::cpp_double_float128>();
#endif
#endif
return boost::report_errors();
}
template<typename FloatType> auto my_zero() -> FloatType&
{
using float_type = FloatType;
static float_type my_val_zero(0);
return my_val_zero;
}