Skip to content

Commit 529f3a7

Browse files
mborlandmborland
authored
Merge pull request #1246 from boostorg/develop
Merge for 1.88
2 parents 0b574c3 + a5c0625 commit 529f3a7

24 files changed

+106
-86
lines changed

Diff for: ‎doc/background/special_tut.qbk

+1-1
Original file line numberDiff line numberDiff line change
@@ -373,7 +373,7 @@ Now we just need to write the test driver program, at it's most basic it looks s
373373
std::cout << "<note>The long double tests have been disabled on this platform "
374374
"either because the long double overloads of the usual math functions are "
375375
"not available at all, or because they are too inaccurate for these tests "
376-
"to pass.</note>" << std::cout;
376+
"to pass.</note>" << std::endl;
377377
#endif
378378
}
379379

Diff for: ‎doc/distributions/cauchy.qbk

+6-11
Original file line numberDiff line numberDiff line change
@@ -105,28 +105,23 @@ In the following table __x0 is the location parameter of the distribution,
105105

106106
[table
107107
[[Function][Implementation Notes]]
108-
[[pdf][Using the relation: ['pdf = 1 / ([pi] * [gamma] * (1 + ((x - __x0) / [gamma])[super 2]) ]]]
108+
[[pdf][Using the relation: ['pdf = 1 / ([pi] * [gamma] * (1 + ((x - __x0) / [gamma])[super 2])) ]]]
109109
[[cdf and its complement][
110110
The cdf is normally given by:
111111

112112
[expression p = 0.5 + atan(x)/[pi]]
113113

114114
But that suffers from cancellation error as x -> -[infin].
115-
So recall that for `x < 0`:
116115

117-
[expression atan(x) = -[pi]/2 - atan(1/x)]
116+
Instead, the mathematically equivalent expression based on the function atan2
117+
is used:
118118

119-
Substituting into the above we get:
119+
[expression p = atan2(1, -x)/[pi]]
120120

121-
[expression p = -atan(1/x) / [pi] ; x < 0]
121+
By symmetry, the complement is
122122

123-
So the procedure is to calculate the cdf for -fabs(x)
124-
using the above formula. Note that to factor in the location and scale
125-
parameters you must substitute (x - __x0) / [gamma] for x in the above.
123+
[expression q = atan2(1, x)/[pi]]
126124

127-
This procedure yields the smaller of /p/ and /q/, so the result
128-
may need subtracting from 1 depending on whether we want the complement
129-
or not, and whether /x/ is less than __x0 or not.
130125
]]
131126
[[quantile][The same procedure is used irrespective of whether we're starting
132127
from the probability or its complement. First the argument /p/ is

Diff for: ‎include/boost/math/ccmath/floor.hpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,7 @@ inline constexpr T floor_pos_impl(T arg) noexcept
3333

3434
T result = 1;
3535

36-
if(result < arg)
36+
if(result <= arg)
3737
{
3838
while(result < arg)
3939
{

Diff for: ‎include/boost/math/ccmath/isinf.hpp

+7-2
Original file line numberDiff line numberDiff line change
@@ -23,8 +23,13 @@ constexpr bool isinf BOOST_MATH_PREVENT_MACRO_SUBSTITUTION(T x) noexcept
2323
if constexpr (std::numeric_limits<T>::is_signed)
2424
{
2525
#if defined(__clang_major__) && __clang_major__ >= 6
26-
#pragma clang diagnostic push
27-
#pragma clang diagnostic ignored "-Wtautological-constant-compare"
26+
# pragma clang diagnostic push
27+
# pragma clang diagnostic ignored "-Wtautological-constant-compare"
28+
# if defined(__has_warning)
29+
# if __has_warning("-Wnan-infinity-disabled")
30+
# pragma clang diagnostic ignored "-Wnan-infinity-disabled"
31+
# endif
32+
# endif
2833
#endif
2934
return x == std::numeric_limits<T>::infinity() || -x == std::numeric_limits<T>::infinity();
3035
#if defined(__clang_major__) && __clang_major__ >= 6

Diff for: ‎include/boost/math/distributions/cauchy.hpp

+5-22
Original file line numberDiff line numberDiff line change
@@ -45,23 +45,11 @@ BOOST_MATH_GPU_ENABLED RealType cdf_imp(const cauchy_distribution<RealType, Poli
4545
//
4646
// This calculates the cdf of the Cauchy distribution and/or its complement.
4747
//
48-
// The usual formula for the Cauchy cdf is:
48+
// This implementation uses the formula
4949
//
50-
// cdf = 0.5 + atan(x)/pi
50+
// cdf = atan2(1, -x)/pi
5151
//
52-
// But that suffers from cancellation error as x -> -INF.
53-
//
54-
// Recall that for x < 0:
55-
//
56-
// atan(x) = -pi/2 - atan(1/x)
57-
//
58-
// Substituting into the above we get:
59-
//
60-
// CDF = -atan(1/x)/pi ; x < 0
61-
//
62-
// So the procedure is to calculate the cdf for -fabs(x)
63-
// using the above formula, and then subtract from 1 when required
64-
// to get the result.
52+
// where x is the standardized (i.e. shifted and scaled) domain variable.
6553
//
6654
BOOST_MATH_STD_USING // for ADL of std functions
6755
constexpr auto function = "boost::math::cdf(cauchy<%1%>&, %1%)";
@@ -99,13 +87,8 @@ BOOST_MATH_GPU_ENABLED RealType cdf_imp(const cauchy_distribution<RealType, Poli
9987
{ // Catches x == NaN
10088
return result;
10189
}
102-
RealType mx = -fabs((x - location) / scale); // scale is > 0
103-
if(mx > -tools::epsilon<RealType>() / 8)
104-
{ // special case first: x extremely close to location.
105-
return static_cast<RealType>(0.5f);
106-
}
107-
result = -atan(1 / mx) / constants::pi<RealType>();
108-
return (((x > location) != complement) ? 1 - result : result);
90+
RealType x_std = static_cast<RealType>((complement) ? 1 : -1)*(x - location) / scale;
91+
return atan2(static_cast<RealType>(1), x_std) / constants::pi<RealType>();
10992
} // cdf
11093

11194
template <class RealType, class Policy>

Diff for: ‎include/boost/math/distributions/fwd.hpp

-2
Original file line numberDiff line numberDiff line change
@@ -11,8 +11,6 @@
1111
#ifndef BOOST_MATH_DISTRIBUTIONS_FWD_HPP
1212
#define BOOST_MATH_DISTRIBUTIONS_FWD_HPP
1313

14-
// 33 distributions at Boost 1.9.1 after adding hyperexpon and arcsine
15-
1614
namespace boost{ namespace math{
1715

1816
template <class RealType, class Policy>

Diff for: ‎include/boost/math/distributions/skew_normal.hpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@
1313
// Azzalini, A. (1985). "A class of distributions which includes the normal ones".
1414
// Scand. J. Statist. 12: 171-178.
1515

16-
#include <boost/math/distributions/fwd.hpp> // TODO add skew_normal distribution to fwd.hpp!
16+
#include <boost/math/distributions/fwd.hpp>
1717
#include <boost/math/special_functions/owens_t.hpp> // Owen's T function
1818
#include <boost/math/distributions/complement.hpp>
1919
#include <boost/math/distributions/normal.hpp>

Diff for: ‎include/boost/math/special_functions/bessel.hpp

+2-2
Original file line numberDiff line numberDiff line change
@@ -96,7 +96,7 @@ BOOST_MATH_GPU_ENABLED inline T sph_bessel_j_small_z_series(unsigned v, T x, con
9696
}
9797

9898
template <class T, class Policy>
99-
BOOST_MATH_GPU_ENABLED T cyl_bessel_j_imp_final(T v, T x, const bessel_no_int_tag& t, const Policy& pol)
99+
BOOST_MATH_GPU_ENABLED T cyl_bessel_j_imp_final(T v, T x, const bessel_no_int_tag&, const Policy& pol)
100100
{
101101
BOOST_MATH_STD_USING
102102

@@ -264,7 +264,7 @@ BOOST_MATH_GPU_ENABLED T cyl_bessel_i_imp(T v, T x, const Policy& pol)
264264
}
265265

266266
template <class T, class Policy>
267-
BOOST_MATH_GPU_ENABLED inline T cyl_bessel_k_imp(T v, T x, const bessel_no_int_tag& /* t */, const Policy& pol)
267+
BOOST_MATH_GPU_ENABLED inline T cyl_bessel_k_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol)
268268
{
269269
constexpr auto function = "boost::math::cyl_bessel_k<%1%>(%1%,%1%)";
270270
BOOST_MATH_STD_USING

Diff for: ‎include/boost/math/special_functions/detail/bessel_k0.hpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -470,7 +470,7 @@ BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T& x, const boost::math::integral_c
470470
BOOST_MATH_BIG_CONSTANT(T, 113, 8.370574966987293592457152146806662562e+03),
471471
BOOST_MATH_BIG_CONSTANT(T, 113, 4.871254714311063594080644835895740323e+01)
472472
};
473-
if(x < tools::log_max_value<T>())
473+
if(-x > tools::log_min_value<T>())
474474
return ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
475475
else
476476
{

Diff for: ‎include/boost/math/special_functions/gamma.hpp

+4-2
Original file line numberDiff line numberDiff line change
@@ -111,11 +111,13 @@ BOOST_MATH_GPU_ENABLED T sinpx(T z)
111111
// tgamma(z), with Lanczos support:
112112
//
113113
template <class T, class Policy, class Lanczos>
114-
BOOST_MATH_GPU_ENABLED T gamma_imp_final(T z, const Policy& pol, const Lanczos&)
114+
BOOST_MATH_GPU_ENABLED T gamma_imp_final(T z, const Policy& pol, const Lanczos& l)
115115
{
116116
BOOST_MATH_STD_USING
117+
118+
(void)l; // Suppresses unused variable warning when BOOST_MATH_INSTRUMENT is not defined
117119

118-
T result = 1;
120+
T result = 1;
119121

120122
#ifdef BOOST_MATH_INSTRUMENT
121123
static bool b = false;

Diff for: ‎include/boost/math/tools/precision.hpp

+9-5
Original file line numberDiff line numberDiff line change
@@ -182,9 +182,13 @@ struct log_limit_traits
182182
{
183183
typedef typename boost::math::conditional<
184184
(boost::math::numeric_limits<T>::radix == 2) &&
185-
(boost::math::numeric_limits<T>::max_exponent == 128
186-
|| boost::math::numeric_limits<T>::max_exponent == 1024
187-
|| boost::math::numeric_limits<T>::max_exponent == 16384),
185+
(
186+
( boost::math::numeric_limits<T>::max_exponent == 128
187+
|| boost::math::numeric_limits<T>::max_exponent == 1024
188+
|| boost::math::numeric_limits<T>::max_exponent == 16384
189+
)
190+
&& (-boost::math::numeric_limits<T>::min_exponent10 + 1 == boost::math::numeric_limits<T>::max_exponent10)
191+
),
188192
boost::math::integral_constant<int, (boost::math::numeric_limits<T>::max_exponent > (boost::math::numeric_limits<int>::max)() ? (boost::math::numeric_limits<int>::max)() : static_cast<int>(boost::math::numeric_limits<T>::max_exponent))>,
189193
boost::math::integral_constant<int, 0>
190194
>::type tag_type;
@@ -206,7 +210,7 @@ struct log_limit_noexcept_traits : public log_limit_noexcept_traits_imp<T, boost
206210
#endif
207211

208212
template <class T>
209-
BOOST_MATH_GPU_ENABLED inline constexpr T log_max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)
213+
BOOST_MATH_GPU_ENABLED inline T log_max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)
210214
{
211215
#ifndef BOOST_MATH_HAS_NVRTC
212216
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
@@ -223,7 +227,7 @@ BOOST_MATH_GPU_ENABLED inline constexpr T log_max_value(BOOST_MATH_EXPLICIT_TEMP
223227
}
224228

225229
template <class T>
226-
BOOST_MATH_GPU_ENABLED inline constexpr T log_min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)
230+
BOOST_MATH_GPU_ENABLED inline T log_min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)
227231
{
228232
#ifndef BOOST_MATH_HAS_NVRTC
229233
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS

Diff for: ‎reporting/performance/test_kn.cpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -48,7 +48,7 @@ static const std::array<std::array<T, 3>, 9> kn_data = { {
4848
{ { SC_(-10.0), SC_(1.0), SC_(1.80713289901029454691597861302340015908245782948536080022119e8) } },
4949
{ { SC_(100.0), SC_(5.0), SC_(7.03986019306167654653386616796116726248616158936088056952477e115) } },
5050
{ { SC_(100.0), SC_(80.0), SC_(8.39287107246490782848985384895907681748152272748337807033319e-12) } },
51-
{ { SC_(-1000.0), SC_(700.0), SC_(6.51561979144735818903553852606383312984409361984128221539405e-31) } },
51+
{ { SC_(-129.0), SC_(200.0), SC_(3.61744436315860678558682169223740584132967454950379795115566e-71) } },
5252
} };
5353

5454
int main()

Diff for: ‎test/ccmath_floor_test.cpp

+1
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@ constexpr void test()
3232
static_assert(boost::math::ccmath::isinf(boost::math::ccmath::floor(-std::numeric_limits<T>::infinity())),
3333
"If x is +- inf, it is returned, unmodified");
3434

35+
static_assert(boost::math::ccmath::floor(T(1.0)) == T(1.0));
3536
static_assert(boost::math::ccmath::floor(T(2)) == T(2));
3637
static_assert(boost::math::ccmath::floor(T(2.4)) == T(2));
3738
static_assert(boost::math::ccmath::floor(T(2.9)) == T(2));

Diff for: ‎test/cubic_roots_test.cpp

+3-3
Original file line numberDiff line numberDiff line change
@@ -99,16 +99,16 @@ template <class Real> void test_zero_coefficients() {
9999

100100
auto roots = cubic_roots(a, b, c, d);
101101
// I could check the condition number here, but this is fine right?
102-
if (!CHECK_ULP_CLOSE(r[0], roots[0], (std::numeric_limits<Real>::digits > 100 ? 120 : 25))) {
102+
if (!CHECK_ULP_CLOSE(r[0], roots[0], (std::numeric_limits<Real>::digits > 100 ? 120 : 60))) {
103103
std::cerr << " Polynomial x^3 + " << b << "x^2 + " << c << "x + "
104104
<< d << " has roots {";
105105
std::cerr << r[0] << ", " << r[1] << ", " << r[2]
106106
<< "}, but the computed roots are {";
107107
std::cerr << roots[0] << ", " << roots[1] << ", " << roots[2]
108108
<< "}\n";
109109
}
110-
CHECK_ULP_CLOSE(r[1], roots[1], 25);
111-
CHECK_ULP_CLOSE(r[2], roots[2], (std::numeric_limits<Real>::digits > 100 ? 120 : 25));
110+
CHECK_ULP_CLOSE(r[1], roots[1], 80);
111+
CHECK_ULP_CLOSE(r[2], roots[2], (std::numeric_limits<Real>::digits > 100 ? 120 : 80));
112112
for (auto root : roots) {
113113
auto res = cubic_root_residual(a, b, c, d, root);
114114
CHECK_LE(abs(res[0]), res[1]);

Diff for: ‎test/linear_regression_test.cpp

+2-2
Original file line numberDiff line numberDiff line change
@@ -223,7 +223,7 @@ void test_scaling_relations()
223223
Real c1_lambda = std::get<1>(temp);
224224
Real Rsquared_lambda = std::get<2>(temp);
225225

226-
CHECK_ULP_CLOSE(lambda*c0, c0_lambda, 50);
226+
CHECK_ULP_CLOSE(lambda*c0, c0_lambda, 70);
227227
CHECK_ULP_CLOSE(lambda*c1, c1_lambda, 30);
228228
CHECK_ULP_CLOSE(Rsquared, Rsquared_lambda, 3);
229229

@@ -241,7 +241,7 @@ void test_scaling_relations()
241241
Real c1_ = std::get<1>(temp);
242242
Real Rsquared_ = std::get<2>(temp);
243243

244-
CHECK_ULP_CLOSE(c0, c0_, 70);
244+
CHECK_ULP_CLOSE(c0, c0_, 100);
245245
CHECK_ULP_CLOSE(c1, c1_*lambda, 50);
246246
CHECK_ULP_CLOSE(Rsquared, Rsquared_, 50);
247247

Diff for: ‎test/quartic_roots_test.cpp

+2-2
Original file line numberDiff line numberDiff line change
@@ -117,8 +117,8 @@ void test_zero_coefficients()
117117

118118
roots = quartic_roots(a, b, c, d, e);
119119
// I could check the condition number here, but this is fine right?
120-
CHECK_ULP_CLOSE(r[0], roots[0], 160);
121-
CHECK_ULP_CLOSE(r[1], roots[1], 260);
120+
CHECK_ULP_CLOSE(r[0], roots[0], 340);
121+
CHECK_ULP_CLOSE(r[1], roots[1], 440);
122122
CHECK_ULP_CLOSE(r[2], roots[2], 220);
123123
CHECK_ULP_CLOSE(r[3], roots[3], 160);
124124
}

Diff for: ‎test/test_bessel_i.hpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -214,7 +214,7 @@ void test_bessel(T, const char* name)
214214
BOOST_CHECK_CLOSE_FRACTION(boost::math::cyl_bessel_i(T(0.5), T(11357)), SC_(7.173138695269929329584326974917488634629578339622112563648e4929), tolerance * mul);
215215
}
216216
#endif
217-
BOOST_IF_CONSTEXPR (std::numeric_limits<T>::max_exponent > 1000)
217+
BOOST_IF_CONSTEXPR (std::numeric_limits<T>::max_exponent10 > 304)
218218
{
219219
BOOST_IF_CONSTEXPR(std::is_floating_point<T>::value == false)
220220
tolerance *= 4; // multiprecision type.

Diff for: ‎test/test_bessel_k.hpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -133,7 +133,7 @@ void test_bessel(T, const char* name)
133133
{{ SC_(-10.0), SC_(1.0), SC_(1.80713289901029454691597861302340015908245782948536080022119e8) }},
134134
{{ SC_(100.0), SC_(5.0), SC_(7.03986019306167654653386616796116726248616158936088056952477e115) }},
135135
{{ SC_(100.0), SC_(80.0), SC_(8.39287107246490782848985384895907681748152272748337807033319e-12) }},
136-
{{ SC_(-1000.0), SC_(700.0), SC_(6.51561979144735818903553852606383312984409361984128221539405e-31) }},
136+
{{ SC_(-129.0), SC_(200.0), SC_(3.61744436315860678558682169223740584132967454950379795115566e-71) }},
137137
}};
138138
static const std::array<std::array<typename table_type<T>::type, 3>, 11> kv_data = {{
139139
{{ SC_(0.5), SC_(0.875), SC_(0.558532231646608646115729767013630967055657943463362504577189) }},

Diff for: ‎test/test_bessel_k_prime.hpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -131,7 +131,7 @@ void test_bessel(T, const char* name)
131131
{{ SC_(-10.0), SC_(1.0), SC_(-1.8171379399979651461891429013401068319174853467388121e9) }},
132132
{{ SC_(100.0), SC_(5.0), SC_(-1.4097486373570936520327835736048715219413065916411893e117) }},
133133
{{ SC_(100.0), SC_(80.0), SC_(-1.34557011017664184003144916855685180771861680634827508e-11) }},
134-
{{ SC_(-1000.0), SC_(700.0), SC_(-1.136342773238774160870536985092768591616106526374957e-30) }},
134+
{{ SC_(-129.0), SC_(200.0), SC_(-4.3110345255133348027545113739271337415489367194240537230455182e-71) }},
135135
}};
136136
static const std::array<std::array<T, 3>, 11> kv_prime_data = {{
137137
{{ SC_(0.5), SC_(0.875), SC_(-0.8776935068732421581818610624499915196588910540138553643355820) }},

Diff for: ‎test/test_cauchy.cpp

+32
Original file line numberDiff line numberDiff line change
@@ -124,6 +124,22 @@ void test_spots(RealType T)
124124
static_cast<RealType>(-10.0)), // x
125125
static_cast<RealType>(0.031725517430553569514977118601302L), // probability.
126126
tolerance); // %
127+
BOOST_CHECK_CLOSE(
128+
::boost::math::cdf(
129+
cauchy_distribution<RealType>(),
130+
static_cast<RealType>(-15000000.0)),
131+
static_cast<RealType>(0.000000021220659078919346664504384865488560725L),
132+
tolerance); // %
133+
BOOST_CHECK_CLOSE(
134+
// Test the CDF at -max_value()/4.
135+
// For an input x of this magnitude, the reference value is 4/|x|/pi.
136+
::boost::math::cdf(
137+
cauchy_distribution<RealType>(),
138+
-boost::math::tools::max_value<RealType>()/4),
139+
static_cast<RealType>(4)
140+
/ boost::math::tools::max_value<RealType>()
141+
/ boost::math::constants::pi<RealType>(),
142+
tolerance); // %
127143

128144
//
129145
// Complements:
@@ -188,6 +204,22 @@ void test_spots(RealType T)
188204
static_cast<RealType>(-10.0))), // x
189205
static_cast<RealType>(0.9682744825694464304850228813987L), // probability.
190206
tolerance); // %
207+
BOOST_CHECK_CLOSE(
208+
::boost::math::cdf(
209+
complement(cauchy_distribution<RealType>(),
210+
static_cast<RealType>(15000000.0))),
211+
static_cast<RealType>(0.000000021220659078919346664504384865488560725L),
212+
tolerance); // %
213+
BOOST_CHECK_CLOSE(
214+
// Test the complemented CDF at max_value()/4.
215+
// For an input x of this magnitude, the reference value is 4/x/pi.
216+
::boost::math::cdf(
217+
complement(cauchy_distribution<RealType>(),
218+
boost::math::tools::max_value<RealType>()/4)),
219+
static_cast<RealType>(4)
220+
/ boost::math::tools::max_value<RealType>()
221+
/ boost::math::constants::pi<RealType>(),
222+
tolerance); // %
191223

192224
//
193225
// Quantiles:

0 commit comments

Comments
 (0)