Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 3 commits into from
May 4, 2023
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
15de9e7
Fix case where b is a negative integer and z is also negative.
jzmaddock May 2, 2023
edb28cc
Add missing file.
jzmaddock May 3, 2023
2fc45e3
Correct #include name.
jzmaddock May 3, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 22 additions & 7 deletions include/boost/math/special_functions/hypergeometric_1F1.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
36 changes: 36 additions & 0 deletions test/hypergeometric_1f1_neg_int.ipp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 changes: 7 additions & 0 deletions test/test_1F1.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
9 changes: 9 additions & 0 deletions test/test_1F1.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand All @@ -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);
}


Expand Down
83 changes: 83 additions & 0 deletions tools/hyp_1f1_neg_int_data.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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;
}