Skip to content

Commit

Permalink
Update contrib/restricted/boost/math to 1.82.0
Browse files Browse the repository at this point in the history
  • Loading branch information
@robot-piglet
robot-piglet committed May 2, 2023
1 parent 0b0e771 commit 9f42a93
Show file tree
Hide file tree
Showing 18 changed files with 1,634 additions and 335 deletions.
2 changes: 1 addition & 1 deletion contrib/restricted/boost/math/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ The integration routines are usable for functions returning complex results - an

### Quaternions and Octonions

Quaternion and Octonians are class templates similar to std::complex.
Quaternion and Octonion are class templates similar to std::complex.

The full documentation is available on [boost.org](http://www.boost.org/doc/libs/release/libs/math).

Expand Down
1 change: 0 additions & 1 deletion contrib/restricted/boost/math/include/boost/math/constants/calculate_constants.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1011,7 +1011,6 @@ template <class T>
template<int N>
inline T constant_plastic<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>)))
{
using std::cbrt;
using std::sqrt;
return (cbrt(9-sqrt(T(69))) + cbrt(9+sqrt(T(69))))/cbrt(T(18));
}
Expand Down
4 changes: 2 additions & 2 deletions contrib/restricted/boost/math/include/boost/math/cstdfloat/cstdfloat_types.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,14 @@
#ifndef BOOST_MATH_CSTDFLOAT_TYPES_2014_01_09_HPP_
#define BOOST_MATH_CSTDFLOAT_TYPES_2014_01_09_HPP_

#include <float.h>
#include <cfloat>
#include <limits>
#include <boost/math/tools/config.hpp>

// This is the beginning of the preamble.

// In this preamble, the preprocessor is used to query certain
// preprocessor definitions from <float.h>. Based on the results
// preprocessor definitions from <cfloat>. Based on the results
// of these queries, an attempt is made to automatically detect
// the presence of built-in floating-point types having specified
// widths. These are *thought* to be conformant with IEEE-754,
Expand Down
7 changes: 4 additions & 3 deletions ...b/restricted/boost/math/include/boost/math/distributions/detail/common_error_handling.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,11 +185,12 @@ inline bool check_non_centrality(
RealType* result,
const Policy& pol)
{
if((ncp < 0) || !(boost::math::isfinite)(ncp))
{ // Assume scale == 0 is NOT valid for any distribution.
static const RealType upper_limit = static_cast<RealType>((std::numeric_limits<long long>::max)()) - boost::math::policies::get_max_root_iterations<Policy>();
if((ncp < 0) || !(boost::math::isfinite)(ncp) || ncp > upper_limit)
{
*result = policies::raise_domain_error<RealType>(
function,
"Non centrality parameter is %1%, but must be > 0 !", ncp, pol);
"Non centrality parameter is %1%, but must be > 0, and a countable value such that x+1 != x", ncp, pol);
return false;
}
return true;
Expand Down
25 changes: 21 additions & 4 deletions contrib/restricted/boost/math/include/boost/math/special_functions/beta.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -373,7 +373,7 @@ T ibeta_power_terms(T a,
if(a < b)
{
T p1 = pow(b2, b / a);
T l3 = a * (log(b1) + log(p1));
T l3 = (b1 != 0) && (p1 != 0) ? (a * (log(b1) + log(p1))) : tools::max_value<T>(); // arbitrary large value if the logs would fail!
if((l3 < tools::log_max_value<T>())
&& (l3 > tools::log_min_value<T>()))
{
Expand Down Expand Up @@ -570,6 +570,9 @@ T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_deriva
T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));

if (!(boost::math::isfinite)(result))
result = 0;

T l1 = log(cgh / bgh) * (b - 0.5f);
T l2 = log(x * cgh / agh) * a;
//
Expand Down Expand Up @@ -960,14 +963,14 @@ T binomial_ccdf(T n, T k, T x, T y, const Policy& pol)
int start = itrunc(n * x);
if(start <= k + 1)
start = itrunc(k + 2);
result = pow(x, start) * pow(y, n - start) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(start), pol);
result = static_cast<T>(pow(x, start) * pow(y, n - start) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(start), pol));
if(result == 0)
{
// OK, starting slightly above the mode didn't work,
// we'll have to sum the terms the old fashioned way:
for(unsigned i = start - 1; i > k; --i)
{
result += pow(x, (int)i) * pow(y, n - i) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(i), pol);
result += static_cast<T>(pow(x, static_cast<int>(i)) * pow(y, n - i) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(i), pol));
}
}
else
Expand Down Expand Up @@ -1017,6 +1020,13 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de

