diff --git a/include/boost/math/special_functions/hypergeometric_1F1.hpp b/include/boost/math/special_functions/hypergeometric_1F1.hpp index 09d51e1e7..f40699629 100644 --- a/include/boost/math/special_functions/hypergeometric_1F1.hpp +++ b/include/boost/math/special_functions/hypergeometric_1F1.hpp @@ -460,13 +460,28 @@ namespace boost { namespace math { namespace detail { return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling); } } - // Let's otherwise make z positive (almost always) - // by Kummer's transformation - // (we also don't transform if z belongs to [-1,0]) - long long scaling = lltrunc(z); - T r = exp(z - scaling) * detail::hypergeometric_1F1_imp(b_minus_a, b, -z, pol, log_scaling); - log_scaling += scaling; - return r; + if ((b < 0) && (floor(b) == b)) + { + // Negative integer b, so a must be a negative integer too. + // Kummer's transformation fails here! + if(a > -50) + return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, function); + // Is there anything better than this?? + return hypergeometric_1F1_imp(a, float_next(b), z, pol, log_scaling); + } + else + { + // Let's otherwise make z positive (almost always) + // by Kummer's transformation + // (we also don't transform if z belongs to [-1,0]) + // Also note that Kummer's transformation fails when b is + // a negative integer, although this seems to be unmentioned + // in the literature... + long long scaling = lltrunc(z); + T r = exp(z - scaling) * detail::hypergeometric_1F1_imp(b_minus_a, b, -z, pol, log_scaling); + log_scaling += scaling; + return r; + } } // // Check for initial divergence: diff --git a/test/hypergeometric_1f1_neg_int.ipp b/test/hypergeometric_1f1_neg_int.ipp new file mode 100644 index 000000000..3cc6debf2 --- /dev/null +++ b/test/hypergeometric_1f1_neg_int.ipp @@ -0,0 +1,36 @@ +// (C) Copyright John Maddock 2023. +// Use, modification and distribution are subject to 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) + +#ifndef SC_ +# define SC_(x) static_cast(BOOST_JOIN(x, L)) +#endif + static const std::array, 24> hypergeometric_1f1_neg_int = {{ + { SC_(-6.7228000000000000000000000000000000000000e+04), SC_(-1.0000400000000000000000000000000000000000e+05), SC_(-1.9377357482910156250000000000000000000000e+01), SC_(2.2003316788772569594717208669486088054864e-06) }, + { SC_(-6.7228000000000000000000000000000000000000e+04), SC_(-1.0000400000000000000000000000000000000000e+05), SC_(1.9377349853515625000000000000000000000000e+01), SC_(4.5409878747153747404956058641767125880287e+05) }, + { SC_(-6.7228000000000000000000000000000000000000e+04), SC_(-6.7231000000000000000000000000000000000000e+04), SC_(-1.0944412231445312500000000000000000000000e+01), SC_(1.7665024210038748546916559691097019867859e-05) }, + { SC_(-6.7228000000000000000000000000000000000000e+04), SC_(-6.7231000000000000000000000000000000000000e+04), SC_(1.0944412231445312500000000000000000000000e+01), SC_(5.6609031983896583682780819362963230062739e+04) }, + { SC_(-9.6040000000000000000000000000000000000000e+03), SC_(-1.0000400000000000000000000000000000000000e+05), SC_(-2.5397367477416992187500000000000000000000e+00), SC_(7.8355865218528803514965713672531279966815e-01) }, + { SC_(-9.6040000000000000000000000000000000000000e+03), SC_(-1.0000400000000000000000000000000000000000e+05), SC_(2.5397357940673828125000000000000000000000e+00), SC_(1.2762213857631427130211900598472782520337e+00) }, + { SC_(-9.6040000000000000000000000000000000000000e+03), SC_(-9.6070000000000000000000000000000000000000e+03), SC_(-1.9508080482482910156250000000000000000000e+00), SC_(1.4224577228142106290591045359570577017827e-01) }, + { SC_(-9.6040000000000000000000000000000000000000e+03), SC_(-9.6070000000000000000000000000000000000000e+03), SC_(1.9508080482482910156250000000000000000000e+00), SC_(7.0300850442564117333542035244241252172210e+00) }, + { SC_(-1.3720000000000000000000000000000000000000e+03), SC_(-1.0000400000000000000000000000000000000000e+05), SC_(-1.6700172424316406250000000000000000000000e+01), SC_(7.9522031935204164278826294658501467907497e-01) }, + { SC_(-1.3720000000000000000000000000000000000000e+03), SC_(-1.0000400000000000000000000000000000000000e+05), SC_(1.6700164794921875000000000000000000000000e+01), SC_(1.2574655536817972236001926227813569573011e+00) }, + { SC_(-1.3720000000000000000000000000000000000000e+03), SC_(-1.3750000000000000000000000000000000000000e+03), SC_(-6.1633415222167968750000000000000000000000e+00), SC_(2.1336434158773801794817467562036858905312e-03) }, + { SC_(-1.3720000000000000000000000000000000000000e+03), SC_(-1.3750000000000000000000000000000000000000e+03), SC_(6.1633396148681640625000000000000000000000e+00), SC_(4.6865277577094696652027472027299505706366e+02) }, + { SC_(-1.9600000000000000000000000000000000000000e+02), SC_(-1.0000400000000000000000000000000000000000e+05), SC_(-1.8115844726562500000000000000000000000000e+01), SC_(9.6511419727638634464802575293834882470584e-01) }, + { SC_(-1.9600000000000000000000000000000000000000e+02), SC_(-1.0000400000000000000000000000000000000000e+05), SC_(1.8115837097167968750000000000000000000000e+01), SC_(1.0361401464672529487190714874053948234800e+00) }, + { SC_(-1.9600000000000000000000000000000000000000e+02), SC_(-1.9900000000000000000000000000000000000000e+02), SC_(-1.2647186279296875000000000000000000000000e+01), SC_(3.8698894687779787281308700716900907184894e-06) }, + { SC_(-1.9600000000000000000000000000000000000000e+02), SC_(-1.9900000000000000000000000000000000000000e+02), SC_(1.2647182464599609375000000000000000000000e+01), SC_(2.5531745685033098733077958740292498372443e+05) }, + { SC_(-2.8000000000000000000000000000000000000000e+01), SC_(-1.0000400000000000000000000000000000000000e+05), SC_(-2.7095403671264648437500000000000000000000e+00), SC_(9.9924163647108768825421740835291862562688e-01) }, + { SC_(-2.8000000000000000000000000000000000000000e+01), SC_(-1.0000400000000000000000000000000000000000e+05), SC_(2.7095394134521484375000000000000000000000e+00), SC_(1.0007589182485115873525832553163941579882e+00) }, + { SC_(-2.8000000000000000000000000000000000000000e+01), SC_(-3.1000000000000000000000000000000000000000e+01), SC_(-4.4206809997558593750000000000000000000000e+00), SC_(1.7967518222315244889828990302685332933976e-02) }, + { SC_(-2.8000000000000000000000000000000000000000e+01), SC_(-3.1000000000000000000000000000000000000000e+01), SC_(4.4206809997558593750000000000000000000000e+00), SC_(5.2555040655952316061197384365654367193379e+01) }, + { SC_(-4.0000000000000000000000000000000000000000e+00), SC_(-1.0000400000000000000000000000000000000000e+05), SC_(-1.6294479370117187500000000000000000000000e+01), SC_(9.9934840617290030222434494769397556335641e-01) }, + { SC_(-4.0000000000000000000000000000000000000000e+00), SC_(-1.0000400000000000000000000000000000000000e+05), SC_(1.6294471740722656250000000000000000000000e+01), SC_(1.0006519121115562043462443592222315867739e+00) }, + { SC_(-4.0000000000000000000000000000000000000000e+00), SC_(-7.0000000000000000000000000000000000000000e+00), SC_(-1.8267517089843750000000000000000000000000e+01), SC_(5.4688612175408948919092374677372195687245e+01) }, + { SC_(-4.0000000000000000000000000000000000000000e+00), SC_(-7.0000000000000000000000000000000000000000e+00), SC_(1.8267517089843750000000000000000000000000e+01), SC_(3.0779092837466836253322973896270501167497e+02) } + }}; +//#undef SC_ + diff --git a/test/test_1F1.cpp b/test/test_1F1.cpp index e26c90421..f3b226efd 100644 --- a/test/test_1F1.cpp +++ b/test/test_1F1.cpp @@ -154,6 +154,13 @@ void expected_results() largest_type, // test type(s) "Bug cases.*", // test data group ".*", 1500000, 430000); // test function + add_expected_result( + ".*", // compiler + ".*", // stdlib + ".*", // platform + largest_type, // test type(s) + ".*negative.*", // test data group + ".*", 200, 100); // test function // // Finish off by printing out the compiler/stdlib/platform names, diff --git a/test/test_1F1.hpp b/test/test_1F1.hpp index 8892319c4..52f1f1b62 100644 --- a/test/test_1F1.hpp +++ b/test/test_1F1.hpp @@ -413,6 +413,14 @@ void test_spots6(T, const char* type_name) } } +template +void test_spots7(T, const char* type_name) +{ +#include "hypergeometric_1f1_neg_int.ipp" + + do_test_1F1(hypergeometric_1f1_neg_int, type_name, "Both parameters negative integers."); +} + template void test_spots(T z, const char* type_name) { @@ -434,6 +442,7 @@ void test_spots(T z, const char* type_name) // if(std::numeric_limits::digits >= std::numeric_limits::digits && std::numeric_limits::digits <= 128) test_spots6(z, type_name); + test_spots7(z, type_name); } diff --git a/tools/hyp_1f1_neg_int_data.cpp b/tools/hyp_1f1_neg_int_data.cpp new file mode 100644 index 000000000..c2f77028d --- /dev/null +++ b/tools/hyp_1f1_neg_int_data.cpp @@ -0,0 +1,83 @@ +// (C) Copyright John Maddock 2023. +// Use, modification and distribution are subject to 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_MATH_MAX_SERIES_ITERATION_POLICY 10000000 +#define BOOST_MATH_USE_MPFR + +#include "mp_t.hpp" +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace boost::math::tools; +using namespace boost::math; +using namespace std; + +struct hypergeometric_1f1_gen +{ + mp_t operator()(mp_t arg1, mp_t arg2, mp_t z) + { + boost::multiprecision::mpfr_float a1(arg1), a2(arg2), a3(z), r; + r = boost::math::hypergeometric_pFq_precision({ a1 }, { a2 }, a3, 50); + return mp_t(r); + } +}; + + +int main(int, char* []) +{ + parameter_info arg1, arg2, arg3; + test_data data; + + std::cout << "Welcome.\n" + "This program will generate spot tests for 1F1 with negative integer arguments:\n"; + + std::string line; + bool cont; + + random_ns::mt19937 rnd; + random_ns::uniform_real_distribution ur_a(0, 20); + + + for (int i = -4; i > -100000; i *= 7) + { + arg1 = make_single_param(mp_t(i)); + arg2 = make_single_param(mp_t(-100004)); + arg3 = make_single_param(mp_t(ur_a(rnd))); + + data.insert(hypergeometric_1f1_gen(), arg1, arg2, arg3); + arg3.z1 = -arg3.z1; + data.insert(hypergeometric_1f1_gen(), arg1, arg2, arg3); + } + + for (int i = -4; i > -100000; i *= 7) + { + arg1 = make_single_param(mp_t(i)); + arg2 = make_single_param(mp_t(mp_t(i) - 3)); + arg3 = make_single_param(mp_t(ur_a(rnd))); + + data.insert(hypergeometric_1f1_gen(), arg1, arg2, arg3); + arg3.z1 = -arg3.z1; + data.insert(hypergeometric_1f1_gen(), arg1, arg2, arg3); + } + + std::cout << "Enter name of test data file [default=hypergeometric_1f1_neg_int.ipp]"; + std::getline(std::cin, line); + boost::algorithm::trim(line); + if (line == "") + line = "hypergeometric_1f1_neg_int.ipp"; + std::ofstream ofs(line.c_str()); + ofs << std::scientific << std::setprecision(40); + write_code(ofs, data, line.c_str()); + + return 0; +} + +