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

Merge pull request #1141 from boostorg/nc_t_improvements

Browse Source 
Prevent passing denormals in calculation.
This commit is contained in:
jzmaddock
2024-06-02 15:25:36 +01:00
committed by GitHub
parent bcd56ee881 20f44d1275
commit 70cdb37759
17 changed files with 12519 additions and 61 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
include/boost/math/concepts/distributions.hpp
View File

@@ -343,7 +343,7 @@ struct DistributionConcept
template <class R, class P>
static void test_extra_members(const boost::math::hypergeometric_distribution<R, P>& d)
{
unsigned u = d.defective();
std::uint64_t u = d.defective();
u = d.sample_count();
u = d.total();
suppress_unused_variable_warning(u);

3
include/boost/math/concepts/std_real_concept.hpp
View File

@@ -90,6 +90,9 @@ public:
std_real_concept& operator=(float c) { m_value = c; return *this; }
std_real_concept& operator=(double c) { m_value = c; return *this; }
std_real_concept& operator=(long double c) { m_value = c; return *this; }
#ifdef BOOST_MATH_USE_FLOAT128
std_real_concept& operator=(BOOST_MATH_FLOAT128_TYPE c) { m_value = c; return *this; }
#endif
// Access:
std_real_concept_base_type value()const{ return m_value; }

28
include/boost/math/distributions/detail/hypergeometric_pdf.hpp
View File

@@ -297,7 +297,7 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
while(data.current_prime <= data.N)
{
std::uint64_t base = data.current_prime;
std::int64_t prime_powers = 0;
std::uint64_t prime_powers = 0;
while(base <= data.N)
{
prime_powers += data.n / base;
@@ -313,7 +313,7 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
}
if(prime_powers)
{
T p = integer_power<T>(static_cast<T>(data.current_prime), prime_powers);
T p = integer_power<T>(static_cast<T>(data.current_prime), static_cast<int>(prime_powers));
if((p > 1) && (tools::max_value<T>() / p < result.value))
{
//
@@ -321,7 +321,7 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
// to sidestep the issue:
//
hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
data.current_prime = prime(++data.prime_index);
data.current_prime = prime(static_cast<unsigned>(++data.prime_index));
return hypergeometric_pdf_prime_loop_imp<T>(data, t);
}
if((p < 1) && (tools::min_value<T>() / p > result.value))
@@ -331,12 +331,12 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
// to sidestep the issue:
//
hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
data.current_prime = prime(++data.prime_index);
data.current_prime = prime(static_cast<unsigned>(++data.prime_index));
return hypergeometric_pdf_prime_loop_imp<T>(data, t);
}
result.value *= p;
}
data.current_prime = prime(++data.prime_index);
data.current_prime = prime(static_cast<unsigned>(++data.prime_index));
}
//
// When we get to here we have run out of prime factors,
@@ -395,18 +395,18 @@ T hypergeometric_pdf_factorial_imp(std::uint64_t x, std::uint64_t r, std::uint64
{
BOOST_MATH_STD_USING
BOOST_MATH_ASSERT(N <= boost::math::max_factorial<T>::value);
T result = boost::math::unchecked_factorial<T>(n);
T result = boost::math::unchecked_factorial<T>(static_cast<unsigned>(n));
T num[3] = {
boost::math::unchecked_factorial<T>(r),
boost::math::unchecked_factorial<T>(N - n),
boost::math::unchecked_factorial<T>(N - r)
boost::math::unchecked_factorial<T>(static_cast<unsigned>(r)),
boost::math::unchecked_factorial<T>(static_cast<unsigned>(N - n)),
boost::math::unchecked_factorial<T>(static_cast<unsigned>(N - r))
};
T denom[5] = {
boost::math::unchecked_factorial<T>(N),
boost::math::unchecked_factorial<T>(x),
boost::math::unchecked_factorial<T>(n - x),
boost::math::unchecked_factorial<T>(r - x),
boost::math::unchecked_factorial<T>(N - n - r + x)
boost::math::unchecked_factorial<T>(static_cast<unsigned>(N)),
boost::math::unchecked_factorial<T>(static_cast<unsigned>(x)),
boost::math::unchecked_factorial<T>(static_cast<unsigned>(n - x)),
boost::math::unchecked_factorial<T>(static_cast<unsigned>(r - x)),
boost::math::unchecked_factorial<T>(static_cast<unsigned>(N - n - r + x))
};
std::size_t i = 0;
std::size_t j = 0;

32
include/boost/math/distributions/detail/hypergeometric_quantile.hpp
View File

@@ -14,7 +14,7 @@
namespace boost{ namespace math{ namespace detail{
template <class T>
inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
{
if((p < cum * fudge_factor) && (x != lbound))
{
@@ -25,7 +25,7 @@ inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned
}
template <class T>
inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned /*lbound*/, unsigned ubound, const policies::discrete_quantile<policies::integer_round_up>&)
inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t /*lbound*/, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_up>&)
{
if((cum < p * fudge_factor) && (x != ubound))
{
@@ -36,7 +36,7 @@ inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned
}
template <class T>
inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
{
if(p >= 0.5)
return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
@@ -44,7 +44,7 @@ inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned
}
template <class T>
inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
{
if(p >= 0.5)
return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
@@ -52,13 +52,13 @@ inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned
}
template <class T>
inline unsigned round_x_from_p(unsigned x, T /*p*/, T /*cum*/, T /*fudge_factor*/, unsigned /*lbound*/, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
inline std::uint64_t round_x_from_p(std::uint64_t x, T /*p*/, T /*cum*/, T /*fudge_factor*/, std::uint64_t /*lbound*/, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
{
return x;
}
template <class T>
inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
{
if((q * fudge_factor > cum) && (x != lbound))
{
@@ -69,7 +69,7 @@ inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned
}
template <class T>
inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned /*lbound*/, unsigned ubound, const policies::discrete_quantile<policies::integer_round_up>&)
inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t /*lbound*/, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_up>&)
{
if((q < cum * fudge_factor) && (x != ubound))
{
@@ -80,7 +80,7 @@ inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned
}
template <class T>
inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
{
if(q < 0.5)
return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
@@ -88,7 +88,7 @@ inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned
}
template <class T>
inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
{
if(q >= 0.5)
return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
@@ -96,13 +96,13 @@ inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned
}
template <class T>
inline unsigned round_x_from_q(unsigned x, T /*q*/, T /*cum*/, T /*fudge_factor*/, unsigned /*lbound*/, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
inline std::uint64_t round_x_from_q(std::uint64_t x, T /*q*/, T /*cum*/, T /*fudge_factor*/, std::uint64_t /*lbound*/, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
{
return x;
}
template <class T, class Policy>
unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned N, const Policy& pol)
std::uint64_t hypergeometric_quantile_imp(T p, T q, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy& pol)
{
#ifdef _MSC_VER
# pragma warning(push)
@@ -113,8 +113,8 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned
BOOST_FPU_EXCEPTION_GUARD
T result;
T fudge_factor = 1 + tools::epsilon<T>() * ((N <= boost::math::prime(boost::math::max_prime - 1)) ? 50 : 2 * N);
unsigned base = static_cast<unsigned>((std::max)(0, static_cast<int>(n + r) - static_cast<int>(N)));
unsigned lim = (std::min)(r, n);
std::uint64_t base = static_cast<std::uint64_t>((std::max)(0, static_cast<int>(n + r) - static_cast<int>(N)));
std::uint64_t lim = (std::min)(r, n);
BOOST_MATH_INSTRUMENT_VARIABLE(p);
BOOST_MATH_INSTRUMENT_VARIABLE(q);
@@ -127,7 +127,7 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned
if(p <= 0.5)
{
unsigned x = base;
std::uint64_t x = base;
result = hypergeometric_pdf<T>(x, r, n, N, pol);
T diff = result;
if (diff == 0)
@@ -175,7 +175,7 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned
}
else
{
unsigned x = lim;
std::uint64_t x = lim;
result = 0;
T diff = hypergeometric_pdf<T>(x, r, n, N, pol);
if (diff == 0)
@@ -225,7 +225,7 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned
}
template <class T, class Policy>
inline unsigned hypergeometric_quantile(T p, T q, unsigned r, unsigned n, unsigned N, const Policy&)
inline std::uint64_t hypergeometric_quantile(T p, T q, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)
{
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T>::type result_type;

21
include/boost/math/distributions/non_central_beta.hpp
View File

@@ -63,7 +63,7 @@ namespace boost
? detail::ibeta_imp(T(a + k), b, x, pol, false, true, &xterm)
: detail::ibeta_imp(b, T(a + k), y, pol, true, true, &xterm);
while (beta * pois == 0)
while (fabs(beta * pois) < tools::min_value<T>())
{
if ((k == 0) || (pois == 0))
return init_val;
@@ -91,7 +91,7 @@ namespace boost
{
T term = beta * pois;
sum += term;
if(((fabs(term/sum) < errtol) && (last_term >= term)) || (term == 0))
if(((fabs(term/sum) < errtol) && (fabs(last_term) >= fabs(term))) || (term == 0))
{
count = k - i;
break;
@@ -106,6 +106,7 @@ namespace boost
last_term = term;
}
last_term = 0;
for(auto i = k + 1; ; ++i)
{
poisf *= l2 / i;
@@ -114,10 +115,11 @@ namespace boost
T term = poisf * betaf;
sum += term;
if((fabs(term/sum) < errtol) || (term == 0))
if(((fabs(term/sum) < errtol) && (fabs(last_term) >= fabs(term))) || (term == 0))
{
break;
}
last_term = term;
if(static_cast<std::uintmax_t>(count + i - k) > max_iter)
{
return policies::raise_evaluation_error("cdf(non_central_beta_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
@@ -554,7 +556,7 @@ namespace boost
ibeta_derivative(a + k, b, x, pol)
: ibeta_derivative(b, a + k, y, pol);
while (beta * pois == 0)
while (fabs(beta * pois) < tools::min_value<T>())
{
if ((k == 0) || (pois == 0))
return 0; // Nothing else we can do!
@@ -579,15 +581,19 @@ namespace boost
// Stable backwards recursion first:
//
std::uintmax_t count = k;
T ratio = 0;
T old_ratio = 0;
for(auto i = k; i >= 0; --i)
{
T term = beta * pois;
sum += term;
if((fabs(term/sum) < errtol) || (term == 0))
ratio = fabs(term / sum);
if(((ratio < errtol) && (ratio < old_ratio)) || (term == 0))
{
count = k - i;
break;
}
ratio = old_ratio;
pois *= i / l2;
if (a + b + i != 1)
@@ -595,6 +601,7 @@ namespace boost
beta *= (a + i - 1) / (x * (a + i + b - 1));
}
}
old_ratio = 0;
for(auto i = k + 1; ; ++i)
{
poisf *= l2 / i;
@@ -602,10 +609,12 @@ namespace boost
T term = poisf * betaf;
sum += term;
if((fabs(term/sum) < errtol) || (term == 0))
ratio = fabs(term / sum);
if(((ratio < errtol) && (ratio < old_ratio)) || (term == 0))
{
break;
}
old_ratio = ratio;
if(static_cast<std::uintmax_t>(count + i - k) > max_iter)
{
return policies::raise_evaluation_error("pdf(non_central_beta_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE

111
include/boost/math/distributions/non_central_t.hpp
View File

@@ -16,6 +16,8 @@
#include <boost/math/distributions/students_t.hpp>
#include <boost/math/distributions/detail/generic_quantile.hpp> // quantile
#include <boost/math/special_functions/trunc.hpp>
#include <boost/math/special_functions/detail/hypergeometric_series.hpp>
#include <boost/math/quadrature/exp_sinh.hpp>
namespace boost
{
@@ -59,7 +61,7 @@ namespace boost
? detail::ibeta_imp(T(k + 1), T(v / 2), x, pol, false, true, &xterm)
: detail::ibeta_imp(T(v / 2), T(k + 1), y, pol, true, true, &xterm);
while (beta * pois == 0)
while (fabs(beta * pois) < tools::min_value<T>())
{
if ((k == 0) || (pois == 0))
return init_val;
@@ -390,6 +392,60 @@ namespace boost
function);
}
template <class T, class Policy>
T non_central_t_pdf_integral(T x, T v, T mu, const Policy& pol)
{
BOOST_MATH_STD_USING
boost::math::quadrature::exp_sinh<T, Policy> integrator;
T integral = pow(v, v / 2) * exp(-v * mu * mu / (2 * (x * x + v)));
if (integral != 0)
{
integral *= integrator.integrate([&x, v, mu](T y)
{
T p;
if (v * log(y) < tools::log_max_value<T>())
p = pow(y, v) * exp(boost::math::pow<2>((y - mu * x / sqrt(x * x + v))) / -2);
else
p = exp(log(y) * v + boost::math::pow<2>((y - mu * x / sqrt(x * x + v))) / -2);
return p;
});
}
integral /= boost::math::constants::root_pi<T>() * boost::math::tgamma(v / 2, pol) * pow(T(2), (v - 1) / 2) * pow(x * x + v, (v + 1) / 2);
return integral;
}
template <class T, class Policy>
T non_central_t_pdf_hypergeometric(T x, T v, T mu, const Policy& pol)
{
BOOST_MATH_STD_USING
long long scale = 0;
const char* function = "non central T PDF";
//
// We only call this routine when we know that the series form of 1F1 is cheap to evaluate,
// so no need to call the whole 1F1 function, just the series will do:
//
T Av = hypergeometric_1F1_generic_series(static_cast<T>((v + 1) / 2), boost::math::constants::half<T>(), static_cast<T>(mu * mu * x * x / (2 * (x * x + v))), pol, scale, function);
Av = ldexp(Av, static_cast<int>(scale));
scale = 0;
T Bv = hypergeometric_1F1_generic_series(static_cast<T>(v / 2 + T(1)), static_cast<T>(T(3) / 2), static_cast<T>(mu * mu * x * x / (2 * (x * x + v))), pol, scale, function);
Bv = ldexp(Bv, static_cast<int>(scale));
Bv *= boost::math::tgamma_delta_ratio(v / 2 + T(1), -constants::half<T>(), pol);
Bv *= boost::math::constants::root_two<T>() * mu * x / sqrt(x * x + v);
T tolerance = tools::root_epsilon<T>() * Av * 4;
Av += Bv;
if (Av < tolerance)
{
// More than half the digits have cancelled, fall back to the integral method:
return non_central_t_pdf_integral(x, v, mu, pol);
}
Av *= exp(-mu * mu / 2) * pow(1 + x * x / v, -(v + 1) / 2) * boost::math::tgamma_delta_ratio(v / 2 + constants::half<T>(), -constants::half<T>(), pol);
Av /= sqrt(v) * boost::math::constants::root_pi<T>();
return Av;
}
template <class T, class Policy>
T non_central_t2_pdf(T n, T delta, T x, T y, const Policy& pol, T init_val)
{
@@ -419,7 +475,7 @@ namespace boost
: ibeta_derivative(n / 2, T(k + 1), y, pol);
BOOST_MATH_INSTRUMENT_VARIABLE(xterm);
while (xterm * pois == 0)
while (fabs(xterm * pois) < tools::min_value<T>())
{
if (k == 0)
return init_val;
@@ -462,15 +518,18 @@ namespace boost
}
}
BOOST_MATH_INSTRUMENT_VARIABLE(sum);
old_ratio = 0;
for(auto i = k + 1; ; ++i)
{
poisf *= d2 / (i + 0.5f);
xtermf *= (x * (n / 2 + i)) / (i);
T term = poisf * xtermf;
sum += term;
if((fabs(term/sum) < errtol) || (term == 0))
T ratio = fabs(term / sum);
if(((ratio < errtol) && (ratio < old_ratio)) || (term == 0))
break;
++count;
old_ratio = ratio;
if(count > max_iter)
{
return policies::raise_evaluation_error("pdf(non_central_t_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE
@@ -489,14 +548,7 @@ namespace boost
normal_distribution<T, Policy> norm(delta, 1);
return pdf(norm, t);
}
//
// Otherwise, for t < 0 we have to use the reflection formula:
if(t < 0)
{
t = -t;
delta = -delta;
}
if(t * t == 0)
if(t * t < tools::epsilon<T>())
{
//
// Handle this as a special case, using the formula
@@ -527,13 +579,30 @@ namespace boost
return pdf(students_t_distribution<T, Policy>(n), t - delta);
}
//
// Figure out if the hypergeometric formula will be efficient or not,
// based on where the summit of the series.
//
T a = (n + 1) / 2;
T x = delta * delta * t * t / (2 * (t * t + n));
T summit = (sqrt(x * (4 * a + x)) + x) / 2;
if (summit < 40)
return non_central_t_pdf_hypergeometric(t, n, delta, pol);
//
// Otherwise, for t < 0 we have to use the reflection formula:
//
if (t < 0)
{
t = -t;
delta = -delta;
}
//
// x and y are the corresponding random
// variables for the noncentral beta distribution,
// with y = 1 - x:
//
T x = t * t / (n + t * t);
x = t * t / (n + t * t);
T y = n / (n + t * t);
T a = 0.5f;
a = 0.5f;
T b = n / 2;
T d2 = delta * delta;
//
@@ -543,12 +612,22 @@ namespace boost
BOOST_MATH_INSTRUMENT_VARIABLE(dt);
T result = non_central_beta_pdf(a, b, d2, x, y, pol);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
T tol = tools::epsilon<T>() * result * 500;
T tol = tools::root_epsilon<T>() * result;
result = non_central_t2_pdf(n, delta, x, y, pol, result);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
if(result <= tol)
result = 0;
result *= dt;
if (result <= tol)
{
// More than half the digits in the result have cancelled,
// Try direct integration... super slow but reliable as far as we can tell...
if (delta < 0)
{
// reflect back:
delta = -delta;
t = -t;
}
result = non_central_t_pdf_integral(t, n, delta, pol);
}
return result;
}

22
include/boost/math/quadrature/detail/exp_sinh_detail.hpp
View File

@@ -198,13 +198,29 @@ auto exp_sinh_detail<Real, Policy>::integrate(const F& f, Real* error, Real* L1,
if (!have_first_j && (I1_last == I1))
{
// No change to the sum, disregard these values on the LHS:
min_abscissa = m_abscissas[1][i];
first_j = i;
if ((i < m_abscissas[1].size() - 1) && (m_abscissas[1][i + 1] > max_abscissa))
{
// The summit is so high, that we found nothing in this row which added to the integral!!
have_first_j = true;
}
else
{
min_abscissa = m_abscissas[1][i];
first_j = i;
}
}
else
have_first_j = true;
}
if (I0 == static_cast<Real>(0))
{
// We failed to find anything, is the integral zero, or have we just not found it yet?
// We'll try one more level, if that still finds nothing then it'll terminate.
min_abscissa = 0;
max_abscissa = boost::math::tools::max_value<Real>();
}
I1 *= half<Real>();
L1_I1 *= half<Real>();
Real err = abs(I0 - I1);
@@ -243,7 +259,7 @@ auto exp_sinh_detail<Real, Policy>::integrate(const F& f, Real* error, Real* L1,
err = abs(I0 - I1);
//std::cout << "Estimate: " << I1 << " Error estimate at level " << i << " = " << err << std::endl;
// Use L1_I1 here to make it work with both complex and real valued integrands:
if (!isfinite(L1_I1))
if (!(boost::math::isfinite)(L1_I1))
{
return static_cast<K>(policies::raise_evaluation_error(function, "The exp_sinh quadrature evaluated your function at a singular point and returned %1%. Please ensure your function evaluates to a finite number over its entire domain.", I1, Policy()));
}

3
include/boost/math/quadrature/exp_sinh.hpp
View File

@@ -91,8 +91,7 @@ auto exp_sinh<Real, Policy>::integrate(const F& f, Real tolerance, Real* error,
static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate";
using std::abs;
if (abs(tolerance) > 1) {
std::string msg = std::string(__FILE__) + ":" + std::to_string(__LINE__) + ":" + std::string(function) + ": The tolerance provided is unusually large; did you confuse it with a domain bound?";
throw std::domain_error(msg);
return policies::raise_domain_error(function, "The tolerance provided (%1%) is unusually large; did you confuse it with a domain bound?", tolerance, Policy());
}
return m_imp->integrate(f, error, L1, function, tolerance, levels);
}

1
include/boost/math/special_functions/jacobi_zeta.hpp
View File

@@ -43,7 +43,6 @@ T jacobi_zeta_imp(T phi, T k, const Policy& pol, T kp)
T result;
T sinp = sin(phi);
T cosp = cos(phi);
T s2 = sinp * sinp;
T c2 = cosp * cosp;
T one_minus_ks2 = kp + c2 - kp * c2;
T k2 = k * k;

2
test/Jamfile.v2
View File

@@ -889,7 +889,7 @@ test-suite new_floats :
test-suite mp :
[ run test_nc_t_quad.cpp pch ../../test/build//boost_unit_test_framework : : : release [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] ]
[ run test_nc_t_quad.cpp pch ../../test/build//boost_unit_test_framework : : : release <toolset>gcc-mingw:<cxxflags>-Wa,-mbig-obj <debug-symbols>off [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] ]
[ run test_polynomial.cpp ../../test/build//boost_unit_test_framework : : : <define>TEST1 : test_polynomial_1 ]
[ run test_polynomial.cpp ../../test/build//boost_unit_test_framework : : : <define>TEST2 : test_polynomial_2 ]
[ run test_polynomial.cpp ../../test/build//boost_unit_test_framework : : : <define>TEST3 : test_polynomial_3 ]

37
test/exp_sinh_quadrature_test.cpp
View File

@@ -573,6 +573,37 @@ void test_complex_exponential_integral_E1(){
}
}
template <class T>
void test_non_central_t()
{
//
// Bug case from the non-central t distribution:
//
using std::pow;
using std::exp;
using std::sqrt;
std::cout << "Testing non-central T PDF integral" << std::endl;
T x = -1.882352352142334;
T v = 77.384613037109375;
T mu = 8.1538467407226562;
T expected = static_cast<T>(4.5098555913703146875364186893655197e+49L);
boost::math::quadrature::exp_sinh<T> integrator;
T err;
T L1;
std::size_t levels;
T integral = integrator.integrate([&x, v, mu](T y)
{
return pow(y, v) * exp(boost::math::pow<2>((y - mu * x / sqrt(x * x + v))) / -2);
},
boost::math::tools::root_epsilon<T>(), &err, &L1, &levels);
T tol = 100 * boost::math::tools::epsilon<T>();
BOOST_CHECK_CLOSE_FRACTION(integral, expected, tol);
}
BOOST_AUTO_TEST_CASE(exp_sinh_quadrature_test)
{
@@ -596,11 +627,13 @@ BOOST_AUTO_TEST_CASE(exp_sinh_quadrature_test)
test_right_limit_infinite<std::float32_t>();
test_nr_examples<std::float32_t>();
test_crc<std::float32_t>();
//test_non_central_t<float32_t>();
#else
test_left_limit_infinite<float>();
test_right_limit_infinite<float>();
test_nr_examples<float>();
test_crc<float>();
//test_non_central_t<float>();
#endif
#endif
@@ -611,11 +644,13 @@ BOOST_AUTO_TEST_CASE(exp_sinh_quadrature_test)
test_right_limit_infinite<std::float64_t>();
test_nr_examples<std::float64_t>();
test_crc<std::float64_t>();
test_non_central_t<std::float64_t>();
#else
test_left_limit_infinite<double>();
test_right_limit_infinite<double>();
test_nr_examples<double>();
test_crc<double>();
test_non_central_t<double>();
#endif
#endif
@@ -626,6 +661,7 @@ BOOST_AUTO_TEST_CASE(exp_sinh_quadrature_test)
test_right_limit_infinite<long double>();
test_nr_examples<long double>();
test_crc<long double>();
test_non_central_t<long double>();
#endif
#endif
#endif
@@ -642,6 +678,7 @@ BOOST_AUTO_TEST_CASE(exp_sinh_quadrature_test)
test_right_limit_infinite<boost::math::concepts::real_concept>();
test_nr_examples<boost::math::concepts::real_concept>();
test_crc<boost::math::concepts::real_concept>();
test_non_central_t<boost::math::concepts::real_concept>();
#endif
#endif
#if defined(TEST6) && defined(BOOST_MATH_RUN_MP_TESTS)

12125
test/nc_t_pdf_data.ipp Normal file
View File

File diff suppressed because it is too large Load Diff

40
test/tanh_sinh_quadrature_test.cpp
View File

@@ -883,6 +883,41 @@ void test_complex()
}
}
template <class T>
void test_non_central_t()
{
//
// Bug case from the non-central t distribution:
//
using std::pow;
using std::exp;
using std::sqrt;
std::cout << "Testing non-central T PDF integral" << std::endl;
T x = -1.882352352142334;
T v = 77.384613037109375;
T mu = 8.1538467407226562;
T expected = static_cast<T>(4.5098555913703146875364186893655197e+49L);
T left = 0;
T right = std::numeric_limits<T>::has_infinity ? std::numeric_limits<T>::infinity() : boost::math::tools::max_value<T>();
boost::math::quadrature::tanh_sinh<T> integrator;
T err;
T L1;
std::size_t levels;
T integral = integrator.integrate([&x, v, mu](T y)
{
return pow(y, v) * exp(boost::math::pow<2>((y - mu * x / sqrt(x * x + v))) / -2);
},
left, right,
boost::math::tools::root_epsilon<T>(), &err, &L1, &levels);
T tol = 100 * boost::math::tools::epsilon<T>();
BOOST_CHECK_CLOSE_FRACTION(integral, expected, tol);
}
BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test)
{
@@ -944,6 +979,7 @@ BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test)
test_early_termination<double>();
test_sf<double>();
test_2_arg<double>();
test_non_central_t<double>();
#endif
#ifdef TEST2A
#ifndef BOOST_MATH_STANDALONE
@@ -967,6 +1003,7 @@ BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test)
test_early_termination<long double>();
test_sf<long double>();
test_2_arg<long double>();
test_non_central_t<long double>();
#endif
#ifdef TEST3A
#ifndef BOOST_MATH_STANDALONE
@@ -991,7 +1028,7 @@ BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test)
test_crc<cpp_bin_float_quad>();
test_sf<cpp_bin_float_quad>();
test_2_arg<cpp_bin_float_quad>();
#endif
#endif
#endif
#ifdef TEST5
#ifdef BOOST_MATH_RUN_MP_TESTS
@@ -1017,6 +1054,7 @@ BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test)
test_early_termination<boost::math::concepts::real_concept>();
test_sf<boost::math::concepts::real_concept>();
test_2_arg<boost::math::concepts::real_concept>();
test_non_central_t<boost::math::concepts::real_concept>();
#endif
#ifdef TEST6A
test_crc<boost::math::concepts::real_concept>();

14
test/test_nc_t.cpp
View File

@@ -113,6 +113,20 @@ void expected_results()
largest_type, // test type(s)
"[^|]*", // test data group
"[^|]*", 400, 100); // test function
add_expected_result(
"[^|]*", // compiler
"[^|]*", // stdlib
"[^|]*", // platform
"double", // test type(s)
"[^|]*PDF", // test data group
"[^|]*", static_cast<std::uintmax_t>(1 / boost::math::tools::root_epsilon<double>()), static_cast<std::uintmax_t>(1 / boost::math::tools::root_epsilon<double>())); // test function
add_expected_result(
"[^|]*", // compiler
"[^|]*", // stdlib
"[^|]*", // platform
"long double", // test type(s)
"[^|]*PDF", // test data group
"[^|]*", static_cast<std::uintmax_t>(1 / boost::math::tools::root_epsilon<long double>()), static_cast<std::uintmax_t>(1 / boost::math::tools::root_epsilon<long double>())); // test function
add_expected_result(
"[^|]*", // compiler
"[^|]*", // stdlib

30
test/test_nc_t.hpp
View File

@@ -353,6 +353,12 @@ T nct_cdf(T df, T nc, T x)
return cdf(boost::math::non_central_t_distribution<T>(df, nc), x);
}
template <class T>
T nct_pdf(T df, T nc, T x)
{
return pdf(boost::math::non_central_t_distribution<T>(df, nc), x);
}
template <class T>
T nct_ccdf(T df, T nc, T x)
{
@@ -399,6 +405,27 @@ void do_test_nc_t(T& data, const char* type_name, const char* test)
#endif
}
template <typename Real, typename T>
void do_test_nc_t_pdf(T& data, const char* type_name, const char* test)
{
typedef Real value_type;
std::cout << "Testing: " << test << std::endl;
value_type(*fp1)(value_type, value_type, value_type) = nct_pdf;
boost::math::tools::test_result<value_type> result;
result = boost::math::tools::test_hetero<Real>(
data,
bind_func<Real>(fp1, 0, 1, 2),
extract_result<Real>(3));
handle_test_result(result, data[result.worst()], result.worst(),
type_name, "non central t PDF", test);
std::cout << std::endl;
}
template <typename Real, typename T>
void quantile_sanity_check(T& data, const char* type_name, const char* test)
{
@@ -522,6 +549,9 @@ void test_accuracy(T, const char* type_name)
#include "nct_asym.ipp"
do_test_nc_t<T>(nct_asym, type_name, "Non Central T (large parameters)");
quantile_sanity_check<T>(nct_asym, type_name, "Non Central T (large parameters)");
#include "nc_t_pdf_data.ipp"
do_test_nc_t_pdf<T>(nc_t_pdf_data, type_name, "Non Central T PDF");
}
}

17
test/test_nc_t_quad.cpp
View File

@@ -40,6 +40,14 @@ using std::endl;
#include <limits>
using std::numeric_limits;
#if defined(__GNUC__) && (__GNUC__ < 13) && (defined(__CYGWIN__) || defined(_WIN32))
//
// We either get an internal compiler error, or worse, the compiler prints nothing at all
// and exits with an error code :(
//
# define DISABLE_THIS_TEST
#endif
void expected_results()
{
@@ -54,6 +62,13 @@ void expected_results()
"cpp_bin_float_quad", // test type(s)
".*large parameters.*", // test data group
"[^|]*", 300000, 100000); // test function
add_expected_result(
"[^|]*", // compiler
"[^|]*", // stdlib
"[^|]*", // platform
"cpp_bin_float_quad", // test type(s)
"[^|]*PDF", // test data group
"[^|]*", static_cast<std::uintmax_t>(1 / boost::math::tools::root_epsilon<boost::multiprecision::cpp_bin_float_quad>()), static_cast<std::uintmax_t>(1 / boost::math::tools::root_epsilon<boost::multiprecision::cpp_bin_float_quad>())); // test function
add_expected_result(
"[^|]*", // compiler
"[^|]*", // stdlib
@@ -76,10 +91,12 @@ BOOST_AUTO_TEST_CASE( test_main )
BOOST_MATH_CONTROL_FP;
// Basic sanity-check spot values.
expected_results();
#ifndef DISABLE_THIS_TEST
#if !BOOST_WORKAROUND(BOOST_MSVC, < 1920)
test_spots(boost::multiprecision::cpp_bin_float_quad(0));
#endif
test_accuracy(boost::multiprecision::cpp_bin_float_quad(0), "cpp_bin_float_quad");
#endif
// double precision tests only:
//test_big_df(boost::multiprecision::cpp_bin_float_quad(0));
} // BOOST_AUTO_TEST_CASE( test_main )

92
tools/non_central_t_pdf_data.cpp Normal file
View File

@@ -0,0 +1,92 @@
// (C) Copyright John Maddock 2024.
// 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/distributions/non_central_t.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 nc_t_pdf_gen
{
mp_t operator()(mp_t v, mp_t mu, mp_t x)
{
//
// Arbitrary smallest value we accept, below this, we will likely underflow to zero
// when computing at double precision:
//
static const mp_t minimum = 1e-270;
mp_t Av = boost::math::hypergeometric_1F1((v + 1) / 2, boost::math::constants::half<mp_t>(), mu * mu * x * x / (2 * (x * x + v)));
mp_t Bv = boost::math::hypergeometric_1F1(v / 2 + mp_t(1), mp_t(3) / 2, mu * mu * x * x / (2 * (x * x + v)));
Bv *= boost::math::tgamma_ratio(v / 2 + mp_t(1), (v + mp_t(1)) / 2);
Bv *= boost::math::constants::root_two<mp_t>() * mu * x / sqrt(x * x + v);
Av += Bv;
Av *= exp(-mu * mu / 2) * pow(1 + x * x / v, -(v + 1) / 2) * boost::math::tgamma_ratio((v + mp_t(1)) / 2, v / 2);
Av /= sqrt(v) * boost::math::constants::root_pi<mp_t>();
if (Av < minimum)
throw std::domain_error("Result too small!");
return Av;
}
};
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 the non central t PDF:\n";
std::string line;
bool cont;
random_ns::mt19937 rnd;
random_ns::uniform_real_distribution<float> ur_a(0, 20);
do
{
if (0 == get_user_parameter_info(arg1, "v"))
return 1;
if (0 == get_user_parameter_info(arg2, "nc"))
return 1;
arg3 = make_periodic_param(mp_t(-20), mp_t(200), 170);
data.insert(nc_t_pdf_gen(), arg1, arg2, arg3);
std::cout << "Any more data [y/n]?";
std::getline(std::cin, line);
boost::algorithm::trim(line);
cont = (line == "y");
} while (cont);
std::cout << "Enter name of test data file [default=nc_t_pdf_data.ipp]";
std::getline(std::cin, line);
boost::algorithm::trim(line);
if (line == "")
line = "nc_t_pdf_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific << std::setprecision(40);
write_code(ofs, data, line.c_str());
return 0;
}