boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-02-13 00:22:23 +00:00
Code Issues Projects Releases Wiki Activity

Hypergeometric 1F1: Tentatively fix more issues.

Browse Source
This commit is contained in:
jzmaddock
2018-11-07 16:23:35 +00:00
parent cd626ffc89
commit bfb204937e
8 changed files with 484 additions and 224 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
include/boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp
View File

@@ -109,7 +109,7 @@
// function for 13_3_7 evaluation
template <class T, class Policy>
inline T hypergeometric_1F1_13_3_7_series(const T& a, const T& b, const T& z, const Policy& pol, const char* function)
inline T hypergeometric_1F1_13_3_7_series(const T& a, const T& b, const T& z, const Policy& pol, const char* function, int& log_scaling)
{
BOOST_MATH_STD_USING
@@ -156,13 +156,13 @@
}
catch (const std::exception&)
{
// The bessel functions can overflow, use a chacked series in that case:
return hypergeometric_1F1_checked_series_impl(a, b, z, pol);
// The bessel functions can overflow, use a checked series in that case:
return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
}
if (!(boost::math::isfinite)(result))
{
// as above:
return hypergeometric_1F1_checked_series_impl(a, b, z, pol);
return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
}
boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_13_3_7_series<%1%>(%1%,%1%,%1%)", max_iter, pol);

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

@@ -206,6 +206,50 @@
return a < -10; // TODO: make dependent on precision
}
template <class T, class Policy>
inline T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling);
template <class T, class Policy>
T hypergeometric_1F1_fwd_on_b_imp(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling)
{
BOOST_MATH_STD_USING // modf, fabs
int integer_part = 1 + itrunc(ceil(fabs(a) - b));
T bk = b + integer_part;
T first = boost::math::detail::hypergeometric_1F1_imp(a, bk, z, pol, log_scaling);
--bk;
int log_scaling2 = 0;
T second = boost::math::detail::hypergeometric_1F1_imp(a, bk, z, pol, log_scaling2);
if (log_scaling2 != log_scaling)
{
second *= exp(log_scaling2 - log_scaling);
}
//std::cout << first * exp(boost::multiprecision::mpfr_float(log_scaling)) << std::endl;
//std::cout << second * exp(boost::multiprecision::mpfr_float(log_scaling)) << std::endl;
if ((fabs(first) > 1.0e3) || (fabs(first) < 1.0e-3))
{
int rescaling = itrunc(floor(log(fabs(first))));
T scale = exp(-rescaling);
first *= scale;
second *= scale;
log_scaling += rescaling;
}
//std::cout << first * exp(boost::multiprecision::mpfr_float(log_scaling)) << std::endl;
//std::cout << second * exp(boost::multiprecision::mpfr_float(log_scaling)) << std::endl;
boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> s(a, bk, z);
return boost::math::tools::solve_recurrence_relation_backward(s, static_cast<unsigned int>(integer_part), first, second);
}
} } } // namespaces
#endif // BOOST_HYPERGEOMETRIC_1F1_RECURRENCE_HPP_

2
include/boost/math/special_functions/detail/hypergeometric_asym.hpp
View File

@@ -56,7 +56,7 @@
log_scaling -= e;
prefix /= s * exp(t - e);
return prefix * boost::math::detail::hypergeometric_2F0_imp(1 - a, b - a, 1 / z, pol, true);
return prefix * boost::math::detail::hypergeometric_2F0_imp(T(1 - a), T(b - a), T(1 / z), pol, true);
}

30
include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp
View File

