- Extreme
+ Extreme
Value Distributions, Theory and Applications Samuel Kotz & Saralees
Nadarajah discuss the relationship of the types of extreme
value distributions.
diff --git a/doc/html/math_toolkit/dist_ref/dists/extreme_dist.html b/doc/html/math_toolkit/dist_ref/dists/extreme_dist.html
index 2cad77f118..261413c36b 100644
--- a/doc/html/math_toolkit/dist_ref/dists/extreme_dist.html
+++ b/doc/html/math_toolkit/dist_ref/dists/extreme_dist.html
@@ -65,7 +65,7 @@
The relationship of the types of extreme value distributions, of which
- this is but one, is discussed by Extreme
+ this is but one, is discussed by Extreme
Value Distributions, Theory and Applications Samuel Kotz & Saralees
Nadarajah.
diff --git a/doc/html/math_toolkit/dist_ref/dists/triangular_dist.html b/doc/html/math_toolkit/dist_ref/dists/triangular_dist.html
index 5fb9b76a71..8c9da8079a 100644
--- a/doc/html/math_toolkit/dist_ref/dists/triangular_dist.html
+++ b/doc/html/math_toolkit/dist_ref/dists/triangular_dist.html
@@ -64,7 +64,7 @@
vaguely known, but, like the uniform
distribution, upper and limits are 'known', but a 'best guess',
the mode or center point, is also added. It has been recommended as a
- proxy
+ proxy
for the beta distribution. The distribution is used in business
decision making and project planning.
diff --git a/doc/html/math_toolkit/dist_ref/dists/weibull_dist.html b/doc/html/math_toolkit/dist_ref/dists/weibull_dist.html
index 06a9970cc7..f94d0800bb 100644
--- a/doc/html/math_toolkit/dist_ref/dists/weibull_dist.html
+++ b/doc/html/math_toolkit/dist_ref/dists/weibull_dist.html
@@ -108,7 +108,7 @@
distribution. When α = 1, the Weibull distribution
reduces to the exponential
distribution. The relationship of the types of extreme value distributions,
- of which the Weibull is but one, is discussed by Extreme
+ of which the Weibull is but one, is discussed by Extreme
Value Distributions, Theory and Applications Samuel Kotz & Saralees
Nadarajah.
diff --git a/doc/html/math_toolkit/refs.html b/doc/html/math_toolkit/refs.html
index f59e301ac8..4dd108f64f 100644
--- a/doc/html/math_toolkit/refs.html
+++ b/doc/html/math_toolkit/refs.html
@@ -78,7 +78,7 @@
by N.A.J. Hastings, Brian Peacock, Merran Evans, ISBN: 0471371246, Wiley 2000.
- Extreme Value
+ Extreme Value
Distributions, Theory and Applications Samuel Kotz & Saralees Nadarajah,
ISBN 978-1-86094-224-2 & 1-86094-224-5 Oct 2000, Chapter 1.2 discusses
the various extreme value distributions.
diff --git a/doc/math.qbk b/doc/math.qbk
index 902db34dcf..f3374a73d6 100644
--- a/doc/math.qbk
+++ b/doc/math.qbk
@@ -615,6 +615,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22.
[include statistics/runs_test.qbk]
[include statistics/ljung_box.qbk]
[include statistics/linear_regression.qbk]
+[include statistics/chatterjee_correlation.qbk]
[endmathpart] [/section:statistics Statistics]
[mathpart vector_functionals Vector Functionals - Norms]
diff --git a/doc/sf/daubechies.qbk b/doc/sf/daubechies.qbk
index 2d280f81b7..e1e7f52539 100644
--- a/doc/sf/daubechies.qbk
+++ b/doc/sf/daubechies.qbk
@@ -127,6 +127,48 @@ The 2 vanishing moment scaling function.
[$../graphs/daubechies_8_scaling.svg]
The 8 vanishing moment scaling function.
+Boost.Math also provides numerical evaluation of the Fourier transform of these functions.
+This is useful in sparse recovery problems where the measurements are taken in the Fourier basis.
+The usage is exhibited below:
+
+ #include
+ using boost::math::fourier_transform_daubechies_scaling;
+ // Evaluate the Fourier transform of the 4-vanishing moment Daubechies scaling function at ω=1.8:
+ std::complex hat_phi = fourier_transform_daubechies_scaling(1.8f);
+
+The Fourier transform convention is unitary with the sign of the imaginary unit being given in Daubechies Ten Lectures.
+In particular, this means that `fourier_transform_daubechies_scaling(0.0)` returns 1/sqrt(2π).
+
+The implementation computes an infinite product of trigonometric polynomials as can be found from recursive application of the identity 𝓕[φ](ω) = m(ω/2)𝓕[φ](ω/2).
+This is neither particularly fast nor accurate, but there appears to be no literature on this extremely useful topic, and hence the naive method must suffice.
+
+[$../graphs/fourier_transform_daubechies.png]
+
+A benchmark can be found in `reporting/performance/fourier_transform_daubechies_performance.cpp`; the results on a ~2021 M1 Macbook pro are presented below:
+
+
+Run on (10 X 24.1212 MHz CPU s)
+CPU Caches:
+ L1 Data 64 KiB (x10)
+ L1 Instruction 128 KiB (x10)
+ L2 Unified 4096 KiB (x5)
+Load Average: 1.33, 1.52, 1.62
+-----------------------------------------------------------
+Benchmark Time
+-----------------------------------------------------------
+FourierTransformDaubechiesScaling 70.3 ns
+FourierTransformDaubechiesScaling 330 ns
+FourierTransformDaubechiesScaling 335 ns
+FourierTransformDaubechiesScaling 364 ns
+FourierTransformDaubechiesScaling 386 ns
+FourierTransformDaubechiesScaling 436 ns
+FourierTransformDaubechiesScaling 447 ns
+FourierTransformDaubechiesScaling 473 ns
+FourierTransformDaubechiesScaling 503 ns
+FourierTransformDaubechiesScaling 554 ns
+
+Due to the low accuracy of this method, `float` precision is arg-promoted to `double`, and hence takes just as long as `double` precision to execute.
+
[heading References]
* Daubechies, Ingrid. ['Ten Lectures on Wavelets.] Vol. 61. Siam, 1992.
diff --git a/example/calculate_fourier_transform_daubechies_constants.cpp b/example/calculate_fourier_transform_daubechies_constants.cpp
new file mode 100644
index 0000000000..d4ad1f3d18
--- /dev/null
+++ b/example/calculate_fourier_transform_daubechies_constants.cpp
@@ -0,0 +1,43 @@
+// (C) Copyright Nick Thompson 2023.
+// 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)
+
+#include
+#include
+#include
+#include
+#include
+
+using std::pow;
+using boost::multiprecision::cpp_bin_float_100;
+using boost::math::filters::daubechies_scaling_filter;
+using boost::math::tools::polynomial;
+using boost::math::constants::half;
+using boost::math::constants::root_two;
+
+template
+std::vector get_constants() {
+ auto h = daubechies_scaling_filter();
+ auto p = polynomial(h.begin(), h.end());
+
+ auto q = polynomial({half(), half()});
+ q = pow(q, N);
+ auto l = p/q;
+ return l.data();
+}
+
+template
+void print_constants(std::vector const & l) {
+ std::cout << std::setprecision(std::numeric_limits::digits10 -10);
+ std::cout << "return std::array{";
+ for (size_t i = 0; i < l.size() - 1; ++i) {
+ std::cout << "BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits::digits, " << l[i]/root_two() << "), ";
+ }
+ std::cout << "BOOST_MATH_BIG_CONSTANT(Real, std::numeric_limits::digits, " << l.back()/root_two() << ")};\n";
+}
+
+int main() {
+ auto constants = get_constants();
+ print_constants(constants);
+}
diff --git a/example/fourier_transform_daubechies_ulp_plot.cpp b/example/fourier_transform_daubechies_ulp_plot.cpp
new file mode 100644
index 0000000000..3cde316708
--- /dev/null
+++ b/example/fourier_transform_daubechies_ulp_plot.cpp
@@ -0,0 +1,60 @@
+// boost-no-inspect
+// (C) Copyright Nick Thompson 2023.
+// 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)
+
+#include
+#include
+
+using boost::math::fourier_transform_daubechies_scaling;
+using boost::math::tools::ulps_plot;
+
+template
+void real_part() {
+ auto phi_real_hi_acc = [](double omega) {
+ auto z = fourier_transform_daubechies_scaling(omega);
+ return z.real();
+ };
+
+ auto phi_real_lo_acc = [](float omega) {
+ auto z = fourier_transform_daubechies_scaling(omega);
+ return z.real();
+ };
+ auto plot = ulps_plot(phi_real_hi_acc, float(0.0), float(100.0), 20000);
+ plot.ulp_envelope(false);
+ plot.add_fn(phi_real_lo_acc);
+ plot.clip(100);
+ plot.title("Accuracy of 𝔑(𝓕[𝜙](ω)) with " + std::to_string(p) + " vanishing moments.");
+ plot.write("real_ft_daub_scaling_" + std::to_string(p) + ".svg");
+
+}
+
+template
+void imaginary_part() {
+ auto phi_imag_hi_acc = [](double omega) {
+ auto z = fourier_transform_daubechies_scaling(omega);
+ return z.imag();
+ };
+
+ auto phi_imag_lo_acc = [](float omega) {
+ auto z = fourier_transform_daubechies_scaling(omega);
+ return z.imag();
+ };
+ auto plot = ulps_plot(phi_imag_hi_acc, float(0.0), float(100.0), 20000);
+ plot.ulp_envelope(false);
+ plot.add_fn(phi_imag_lo_acc);
+ plot.clip(100);
+ plot.title("Accuracy of 𝕴(𝓕[𝜙](ω)) with " + std::to_string(p) + " vanishing moments.");
+ plot.write("imag_ft_daub_scaling_" + std::to_string(p) + ".svg");
+
+}
+
+
+int main() {
+ real_part<3>();
+ imaginary_part<3>();
+ real_part<6>();
+ imaginary_part<6>();
+ return 0;
+}
diff --git a/example/lambert_w_example.cpp b/example/lambert_w_example.cpp
index 5c862a3b76..f4e1880d00 100644
--- a/example/lambert_w_example.cpp
+++ b/example/lambert_w_example.cpp
@@ -12,8 +12,6 @@
#include // for BOOST_PLATFORM, BOOST_COMPILER, BOOST_STDLIB ...
#include // for BOOST_MSVC versions.
-#include
-#include // boost::exception
#include // For exp_minus_one == 3.67879441171442321595523770161460867e-01.
#define BOOST_MATH_INSTRUMENT_LAMBERT_W // #define only for diagnostic output.
diff --git a/example/lambert_w_precision_example.cpp b/example/lambert_w_precision_example.cpp
index e95f3b1931..18120c2b4e 100644
--- a/example/lambert_w_precision_example.cpp
+++ b/example/lambert_w_precision_example.cpp
@@ -10,8 +10,6 @@
#include // for BOOST_PLATFORM, BOOST_COMPILER, BOOST_STDLIB ...
#include // for BOOST_MSVC versions.
-#include
-#include // boost::exception
#include // For exp_minus_one == 3.67879441171442321595523770161460867e-01.
#include
#include // for float_distance.
diff --git a/example/lambert_w_simple_examples.cpp b/example/lambert_w_simple_examples.cpp
index 3c422e0626..443e22fb42 100644
--- a/example/lambert_w_simple_examples.cpp
+++ b/example/lambert_w_simple_examples.cpp
@@ -20,8 +20,6 @@
#include // for BOOST_PLATFORM, BOOST_COMPILER, BOOST_STDLIB ...
#include // for BOOST_MSVC versions.
-#include
-#include // boost::exception
#include // For exp_minus_one == 3.67879441171442321595523770161460867e-01.
#include
diff --git a/include/boost/math/ccmath/abs.hpp b/include/boost/math/ccmath/abs.hpp
index fa798eba12..84e01fba01 100644
--- a/include/boost/math/ccmath/abs.hpp
+++ b/include/boost/math/ccmath/abs.hpp
@@ -31,7 +31,7 @@ namespace detail {
template
constexpr T abs_impl(T x) noexcept
{
- if (boost::math::ccmath::isnan(x))
+ if ((boost::math::ccmath::isnan)(x))
{
return std::numeric_limits::quiet_NaN();
}
diff --git a/include/boost/math/ccmath/fpclassify.hpp b/include/boost/math/ccmath/fpclassify.hpp
index 2750af6d3e..4faba80310 100644
--- a/include/boost/math/ccmath/fpclassify.hpp
+++ b/include/boost/math/ccmath/fpclassify.hpp
@@ -10,6 +10,7 @@
#include
#include
#include
+#include
#include
#include
#include
@@ -18,19 +19,19 @@
namespace boost::math::ccmath {
template , bool> = true>
-inline constexpr int fpclassify(T x)
+inline constexpr int fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(T x)
{
if(BOOST_MATH_IS_CONSTANT_EVALUATED(x))
{
- return boost::math::ccmath::isnan(x) ? FP_NAN :
- boost::math::ccmath::isinf(x) ? FP_INFINITE :
+ return (boost::math::ccmath::isnan)(x) ? FP_NAN :
+ (boost::math::ccmath::isinf)(x) ? FP_INFINITE :
boost::math::ccmath::abs(x) == T(0) ? FP_ZERO :
boost::math::ccmath::abs(x) > 0 && boost::math::ccmath::abs(x) < (std::numeric_limits::min)() ? FP_SUBNORMAL : FP_NORMAL;
}
else
{
- using std::fpclassify;
- return fpclassify(x);
+ using boost::math::fpclassify;
+ return (fpclassify)(x);
}
}
diff --git a/include/boost/math/ccmath/isinf.hpp b/include/boost/math/ccmath/isinf.hpp
index 706b9605c2..7b942dcf72 100644
--- a/include/boost/math/ccmath/isinf.hpp
+++ b/include/boost/math/ccmath/isinf.hpp
@@ -10,6 +10,7 @@
#include
#include
#include
+#include
#include
#ifndef BOOST_MATH_STANDALONE
@@ -22,7 +23,7 @@
namespace boost::math::ccmath {
template
-constexpr bool isinf(T x) noexcept
+constexpr bool isinf BOOST_PREVENT_MACRO_SUBSTITUTION(T x) noexcept
{
if(BOOST_MATH_IS_CONSTANT_EVALUATED(x))
{
@@ -37,15 +38,15 @@ constexpr bool isinf(T x) noexcept
}
else
{
- using std::isinf;
+ using boost::math::isinf;
if constexpr (!std::is_integral_v)
{
- return isinf(x);
+ return (isinf)(x);
}
else
{
- return isinf(static_cast(x));
+ return (isinf)(static_cast(x));
}
}
}
diff --git a/include/boost/math/ccmath/isnan.hpp b/include/boost/math/ccmath/isnan.hpp
index ca915a4642..a5e121a97f 100644
--- a/include/boost/math/ccmath/isnan.hpp
+++ b/include/boost/math/ccmath/isnan.hpp
@@ -9,6 +9,7 @@
#include
#include
#include
+#include
#include
#ifndef BOOST_MATH_STANDALONE
@@ -21,7 +22,7 @@
namespace boost::math::ccmath {
template
-inline constexpr bool isnan(T x)
+inline constexpr bool isnan BOOST_PREVENT_MACRO_SUBSTITUTION(T x)
{
if(BOOST_MATH_IS_CONSTANT_EVALUATED(x))
{
@@ -29,15 +30,15 @@ inline constexpr bool isnan(T x)
}
else
{
- using std::isnan;
+ using boost::math::isnan;
if constexpr (!std::is_integral_v)
{
- return isnan(x);
+ return (isnan)(x);
}
else
{
- return isnan(static_cast(x));
+ return (isnan)(static_cast(x));
}
}
}
diff --git a/include/boost/math/ccmath/ldexp.hpp b/include/boost/math/ccmath/ldexp.hpp
index 92277706fb..31730b265f 100644
--- a/include/boost/math/ccmath/ldexp.hpp
+++ b/include/boost/math/ccmath/ldexp.hpp
@@ -45,8 +45,8 @@ inline constexpr Real ldexp(Real arg, int exp) noexcept
if(BOOST_MATH_IS_CONSTANT_EVALUATED(arg))
{
return boost::math::ccmath::abs(arg) == Real(0) ? arg :
- boost::math::ccmath::isinf(arg) ? arg :
- boost::math::ccmath::isnan(arg) ? arg :
+ (boost::math::ccmath::isinf)(arg) ? arg :
+ (boost::math::ccmath::isnan)(arg) ? arg :
boost::math::ccmath::detail::ldexp_impl(arg, exp);
}
else
diff --git a/include/boost/math/concepts/real_concept.hpp b/include/boost/math/concepts/real_concept.hpp
index 446bfa8d5a..56e4342e6b 100644
--- a/include/boost/math/concepts/real_concept.hpp
+++ b/include/boost/math/concepts/real_concept.hpp
@@ -45,6 +45,10 @@
# include
#endif
+#if __has_include()
+# include
+#endif
+
namespace boost{ namespace math{
namespace concepts
@@ -79,6 +83,12 @@ class real_concept
#ifdef BOOST_MATH_USE_FLOAT128
real_concept(BOOST_MATH_FLOAT128_TYPE c) : m_value(c){}
#endif
+#ifdef __STDCPP_FLOAT32_T__
+ real_concept(std::float32_t c) : m_value(static_cast(c)){}
+#endif
+#ifdef __STDCPP_FLOAT64_T__
+ real_concept(std::float64_t c) : m_value(static_cast(c)){}
+#endif
// Assignment:
real_concept& operator=(char c) { m_value = c; return *this; }
@@ -96,6 +106,12 @@ class real_concept
real_concept& operator=(float c) { m_value = c; return *this; }
real_concept& operator=(double c) { m_value = c; return *this; }
real_concept& operator=(long double c) { m_value = c; return *this; }
+ #ifdef __STDCPP_FLOAT32_T__
+ real_concept& operator=(std::float32_t c) { m_value = c; return *this; }
+ #endif
+ #ifdef __STDCPP_FLOAT64_T__
+ real_concept& operator=(std::float64_t c) { m_value = c; return *this; }
+ #endif
// Access:
real_concept_base_type value()const{ return m_value; }
diff --git a/include/boost/math/cstdfloat/cstdfloat_limits.hpp b/include/boost/math/cstdfloat/cstdfloat_limits.hpp
index 43265424d9..9661ab6526 100644
--- a/include/boost/math/cstdfloat/cstdfloat_limits.hpp
+++ b/include/boost/math/cstdfloat/cstdfloat_limits.hpp
@@ -24,7 +24,7 @@
#pragma GCC system_header
#endif
- #if defined(BOOST_CSTDFLOAT_HAS_INTERNAL_FLOAT128_T) && defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT)
+ #if defined(BOOST_CSTDFLOAT_HAS_INTERNAL_FLOAT128_T) && defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT) && (!defined(__GNUC__) || (defined(__GNUC__) && __GNUC__ < 14))
#include
#include
diff --git a/include/boost/math/differentiation/autodiff.hpp b/include/boost/math/differentiation/autodiff.hpp
index 1286326c2e..7cda196852 100644
--- a/include/boost/math/differentiation/autodiff.hpp
+++ b/include/boost/math/differentiation/autodiff.hpp
@@ -1491,7 +1491,7 @@ fvar sqrt(fvar const& cr) {
BOOST_IF_CONSTEXPR (order == 0)
return fvar(*derivatives);
else {
- root_type numerator = 0.5;
+ root_type numerator = root_type(0.5);
root_type powers = 1;
#ifndef BOOST_NO_CXX17_IF_CONSTEXPR
derivatives[1] = numerator / *derivatives;
diff --git a/include/boost/math/differentiation/finite_difference.hpp b/include/boost/math/differentiation/finite_difference.hpp
index 5e4433e000..1375bac7ad 100644
--- a/include/boost/math/differentiation/finite_difference.hpp
+++ b/include/boost/math/differentiation/finite_difference.hpp
@@ -153,7 +153,7 @@ namespace detail {
const Real eps = (numeric_limits::epsilon)();
// Error bound ~eps^4/5
- Real h = pow(11.25*eps, static_cast(1) / static_cast(5));
+ Real h = pow(Real(11.25)*eps, static_cast(1) / static_cast(5));
h = detail::make_xph_representable(x, h);
Real ymth = f(x - 2 * h);
Real yth = f(x + 2 * h);
@@ -222,7 +222,7 @@ namespace detail {
// Mathematica code to get the error:
// Series[(f[x+h]-f[x-h])*(4/5) + (1/5)*(f[x-2*h] - f[x+2*h]) + (4/105)*(f[x+3*h] - f[x-3*h]) + (1/280)*(f[x-4*h] - f[x+4*h]), {h, 0, 9}]
// If we used Kahan summation, we could get the max error down to h^8|f^(9)(x)|/630 + |f(x)|eps/h.
- Real h = pow(551.25*eps, static_cast(1) / static_cast(9));
+ Real h = pow(Real(551.25)*eps, static_cast(1) / static_cast(9));
h = detail::make_xph_representable(x, h);
Real yh = f(x + h);
diff --git a/include/boost/math/differentiation/lanczos_smoothing.hpp b/include/boost/math/differentiation/lanczos_smoothing.hpp
index 1f18c2f6d5..9406b302d7 100644
--- a/include/boost/math/differentiation/lanczos_smoothing.hpp
+++ b/include/boost/math/differentiation/lanczos_smoothing.hpp
@@ -10,6 +10,7 @@
#include // to nan initialize
#include
#include
+#include
#include
#include
#include
diff --git a/include/boost/math/distributions/beta.hpp b/include/boost/math/distributions/beta.hpp
index 1de399ead2..6c17ffa1a2 100644
--- a/include/boost/math/distributions/beta.hpp
+++ b/include/boost/math/distributions/beta.hpp
@@ -2,6 +2,7 @@
// Copyright John Maddock 2006.
// Copyright Paul A. Bristow 2006.
+// Copyright Matt Borland 2023.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
@@ -241,7 +242,7 @@ namespace boost
{
return result;
}
- return ibeta_inva(beta, x, probability, Policy());
+ return static_cast(ibeta_inva(beta, x, probability, Policy()));
} // RealType find_alpha(beta, a, probability)
static RealType find_beta(
@@ -264,7 +265,7 @@ namespace boost
{
return result;
}
- return ibeta_invb(alpha, x, probability, Policy());
+ return static_cast(ibeta_invb(alpha, x, probability, Policy()));
} // RealType find_beta(alpha, x, probability)
private:
@@ -390,13 +391,40 @@ namespace boost
}
using boost::math::beta;
- // Corner case: check_x ensures x element of [0, 1], but PDF is 0 for x = 0 and x = 1. PDF EQN:
+ // Corner cases: check_x ensures x element of [0, 1], but PDF is 0 for x = 0 and x = 1. PDF EQN:
// https://wikimedia.org/api/rest_v1/media/math/render/svg/125fdaa41844a8703d1a8610ac00fbf3edacc8e7
- if(x == 0 || x == 1)
+ if(x == 0)
{
- return RealType(0);
+ if (a == 1)
+ {
+ return static_cast(1 / beta(a, b));
+ }
+ else if (a < 1)
+ {
+ policies::raise_overflow_error(function, nullptr, Policy());
+ }
+ else
+ {
+ return RealType(0);
+ }
}
- return ibeta_derivative(a, b, x, Policy());
+ else if (x == 1)
+ {
+ if (b == 1)
+ {
+ return static_cast(1 / beta(a, b));
+ }
+ else if (b < 1)
+ {
+ policies::raise_overflow_error(function, nullptr, Policy());
+ }
+ else
+ {
+ return RealType(0);
+ }
+ }
+
+ return static_cast(ibeta_derivative(a, b, x, Policy()));
} // pdf
template
@@ -427,7 +455,7 @@ namespace boost
{
return 1;
}
- return ibeta(a, b, x, Policy());
+ return static_cast(ibeta(a, b, x, Policy()));
} // beta cdf
template
@@ -454,16 +482,16 @@ namespace boost
}
if (x == 0)
{
- return 1;
+ return RealType(1);
}
else if (x == 1)
{
- return 0;
+ return RealType(0);
}
// Calculate cdf beta using the incomplete beta function.
// Use of ibeta here prevents cancellation errors in calculating
// 1 - x if x is very small, perhaps smaller than machine epsilon.
- return ibetac(a, b, x, Policy());
+ return static_cast(ibetac(a, b, x, Policy()));
} // beta cdf
template
@@ -492,13 +520,13 @@ namespace boost
// Special cases:
if (p == 0)
{
- return 0;
+ return RealType(0);
}
if (p == 1)
{
- return 1;
+ return RealType(1);
}
- return ibeta_inv(a, b, p, static_cast(nullptr), Policy());
+ return static_cast(ibeta_inv(a, b, p, static_cast(nullptr), Policy()));
} // quantile
template
@@ -528,14 +556,14 @@ namespace boost
// Special cases:
if(q == 1)
{
- return 0;
+ return RealType(0);
}
if(q == 0)
{
- return 1;
+ return RealType(1);
}
- return ibetac_inv(a, b, q, static_cast(nullptr), Policy());
+ return static_cast(ibetac_inv(a, b, q, static_cast(nullptr), Policy()));
} // Quantile Complement
} // namespace math
diff --git a/include/boost/math/distributions/cauchy.hpp b/include/boost/math/distributions/cauchy.hpp
index 8588686c8f..d914cca77e 100644
--- a/include/boost/math/distributions/cauchy.hpp
+++ b/include/boost/math/distributions/cauchy.hpp
@@ -81,7 +81,7 @@ RealType cdf_imp(const cauchy_distribution& dist, const RealTy
RealType mx = -fabs((x - location) / scale); // scale is > 0
if(mx > -tools::epsilon() / 8)
{ // special case first: x extremely close to location.
- return 0.5;
+ return static_cast(0.5f);
}
result = -atan(1 / mx) / constants::pi();
return (((x > location) != complement) ? 1 - result : result);
diff --git a/include/boost/math/distributions/detail/derived_accessors.hpp b/include/boost/math/distributions/detail/derived_accessors.hpp
index d278e78776..eb76409a1c 100644
--- a/include/boost/math/distributions/detail/derived_accessors.hpp
+++ b/include/boost/math/distributions/detail/derived_accessors.hpp
@@ -122,6 +122,13 @@ inline typename Distribution::value_type cdf(const Distribution& dist, const Rea
typedef typename Distribution::value_type value_type;
return cdf(dist, static_cast(x));
}
+template
+inline typename Distribution::value_type logcdf(const Distribution& dist, const Realtype& x)
+{
+ using std::log;
+ using value_type = typename Distribution::value_type;
+ return log(cdf(dist, static_cast(x)));
+}
template
inline typename Distribution::value_type quantile(const Distribution& dist, const RealType& x)
{
@@ -143,6 +150,14 @@ inline typename Distribution::value_type cdf(const complemented2_type(c.param)));
}
+template
+inline typename Distribution::value_type logcdf(const complemented2_type& c)
+{
+ using std::log;
+ typedef typename Distribution::value_type value_type;
+ return log(cdf(complement(c.dist, static_cast(c.param))));
+}
+
template
inline typename Distribution::value_type quantile(const complemented2_type& c)
{
diff --git a/include/boost/math/distributions/detail/hypergeometric_cdf.hpp b/include/boost/math/distributions/detail/hypergeometric_cdf.hpp
index 067be3e525..60e29bc10c 100644
--- a/include/boost/math/distributions/detail/hypergeometric_cdf.hpp
+++ b/include/boost/math/distributions/detail/hypergeometric_cdf.hpp
@@ -10,11 +10,12 @@
#include
#include
+#include
namespace boost{ namespace math{ namespace detail{
template
- T hypergeometric_cdf_imp(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, const Policy& pol)
+ T hypergeometric_cdf_imp(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, bool invert, const Policy& pol)
{
#ifdef _MSC_VER
# pragma warning(push)
@@ -27,7 +28,7 @@ namespace boost{ namespace math{ namespace detail{
{
result = hypergeometric_pdf(x, r, n, N, pol);
T diff = result;
- unsigned lower_limit = static_cast((std::max)(0, static_cast(n + r) - static_cast(N)));
+ const auto lower_limit = static_cast((std::max)(INT64_C(0), static_cast(n + r) - static_cast(N)));
while(diff > (invert ? T(1) : result) * tools::epsilon())
{
diff = T(x) * T((N + x) - n - r) * diff / (T(1 + n - x) * T(1 + r - x));
@@ -43,7 +44,7 @@ namespace boost{ namespace math{ namespace detail{
else
{
invert = !invert;
- unsigned upper_limit = (std::min)(r, n);
+ const auto upper_limit = (std::min)(r, n);
if(x != upper_limit)
{
++x;
@@ -69,7 +70,7 @@ namespace boost{ namespace math{ namespace detail{
}
template
- inline T hypergeometric_cdf(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, const Policy&)
+ inline T hypergeometric_cdf(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, bool invert, const Policy&)
{
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args::type result_type;
diff --git a/include/boost/math/distributions/detail/hypergeometric_pdf.hpp b/include/boost/math/distributions/detail/hypergeometric_pdf.hpp
index 220c8d8e0c..9eeef270d8 100644
--- a/include/boost/math/distributions/detail/hypergeometric_pdf.hpp
+++ b/include/boost/math/distributions/detail/hypergeometric_pdf.hpp
@@ -15,6 +15,8 @@
#include
#include
#include
+#include
+#include
#ifdef BOOST_MATH_INSTRUMENT
#include
@@ -40,7 +42,7 @@ template
struct sort_functor
{
sort_functor(const T* exponents) : m_exponents(exponents){}
- bool operator()(int i, int j)
+ bool operator()(std::size_t i, std::size_t j)
{
return m_exponents[i] > m_exponents[j];
}
@@ -49,7 +51,7 @@ struct sort_functor
};
template
-T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n, unsigned N, const Lanczos&, const Policy&)
+T hypergeometric_pdf_lanczos_imp(T /*dummy*/, std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Lanczos&, const Policy&)
{
BOOST_MATH_STD_USING
@@ -137,7 +139,7 @@ T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n
//
// Combine equal powers:
//
- int j = 8;
+ std::size_t j = 8;
while(exponents[sorted_indexes[j]] == 0) --j;
while(j)
{
@@ -184,7 +186,7 @@ T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n
result = pow(bases[sorted_indexes[0]] * exp(static_cast(base_e_factors[sorted_indexes[0]])), exponents[sorted_indexes[0]]);
}
BOOST_MATH_INSTRUMENT_VARIABLE(result);
- for(unsigned i = 1; (i < 9) && (exponents[sorted_indexes[i]] > 0); ++i)
+ for(std::size_t i = 1; (i < 9) && (exponents[sorted_indexes[i]] > 0); ++i)
{
BOOST_FPU_EXCEPTION_GUARD
if(result < tools::min_value())
@@ -215,7 +217,7 @@ T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n
}
template
-T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n, unsigned N, const boost::math::lanczos::undefined_lanczos&, const Policy& pol)
+T hypergeometric_pdf_lanczos_imp(T /*dummy*/, std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const boost::math::lanczos::undefined_lanczos&, const Policy& pol)
{
BOOST_MATH_STD_USING
return exp(
@@ -260,7 +262,7 @@ inline T integer_power(const T& x, int ex)
#ifdef __SUNPRO_CC
return pow(x, T(ex));
#else
- return pow(x, ex);
+ return static_cast(pow(x, ex));
#endif
}
template
@@ -277,12 +279,12 @@ struct hypergeometric_pdf_prime_loop_result_entry
struct hypergeometric_pdf_prime_loop_data
{
- const unsigned x;
- const unsigned r;
- const unsigned n;
- const unsigned N;
- unsigned prime_index;
- unsigned current_prime;
+ const std::uint64_t x;
+ const std::uint64_t r;
+ const std::uint64_t n;
+ const std::uint64_t N;
+ std::size_t prime_index;
+ std::uint64_t current_prime;
};
#ifdef _MSC_VER
@@ -294,8 +296,8 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
{
while(data.current_prime <= data.N)
{
- unsigned base = data.current_prime;
- int prime_powers = 0;
+ std::uint64_t base = data.current_prime;
+ std::int64_t prime_powers = 0;
while(base <= data.N)
{
prime_powers += data.n / base;
@@ -381,7 +383,7 @@ T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hy
}
template
-inline T hypergeometric_pdf_prime_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
+inline T hypergeometric_pdf_prime_imp(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)
{
hypergeometric_pdf_prime_loop_result_entry result = { 1, 0 };
hypergeometric_pdf_prime_loop_data data = { x, r, n, N, 0, prime(0) };
@@ -389,7 +391,7 @@ inline T hypergeometric_pdf_prime_imp(unsigned x, unsigned r, unsigned n, unsign
}
template
-T hypergeometric_pdf_factorial_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
+T hypergeometric_pdf_factorial_imp(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)
{
BOOST_MATH_STD_USING
BOOST_MATH_ASSERT(N <= boost::math::max_factorial::value);
@@ -406,8 +408,8 @@ T hypergeometric_pdf_factorial_imp(unsigned x, unsigned r, unsigned n, unsigned
boost::math::unchecked_factorial(r - x),
boost::math::unchecked_factorial(N - n - r + x)
};
- int i = 0;
- int j = 0;
+ std::size_t i = 0;
+ std::size_t j = 0;
while((i < 3) || (j < 5))
{
while((j < 5) && ((result >= 1) || (i >= 3)))
@@ -427,7 +429,7 @@ T hypergeometric_pdf_factorial_imp(unsigned x, unsigned r, unsigned n, unsigned
template
inline typename tools::promote_args::type
- hypergeometric_pdf(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&)
+ hypergeometric_pdf(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)
{
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args::type result_type;
diff --git a/include/boost/math/distributions/exponential.hpp b/include/boost/math/distributions/exponential.hpp
index 5214575a64..164e01f205 100644
--- a/include/boost/math/distributions/exponential.hpp
+++ b/include/boost/math/distributions/exponential.hpp
@@ -163,6 +163,24 @@ inline RealType cdf(const exponential_distribution& dist, cons
return result;
} // cdf
+template
+inline RealType logcdf(const exponential_distribution& dist, const RealType& x)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ static const char* function = "boost::math::logcdf(const exponential_distribution<%1%>&, %1%)";
+
+ RealType result = 0;
+ RealType lambda = dist.lambda();
+ if(0 == detail::verify_lambda(function, lambda, &result, Policy()))
+ return result;
+ if(0 == detail::verify_exp_x(function, x, &result, Policy()))
+ return result;
+ result = boost::math::log1p(-exp(-x * lambda), Policy());
+
+ return result;
+} // cdf
+
template
inline RealType quantile(const exponential_distribution& dist, const RealType& p)
{
@@ -207,6 +225,27 @@ inline RealType cdf(const complemented2_type
+inline RealType logcdf(const complemented2_type, RealType>& c)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ static const char* function = "boost::math::logcdf(const exponential_distribution<%1%>&, %1%)";
+
+ RealType result = 0;
+ RealType lambda = c.dist.lambda();
+ if(0 == detail::verify_lambda(function, lambda, &result, Policy()))
+ return result;
+ if(0 == detail::verify_exp_x(function, c.param, &result, Policy()))
+ return result;
+ // Workaround for VC11/12 bug:
+ if (c.param >= tools::max_value())
+ return 0;
+ result = -c.param * lambda;
+
+ return result;
+}
+
template
inline RealType quantile(const complemented2_type, RealType>& c)
{
diff --git a/include/boost/math/distributions/extreme_value.hpp b/include/boost/math/distributions/extreme_value.hpp
index b2e8bea966..1bde2743c0 100644
--- a/include/boost/math/distributions/extreme_value.hpp
+++ b/include/boost/math/distributions/extreme_value.hpp
@@ -174,6 +174,32 @@ inline RealType cdf(const extreme_value_distribution& dist, co
return result;
} // cdf
+template
+inline RealType logcdf(const extreme_value_distribution& dist, const RealType& x)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ static const char* function = "boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)";
+
+ if((boost::math::isinf)(x))
+ return x < 0 ? 0.0f : 1.0f;
+ RealType a = dist.location();
+ RealType b = dist.scale();
+ RealType result = 0;
+ if(0 == detail::verify_scale_b(function, b, &result, Policy()))
+ return result;
+ if(0 == detail::check_finite(function, a, &result, Policy()))
+ return result;
+ if(0 == detail::check_finite(function, a, &result, Policy()))
+ return result;
+ if(0 == detail::check_x("boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)", x, &result, Policy()))
+ return result;
+
+ result = -exp((a-x)/b);
+
+ return result;
+} // logcdf
+
template
RealType quantile(const extreme_value_distribution& dist, const RealType& p)
{
@@ -225,6 +251,30 @@ inline RealType cdf(const complemented2_type
+inline RealType logcdf(const complemented2_type, RealType>& c)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ static const char* function = "boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)";
+
+ if((boost::math::isinf)(c.param))
+ return c.param < 0 ? 1.0f : 0.0f;
+ RealType a = c.dist.location();
+ RealType b = c.dist.scale();
+ RealType result = 0;
+ if(0 == detail::verify_scale_b(function, b, &result, Policy()))
+ return result;
+ if(0 == detail::check_finite(function, a, &result, Policy()))
+ return result;
+ if(0 == detail::check_x(function, c.param, &result, Policy()))
+ return result;
+
+ result = log1p(-exp(-exp((a-c.param)/b)), Policy());
+
+ return result;
+}
+
template
RealType quantile(const complemented2_type, RealType>& c)
{
diff --git a/include/boost/math/distributions/fisher_f.hpp b/include/boost/math/distributions/fisher_f.hpp
index cec79eccd6..e22cdf50ae 100644
--- a/include/boost/math/distributions/fisher_f.hpp
+++ b/include/boost/math/distributions/fisher_f.hpp
@@ -306,10 +306,10 @@ inline RealType mode(const fisher_f_distribution& dist)
&& detail::check_df(
function, df2, &error_result, Policy()))
return error_result;
- if(df2 <= 2)
+ if(df1 <= 2)
{
return policies::raise_domain_error(
- function, "Second degree of freedom was %1% but must be > 2 in order for the distribution to have a mode.", df2, Policy());
+ function, "First degree of freedom was %1% but must be > 2 in order for the distribution to have a mode.", df1, Policy());
}
return df2 * (df1 - 2) / (df1 * (df2 + 2));
}
diff --git a/include/boost/math/distributions/geometric.hpp b/include/boost/math/distributions/geometric.hpp
index b091bf1c63..baff0e249f 100644
--- a/include/boost/math/distributions/geometric.hpp
+++ b/include/boost/math/distributions/geometric.hpp
@@ -43,9 +43,11 @@
#include // isnan.
#include // for root finding.
#include
+#include
#include // using std::numeric_limits;
#include
+#include
#if defined (BOOST_MSVC)
# pragma warning(push)
@@ -378,9 +380,39 @@ namespace boost
return probability;
} // cdf Cumulative Distribution Function geometric.
- template
- inline RealType cdf(const complemented2_type, RealType>& c)
- { // Complemented Cumulative Distribution Function geometric.
+ template
+ inline RealType logcdf(const geometric_distribution& dist, const RealType& k)
+ { // Cumulative Distribution Function of geometric.
+ using std::pow;
+ static const char* function = "boost::math::logcdf(const geometric_distribution<%1%>&, %1%)";
+
+ // k argument may be integral, signed, or unsigned, or floating point.
+ // If necessary, it has already been promoted from an integral type.
+ RealType p = dist.success_fraction();
+ // Error check:
+ RealType result = 0;
+ if(false == geometric_detail::check_dist_and_k(
+ function,
+ p,
+ k,
+ &result, Policy()))
+ {
+ return -std::numeric_limits::infinity();
+ }
+ if(k == 0)
+ {
+ return log(p); // success_fraction
+ }
+ //RealType q = 1 - p; // Bad for small p
+ //RealType probability = 1 - std::pow(q, k+1);
+
+ RealType z = boost::math::log1p(-p, Policy()) * (k + 1);
+ return log1p(-exp(z), Policy());
+ } // logcdf Cumulative Distribution Function geometric.
+
+ template
+ inline RealType cdf(const complemented2_type, RealType>& c)
+ { // Complemented Cumulative Distribution Function geometric.
BOOST_MATH_STD_USING
static const char* function = "boost::math::cdf(const geometric_distribution<%1%>&, %1%)";
// k argument may be integral, signed, or unsigned, or floating point.
@@ -403,6 +435,30 @@ namespace boost
return probability;
} // cdf Complemented Cumulative Distribution Function geometric.
+ template
+ inline RealType logcdf(const complemented2_type, RealType>& c)
+ { // Complemented Cumulative Distribution Function geometric.
+ BOOST_MATH_STD_USING
+ static const char* function = "boost::math::logcdf(const geometric_distribution<%1%>&, %1%)";
+ // k argument may be integral, signed, or unsigned, or floating point.
+ // If necessary, it has already been promoted from an integral type.
+ RealType const& k = c.param;
+ geometric_distribution const& dist = c.dist;
+ RealType p = dist.success_fraction();
+ // Error check:
+ RealType result = 0;
+ if(false == geometric_detail::check_dist_and_k(
+ function,
+ p,
+ k,
+ &result, Policy()))
+ {
+ return -std::numeric_limits::infinity();
+ }
+
+ return boost::math::log1p(-p, Policy()) * (k+1);
+ } // logcdf Complemented Cumulative Distribution Function geometric.
+
template
inline RealType quantile(const geometric_distribution& dist, const RealType& x)
{ // Quantile, percentile/100 or Percent Point geometric function.
diff --git a/include/boost/math/distributions/hypergeometric.hpp b/include/boost/math/distributions/hypergeometric.hpp
index 58b8d42727..d7d1b6a069 100644
--- a/include/boost/math/distributions/hypergeometric.hpp
+++ b/include/boost/math/distributions/hypergeometric.hpp
@@ -16,6 +16,7 @@
#include
#include
#include
+#include
namespace boost { namespace math {
@@ -26,7 +27,7 @@ namespace boost { namespace math {
typedef RealType value_type;
typedef Policy policy_type;
- hypergeometric_distribution(unsigned r, unsigned n, unsigned N) // Constructor. r=defective/failures/success, n=trials/draws, N=total population.
+ hypergeometric_distribution(std::uint64_t r, std::uint64_t n, std::uint64_t N) // Constructor. r=defective/failures/success, n=trials/draws, N=total population.
: m_n(n), m_N(N), m_r(r)
{
static const char* function = "boost::math::hypergeometric_distribution<%1%>::hypergeometric_distribution";
@@ -34,17 +35,17 @@ namespace boost { namespace math {
check_params(function, &ret);
}
// Accessor functions.
- unsigned total()const
+ std::uint64_t total() const
{
return m_N;
}
- unsigned defective()const // successes/failures/events
+ std::uint64_t defective() const // successes/failures/events
{
return m_r;
}
- unsigned sample_count()const
+ std::uint64_t sample_count()const
{
return m_n;
}
@@ -65,9 +66,9 @@ namespace boost { namespace math {
}
return true;
}
- bool check_x(unsigned x, const char* function, RealType* result)const
+ bool check_x(std::uint64_t x, const char* function, RealType* result)const
{
- if(x < static_cast((std::max)(0, static_cast(m_n + m_r) - static_cast(m_N))))
+ if(x < static_cast((std::max)(INT64_C(0), static_cast(m_n + m_r) - static_cast(m_N))))
{
*result = boost::math::policies::raise_domain_error(
function, "Random variable out of range: must be > 0 and > m + r - N but got %1%", static_cast(x), Policy());
@@ -84,40 +85,40 @@ namespace boost { namespace math {
private:
// Data members:
- unsigned m_n; // number of items picked or drawn.
- unsigned m_N; // number of "total" items.
- unsigned m_r; // number of "defective/successes/failures/events items.
+ std::uint64_t m_n; // number of items picked or drawn.
+ std::uint64_t m_N; // number of "total" items.
+ std::uint64_t m_r; // number of "defective/successes/failures/events items.
}; // class hypergeometric_distribution
typedef hypergeometric_distribution hypergeometric;
template
- inline const std::pair range(const hypergeometric_distribution& dist)
+ inline const std::pair range(const hypergeometric_distribution& dist)
{ // Range of permissible values for random variable x.
#ifdef _MSC_VER
# pragma warning(push)
# pragma warning(disable:4267)
#endif
- unsigned r = dist.defective();
- unsigned n = dist.sample_count();
- unsigned N = dist.total();
- unsigned l = static_cast((std::max)(0, static_cast(n + r) - static_cast(N)));
- unsigned u = (std::min)(r, n);
- return std::pair(l, u);
+ const auto r = dist.defective();
+ const auto n = dist.sample_count();
+ const auto N = dist.total();
+ const auto l = static_cast((std::max)(INT64_C(0), static_cast