BOOST_MATH_ASSERT((p_derivative == 0) || normalised);

if(!(boost::math::isfinite)(a))
return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
if(!(boost::math::isfinite)(b))
return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
if(!(boost::math::isfinite)(x))
return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol);

if(p_derivative)
*p_derivative = -1; // value not set.

Expand Down Expand Up @@ -1411,6 +1421,13 @@ T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
//
// start with the usual error checks:
//
if (!(boost::math::isfinite)(a))
return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
if (!(boost::math::isfinite)(b))
return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
if (!(boost::math::isfinite)(x))
return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol);

if(a <= 0)
return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol);
if(b <= 0)
Expand Down Expand Up @@ -1439,7 +1456,7 @@ T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
return f1;
}
//
// Some forwarding functions that dis-ambiguate the third argument type:
// Some forwarding functions that disambiguate the third argument type:
//
template <class RT1, class RT2, class Policy>
inline typename tools::promote_args<RT1, RT2>::type
Expand Down
2 changes: 1 addition & 1 deletion contrib/restricted/boost/math/include/boost/math/special_functions/detail/fp_traits.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -493,7 +493,7 @@ struct fp_traits_non_native<long double, extended_double_precision>
// size_to_precision is a type switch for converting a C++ floating point type
// to the corresponding precision type.

template<int n, bool fp> struct size_to_precision
template<size_t n, bool fp> struct size_to_precision
{
typedef unknown_precision type;
};
Expand Down
82 changes: 69 additions & 13 deletions contrib/restricted/boost/math/include/boost/math/special_functions/detail/ibeta_inverse.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,9 +118,17 @@ T temme_method_1_ibeta_inverse(T a, T b, T z, const Policy& pol)
x = 0.5;
else
x = (1 + eta * sqrt((1 + c) / eta_2)) / 2;

BOOST_MATH_ASSERT(x >= 0);
BOOST_MATH_ASSERT(x <= 1);
//
// These are post-conditions of the method, but the addition above
// may result in us being out by 1ulp either side of the boundary,
// so just check that we're in bounds and adjust as needed.
// See https://github.com/boostorg/math/issues/961
//
if (x < 0)
x = 0;
else if (x > 1)
x = 1;

BOOST_MATH_ASSERT(eta * (x - 0.5) >= 0);
#ifdef BOOST_INSTRUMENT
std::cout << "Estimating x with Temme method 1: " << x << std::endl;
Expand Down Expand Up @@ -368,11 +376,11 @@ T temme_method_3_ibeta_inverse(T a, T b, T p, T q, const Policy& pol)
T e2 = (28 * w_4 + 131 * w_3 + 402 * w_2 + 581 * w + 208) * (w - 1) / (1620 * w1 * w_3);
e2 -= (35 * w_6 - 154 * w_5 - 623 * w_4 - 1636 * w_3 - 3983 * w_2 - 3514 * w - 925) * d / (12960 * w1_2 * w_4);
e2 -= (2132 * w_7 + 7915 * w_6 + 16821 * w_5 + 35066 * w_4 + 87490 * w_3 + 141183 * w_2 + 95993 * w + 21640) * d_2 / (816480 * w_5 * w1_3);
e2 -= (11053 * w_8 + 53308 * w_7 + 117010 * w_6 + 163924 * w_5 + 116188 * w_4 - 258428 * w_3 - 677042 * w_2 - 481940 * w - 105497) * d_3 / (14696640 * w1_4 * w_6);
e2 -= (11053 * w_8 + 53308 * w_7 + 117010 * w_6 + 163924 * w_5 + 116188 * w_4 - 258428 * w_3 - 677042 * w_2 - 481940 * w - 105497) * d_3 / (T(14696640) * w1_4 * w_6);

