From 90779e715f8126dd79782efbcba1c6150c50c166 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sun, 11 Feb 2024 14:21:24 +0100 Subject: [PATCH 1/6] Add several tgamma() edge cases --- test/Jamfile.v2 | 1 + test/test_gamma_edge.cpp | 168 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 169 insertions(+) create mode 100644 test/test_gamma_edge.cpp diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 685bc15ccd..3344396ba9 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -225,6 +225,7 @@ test-suite special_fun : [ run test_expint.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] [ run test_factorials.cpp pch ../../test/build//boost_unit_test_framework ] [ run test_gamma.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] + [ run test_gamma_edge.cpp ] [ run test_gamma_mp.cpp ../../test/build//boost_unit_test_framework : : : release TEST=1 [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] : test_gamma_mp_1 ] [ run test_gamma_mp.cpp ../../test/build//boost_unit_test_framework : : : release TEST=2 [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] : test_gamma_mp_2 ] [ run test_gamma_mp.cpp ../../test/build//boost_unit_test_framework : : : release TEST=3 [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] : test_gamma_mp_3 ] diff --git a/test/test_gamma_edge.cpp b/test/test_gamma_edge.cpp new file mode 100644 index 0000000000..6f96484301 --- /dev/null +++ b/test/test_gamma_edge.cpp @@ -0,0 +1,168 @@ +// Copyright 2024 Christopher Kormanyos +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include +#include +#include + +namespace local +{ + template + 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(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 + // 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::et_off>; + + static_assert( (std::numeric_limits::digits10 >= 248) + && (std::numeric_limits::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(UINT8_C(10))>; + + const local_data_array_type ctrl_data = + {{ + static_cast("9.9999999999999999999999999999999999999999999999999999999999999999999999999999999999942278433509846713939348790991759756895784066406007640119423276511513227322233532906404199270358122304076394840169355222697867950624167125421902178604121309579597843800360126781E+83"), + static_cast("4.9999999999999999999999999999999999999999999999999999999999999999999999999999999999942278433509846713939348790991759756895784066406007640119423276511513227322233532906503104869890919559615934405319418693491786302696988534466221756972734972784545720974832491521E+83"), + static_cast("3.3333333333333333333333333333333333333333333333333333333333333333333333333333333333275611766843180047272682124325093090229117399739340973452756609844846560655566866239935343802757050148488807303802815497619037988103143276843874668674681969322826931482456693779E+83"), + static_cast("2.4999999999999999999999999999999999999999999999999999999999999999999999999999999999942278433509846713939348790991759756895784066406007640119423276511513227322233532906700916068956514070695013535619545635079623006842631352554860913709962299194441475323232733555E+83"), + static_cast("1.9999999999999999999999999999999999999999999999999999999999999999999999999999999999942278433509846713939348790991759756895784066406007640119423276511513227322233532906799821668489311326234553100769609105873541358915452761599180492078575962399389352497160610849E+83"), + static_cast("1.6666666666666666666666666666666666666666666666666666666666666666666666666666666666608945100176513380606015457658426423562450733072674306786089943178179893988900199573565393934688775248440759332586339243334126377654940837310166737113856292271003896337573658994E+83"), + static_cast("1.4285714285714285714285714285714285714285714285714285714285714285714285714285714285656564147795560999653634505277474042610069780691721925833708990797227513036519247192711918581840620123027917945355450333175663777346809865402105363101517574523570821130186163706E+83"), + static_cast("1.2499999999999999999999999999999999999999999999999999999999999999999999999999999999942278433509846713939348790991759756895784066406007640119423276511513227322233532907096538467087703092853171796219799518255296415133916988732139227184416952014232984017855267840E+83"), + static_cast("1.1111111111111111111111111111111111111111111111111111111111111111111111111111111111053389544620957825050459902102870868006895177517118751230534387622624338433344644018306555177731611459503822472480974100160325878317849508887569916664141726330291972302168272984E+83"), + static_cast("9.9999999999999999999999999999999999999999999999999999999999999999999999999999999999422784335098467139393487909917597568957840664060076401194232765115132273222335329072943496661532976039322509265199264598431331192795598068207783839216442784241287383640775600912E+82"), + }}; + + unsigned index = 1U; + + const local_float_type little { "1E-84" }; + const local_float_type my_tol { std::numeric_limits::epsilon() * 256 }; + + for(const auto& ctrl : ctrl_data) + { + const auto x_small = static_cast(static_cast(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; + } + } + } + + auto tgamma_undefined_lanczos_known_error() -> void + { + // This test is intended to hit the lines: + + // template + // 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::et_off>; + + { + local_float_type zero { 0 }; + + bool domain_error_is_ok { false }; + + try + { + boost::math::tgamma(zero); + } + catch(std::domain_error& err) + { + static_cast(err.what()); + + domain_error_is_ok = true; + } + + BOOST_TEST(domain_error_is_ok); + } + + { + local_float_type my_nan = std::numeric_limits::quiet_NaN(); + + bool domain_error_is_ok { false }; + + try + { + boost::math::tgamma(my_nan); + } + catch(std::domain_error& err) + { + static_cast(err.what()); + + domain_error_is_ok = true; + } + + BOOST_TEST(domain_error_is_ok); + } + + { + local_float_type my_inf = -std::numeric_limits::infinity(); + + bool domain_error_is_ok { false }; + + try + { + boost::math::tgamma(my_inf); + } + catch(std::domain_error& err) + { + static_cast(err.what()); + + domain_error_is_ok = true; + } + + BOOST_TEST(domain_error_is_ok); + } + } +} + +auto main() -> int +{ + local::tgamma_under_cbrt_epsilon(); + local::tgamma_undefined_lanczos_known_error(); + + const auto result_is_ok = (boost::report_errors() == 0); + + return (result_is_ok ? 0 : -1); +} From fd917982dfbf9e05d55f4ee2c140fa42ec50cb4c Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sun, 11 Feb 2024 15:20:52 +0100 Subject: [PATCH 2/6] Add a few corresponding lgamma() tests --- test/test_gamma_edge.cpp | 106 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 106 insertions(+) diff --git a/test/test_gamma_edge.cpp b/test/test_gamma_edge.cpp index 6f96484301..1bab92b1e3 100644 --- a/test/test_gamma_edge.cpp +++ b/test/test_gamma_edge.cpp @@ -155,12 +155,118 @@ namespace local 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::et_off>; + + static_assert( (std::numeric_limits::digits10 >= 248) + && (std::numeric_limits::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::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 + // 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::et_off>; + + { + local_float_type zero { 0 }; + + bool domain_error_is_ok { false }; + + try + { + boost::math::lgamma(zero); + } + catch(std::domain_error& err) + { + static_cast(err.what()); + + domain_error_is_ok = true; + } + + BOOST_TEST(domain_error_is_ok); + } + + { + local_float_type my_nan = std::numeric_limits::quiet_NaN(); + + bool domain_error_is_ok { false }; + + try + { + boost::math::lgamma(my_nan); + } + catch(std::domain_error& err) + { + static_cast(err.what()); + + domain_error_is_ok = true; + } + + BOOST_TEST(domain_error_is_ok); + } + + { + local_float_type my_inf = -std::numeric_limits::infinity(); + + bool overflow_error_is_ok { false }; + + try + { + boost::math::lgamma(my_inf); + } + catch(std::overflow_error& err) + { + static_cast(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); From 9aabd9a07214079aa65ebaf9ce33fc2acec7022f Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sun, 11 Feb 2024 17:35:29 +0100 Subject: [PATCH 3/6] Cover the test line(s) themselves --- test/test_gamma_edge.cpp | 56 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 55 insertions(+), 1 deletion(-) diff --git a/test/test_gamma_edge.cpp b/test/test_gamma_edge.cpp index 1bab92b1e3..bb29b6b474 100644 --- a/test/test_gamma_edge.cpp +++ b/test/test_gamma_edge.cpp @@ -6,8 +6,13 @@ #include #include +#include + namespace local { + std::mt19937 eng(static_cast(UINT8_C(42))); + std::uniform_int_distribution dst_one(1, 1); + template auto is_close_fraction(const NumericType& a, const NumericType& b, @@ -82,7 +87,7 @@ namespace local if(!result_tgamma_x_small_is_ok) { - break; + break; // LCOV_EXCL_LINE } } } @@ -99,8 +104,32 @@ namespace local using local_float_type = boost::multiprecision::number, boost::multiprecision::et_off>; { + const local_float_type my_tol { std::numeric_limits::epsilon() * 256 }; + + for(auto index = static_cast(UINT8_C(0)); index < static_cast(UINT8_C(3)); ++index) + { + static_cast(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(UINT8_C(0)); index < static_cast(UINT8_C(3)); ++index) + { + static_cast(index); + local_float_type zero { 0 }; + zero *= dst_one(eng); + bool domain_error_is_ok { false }; try @@ -117,9 +146,14 @@ namespace local BOOST_TEST(domain_error_is_ok); } + for(auto index = static_cast(UINT8_C(0)); index < static_cast(UINT8_C(3)); ++index) { + static_cast(index); + local_float_type my_nan = std::numeric_limits::quiet_NaN(); + my_nan *= dst_one(eng); + bool domain_error_is_ok { false }; try @@ -136,9 +170,14 @@ namespace local BOOST_TEST(domain_error_is_ok); } + for(auto index = static_cast(UINT8_C(0)); index < static_cast(UINT8_C(3)); ++index) { + static_cast(index); + local_float_type my_inf = -std::numeric_limits::infinity(); + my_inf *= dst_one(eng); + bool domain_error_is_ok { false }; try @@ -202,9 +241,14 @@ namespace local using local_float_type = boost::multiprecision::number, boost::multiprecision::et_off>; + for(auto index = static_cast(UINT8_C(0)); index < static_cast(UINT8_C(3)); ++index) { + static_cast(index); + local_float_type zero { 0 }; + zero *= dst_one(eng); + bool domain_error_is_ok { false }; try @@ -221,9 +265,14 @@ namespace local BOOST_TEST(domain_error_is_ok); } + for(auto index = static_cast(UINT8_C(0)); index < static_cast(UINT8_C(3)); ++index) { + static_cast(index); + local_float_type my_nan = std::numeric_limits::quiet_NaN(); + my_nan *= dst_one(eng); + bool domain_error_is_ok { false }; try @@ -240,9 +289,14 @@ namespace local BOOST_TEST(domain_error_is_ok); } + for(auto index = static_cast(UINT8_C(0)); index < static_cast(UINT8_C(3)); ++index) { + static_cast(index); + local_float_type my_inf = -std::numeric_limits::infinity(); + my_inf *= dst_one(eng); + bool overflow_error_is_ok { false }; try From 6ac438c8886df2f555b53b1a6887f82daf7bc179 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sun, 11 Feb 2024 17:48:22 +0100 Subject: [PATCH 4/6] Pick up a few more coverage lines --- include/boost/math/special_functions/gamma.hpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 4af2ebdcfc..f1682b15b7 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1,7 +1,7 @@ // Copyright John Maddock 2006-7, 2013-20. // Copyright Paul A. Bristow 2007, 2013-14. // Copyright Nikhar Agrawal 2013-14 -// Copyright Christopher Kormanyos 2013-14, 2020 +// Copyright Christopher Kormanyos 2013-14, 2020, 2024 // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. (See accompanying file @@ -813,11 +813,13 @@ T full_igamma_prefix(T a, T z, const Policy& pol) { BOOST_MATH_STD_USING - T prefix; if (z > tools::max_value()) return 0; + T alz = a * log(z); + T prefix { }; + if(z >= 1) { if((alz < tools::log_max_value()) && (-z > tools::log_min_value())) @@ -1021,8 +1023,9 @@ inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool // // Compute the full upper fraction (Q) when a is very small: // - T result; - result = boost::math::tgamma1pm1(a, pol); + + T result { boost::math::tgamma1pm1(a, pol); } + if(pgam) *pgam = (result + 1) / a; T p = boost::math::powm1(x, a, pol); From 43c7cd73811ea53c5e21433b126bba4861a530a4 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sun, 11 Feb 2024 17:58:45 +0100 Subject: [PATCH 5/6] Repair typo semicolon in proper place --- include/boost/math/special_functions/gamma.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index f1682b15b7..36e5cc3405 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1024,7 +1024,7 @@ inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool // Compute the full upper fraction (Q) when a is very small: // - T result { boost::math::tgamma1pm1(a, pol); } + T result { boost::math::tgamma1pm1(a, pol) }; if(pgam) *pgam = (result + 1) / a; From bfe1c190ad52ad21d78704bf817c556f88c610bb Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sun, 11 Feb 2024 18:18:19 +0100 Subject: [PATCH 6/6] Update checkout/cache actions CodeCov CI --- .github/workflows/codecov.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/codecov.yml b/.github/workflows/codecov.yml index 78edaff7c0..7978200e4a 100644 --- a/.github/workflows/codecov.yml +++ b/.github/workflows/codecov.yml @@ -88,13 +88,13 @@ jobs: fi git config --global pack.threads 0 - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 with: # For coverage builds fetch the whole history, else only 1 commit using a 'fake ternary' fetch-depth: ${{ matrix.coverage && '0' || '1' }} - name: Cache ccache - uses: actions/cache@v3 + uses: actions/cache@v4 if: env.B2_USE_CCACHE with: path: ~/.ccache @@ -102,7 +102,7 @@ jobs: restore-keys: ${{matrix.os}}-${{matrix.container}}-${{matrix.compiler}}- - name: Fetch Boost.CI - uses: actions/checkout@v3 + uses: actions/checkout@v4 with: repository: boostorg/boost-ci ref: master