boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Fix case where b is a negative integer and z is also negative. (#983)

Browse Source 
* Fix case where b is a negative integer and z is also negative.
Add tests etc.
Fixes: https://github.com/boostorg/math/issues/982.
This commit is contained in:
jzmaddock
2023-05-04 09:17:16 +01:00
committed by GitHub
parent 298a243ccd
commit d5960de3db
5 changed files with 157 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files

29
include/boost/math/special_functions/hypergeometric_1F1.hpp
View File

@@ -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<T>(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<T>(b_minus_a, b, -z, pol, log_scaling);
log_scaling += scaling;
return r;
}
}
//
// Check for initial divergence:

36
test/hypergeometric_1f1_neg_int.ipp Normal file
View File

@@ -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<T>(BOOST_JOIN(x, L))
#endif
static const std::array<std::array<T, 4>, 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_

7
test/test_1F1.cpp
View File

@@ -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,

9
test/test_1F1.hpp
View File

@@ -413,6 +413,14 @@ void test_spots6(T, const char* type_name)
}
}
template <class T>
void test_spots7(T, const char* type_name)
{
#include "hypergeometric_1f1_neg_int.ipp"
do_test_1F1<T>(hypergeometric_1f1_neg_int, type_name, "Both parameters negative integers.");
}
template <class T>
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<T>::digits >= std::numeric_limits<double>::digits && std::numeric_limits<T>::digits <= 128)
test_spots6(z, type_name);
test_spots7(z, type_name);
}

83
tools/hyp_1f1_neg_int_data.cpp Normal file
View File

@@ -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 <boost/math/special_functions/hypergeometric_1F1.hpp>
#include <boost/math/special_functions/hypergeometric_pFq.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/lexical_cast.hpp>
#include <fstream>
#include <map>
#include <boost/math/tools/test_data.hpp>
#include <boost/random.hpp>
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<mp_t> arg1, arg2, arg3;
test_data<mp_t> 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<float> 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;
}