From da9d77ef2ef5081cc289f8ea746db63eb5e6c9c2 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sat, 20 Jan 2018 17:48:51 +0000 Subject: [PATCH] Hypergeometric: Begin fixing errors from large params in 1F1 --- .../detail/hypergeometric_1f1_bessel.hpp | 46 ++++++++++- .../detail/hypergeometric_rational.hpp | 2 +- .../special_functions/hypergeometric_1F1.hpp | 45 ++++++++--- .../boost/math/tools/test_data.hpp | 31 +++++++ test/test_1f1.cpp | 8 +- test/test_1f1.hpp | 9 +++ tools/hyp_1f1_data.cpp | 80 ++++++++++++++++--- 7 files changed, 191 insertions(+), 30 deletions(-) diff --git a/include/boost/math/special_functions/detail/hypergeometric_1f1_bessel.hpp b/include/boost/math/special_functions/detail/hypergeometric_1f1_bessel.hpp index b28749112..338386b3e 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_1f1_bessel.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_1f1_bessel.hpp @@ -107,13 +107,40 @@ // function for 13_3_7 evaluation template - inline T hypergeometric_1f1_13_3_7_series(const T& a, const T& b, const T& z, const Policy& pol) + inline T hypergeometric_1f1_13_3_7_series(const T& a, const T& b, const T& z, const Policy& pol, const char* function) { BOOST_MATH_STD_USING const T sqrt_bz_div_2_minus_az = sqrt(((b * z) / 2) - (a * z)); - const T prefix = ((boost::math::tgamma(b, pol) * sqrt_bz_div_2_minus_az) / - pow(sqrt_bz_div_2_minus_az, b)) * exp(z / 2); + + bool use_logs = false; + if ((b > max_factorial::value) && (b * b > tools::log_max_value())) + use_logs = true; + if(z / 2 > tools::log_max_value()) + use_logs = true; + + T prefix; + int sign = 1; + try + { + prefix = boost::math::tgamma(b, pol); + } + catch (std::overflow_error const&) + { + use_logs = true; + } + if(!use_logs && !(boost::math::isfinite)(prefix)) + use_logs = true; + + if (!use_logs) + { + prefix *= (sqrt_bz_div_2_minus_az / pow(sqrt_bz_div_2_minus_az, b)) * exp(z / 2); + } + else + { + prefix = boost::math::lgamma(b, &sign, pol) + (b - 1) * log(sqrt_bz_div_2_minus_az); + prefix += z / 2; + } detail::hypergeometric_1f1_13_3_7_series_term s(a, b, z); boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations(); @@ -124,7 +151,18 @@ T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter); #endif boost::math::policies::check_series_iterations("boost::math::hypergeometric_1f1_13_3_7_series<%1%>(%1%,%1%,%1%)", max_iter, pol); - return prefix * result; + + if (use_logs) + { + result = log(result) + prefix; + if (result > tools::log_max_value()) + return sign * policies::raise_overflow_error(function, 0, pol); + result = sign * exp(result); + } + else + result *= prefix; + + return result; } // term class of Abramowitz & Stegun 13_3_8 formula diff --git a/include/boost/math/special_functions/detail/hypergeometric_rational.hpp b/include/boost/math/special_functions/detail/hypergeometric_rational.hpp index 4b0155683..da2e26a52 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_rational.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_rational.hpp @@ -67,7 +67,7 @@ result = a4 / b4; // condition for interruption - if ((fabs(result) * boost::math::tools::epsilon()) > fabs(result - prev_result)) + if ((fabs(result) * boost::math::tools::epsilon()) > fabs(result - prev_result) / fabs(result)) break; b1 = b2; b2 = b3; b3 = b4; diff --git a/include/boost/math/special_functions/hypergeometric_1F1.hpp b/include/boost/math/special_functions/hypergeometric_1F1.hpp index 87d942b7a..8ef88ad36 100644 --- a/include/boost/math/special_functions/hypergeometric_1F1.hpp +++ b/include/boost/math/special_functions/hypergeometric_1F1.hpp @@ -19,6 +19,7 @@ #include #include #include +#include namespace boost { namespace math { namespace detail { @@ -94,7 +95,7 @@ namespace boost { namespace math { namespace detail { } } - if (fabs(b) >= fabs(100 * z)) // TODO: extend to multuiprecision + if ((fabs(a * z / b) < 3.5) && (fabs(z * 100) < fabs(b))) return detail::hypergeometric_1f1_rational(a, b, z, pol); if (z < -1) @@ -105,7 +106,17 @@ namespace boost { namespace math { namespace detail { // Let's otherwise make z positive (almost always) // by Kummer's transformation // (we also don't transform if z belongs to [-1,0]) - return exp(z) * detail::hypergeometric_1f1_imp(b_minus_a, b, -z, pol); + if (z > -tools::log_max_value() / 4) + return exp(z) * detail::hypergeometric_1f1_imp(b_minus_a, b, -z, pol); + else + { + // + // z is so large that we can't easily apply Kummer's relation, + // use a scaled series when appropriate: + // + if (-z * (b - a) * b > 0) + return hypergeometric_1f1_scaled_series(T(b - a), b, T(-z), pol, function); + } } // // Check for initial divergence: @@ -116,11 +127,25 @@ namespace boost { namespace math { namespace detail { // then it's because both a and b are negative, so check for later // divergence as well: // - if ((a < 0) && (b < 0) && (b > a)) + if (!series_is_divergent && (a < 0) && (b < 0) && (b > a)) { - T n = -floor(b); - series_is_divergent = (a + n) * z / ((b + n) * n) < -1; + // + // We need to exclude situations where we're over the initial "hump" + // in the series terms (ie series has already converged by the time + // b crosses the origin: + // + T fa = fabs(a); + T fb = fabs(b); + T convergence_point = sqrt((a - 1) * (a - b)) - a; + if (-b < convergence_point) + { + T n = -floor(b); + series_is_divergent = (a + n) * z / ((b + n) * n) < -1; + } } + if (series_is_divergent && (b < -1) && (b > -5)) + series_is_divergent = false; // don't bother with divergence, series will be OK + // // Test for alternating series due to negative a, // in particular, see if the series is initially divergent @@ -128,17 +153,17 @@ namespace boost { namespace math { namespace detail { // if (series_is_divergent) { - if((a < 0) && (b > 0)) - return detail::hypergeometric_1f1_backward_recurrence_for_negative_a(a, b, z, pol); + if((a < 0) && (b > -1)) + return detail::hypergeometric_1f1_backward_recurrence_for_negative_a(a, b, z, pol, function); if (detail::hypergeometric_1f1_is_a_small_enough(a)) { // TODO: this part has to be researched deeper const bool b_is_negative_and_greater_than_z = b < 0 ? (fabs(b) > fabs(z) ? 1 : 0) : 0; if ((a == ceil(a)) && !b_is_negative_and_greater_than_z) - return detail::hypergeometric_1f1_backward_recurrence_for_negative_a(a, b, z, pol); - else if ((2 * (z * (b - (2 * a)))) > 0) // TODO: see when this methd is bad in opposite to usual taylor - return detail::hypergeometric_1f1_13_3_7_series(a, b, z, pol); + return detail::hypergeometric_1f1_backward_recurrence_for_negative_a(a, b, z, pol, function); + else if ((z > 0) && (2 * (z * (b - (2 * a)))) > 0) // TODO: see when this methd is bad in opposite to usual taylor + return detail::hypergeometric_1f1_13_3_7_series(a, b, z, pol, function); else if (b < a) return detail::hypergeometric_1f1_backward_recurrence_for_negative_b(a, b, z, pol); } diff --git a/include_private/boost/math/tools/test_data.hpp b/include_private/boost/math/tools/test_data.hpp index 00b17403b..d6effef4a 100644 --- a/include_private/boost/math/tools/test_data.hpp +++ b/include_private/boost/math/tools/test_data.hpp @@ -52,6 +52,8 @@ enum parameter_type random_in_range = 0, periodic_in_range = 1, power_series = 2, + single_value = 3, + plus_minus_value = 4, dummy_param = 0x80 }; @@ -75,6 +77,10 @@ parameter_type& operator |= (parameter_type& a, parameter_type b) // If type == power_series then // n1 and n2 are the endpoints of the exponents (closed range) and z1 is the basis. // +// If type == single_value then z1 contains the single value to add. +// +// If type == plus_minus_value then test at +-z1 +// // If type & dummy_param then this data is ignored and not stored in the output, it // is passed to the generator function however which can do with it as it sees fit. // @@ -107,6 +113,20 @@ inline parameter_info make_power_param(T basis, int start_exponent, int end_e return result; } +template +inline parameter_info make_single_param(T val) +{ + parameter_info result = { single_value, val }; + return result; +} + +template +inline parameter_info make_plus_minus_param(T val) +{ + parameter_info result = { plus_minus_value, val }; + return result; +} + namespace detail{ template @@ -420,6 +440,17 @@ void test_data::create_test_points(std::set& points, const parameter_info< } } break; + case single_value: + { + points.insert(truncate_to_float(real_cast(arg1.z1))); + break; + } + case plus_minus_value: + { + points.insert(truncate_to_float(real_cast(arg1.z1))); + points.insert(truncate_to_float(-real_cast(arg1.z1))); + break; + } default: BOOST_ASSERT(0 == "Invalid parameter_info object"); // Assert will fail if get here. diff --git a/test/test_1f1.cpp b/test/test_1f1.cpp index 2930b4ce7..304a3f98a 100644 --- a/test/test_1f1.cpp +++ b/test/test_1f1.cpp @@ -64,13 +64,13 @@ BOOST_AUTO_TEST_CASE( test_main ) #endif test_spots(0.0, "double"); #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS - test_spots(0.0L, "long double"); + //test_spots(0.0L, "long double"); #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS - test_spots(boost::math::concepts::real_concept(0.1), "real_concept"); + //test_spots(boost::math::concepts::real_concept(0.1), "real_concept"); #endif #endif - test_spots(boost::multiprecision::cpp_bin_float_quad(), "cpp_bin_float_quad"); - test_spots(dec_40(), "dec_40"); + //test_spots(boost::multiprecision::cpp_bin_float_quad(), "cpp_bin_float_quad"); + //test_spots(dec_40(), "dec_40"); } diff --git a/test/test_1f1.hpp b/test/test_1f1.hpp index 10e880b71..07660a9d9 100644 --- a/test/test_1f1.hpp +++ b/test/test_1f1.hpp @@ -68,9 +68,18 @@ void test_spots1(T, const char* type_name) do_test_2F0(hypergeometric_1f1_small_random, type_name, "Small random values"); } +template +void test_spots2(T, const char* type_name) +{ +#include "hypergeometric_1f1_big.ipp" + + do_test_2F0(hypergeometric_1f1_big, type_name, "Large random values"); +} + template void test_spots(T z, const char* type_name) { test_spots1(z, type_name); + test_spots2(z, type_name); } diff --git a/tools/hyp_1f1_data.cpp b/tools/hyp_1f1_data.cpp index b4f7d9d2a..0205b42b3 100644 --- a/tools/hyp_1f1_data.cpp +++ b/tools/hyp_1f1_data.cpp @@ -3,6 +3,8 @@ // 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_MATH_MAX_SERIES_ITERATION_POLICY 10000000 + #include #include #include @@ -27,16 +29,27 @@ struct hypergeometric_1f1_gen } }; -struct hypergeometric_2f0_gen_spec1 +struct hypergeometric_1f1_gen_2 { - boost::math::tuple operator()(mp_t a1, mp_t z) + mp_t operator()(mp_t a1, mp_t a2, mp_t z) { - mp_t result = boost::math::detail::hypergeometric_2f0_generic_series(a1, a1 + 0.5, z, boost::math::policies::policy<>()); - std::cout << a1 << " " << a1 + 0.5 << " " << z << " " << result << std::endl; - return boost::math::make_tuple(a1, a1 + 0.5, z, result); + mp_t result = boost::math::detail::hypergeometric_1f1_generic_series(a1, a2, z, boost::math::policies::policy<>()); + std::cout << a1 << " " << a2 << " " << z << " " << result << std::endl; + if (fabs(result) > (std::numeric_limits::max)()) + { + std::cout << "Discarding result as too large\n"; + throw std::domain_error(""); + } + if (static_cast(result) == 1) + { + std::cout << "Discarding result as unity\n"; + throw std::domain_error(""); // uninteresting result. + } + return result; } }; + int main(int, char* []) { parameter_info arg1, arg2, arg3; @@ -48,12 +61,57 @@ int main(int, char* []) std::string line; bool cont; -#if 0 - arg1 = make_periodic_param(mp_t(-20), mp_t(-1), 19); - arg2 = make_random_param(mp_t(-5), mp_t(5), 8); - arg1.type |= dummy_param; - arg2.type |= dummy_param; - data.insert(hypergeometric_2f0_gen_spec1(), arg1, arg2); +#if 1 + std::vector v; + random_ns::mt19937 rnd; + random_ns::uniform_real_distribution ur_a(0, 1); + + mp_t p = ur_a(rnd); + p *= 1e6; + v.push_back(p); + v.push_back(-p); + p = ur_a(rnd); + p *= 1e5; + v.push_back(p); + v.push_back(-p); + p = ur_a(rnd); + p *= 1e4; + v.push_back(p); + v.push_back(-p); + p = ur_a(rnd); + p *= 1e3; + v.push_back(p); + v.push_back(-p); + p = ur_a(rnd); + p *= 1e2; + v.push_back(p); + v.push_back(-p); + p = ur_a(rnd); + p *= 1e-5; + v.push_back(p); + v.push_back(-p); + p = ur_a(rnd); + p *= 1e-12; + v.push_back(p); + v.push_back(-p); + p = ur_a(rnd); + p *= 1e-30; + v.push_back(p); + v.push_back(-p); + + for (unsigned i = 0; i < v.size(); ++i) + { + for (unsigned j = 0; j < v.size(); ++j) + { + for (unsigned k = 0; k < v.size(); ++k) + { + arg1 = make_single_param(v[i]); + arg2 = make_single_param(v[j] * 3 / 2); + arg3 = make_single_param(v[k] * 5 / 4); + data.insert(hypergeometric_1f1_gen_2(), arg1, arg2, arg3); + } + } + } #else