T e3 = -((3592 * w_7 + 8375 * w_6 - 1323 * w_5 - 29198 * w_4 - 89578 * w_3 - 154413 * w_2 - 116063 * w - 29632) * (w - 1)) / (816480 * w_5 * w1_2);
e3 -= (442043 * w_9 + 2054169 * w_8 + 3803094 * w_7 + 3470754 * w_6 + 2141568 * w_5 - 2393568 * w_4 - 19904934 * w_3 - 34714674 * w_2 - 23128299 * w - 5253353) * d / (146966400 * w_6 * w1_3);
e3 -= (116932 * w_10 + 819281 * w_9 + 2378172 * w_8 + 4341330 * w_7 + 6806004 * w_6 + 10622748 * w_5 + 18739500 * w_4 + 30651894 * w_3 + 30869976 * w_2 + 15431867 * w + 2919016) * d_2 / (146966400 * w1_4 * w_7);
e3 -= (442043 * w_9 + T(2054169) * w_8 + T(3803094) * w_7 + T(3470754) * w_6 + T(2141568) * w_5 - T(2393568) * w_4 - T(19904934) * w_3 - T(34714674) * w_2 - T(23128299) * w - T(5253353)) * d / (T(146966400) * w_6 * w1_3);
e3 -= (116932 * w_10 + 819281 * w_9 + T(2378172) * w_8 + T(4341330) * w_7 + T(6806004) * w_6 + T(10622748) * w_5 + T(18739500) * w_4 + T(30651894) * w_3 + T(30869976) * w_2 + T(15431867) * w + T(2919016)) * d_2 / (T(146966400) * w1_4 * w_7);
//
// Combine eta0 and the error terms to compute eta (Second equation p155):
//
Expand Down Expand Up @@ -631,8 +639,27 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
// Try and compute the easy way first:
//
T bet = 0;
if(b < 2)
bet = boost::math::beta(a, b, pol);
if (b < 2)
{
#ifndef BOOST_NO_EXCEPTIONS
try
#endif
{
bet = boost::math::beta(a, b, pol);

typedef typename Policy::overflow_error_type overflow_type;

BOOST_IF_CONSTEXPR(overflow_type::value != boost::math::policies::throw_on_error)
if(bet > tools::max_value<T>())
bet = tools::max_value<T>();
}
#ifndef BOOST_NO_EXCEPTIONS
catch (const std::overflow_error&)
{
bet = tools::max_value<T>();
}
#endif
}
if(bet != 0)
{
y = pow(b * q * bet, 1/b);
Expand Down Expand Up @@ -673,7 +700,7 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
invert = !invert;
xs = 1 - xs;
}
if (a < tools::min_value<T>())
if ((a < tools::min_value<T>()) && (b > tools::min_value<T>()))
{
if (py)
{
Expand Down Expand Up @@ -702,6 +729,8 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
if (overflow || !(boost::math::isfinite)(bet))
{
xg = exp((boost::math::lgamma(a + 1, pol) + boost::math::lgamma(b, pol) - boost::math::lgamma(a + b, pol) + log(p)) / a);
if (xg > 2 / tools::epsilon<T>())
xg = 2 / tools::epsilon<T>();
}
else
xg = pow(a * p * bet, 1/a);
Expand Down Expand Up @@ -801,15 +830,41 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
}
if(pow(p, 1/a) < 0.5)
{
x = pow(p * a * boost::math::beta(a, b, pol), 1 / a);
#ifndef BOOST_NO_EXCEPTIONS
try
{
#endif
x = pow(p * a * boost::math::beta(a, b, pol), 1 / a);
if ((x > 1) || !(boost::math::isfinite)(x))
x = 1;
#ifndef BOOST_NO_EXCEPTIONS
}
catch (const std::overflow_error&)
{
x = 1;
}
#endif
if(x == 0)
x = boost::math::tools::min_value<T>();
y = 1 - x;
}
else /*if(pow(q, 1/b) < 0.1)*/
{
// model a distorted quarter circle:
y = pow(1 - pow(p, b * boost::math::beta(a, b, pol)), 1/b);
#ifndef BOOST_NO_EXCEPTIONS
try
{
#endif
y = pow(1 - pow(p, b * boost::math::beta(a, b, pol)), 1/b);
if ((y > 1) || !(boost::math::isfinite)(y))
y = 1;
#ifndef BOOST_NO_EXCEPTIONS
}
catch (const std::overflow_error&)
{
y = 1;
}
#endif
if(y == 0)
y = boost::math::tools::min_value<T>();
x = 1 - y;
Expand Down Expand Up @@ -854,6 +909,8 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
if(x < lower)
x = lower;
}
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
std::uintmax_t max_iter_used = 0;
//
// Figure out how many digits to iterate towards:
//
Expand All @@ -876,10 +933,9 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
// Now iterate, we can use either p or q as the target here
// depending on which is smaller:
//
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
x = boost::math::tools::halley_iterate(
boost::math::detail::ibeta_roots<T, Policy>(a, b, (p < q ? p : q), (p < q ? false : true)), x, lower, upper, digits, max_iter);
policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter, pol);
policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter + max_iter_used, pol);
//
// We don't really want these asserts here, but they are useful for sanity
// checking that we have the limits right, uncomment if you suspect bugs *only*.
Expand Down
Loading

0 comments on commit 9f42a93

Please sign in to comment.