boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-26 06:42:12 +00:00
Code Issues Projects Releases Wiki Activity

1F1: Apply backwards recurrence relations for GammaP in large a,b,z approximation.

Browse Source 
[CI SKIP]
This commit is contained in:
jzmaddock
2019-01-31 19:06:17 +00:00
parent 30af3aebd2
commit a187b714e9
3 changed files with 27 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
include/boost/math/special_functions/detail/hypergeometric_1F1_large_abz.hpp
View File

@@ -19,26 +19,46 @@
template <class T, class Policy>
struct hypergeometric_1F1_igamma_series
{
enum{ cache_size = 64 };
typedef T result_type;
hypergeometric_1F1_igamma_series(const T& alpha, const T& delta, const T& x, const Policy& pol)
: delta_poch(-delta), alpha_poch(alpha), x(x), k(0), pol(pol)
: delta_poch(-delta), alpha_poch(alpha), x(x), k(0), cache_offset(0), pol(pol)
{
BOOST_MATH_STD_USING
T log_term = log(x) * -alpha;
log_scaling = itrunc(log_term - 3 - boost::math::tools::log_min_value<T>() / 50);
term = exp(log_term - log_scaling);
refill_cache();
}
T operator()()
{
T result = term * boost::math::gamma_p(alpha_poch, x, pol);
if (k - cache_offset >= cache_size)
{
cache_offset += cache_size;
refill_cache();
}
T result = term * gamma_cache[k - cache_offset];
term *= delta_poch * alpha_poch / (++k * x);
delta_poch += 1;
alpha_poch += 1;
return result;
}
void refill_cache()
{
typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
gamma_cache[cache_size - 1] = boost::math::gamma_p(alpha_poch + cache_size - 1, x, pol);
for (int i = cache_size - 1; i > 0; --i)
{
gamma_cache[i - 1] = gamma_cache[i] >= 1 ? 1 : gamma_cache[i] + regularised_gamma_prefix(alpha_poch + i - 1, x, pol, lanczos_type()) / (alpha_poch + i - 1);
}
}
T delta_poch, alpha_poch, x, term;
T gamma_cache[cache_size];
int k;
int log_scaling;
int cache_offset;
Policy pol;
};

1
include/boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp
View File

@@ -237,7 +237,6 @@
BOOST_ASSERT(leading_a_shift > 1);
BOOST_ASSERT(a_b_shift > 1);
boost::uintmax_t max_iter = 200;
T first, second;
first = boost::math::detail::hypergeometric_1F1_imp(T(a + a_shift), T(b + b_shift), z, pol, log_scaling);
//

6
test/test_1F1.hpp
View File

@@ -25,6 +25,10 @@
#include <boost/math/special_functions/hypergeometric_1F1.hpp>
#include <boost/math/quadrature/exp_sinh.hpp>
#ifdef BOOST_MSVC
#pragma warning(disable:4127)
#endif
template <class Real, class T>
void do_test_1F1(const T& data, const char* type_name, const char* test_name)
{
@@ -96,7 +100,7 @@ void test_spots4(T, const char* type_name)
template <class T>
void test_spots5(T, const char* type_name)
{
std::cout << "Testing special cases" << std::endl;
std::cout << "Testing special cases for type " << type_name << std::endl;
//
// Special cases:
//