Skip to content

Commit 8adfc2f

Browse files
robot-pigletrobot-piglet
committed
Update contrib/restricted/boost/math to 1.86.0
26ffe6aa485d878318cbc7d5a6d850e6c1121dd1
1 parent 9424b3c commit 8adfc2f

File tree

66 files changed

+1520
-641
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

66 files changed

+1520
-641
lines changed

‎contrib/restricted/boost/math/include/boost/math/distributions/detail/hypergeometric_pdf.hpp

Lines changed: 14 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -297,7 +297,7 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
297297
while(data.current_prime <= data.N)
298298
{
299299
std::uint64_t base = data.current_prime;
300-
std::int64_t prime_powers = 0;
300+
std::uint64_t prime_powers = 0;
301301
while(base <= data.N)
302302
{
303303
prime_powers += data.n / base;
@@ -313,15 +313,15 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
313313
}
314314
if(prime_powers)
315315
{
316-
T p = integer_power<T>(static_cast<T>(data.current_prime), prime_powers);
316+
T p = integer_power<T>(static_cast<T>(data.current_prime), static_cast<int>(prime_powers));
317317
if((p > 1) && (tools::max_value<T>() / p < result.value))
318318
{
319319
//
320320
// The next calculation would overflow, use recursion
321321
// to sidestep the issue:
322322
//
323323
hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
324-
data.current_prime = prime(++data.prime_index);
324+
data.current_prime = prime(static_cast<unsigned>(++data.prime_index));
325325
return hypergeometric_pdf_prime_loop_imp<T>(data, t);
326326
}
327327
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
331331
// to sidestep the issue:
332332
//
333333
hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
334-
data.current_prime = prime(++data.prime_index);
334+
data.current_prime = prime(static_cast<unsigned>(++data.prime_index));
335335
return hypergeometric_pdf_prime_loop_imp<T>(data, t);
336336
}
337337
result.value *= p;
338338
}
339-
data.current_prime = prime(++data.prime_index);
339+
data.current_prime = prime(static_cast<unsigned>(++data.prime_index));
340340
}
341341
//
342342
// 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
395395
{
396396
BOOST_MATH_STD_USING
397397
BOOST_MATH_ASSERT(N <= boost::math::max_factorial<T>::value);
398-
T result = boost::math::unchecked_factorial<T>(n);
398+
T result = boost::math::unchecked_factorial<T>(static_cast<unsigned>(n));
399399
T num[3] = {
400-
boost::math::unchecked_factorial<T>(r),
401-
boost::math::unchecked_factorial<T>(N - n),
402-
boost::math::unchecked_factorial<T>(N - r)
400+
boost::math::unchecked_factorial<T>(static_cast<unsigned>(r)),
401+
boost::math::unchecked_factorial<T>(static_cast<unsigned>(N - n)),
402+
boost::math::unchecked_factorial<T>(static_cast<unsigned>(N - r))
403403
};
404404
T denom[5] = {
405-
boost::math::unchecked_factorial<T>(N),
406-
boost::math::unchecked_factorial<T>(x),
407-
boost::math::unchecked_factorial<T>(n - x),
408-
boost::math::unchecked_factorial<T>(r - x),
409-
boost::math::unchecked_factorial<T>(N - n - r + x)
405+
boost::math::unchecked_factorial<T>(static_cast<unsigned>(N)),
406+
boost::math::unchecked_factorial<T>(static_cast<unsigned>(x)),
407+
boost::math::unchecked_factorial<T>(static_cast<unsigned>(n - x)),
408+
boost::math::unchecked_factorial<T>(static_cast<unsigned>(r - x)),
409+
boost::math::unchecked_factorial<T>(static_cast<unsigned>(N - n - r + x))
410410
};
411411
std::size_t i = 0;
412412
std::size_t j = 0;

‎contrib/restricted/boost/math/include/boost/math/distributions/detail/hypergeometric_quantile.hpp

Lines changed: 16 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@
1414
namespace boost{ namespace math{ namespace detail{
1515

1616
template <class T>
17-
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>&)
17+
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>&)
1818
{
1919
if((p < cum * fudge_factor) && (x != lbound))
2020
{
@@ -25,7 +25,7 @@ inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned
2525
}
2626

2727
template <class T>
28-
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>&)
28+
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>&)
2929
{
3030
if((cum < p * fudge_factor) && (x != ubound))
3131
{
@@ -36,29 +36,29 @@ inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned
3636
}
3737

3838
template <class T>
39-
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>&)
39+
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>&)
4040
{
4141
if(p >= 0.5)
4242
return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
4343
return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
4444
}
4545

