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

1F1: Allow log_pochhammer to work when the argument is a negative integer as long as it doesn't cross the origin. Remove special series code for b < 0 as this now always goes to the checked implementation.

Browse Source
This commit is contained in:
jzmaddock
2019-04-14 09:50:35 +01:00
parent b265ce14cf
commit f520bc9dff
1 changed files with 7 additions and 86 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -231,6 +231,13 @@
else
#endif
{
if (z + n < 0)
{
T r = log_pochhammer(T(-z - n + 1), n, pol, s);
if (s)
*s *= (n & 1 ? -1 : 1);
return r;
}
int s1, s2;
T r = boost::math::lgamma(T(z + n), &s1, pol) - boost::math::lgamma(z, &s2, pol);
if(s)
@@ -286,92 +293,6 @@
diff = fabs(term / sum);
} while ((diff > boost::math::policies::get_epsilon<T, Policy>()) || (fabs(term_m1) < fabs(term)) || (small_a && n < 10));
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;
}