Skip to content

Commit

Permalink
Merge pull request #1141 from boostorg/nc_t_improvements
Browse files Browse the repository at this point in the history 
Prevent passing denormals in calculation.
  • Loading branch information
@jzmaddock
jzmaddock authored Jun 2, 2024
2 parents bcd56ee + 20f44d1 commit 70cdb37
Show file tree
Hide file tree
Showing 17 changed files with 12,519 additions and 61 deletions.
2 changes: 1 addition & 1 deletion include/boost/math/concepts/distributions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
3 changes: 3 additions & 0 deletions include/boost/math/concepts/std_real_concept.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,9 @@ class std_real_concept
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; }
Expand Down
28 changes: 14 additions & 14 deletions include/boost/math/distributions/detail/hypergeometric_pdf.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -313,15 +313,15 @@ 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))
{
//
// The next calculation would overflow, use recursion
// 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))
Expand All @@ -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,
Expand Down Expand Up @@ -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;
Expand Down
32 changes: 16 additions & 16 deletions include/boost/math/distributions/detail/hypergeometric_quantile.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
{
Expand All @@ -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))
{
Expand All @@ -36,29 +36,29 @@ 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>());
return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
}

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>());
return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
}

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))
{
Expand All @@ -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))
{
Expand All @@ -80,29 +80,29 @@ 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>());
return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
}

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>());
return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
}

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)
Expand All @@ -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);
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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;
Expand Down
21 changes: 15 additions & 6 deletions include/boost/math/distributions/non_central_beta.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -106,6 +106,7 @@ namespace boost

last_term = term;
}
last_term = 0;
for(auto i = k + 1; ; ++i)
{
poisf *= l2 / i;
Expand All @@ -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
Expand Down Expand Up @@ -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!
Expand All @@ -579,33 +581,40 @@ 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)
{
beta *= (a + i - 1) / (x * (a + i + b - 1));
}
}
old_ratio = 0;
for(auto i = k + 1; ; ++i)
{
poisf *= l2 / i;
betaf *= x * (a + b + i - 1) / (a + i - 1);

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
Expand Down
Loading

0 comments on commit 70cdb37

Please sign in to comment.