From 1e61cea5df00a320dc730d645ee339338c4e8672 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Mon, 14 Dec 2020 07:21:52 +0100 Subject: [PATCH] Improve fix, no pch, test spec fun locally --- .../boost/math/special_functions/gamma.hpp | 19 ++- test/Jamfile.v2 | 2 +- test/test_tgamma_for_issue396.cpp | 159 +++++++++--------- 3 files changed, 95 insertions(+), 85 deletions(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 8cc4b390b..c1b92f139 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -366,9 +366,9 @@ std::size_t highest_bernoulli_index() template int minimum_argument_for_bernoulli_recursion() { - const float digits10_of_type = (std::numeric_limits::is_specialized - ? static_cast(std::numeric_limits::digits10) - : static_cast(boost::math::tools::digits() * 0.301F)); + BOOST_CONSTEXPR_OR_CONST float digits10_of_type = + (std::numeric_limits::is_specialized ? (float) std::numeric_limits::digits10 + : (float) boost::math::tools::digits() * 0.301F); float min_arg_as_float = (float) digits10_of_type * 1.7F; @@ -377,14 +377,17 @@ int minimum_argument_for_bernoulli_recursion() // The following code sequence has been modified // within the context of issue 396. - // The calculation of the test-variable limit has now been - // made more clearly safe regarding overflow/underflow dangers. + // The calculation of the test-variable limit has now + // been protected against overflow/underflow dangers. + + // The previous line looked like this and did, in fact, + // underflow ldexp when using certain multiprecision types. - // The previous line looked like this (and did underflow ldexp for multiprecision types): // const float limit = std::ceil(std::pow(1.0f / std::ldexp(1.0f, 1-boost::math::tools::digits()), 1.0f / 20.0f)); - // The new safe versoin of the limit check is now here. - const float limit = std::exp(((digits10_of_type / 0.301F) * std::log(2.0F)) / 20.0F); + // The new safe version of the limit check is now here. + const float digits2_of_type = (digits10_of_type + 1.0F) / 0.301F; + const float limit = std::exp((digits2_of_type * std::log(2.0F)) / 20.0F); min_arg_as_float = (std::min)(min_arg_as_float, limit); } diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 94f206322..1f6edfd11 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -487,7 +487,7 @@ test-suite special_fun : [ run test_round.cpp pch ../../test/build//boost_unit_test_framework ] [ run test_spherical_harmonic.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] [ run test_sign.cpp ../../test/build//boost_unit_test_framework ] - [ run test_tgamma_for_issue396.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] + [ run test_tgamma_for_issue396.cpp ../../test/build//boost_unit_test_framework ] [ run test_tgamma_ratio.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] [ run test_trig.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] [ run test_zeta.cpp ../../test/build//boost_unit_test_framework test_instances//test_instances pch_light ] diff --git a/test/test_tgamma_for_issue396.cpp b/test/test_tgamma_for_issue396.cpp index e87c0cfa7..579126550 100644 --- a/test/test_tgamma_for_issue396.cpp +++ b/test/test_tgamma_for_issue396.cpp @@ -1,26 +1,67 @@ /////////////////////////////////////////////////////////////////// -// Copyright John Maddock 2020. // Copyright Christopher Kormanyos 2020. +// Copyright John Maddock 2020. // Distributed under the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt // or copy at http://www.boost.org/LICENSE_1_0.txt) // +#define BOOST_TEST_MAIN +#include // Boost.Test + #include #include -#include "boost/math/distributions/beta.hpp" +#include #include #include #include -#define BOOST_TEST_MODULE test_tgamma_for_issue396 -#include - // In issue 396, a bug regarding overflow was traced back to the tgamma function. // This test file tests the fix in 433 and its corresponding original bug report code. namespace local { +namespace detail { + +template +bool test_tgamma_for_issue396_value_checker() +{ + using floating_point_type = BigFloatType; + + // Table[N[Gamma[(1/2) + (10^n)], 503], {n, 0, 3, 1}] + + const boost::array control = + {{ + floating_point_type("0.88622692545275801364908374167057259139877472806119356410690389492645564229551609068747532836927233270811341181214128533311807643286221130126254685480139353423101884932655256142496258651447541311446604768963398140008731950767573986025835009509261700929272348724745632015696088776295310820270966625045319920380686673873757671683399489468292591820439772558258086938002953369671589566640492742312409245102732742609780662578082373375752136938052805399806355360503018602224183618264830685404716174941583421211"), + floating_point_type("1.1332783889487855673345741655888924755602983082751597766087234145294833900560041537176305387276072906583502717008932373348895801731780765775979953796646009714415152490764416630481375706606053932396039541459764525989187023837695167161085523804417015113740063535865261183579508922972990386756543208549178543857406373798865630303794109491220205170302558277398183764099268751365861892723863412249690833216320407918186480305202146014474770321625907339955121137559264239090240758401696425720048012081453338360E6"), + floating_point_type("9.3209631040827166083491098091419104379064970381623611540161175194120765977611623552218076053836060223609993676387199220631835256331102029826429784793420637988460945604451237342972023988743201341318701614328454618664952897316247603329530308777063116667275003586843755354841307657702809317290363831151480295446074722690100652644579131609996151999119113967501099655433566352849645431012667388627160383486515144610582794470005796689975604764040892168183647321540427819244511610500074895473959438490375652158E156"), + floating_point_type("1.2723011956950554641822441803774445695066347098655278283939929838804808618389143636393314317333622154343715992535881414698586440455330620652019981627229614973177953241634213768203151670660953863412381880742653187501307209325406338924004280546485392703623101051957976224599412003938216329590158926122017907280168159527761842471509358725974702333390709735919152262756462872191402491961250987725812831155116532550035967994387094267848607390288008530653715254376729558412833771092612838971719786622446726968E2566") + }}; + + boost::uint32_t ten_pow_n = UINT32_C(1); + + const floating_point_type tol = std::numeric_limits::epsilon() * UINT32_C(5000); + + bool result_is_ok = true; + + for(typename boost::array::size_type i = 0U; i < control.size(); ++i) + { + const floating_point_type g = boost::math::tgamma(boost::math::constants::half() + ten_pow_n); + + ten_pow_n *= UINT32_C(10); + + using std::fabs; + + const floating_point_type closeness = fabs(1 - (g / control[i])); + + result_is_ok &= (closeness < tol); + } + + return result_is_ok; +} + +} // namespace detail + bool test_cdf____for_issue396_bug_report___() { using b_dist_type = boost::math::beta_distribution; @@ -40,43 +81,6 @@ bool test_cdf____for_issue396_bug_report___() return result_is_ok; } -template -bool test_tgamma_for_issue396() -{ - using floating_point_type = BigFloatType; - - // Table[N[Gamma[(1/2) + (10^n)], 503], {n, 0, 3, 1}] - - const boost::array control = - {{ - floating_point_type("0.88622692545275801364908374167057259139877472806119356410690389492645564229551609068747532836927233270811341181214128533311807643286221130126254685480139353423101884932655256142496258651447541311446604768963398140008731950767573986025835009509261700929272348724745632015696088776295310820270966625045319920380686673873757671683399489468292591820439772558258086938002953369671589566640492742312409245102732742609780662578082373375752136938052805399806355360503018602224183618264830685404716174941583421211"), - floating_point_type("1.1332783889487855673345741655888924755602983082751597766087234145294833900560041537176305387276072906583502717008932373348895801731780765775979953796646009714415152490764416630481375706606053932396039541459764525989187023837695167161085523804417015113740063535865261183579508922972990386756543208549178543857406373798865630303794109491220205170302558277398183764099268751365861892723863412249690833216320407918186480305202146014474770321625907339955121137559264239090240758401696425720048012081453338360E6"), - floating_point_type("9.3209631040827166083491098091419104379064970381623611540161175194120765977611623552218076053836060223609993676387199220631835256331102029826429784793420637988460945604451237342972023988743201341318701614328454618664952897316247603329530308777063116667275003586843755354841307657702809317290363831151480295446074722690100652644579131609996151999119113967501099655433566352849645431012667388627160383486515144610582794470005796689975604764040892168183647321540427819244511610500074895473959438490375652158E156"), - floating_point_type("1.2723011956950554641822441803774445695066347098655278283939929838804808618389143636393314317333622154343715992535881414698586440455330620652019981627229614973177953241634213768203151670660953863412381880742653187501307209325406338924004280546485392703623101051957976224599412003938216329590158926122017907280168159527761842471509358725974702333390709735919152262756462872191402491961250987725812831155116532550035967994387094267848607390288008530653715254376729558412833771092612838971719786622446726968E2566") - }}; - - boost::uint32_t ten_pow_n = UINT32_C(1); - - const floating_point_type tol = std::numeric_limits::epsilon() * UINT32_C(100000); - - bool result_is_ok = true; - - for(typename boost::array::size_type i = 0U; i < control.size(); ++i) - { - const floating_point_type g = boost::math::tgamma(boost::math::constants::half() + ten_pow_n); - - ten_pow_n *= UINT32_C(10); - - using std::fabs; - - const floating_point_type closeness = fabs(1 - (g / control[i])); - - result_is_ok &= (closeness < tol); - } - - return result_is_ok; -} - bool test_tgamma_for_issue396_cpp_dec_float() { using cpp_dec_float_type_020 = boost::multiprecision::number, boost::multiprecision::et_off>; @@ -88,23 +92,23 @@ bool test_tgamma_for_issue396_cpp_dec_float() using cpp_dec_float_type_101 = boost::multiprecision::number, boost::multiprecision::et_off>; using cpp_dec_float_type_501 = boost::multiprecision::number, boost::multiprecision::et_off>; - const bool cpp_dec_float_020_is_ok = test_tgamma_for_issue396(); - const bool cpp_dec_float_030_is_ok = test_tgamma_for_issue396(); - const bool cpp_dec_float_040_is_ok = test_tgamma_for_issue396(); - const bool cpp_dec_float_049_is_ok = test_tgamma_for_issue396(); - const bool cpp_dec_float_050_is_ok = test_tgamma_for_issue396(); - const bool cpp_dec_float_051_is_ok = test_tgamma_for_issue396(); - const bool cpp_dec_float_101_is_ok = test_tgamma_for_issue396(); - const bool cpp_dec_float_501_is_ok = test_tgamma_for_issue396(); + const bool cpp_dec_float_020_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_dec_float_030_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_dec_float_040_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_dec_float_049_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_dec_float_050_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_dec_float_051_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_dec_float_101_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_dec_float_501_is_ok = detail::test_tgamma_for_issue396_value_checker(); - const bool cpp_dec_float_is_ok = cpp_dec_float_020_is_ok - && cpp_dec_float_030_is_ok - && cpp_dec_float_040_is_ok - && cpp_dec_float_049_is_ok - && cpp_dec_float_050_is_ok - && cpp_dec_float_051_is_ok - && cpp_dec_float_101_is_ok - && cpp_dec_float_501_is_ok; + const bool cpp_dec_float_is_ok = ( cpp_dec_float_020_is_ok + && cpp_dec_float_030_is_ok + && cpp_dec_float_040_is_ok + && cpp_dec_float_049_is_ok + && cpp_dec_float_050_is_ok + && cpp_dec_float_051_is_ok + && cpp_dec_float_101_is_ok + && cpp_dec_float_501_is_ok); return cpp_dec_float_is_ok; } @@ -120,23 +124,23 @@ bool test_tgamma_for_issue396_cpp_bin_float() using cpp_bin_float_type_101 = boost::multiprecision::number, boost::multiprecision::et_off>; using cpp_bin_float_type_501 = boost::multiprecision::number, boost::multiprecision::et_off>; - const bool cpp_bin_float_020_is_ok = test_tgamma_for_issue396(); - const bool cpp_bin_float_030_is_ok = test_tgamma_for_issue396(); - const bool cpp_bin_float_040_is_ok = test_tgamma_for_issue396(); - const bool cpp_bin_float_049_is_ok = test_tgamma_for_issue396(); - const bool cpp_bin_float_050_is_ok = test_tgamma_for_issue396(); - const bool cpp_bin_float_051_is_ok = test_tgamma_for_issue396(); - const bool cpp_bin_float_101_is_ok = test_tgamma_for_issue396(); - const bool cpp_bin_float_501_is_ok = test_tgamma_for_issue396(); + const bool cpp_bin_float_020_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_bin_float_030_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_bin_float_040_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_bin_float_049_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_bin_float_050_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_bin_float_051_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_bin_float_101_is_ok = detail::test_tgamma_for_issue396_value_checker(); + const bool cpp_bin_float_501_is_ok = detail::test_tgamma_for_issue396_value_checker(); - const bool cpp_bin_float_is_ok = cpp_bin_float_020_is_ok - && cpp_bin_float_030_is_ok - && cpp_bin_float_040_is_ok - && cpp_bin_float_049_is_ok - && cpp_bin_float_050_is_ok - && cpp_bin_float_051_is_ok - && cpp_bin_float_101_is_ok - && cpp_bin_float_501_is_ok; + const bool cpp_bin_float_is_ok = ( cpp_bin_float_020_is_ok + && cpp_bin_float_030_is_ok + && cpp_bin_float_040_is_ok + && cpp_bin_float_049_is_ok + && cpp_bin_float_050_is_ok + && cpp_bin_float_051_is_ok + && cpp_bin_float_101_is_ok + && cpp_bin_float_501_is_ok); return cpp_bin_float_is_ok; } @@ -145,18 +149,21 @@ bool test_tgamma_for_issue396_cpp_bin_float() BOOST_AUTO_TEST_CASE(test_cdf____for_issue396_bug_report____tag) { - const bool cpp_bug_report_is_ok = local::test_cdf____for_issue396_bug_report___(); - BOOST_CHECK(cpp_bug_report_is_ok); + const bool cpp_bug_check_is_ok = local::test_cdf____for_issue396_bug_report___(); + + BOOST_CHECK(cpp_bug_check_is_ok); } BOOST_AUTO_TEST_CASE(test_tgamma_for_issue396_cpp_dec_float_tag) { const bool cpp_dec_float_is_ok = local::test_tgamma_for_issue396_cpp_dec_float(); + BOOST_CHECK(cpp_dec_float_is_ok); } BOOST_AUTO_TEST_CASE(test_tgamma_for_issue396_cpp_bin_float_tag) { const bool cpp_bin_float_is_ok = local::test_tgamma_for_issue396_cpp_bin_float(); + BOOST_CHECK(cpp_bin_float_is_ok); }