@@ -12,7 +12,7 @@
namespace boost { namespace math { namespace detail {
template <class Seq, class Real, class Policy, class Terminal>
std::pair<Real, Real> hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, const Terminal& termination)
std::pair<Real, Real> hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, const Terminal& termination, int& log_scale)
{
using std::abs;
Real result = 1;
@@ -20,6 +20,10 @@
Real term = 1;
Real tol = boost::math::policies::get_epsilon<Real, Policy>();
boost::uintmax_t k = 0;
int log_scaling_factor = boost::math::itrunc(boost::math::tools::log_max_value<Real>()) - 2;
Real scaling_factor = exp(log_scaling_factor);
Real upper_limit(sqrt(boost::math::tools::max_value<Real>()));
Real lower_limit(1 / upper_limit);
while (!termination(k))
{
@@ -47,6 +51,22 @@
term /= k;
result += term;
abs_result += abs(term);
if (fabs(result) >= upper_limit)
{
result /= scaling_factor;
abs_result /= scaling_factor;
term /= scaling_factor;
log_scale += log_scaling_factor;
}
if (fabs(result) < lower_limit)
{
result *= scaling_factor;
abs_result *= scaling_factor;
term *= scaling_factor;
log_scale -= log_scaling_factor;
}
if (abs(result * tol) > abs(term))
break;
if (abs_result * tol > abs(result))
@@ -69,11 +89,11 @@
};
template <class Seq, class Real, class Policy>
Real hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol)
Real hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, int& log_scale)
{
using std::abs;
iteration_terminator term(boost::math::policies::get_max_series_iterations<Policy>());
std::pair<Real, Real> result = hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, term);
std::pair<Real, Real> result = hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, term, log_scale);
//
// Check to see how many digits we've lost, if it's more than half, raise an evaluation error -
// this is an entirely arbitrary cut off, but not unreasonable.
@@ -86,11 +106,11 @@
}
template <class Real, class Policy>
inline Real hypergeometric_1F1_checked_series_impl(const Real& a, const Real& b, const Real& z, const Policy& pol)
inline Real hypergeometric_1F1_checked_series_impl(const Real& a, const Real& b, const Real& z, const Policy& pol, int& log_scale)
{
boost::array<Real, 1> aj = { a };
boost::array<Real, 1> bj = { b };
return hypergeometric_pFq_checked_series_impl(aj, bj, z, pol);
return hypergeometric_pFq_checked_series_impl(aj, bj, z, pol, log_scale);
}
} } } // namespaces

118
include/boost/math/special_functions/detail/hypergeometric_series.hpp
View File

@@ -208,6 +208,31 @@
return detail::sum_pFq_series(s, pol);
}
template <class T, class Policy>
inline T log_pochhammer(T z, unsigned n, const Policy pol, int* s = 0)
{
if (z < 0)
{
if (n < -z)
{
if(s)
*s = (n & 1 ? -1 : 1);
return log_pochhammer(T(-z - n + 1), n, pol);
}
else
{
int cross = itrunc(ceil(-z));
return log_pochhammer(T(-z - cross + 1), cross, pol, s) + log_pochhammer(T(cross + z), n - cross, pol);
}
}
else
{
if(s)
*s = 1;
return lgamma(z + n, pol) - lgamma(z, pol);
}
}
template <class T, class Policy>
inline T hypergeometric_1F1_generic_series(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling, const char* function)
{
@@ -217,6 +242,8 @@
unsigned n = 0;
int log_scaling_factor = boost::math::itrunc(boost::math::tools::log_max_value<T>()) - 2;
T scaling_factor = exp(log_scaling_factor);
T term_m1 = 0;
int local_scaling = 0;
do
{
@@ -226,19 +253,108 @@
sum /= scaling_factor;
term /= scaling_factor;
log_scaling += log_scaling_factor;
local_scaling += log_scaling_factor;
}
if (fabs(sum) < lower_limit)
{
sum *= scaling_factor;
term *= scaling_factor;
log_scaling -= log_scaling_factor;
local_scaling -= log_scaling_factor;
}
term_m1 = term;
term *= (((a + n) / ((b + n) * (n + 1))) * z);
if (n > boost::math::policies::get_max_series_iterations<Policy>())
return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);
++n;
diff = fabs(term / sum);
} while (diff > boost::math::policies::get_epsilon<T, Policy>());
} while ((diff > boost::math::policies::get_epsilon<T, Policy>()) || (fabs(term_m1) < fabs(term)));
if (b + n < 0)
{
if ((a < 0) && (floor(a) == a) && (a > b))
return sum; // b will never cross the origin!
//
// b hasn't crossed the origin yet and the series may spring back into life at that point
// so we need to jump forward to that term and then evaluate forwards and backwards from there:
//
unsigned s = itrunc(-b);
unsigned backstop = n;
int s1, s2;
T t1 = log_pochhammer(a, s, pol, &s1);
T t2 = log_pochhammer(b, s, pol, &s2);
term = t1 - t2;
t1 = lgamma(T(s + 1), pol);
term -= t1;
term += s * log(fabs(z));
if(z < 0)
s1 *= (s & 1 ? -1 : 1);
term -= local_scaling;
if (term > -tools::log_max_value<T>())
{
if (term > 10)
{
int scale = itrunc(floor(term));
term -= scale;
log_scaling += scale;
sum *= exp(T(-scale));
}
term = s1 * s2 * exp(term);
n = s;
T term0 = term;
do
{
sum += term;
if (fabs(sum) >= upper_limit)
{
sum /= scaling_factor;
term /= scaling_factor;
log_scaling += log_scaling_factor;
term0 /= scaling_factor;
}
if (fabs(sum) < lower_limit)
{
sum *= scaling_factor;
term *= scaling_factor;
log_scaling -= log_scaling_factor;
term0 *= scaling_factor;
}
term_m1 = term;
term *= (((a + n) / ((b + n) * (n + 1))) * z);
//if (n > boost::math::policies::get_max_series_iterations<Policy>())
// return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);
++n;
diff = fabs(term / sum);
} while ((diff > boost::math::policies::get_epsilon<T, Policy>()) || (fabs(term) > fabs(term_m1)));
//
// Now go backwards as well:
//
n = s;
term = term0;
do
{
--n;
if (n == backstop)
break;
term_m1 = term;
term /= (((a + n) / ((b + n) * (n + 1))) * z);
sum += term;
if (fabs(sum) >= upper_limit)
{
sum /= scaling_factor;
term /= scaling_factor;
log_scaling += log_scaling_factor;
}
if (fabs(sum) < lower_limit)
{
sum *= scaling_factor;
term *= scaling_factor;
log_scaling -= log_scaling_factor;
}
diff = fabs(term / sum);
} while ((diff > boost::math::policies::get_epsilon<T, Policy>()) || (fabs(term) > fabs(term_m1)));
}
}
return sum;
}

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

