mirror of
https://github.com/boostorg/math.git
synced 2025-05-11 21:33:52 +00:00
329 lines
13 KiB
C++
329 lines
13 KiB
C++
// Copyright 2024 Christopher Kormanyos
|
|
// Distributed under the Boost Software License, Version 1.0.
|
|
// https://www.boost.org/LICENSE_1_0.txt
|
|
|
|
#include <boost/math/special_functions/gamma.hpp>
|
|
#include <boost/core/lightweight_test.hpp>
|
|
#include <boost/multiprecision/cpp_dec_float.hpp>
|
|
|
|
#include <random>
|
|
|
|
namespace local
|
|
{
|
|
std::mt19937 eng(static_cast<typename std::mt19937::result_type>(UINT8_C(42)));
|
|
std::uniform_int_distribution<int> dst_one(1, 1);
|
|
|
|
template<typename NumericType>
|
|
auto is_close_fraction(const NumericType& a,
|
|
const NumericType& b,
|
|
const NumericType& tol) noexcept -> bool
|
|
{
|
|
using std::fabs;
|
|
|
|
auto result_is_ok = bool { };
|
|
|
|
if(b == static_cast<NumericType>(0))
|
|
{
|
|
result_is_ok = (fabs(a - b) < tol);
|
|
}
|
|
else
|
|
{
|
|
const auto delta = fabs(1 - (a / b));
|
|
|
|
result_is_ok = (delta < tol);
|
|
}
|
|
|
|
return result_is_ok;
|
|
}
|
|
|
|
auto tgamma_under_cbrt_epsilon() -> void
|
|
{
|
|
// This test is intended to hit the lines:
|
|
|
|
// template <class T, class Policy>
|
|
// T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&)
|
|
// ...
|
|
// ...near the comment:
|
|
// Special case for ultra-small z:
|
|
|
|
using local_float_type = boost::multiprecision::number<boost::multiprecision::cpp_dec_float<250>, boost::multiprecision::et_off>;
|
|
|
|
static_assert( (std::numeric_limits<local_float_type>::digits10 >= 248)
|
|
&& (std::numeric_limits<local_float_type>::digits10 <= 252), "Error: Multiprecision wrong number of digits");
|
|
|
|
// Table[N[Gamma[n (10^-84)], 260], {n, 1, 10, 1}]
|
|
using local_data_array_type = std::array<local_float_type, static_cast<std::size_t>(UINT8_C(10))>;
|
|
|
|
const local_data_array_type ctrl_data =
|
|
{{
|
|
static_cast<local_float_type>("9.9999999999999999999999999999999999999999999999999999999999999999999999999999999999942278433509846713939348790991759756895784066406007640119423276511513227322233532906404199270358122304076394840169355222697867950624167125421902178604121309579597843800360126781E+83"),
|
|
static_cast<local_float_type>("4.9999999999999999999999999999999999999999999999999999999999999999999999999999999999942278433509846713939348790991759756895784066406007640119423276511513227322233532906503104869890919559615934405319418693491786302696988534466221756972734972784545720974832491521E+83"),
|
|
static_cast<local_float_type>("3.3333333333333333333333333333333333333333333333333333333333333333333333333333333333275611766843180047272682124325093090229117399739340973452756609844846560655566866239935343802757050148488807303802815497619037988103143276843874668674681969322826931482456693779E+83"),
|
|
static_cast<local_float_type>("2.4999999999999999999999999999999999999999999999999999999999999999999999999999999999942278433509846713939348790991759756895784066406007640119423276511513227322233532906700916068956514070695013535619545635079623006842631352554860913709962299194441475323232733555E+83"),
|
|
static_cast<local_float_type>("1.9999999999999999999999999999999999999999999999999999999999999999999999999999999999942278433509846713939348790991759756895784066406007640119423276511513227322233532906799821668489311326234553100769609105873541358915452761599180492078575962399389352497160610849E+83"),
|
|
static_cast<local_float_type>("1.6666666666666666666666666666666666666666666666666666666666666666666666666666666666608945100176513380606015457658426423562450733072674306786089943178179893988900199573565393934688775248440759332586339243334126377654940837310166737113856292271003896337573658994E+83"),
|
|
static_cast<local_float_type>("1.4285714285714285714285714285714285714285714285714285714285714285714285714285714285656564147795560999653634505277474042610069780691721925833708990797227513036519247192711918581840620123027917945355450333175663777346809865402105363101517574523570821130186163706E+83"),
|
|
static_cast<local_float_type>("1.2499999999999999999999999999999999999999999999999999999999999999999999999999999999942278433509846713939348790991759756895784066406007640119423276511513227322233532907096538467087703092853171796219799518255296415133916988732139227184416952014232984017855267840E+83"),
|
|
static_cast<local_float_type>("1.1111111111111111111111111111111111111111111111111111111111111111111111111111111111053389544620957825050459902102870868006895177517118751230534387622624338433344644018306555177731611459503822472480974100160325878317849508887569916664141726330291972302168272984E+83"),
|
|
static_cast<local_float_type>("9.9999999999999999999999999999999999999999999999999999999999999999999999999999999999422784335098467139393487909917597568957840664060076401194232765115132273222335329072943496661532976039322509265199264598431331192795598068207783839216442784241287383640775600912E+82"),
|
|
}};
|
|
|
|
unsigned index = 1U;
|
|
|
|
const local_float_type little { "1E-84" };
|
|
const local_float_type my_tol { std::numeric_limits<local_float_type>::epsilon() * 256 };
|
|
|
|
for(const auto& ctrl : ctrl_data)
|
|
{
|
|
const auto x_small = static_cast<local_float_type>(static_cast<local_float_type>(index) * little);
|
|
|
|
++index;
|
|
|
|
const auto g_val = boost::math::tgamma(x_small);
|
|
|
|
const auto result_tgamma_x_small_is_ok = is_close_fraction(g_val, ctrl, my_tol);
|
|
|
|
BOOST_TEST(result_tgamma_x_small_is_ok);
|
|
|
|
if(!result_tgamma_x_small_is_ok)
|
|
{
|
|
break; // LCOV_EXCL_LINE
|
|
}
|
|
}
|
|
}
|
|
|
|
auto tgamma_undefined_lanczos_known_error() -> void
|
|
{
|
|
// This test is intended to hit the lines:
|
|
|
|
// template <class T, class Policy>
|
|
// T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&)
|
|
// ...
|
|
// ...for edge cases that raise errors such as domain error.
|
|
|
|
using local_float_type = boost::multiprecision::number<boost::multiprecision::cpp_dec_float<250>, boost::multiprecision::et_off>;
|
|
|
|
{
|
|
const local_float_type my_tol { std::numeric_limits<local_float_type>::epsilon() * 256 };
|
|
|
|
for(auto index = static_cast<unsigned>(UINT8_C(0)); index < static_cast<unsigned>(UINT8_C(3)); ++index)
|
|
{
|
|
static_cast<void>(index);
|
|
|
|
const local_float_type zero_ctrl { 0 };
|
|
|
|
local_float_type zero { 0 };
|
|
|
|
zero *= dst_one(eng);
|
|
|
|
const auto result_zero_is_ok = is_close_fraction(zero, zero_ctrl, my_tol);
|
|
|
|
BOOST_TEST(result_zero_is_ok);
|
|
}
|
|
}
|
|
|
|
for(auto index = static_cast<unsigned>(UINT8_C(0)); index < static_cast<unsigned>(UINT8_C(3)); ++index)
|
|
{
|
|
static_cast<void>(index);
|
|
|
|
local_float_type zero { 0 };
|
|
|
|
zero *= dst_one(eng);
|
|
|
|
bool domain_error_is_ok { false };
|
|
|
|
try
|
|
{
|
|
boost::math::tgamma(zero);
|
|
}
|
|
catch(std::domain_error& err)
|
|
{
|
|
static_cast<void>(err.what());
|
|
|
|
domain_error_is_ok = true;
|
|
}
|
|
|
|
BOOST_TEST(domain_error_is_ok);
|
|
}
|
|
|
|
for(auto index = static_cast<unsigned>(UINT8_C(0)); index < static_cast<unsigned>(UINT8_C(3)); ++index)
|
|
{
|
|
static_cast<void>(index);
|
|
|
|
local_float_type my_nan = std::numeric_limits<local_float_type>::quiet_NaN();
|
|
|
|
my_nan *= dst_one(eng);
|
|
|
|
bool domain_error_is_ok { false };
|
|
|
|
try
|
|
{
|
|
boost::math::tgamma(my_nan);
|
|
}
|
|
catch(std::domain_error& err)
|
|
{
|
|
static_cast<void>(err.what());
|
|
|
|
domain_error_is_ok = true;
|
|
}
|
|
|
|
BOOST_TEST(domain_error_is_ok);
|
|
}
|
|
|
|
for(auto index = static_cast<unsigned>(UINT8_C(0)); index < static_cast<unsigned>(UINT8_C(3)); ++index)
|
|
{
|
|
static_cast<void>(index);
|
|
|
|
local_float_type my_inf = -std::numeric_limits<local_float_type>::infinity();
|
|
|
|
my_inf *= dst_one(eng);
|
|
|
|
bool domain_error_is_ok { false };
|
|
|
|
try
|
|
{
|
|
boost::math::tgamma(my_inf);
|
|
}
|
|
catch(std::domain_error& err)
|
|
{
|
|
static_cast<void>(err.what());
|
|
|
|
domain_error_is_ok = true;
|
|
}
|
|
|
|
BOOST_TEST(domain_error_is_ok);
|
|
}
|
|
}
|
|
|
|
auto lgamma_big_asymp() -> void
|
|
{
|
|
// This test is intended to hit the asymptotic log-gamma expansion for multiprecision.
|
|
|
|
using local_float_type = boost::multiprecision::number<boost::multiprecision::cpp_dec_float<250>, boost::multiprecision::et_off>;
|
|
|
|
static_assert( (std::numeric_limits<local_float_type>::digits10 >= 248)
|
|
&& (std::numeric_limits<local_float_type>::digits10 <= 252), "Error: Multiprecision wrong number of digits");
|
|
|
|
local_float_type big_arg_numerator { 1234567L };
|
|
|
|
// Table[N[Log[Gamma[(1234567 + n)/1000]], 260], {n, 0, 3, 1}]
|
|
const local_float_type ctrl0 { "7551.0278099842760398085493506933061185258592164059260052791257174648102458654516760859347475429811747227042884941464597963128452844941163716092798494933305452087249880911022309522317482008162381529082884245980549740815352929296384544778543502768128060636123031" };
|
|
const local_float_type ctrl1 { "7551.0349280552065308610373629214633349814110368633190642156085097598877230874250481117271260334496206128158535271616589000730715715804390525860149840442193710637326207809853649225510544815601053550751028151966244578864039961973124357117676870769851159530881598" };
|
|
const local_float_type ctrl2 { "7551.0420461269473499872464481408311466395059218944034402487444941839649424484440812974771163653823079980811183096059838617975730358237216549153177603277276543261570202589028877022053234067330484829335430150182956719650949866427225508960918047040961544342635607" };
|
|
const local_float_type ctrl3 { "7551.0491641994984965305455873077287163417673805051480753266140277293097406691212685736787897808049319210175552316063915681371903149757892763598947756743285838029800201937242960777121582844482027617050641997959480250230315703285501155035159657216828260204447321" };
|
|
|
|
const local_float_type my_tol { std::numeric_limits<local_float_type>::epsilon() * 256 };
|
|
|
|
const local_float_type lg_big0 { boost::math::lgamma(big_arg_numerator / 1000) }; ++big_arg_numerator;
|
|
const local_float_type lg_big1 { boost::math::lgamma(big_arg_numerator / 1000) }; ++big_arg_numerator;
|
|
const local_float_type lg_big2 { boost::math::lgamma(big_arg_numerator / 1000) }; ++big_arg_numerator;
|
|
const local_float_type lg_big3 { boost::math::lgamma(big_arg_numerator / 1000) };
|
|
|
|
const auto result_lgamma_big0_is_ok = is_close_fraction(lg_big0, ctrl0, my_tol);
|
|
const auto result_lgamma_big1_is_ok = is_close_fraction(lg_big1, ctrl1, my_tol);
|
|
const auto result_lgamma_big2_is_ok = is_close_fraction(lg_big2, ctrl2, my_tol);
|
|
const auto result_lgamma_big3_is_ok = is_close_fraction(lg_big3, ctrl3, my_tol);
|
|
|
|
BOOST_TEST(result_lgamma_big0_is_ok);
|
|
BOOST_TEST(result_lgamma_big1_is_ok);
|
|
BOOST_TEST(result_lgamma_big2_is_ok);
|
|
BOOST_TEST(result_lgamma_big3_is_ok);
|
|
}
|
|
|
|
auto lgamma_undefined_lanczos_known_error() -> void
|
|
{
|
|
// This test is intended to hit the lines:
|
|
|
|
// template <class T, class Policy>
|
|
// T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign) // ...
|
|
// ...
|
|
// ...for edge cases that raise errors such as domain error.
|
|
|
|
using local_float_type = boost::multiprecision::number<boost::multiprecision::cpp_dec_float<250>, boost::multiprecision::et_off>;
|
|
|
|
for(auto index = static_cast<unsigned>(UINT8_C(0)); index < static_cast<unsigned>(UINT8_C(3)); ++index)
|
|
{
|
|
static_cast<void>(index);
|
|
|
|
local_float_type zero { 0 };
|
|
|
|
zero *= dst_one(eng);
|
|
|
|
bool domain_error_is_ok { false };
|
|
|
|
try
|
|
{
|
|
boost::math::lgamma(zero);
|
|
}
|
|
catch(std::domain_error& err)
|
|
{
|
|
static_cast<void>(err.what());
|
|
|
|
domain_error_is_ok = true;
|
|
}
|
|
|
|
BOOST_TEST(domain_error_is_ok);
|
|
}
|
|
|
|
for(auto index = static_cast<unsigned>(UINT8_C(0)); index < static_cast<unsigned>(UINT8_C(3)); ++index)
|
|
{
|
|
static_cast<void>(index);
|
|
|
|
local_float_type my_nan = std::numeric_limits<local_float_type>::quiet_NaN();
|
|
|
|
my_nan *= dst_one(eng);
|
|
|
|
bool domain_error_is_ok { false };
|
|
|
|
try
|
|
{
|
|
boost::math::lgamma(my_nan);
|
|
}
|
|
catch(std::domain_error& err)
|
|
{
|
|
static_cast<void>(err.what());
|
|
|
|
domain_error_is_ok = true;
|
|
}
|
|
|
|
BOOST_TEST(domain_error_is_ok);
|
|
}
|
|
|
|
for(auto index = static_cast<unsigned>(UINT8_C(0)); index < static_cast<unsigned>(UINT8_C(3)); ++index)
|
|
{
|
|
static_cast<void>(index);
|
|
|
|
local_float_type my_inf = -std::numeric_limits<local_float_type>::infinity();
|
|
|
|
my_inf *= dst_one(eng);
|
|
|
|
bool overflow_error_is_ok { false };
|
|
|
|
try
|
|
{
|
|
boost::math::lgamma(my_inf);
|
|
}
|
|
catch(std::overflow_error& err)
|
|
{
|
|
static_cast<void>(err.what());
|
|
|
|
overflow_error_is_ok = true;
|
|
}
|
|
|
|
BOOST_TEST(overflow_error_is_ok);
|
|
}
|
|
}
|
|
}
|
|
|
|
auto main() -> int
|
|
{
|
|
local::tgamma_under_cbrt_epsilon();
|
|
local::tgamma_undefined_lanczos_known_error();
|
|
local::lgamma_big_asymp();
|
|
local::lgamma_undefined_lanczos_known_error();
|
|
|
|
const auto result_is_ok = (boost::report_errors() == 0);
|
|
|
|
return (result_is_ok ? 0 : -1);
|
|
}
|