4646
template <class T>
47-
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>&)
47+
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>&)
4848
{
4949
if(p >= 0.5)
5050
return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
5151
return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
5252
}
5353

5454
template <class T>
55-
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>&)
55+
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>&)
5656
{
5757
return x;
5858
}
5959

6060
template <class T>
61-
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>&)
61+
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>&)
6262
{
6363
if((q * fudge_factor > cum) && (x != lbound))
6464
{
@@ -69,7 +69,7 @@ inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned
6969
}
7070

7171
template <class T>
72-
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>&)
72+
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>&)
7373
{
7474
if((q < cum * fudge_factor) && (x != ubound))
7575
{
@@ -80,29 +80,29 @@ inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned
8080
}
8181

8282
template <class T>
83-
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>&)
83+
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>&)
8484
{
8585
if(q < 0.5)
8686
return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
8787
return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
8888
}
8989

9090
template <class T>
91-
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>&)
91+
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>&)
9292
{
9393
if(q >= 0.5)
9494
return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
9595
return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
9696
}
9797

9898
template <class T>
99-
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>&)
99+
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>&)
100100
{
101101
return x;
102102
}
103103

104104
template <class T, class Policy>
105-
unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned N, const Policy& pol)
105+
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)
106106
{
107107
#ifdef _MSC_VER
108108
# pragma warning(push)
@@ -113,8 +113,8 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned
113113
BOOST_FPU_EXCEPTION_GUARD
114114
T result;
115115
T fudge_factor = 1 + tools::epsilon<T>() * ((N <= boost::math::prime(boost::math::max_prime - 1)) ? 50 : 2 * N);
116-
unsigned base = static_cast<unsigned>((std::max)(0, static_cast<int>(n + r) - static_cast<int>(N)));
117-
unsigned lim = (std::min)(r, n);
116+
std::uint64_t base = static_cast<std::uint64_t>((std::max)(0, static_cast<int>(n + r) - static_cast<int>(N)));
117+
std::uint64_t lim = (std::min)(r, n);
118118

119119
BOOST_MATH_INSTRUMENT_VARIABLE(p);
120120
BOOST_MATH_INSTRUMENT_VARIABLE(q);
@@ -127,7 +127,7 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned
127127

128128
if(p <= 0.5)
129129
{
130-
unsigned x = base;
130+
std::uint64_t x = base;
131131
result = hypergeometric_pdf<T>(x, r, n, N, pol);
132132
T diff = result;
133133
if (diff == 0)
@@ -175,7 +175,7 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned
175175
}
176176
else
177177
{
178-
unsigned x = lim;
178+
std::uint64_t x = lim;
179179
result = 0;
180180
T diff = hypergeometric_pdf<T>(x, r, n, N, pol);
181181
if (diff == 0)
@@ -225,7 +225,7 @@ unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned
225225
}
226226

227227
template <class T, class Policy>
228-
inline unsigned hypergeometric_quantile(T p, T q, unsigned r, unsigned n, unsigned N, const Policy&)
228+
inline std::uint64_t hypergeometric_quantile(T p, T q, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)
229229
{
230230
BOOST_FPU_EXCEPTION_GUARD
231231
typedef typename tools::promote_args<T>::type result_type;

‎contrib/restricted/boost/math/include/boost/math/distributions/geometric.hpp

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -384,6 +384,8 @@ namespace boost
384384
inline RealType logcdf(const geometric_distribution<RealType, Policy>& dist, const RealType& k)
385385
{ // Cumulative Distribution Function of geometric.
386386
using std::pow;
387+
using std::log;
388+
using std::exp;
387389
static const char* function = "boost::math::logcdf(const geometric_distribution<%1%>&, %1%)";
388390

389391
// k argument may be integral, signed, or unsigned, or floating point.

‎contrib/restricted/boost/math/include/boost/math/distributions/non_central_beta.hpp

Lines changed: 42 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,17 @@ namespace boost
6363
? detail::ibeta_imp(T(a + k), b, x, pol, false, true, &xterm)
6464
: detail::ibeta_imp(b, T(a + k), y, pol, true, true, &xterm);
6565

66+
while (fabs(beta * pois) < tools::min_value<T>())
67+
{
68+
if ((k == 0) || (pois == 0))
69+
return init_val;
70+
k /= 2;
71+
pois = gamma_p_derivative(T(k + 1), l2, pol);
72+
beta = x < y
73+
? detail::ibeta_imp(T(a + k), b, x, pol, false, true, &xterm)
74+
: detail::ibeta_imp(b, T(a + k), y, pol, true, true, &xterm);
75+
}
76+
6677
xterm *= y / (a + b + k - 1);
6778
T poisf(pois), betaf(beta), xtermf(xterm);
6879
T sum = init_val;
@@ -80,7 +91,7 @@ namespace boost
8091
{
8192
T term = beta * pois;
8293
sum += term;
83-
if(((fabs(term/sum) < errtol) && (last_term >= term)) || (term == 0))
94+
if(((fabs(term/sum) < errtol) && (fabs(last_term) >= fabs(term))) || (term == 0))
8495
{
8596
count = k - i;
8697
break;
@@ -95,6 +106,7 @@ namespace boost
95106

96107
last_term = term;
97108
}
109+
last_term = 0;
98110
for(auto i = k + 1; ; ++i)
99111
{
100112
poisf *= l2 / i;
@@ -103,10 +115,11 @@ namespace boost
103115

104116
T term = poisf * betaf;
105117
sum += term;
106-
if((fabs(term/sum) < errtol) || (term == 0))
118+
if(((fabs(term/sum) < errtol) && (fabs(last_term) >= fabs(term))) || (term == 0))
107119
{
108120
break;
109121
}
122+
last_term = term;
110123
if(static_cast<std::uintmax_t>(count + i - k) > max_iter)
111124
{
112125
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
@@ -542,6 +555,24 @@ namespace boost
542555
T beta = x < y ?
543556
ibeta_derivative(a + k, b, x, pol)
544557
: ibeta_derivative(b, a + k, y, pol);
558+
559+
while (fabs(beta * pois) < tools::min_value<T>())
560+
{
561+
if ((k == 0) || (pois == 0))
562+
return 0; // Nothing else we can do!
563+
//
564+
// We only get here when a+k and b are large and x is small,
565+
// in that case reduce k (bisect) until both terms are finite:
566+
//
567+
k /= 2;
568+
pois = gamma_p_derivative(T(k + 1), l2, pol);
569+
// Starting beta term:
570+
beta = x < y ?
571+
ibeta_derivative(a + k, b, x, pol)
572+
: ibeta_derivative(b, a + k, y, pol);
573+
}
574+
575+
545576
T sum = 0;
546577
T poisf(pois);
547578
T betaf(beta);
@@ -550,33 +581,40 @@ namespace boost
550581
// Stable backwards recursion first:
551582
//
552583
std::uintmax_t count = k;
584+
T ratio = 0;
585+
T old_ratio = 0;
553586
for(auto i = k; i >= 0; --i)
554587
{
555588
T term = beta * pois;
556589
sum += term;
557-
if((fabs(term/sum) < errtol) || (term == 0))
590+
ratio = fabs(term / sum);
591+
if(((ratio < errtol) && (ratio < old_ratio)) || (term == 0))
558592
{
559593
count = k - i;
560594
break;
561595
}
596+
ratio = old_ratio;
562597
pois *= i / l2;
563598

564599
if (a + b + i != 1)
565600
{
566601
beta *= (a + i - 1) / (x * (a + i + b - 1));
567602
}
568603
}
604+
old_ratio = 0;
569605
for(auto i = k + 1; ; ++i)
570606
{
571607
poisf *= l2 / i;
572608
betaf *= x * (a + b + i - 1) / (a + i - 1);
573609

574610
T term = poisf * betaf;
575611
sum += term;
576-
if((fabs(term/sum) < errtol) || (term == 0))
612+
ratio = fabs(term / sum);
613+
if(((ratio < errtol) && (ratio < old_ratio)) || (term == 0))
577614
{
578615
break;
579616
}
617+
old_ratio = ratio;
580618
if(static_cast<std::uintmax_t>(count + i - k) > max_iter)
581619
{
582620
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

0 commit comments

Comments
 (0)