@@ -39,7 +39,7 @@ namespace boost { namespace math { namespace detail {
}
template <class T, class Policy>
inline T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling)
T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling)
{
BOOST_MATH_STD_USING // exp, fabs, sqrt
@@ -83,9 +83,11 @@ namespace boost { namespace math { namespace detail {
return exp(z);
if (b - a == b)
return 1;
// asymptotic expansion
// check region
//
// Asymptotic expansion for large z
// TODO: check region for higher precision types.
// Use recurrence relations to move to this region when a and b are also large.
//
if (detail::hypergeometric_1F1_asym_region(a, b, z, pol))
{
return hypergeometric_1F1_asym_large_z_series(a, b, z, pol, log_scaling);
@@ -115,13 +117,19 @@ namespace boost { namespace math { namespace detail {
// Without this we get into an area where the series don't converge if b - a ~ b
return hypergeometric_1F1_backward_recurrence_for_negative_a(a, b, z, pol, function);
}
if ((a > 0) && (b + 1 < a))
{
// Moving to a larger b value will allow us to apply Jummer's relation below:
return hypergeometric_1F1_fwd_on_b_imp(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])
int scaling = itrunc(z);
T r = exp(z - scaling) * detail::hypergeometric_1F1_imp<T>(b_minus_a, b, -z, pol, log_scaling);
log_scaling += scaling;
return exp(z - scaling) * detail::hypergeometric_1F1_imp<T>(b_minus_a, b, -z, pol, log_scaling);
return r;
}
//
// Check for initial divergence:
@@ -168,13 +176,13 @@ namespace boost { namespace math { namespace detail {
if ((a == ceil(a)) && !b_is_negative_and_greater_than_z)
return detail::hypergeometric_1F1_backward_recurrence_for_negative_a(a, b, z, pol, function);
else if ((z > 0) && (2 * (z * (b - (2 * a))) > 0) && (b > -100)) // 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);
return detail::hypergeometric_1F1_13_3_7_series(a, b, z, pol, function, log_scaling);
else if (b < a)
return detail::hypergeometric_1F1_backward_recurrence_for_negative_b(a, b, z, pol);
}
// If we get here, then we've run out of methods to try, use the checked series which will
// raise an error if the result is garbage:
return hypergeometric_1F1_checked_series_impl(a, b, z, pol);
return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
}
return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, function);