boostorg · mborland · Feb 24, 2025 · Nov 28, 2024 · Dec 6, 2024 · Dec 6, 2024
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -5,7 +5,7 @@
 
 cmake_minimum_required(VERSION 3.5...3.16)
 
-project(boost_math VERSION 1.87.0 LANGUAGES CXX)
+project(boost_math VERSION 1.88.0 LANGUAGES CXX)
 
 add_library(boost_math INTERFACE)
 

diff --git a/doc/background/special_tut.qbk b/doc/background/special_tut.qbk
@@ -373,7 +373,7 @@ Now we just need to write the test driver program, at it's most basic it looks s
       std::cout << "<note>The long double tests have been disabled on this platform "
          "either because the long double overloads of the usual math functions are "
          "not available at all, or because they are too inaccurate for these tests "
-         "to pass.</note>" << std::cout;
+         "to pass.</note>" << std::endl;
    #endif
    }
 

diff --git a/doc/distributions/cauchy.qbk b/doc/distributions/cauchy.qbk
@@ -105,28 +105,23 @@ In the following table __x0 is the location parameter of the distribution,
 
 [table
 [[Function][Implementation Notes]]
-[[pdf][Using the relation: ['pdf = 1 / ([pi] * [gamma] * (1 + ((x - __x0) / [gamma])[super 2]) ]]]
+[[pdf][Using the relation: ['pdf = 1 / ([pi] * [gamma] * (1 + ((x - __x0) / [gamma])[super 2])) ]]]
 [[cdf and its complement][
 The cdf is normally given by:
 
 [expression p = 0.5 + atan(x)/[pi]]
 
 But that suffers from cancellation error as x -> -[infin].
-So recall that for `x < 0`:
 
-[expression atan(x) = -[pi]/2 - atan(1/x)]
+Instead, the mathematically equivalent expression based on the function atan2
+is used:
 
-Substituting into the above we get:
+[expression p = atan2(1, -x)/[pi]]
 
-[expression p = -atan(1/x) / [pi]  ; x < 0]
+By symmetry, the complement is
 
-So the procedure is to calculate the cdf for -fabs(x)
-using the above formula.  Note that to factor in the location and scale
-parameters you must substitute (x - __x0) / [gamma] for x in the above.
+[expression q = atan2(1, x)/[pi]]
 
-This procedure yields the smaller of /p/ and /q/, so the result
-may need subtracting from 1 depending on whether we want the complement
-or not, and whether /x/ is less than __x0 or not.
 ]]
 [[quantile][The same procedure is used irrespective of whether we're starting
             from the probability or its complement.  First the argument /p/ is

diff --git a/include/boost/math/ccmath/floor.hpp b/include/boost/math/ccmath/floor.hpp
@@ -33,7 +33,7 @@ inline constexpr T floor_pos_impl(T arg) noexcept
 
     T result = 1;
 
-    if(result < arg)
+    if(result <= arg)
     {
         while(result < arg)
         {

diff --git a/include/boost/math/ccmath/isinf.hpp b/include/boost/math/ccmath/isinf.hpp
@@ -23,8 +23,13 @@ constexpr bool isinf BOOST_MATH_PREVENT_MACRO_SUBSTITUTION(T x) noexcept
         if constexpr (std::numeric_limits<T>::is_signed)
         {
 #if defined(__clang_major__) && __clang_major__ >= 6
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wtautological-constant-compare"
+#  pragma clang diagnostic push
+#  pragma clang diagnostic ignored "-Wtautological-constant-compare"
+#  if defined(__has_warning)
+#    if __has_warning("-Wnan-infinity-disabled")
+#      pragma clang diagnostic ignored "-Wnan-infinity-disabled"
+#    endif
+#  endif
 #endif
             return x == std::numeric_limits<T>::infinity() || -x == std::numeric_limits<T>::infinity();
 #if defined(__clang_major__) && __clang_major__ >= 6

diff --git a/include/boost/math/distributions/cauchy.hpp b/include/boost/math/distributions/cauchy.hpp
@@ -45,23 +45,11 @@ BOOST_MATH_GPU_ENABLED RealType cdf_imp(const cauchy_distribution<RealType, Poli
    //
    // This calculates the cdf of the Cauchy distribution and/or its complement.
    //
-   // The usual formula for the Cauchy cdf is:
+   // This implementation uses the formula
    //
-   // cdf = 0.5 + atan(x)/pi
+   //     cdf = atan2(1, -x)/pi
    //
-   // But that suffers from cancellation error as x -> -INF.
-   //
-   // Recall that for x < 0:
-   //
-   // atan(x) = -pi/2 - atan(1/x)
-   //
-   // Substituting into the above we get:
-   //
-   // CDF = -atan(1/x)/pi  ; x < 0
-   //
-   // So the procedure is to calculate the cdf for -fabs(x)
-   // using the above formula, and then subtract from 1 when required
-   // to get the result.
+   // where x is the standardized (i.e. shifted and scaled) domain variable.
    //
    BOOST_MATH_STD_USING // for ADL of std functions
    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
    { // Catches x == NaN
       return result;
    }
-   RealType mx = -fabs((x - location) / scale); // scale is > 0
-   if(mx > -tools::epsilon<RealType>() / 8)
-   {  // special case first: x extremely close to location.
-      return static_cast<RealType>(0.5f);
-   }
-   result = -atan(1 / mx) / constants::pi<RealType>();
-   return (((x > location) != complement) ? 1 - result : result);
+   RealType x_std = static_cast<RealType>((complement) ? 1 : -1)*(x - location) / scale;
+   return atan2(static_cast<RealType>(1), x_std) / constants::pi<RealType>();
 } // cdf
 
 template <class RealType, class Policy>

diff --git a/include/boost/math/distributions/fwd.hpp b/include/boost/math/distributions/fwd.hpp
@@ -11,8 +11,6 @@
 #ifndef BOOST_MATH_DISTRIBUTIONS_FWD_HPP
 #define BOOST_MATH_DISTRIBUTIONS_FWD_HPP
 
-// 33 distributions at Boost 1.9.1 after adding hyperexpon and arcsine
-
 namespace boost{ namespace math{
 
 template <class RealType, class Policy>

diff --git a/include/boost/math/distributions/skew_normal.hpp b/include/boost/math/distributions/skew_normal.hpp
@@ -13,7 +13,7 @@
 // Azzalini, A. (1985). "A class of distributions which includes the normal ones".
 // Scand. J. Statist. 12: 171-178.
 
-#include <boost/math/distributions/fwd.hpp> // TODO add skew_normal distribution to fwd.hpp!
+#include <boost/math/distributions/fwd.hpp>
 #include <boost/math/special_functions/owens_t.hpp> // Owen's T function
 #include <boost/math/distributions/complement.hpp>
 #include <boost/math/distributions/normal.hpp>

diff --git a/include/boost/math/special_functions/bessel.hpp b/include/boost/math/special_functions/bessel.hpp
@@ -96,7 +96,7 @@ BOOST_MATH_GPU_ENABLED inline T sph_bessel_j_small_z_series(unsigned v, T x, con
 }
 
 template <class T, class Policy>
-BOOST_MATH_GPU_ENABLED T cyl_bessel_j_imp_final(T v, T x, const bessel_no_int_tag& t, const Policy& pol)
+BOOST_MATH_GPU_ENABLED T cyl_bessel_j_imp_final(T v, T x, const bessel_no_int_tag&, const Policy& pol)
 {
    BOOST_MATH_STD_USING
 
@@ -264,7 +264,7 @@ BOOST_MATH_GPU_ENABLED T cyl_bessel_i_imp(T v, T x, const Policy& pol)
 }
 
 template <class T, class Policy>
-BOOST_MATH_GPU_ENABLED inline T cyl_bessel_k_imp(T v, T x, const bessel_no_int_tag& /* t */, const Policy& pol)
+BOOST_MATH_GPU_ENABLED inline T cyl_bessel_k_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol)
 {
    constexpr auto function = "boost::math::cyl_bessel_k<%1%>(%1%,%1%)";
    BOOST_MATH_STD_USING

diff --git a/include/boost/math/special_functions/detail/bessel_k0.hpp b/include/boost/math/special_functions/detail/bessel_k0.hpp
@@ -470,7 +470,7 @@ BOOST_MATH_GPU_ENABLED T bessel_k0_imp(const T& x, const boost::math::integral_c
             BOOST_MATH_BIG_CONSTANT(T, 113, 8.370574966987293592457152146806662562e+03),
             BOOST_MATH_BIG_CONSTANT(T, 113, 4.871254714311063594080644835895740323e+01)
          };
-         if(x < tools::log_max_value<T>())
+         if(-x > tools::log_min_value<T>())
             return  ((tools::evaluate_rational(P, Q, T(1 / x)) + Y) * exp(-x) / sqrt(x));
          else
          {

diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp
@@ -111,11 +111,13 @@ BOOST_MATH_GPU_ENABLED T sinpx(T z)
 // tgamma(z), with Lanczos support:
 //
 template <class T, class Policy, class Lanczos>
-BOOST_MATH_GPU_ENABLED T gamma_imp_final(T z, const Policy& pol, const Lanczos&)
+BOOST_MATH_GPU_ENABLED T gamma_imp_final(T z, const Policy& pol, const Lanczos& l)
 {
    BOOST_MATH_STD_USING
+
+   (void)l; // Suppresses unused variable warning when BOOST_MATH_INSTRUMENT is not defined
 
-   T result = 1;
+   T result = 1;  
 
 #ifdef BOOST_MATH_INSTRUMENT
    static bool b = false;

diff --git a/include/boost/math/tools/precision.hpp b/include/boost/math/tools/precision.hpp
@@ -182,9 +182,13 @@ struct log_limit_traits
 {
    typedef typename boost::math::conditional<
       (boost::math::numeric_limits<T>::radix == 2) &&
-      (boost::math::numeric_limits<T>::max_exponent == 128
-         || boost::math::numeric_limits<T>::max_exponent == 1024
-         || boost::math::numeric_limits<T>::max_exponent == 16384),
+      (
+         (     boost::math::numeric_limits<T>::max_exponent == 128
+            || boost::math::numeric_limits<T>::max_exponent == 1024
+            || boost::math::numeric_limits<T>::max_exponent == 16384
+         )
+         && (-boost::math::numeric_limits<T>::min_exponent10 + 1 == boost::math::numeric_limits<T>::max_exponent10)
+      ),
       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))>,
       boost::math::integral_constant<int, 0>
    >::type tag_type;
@@ -206,7 +210,7 @@ struct log_limit_noexcept_traits : public log_limit_noexcept_traits_imp<T, boost
 #endif
 
 template <class T>
-BOOST_MATH_GPU_ENABLED inline constexpr T log_max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)
+BOOST_MATH_GPU_ENABLED inline T log_max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)
 {
 #ifndef BOOST_MATH_HAS_NVRTC
    #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
 }
 
 template <class T>
-BOOST_MATH_GPU_ENABLED inline constexpr T log_min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)
+BOOST_MATH_GPU_ENABLED inline T log_min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) noexcept(detail::log_limit_noexcept_traits<T>::value)
 {
 #ifndef BOOST_MATH_HAS_NVRTC
    #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS

diff --git a/reporting/performance/test_kn.cpp b/reporting/performance/test_kn.cpp
@@ -48,7 +48,7 @@ static const std::array<std::array<T, 3>, 9> kn_data = { {
    { { SC_(-10.0), SC_(1.0), SC_(1.80713289901029454691597861302340015908245782948536080022119e8) } },
    { { SC_(100.0), SC_(5.0), SC_(7.03986019306167654653386616796116726248616158936088056952477e115) } },
    { { SC_(100.0), SC_(80.0), SC_(8.39287107246490782848985384895907681748152272748337807033319e-12) } },
-   { { SC_(-1000.0), SC_(700.0), SC_(6.51561979144735818903553852606383312984409361984128221539405e-31) } },
+   { { SC_(-129.0), SC_(200.0), SC_(3.61744436315860678558682169223740584132967454950379795115566e-71) } },
    } };
 
 int main()

diff --git a/test/ccmath_floor_test.cpp b/test/ccmath_floor_test.cpp
@@ -32,6 +32,7 @@ constexpr void test()
     static_assert(boost::math::ccmath::isinf(boost::math::ccmath::floor(-std::numeric_limits<T>::infinity())), 
                   "If x is +- inf, it is returned, unmodified");
 
+    static_assert(boost::math::ccmath::floor(T(1.0)) == T(1.0));
     static_assert(boost::math::ccmath::floor(T(2)) == T(2));
     static_assert(boost::math::ccmath::floor(T(2.4)) == T(2));
     static_assert(boost::math::ccmath::floor(T(2.9)) == T(2));

diff --git a/test/cubic_roots_test.cpp b/test/cubic_roots_test.cpp
@@ -99,16 +99,16 @@ template <class Real> void test_zero_coefficients() {
 
         auto roots = cubic_roots(a, b, c, d);
         // I could check the condition number here, but this is fine right?
-        if (!CHECK_ULP_CLOSE(r[0], roots[0], (std::numeric_limits<Real>::digits > 100 ? 120 : 25))) {
+        if (!CHECK_ULP_CLOSE(r[0], roots[0], (std::numeric_limits<Real>::digits > 100 ? 120 : 60))) {
             std::cerr << "  Polynomial x^3 + " << b << "x^2 + " << c << "x + "
                       << d << " has roots {";
             std::cerr << r[0] << ", " << r[1] << ", " << r[2]
                       << "}, but the computed roots are {";
             std::cerr << roots[0] << ", " << roots[1] << ", " << roots[2]
                       << "}\n";
         }
-        CHECK_ULP_CLOSE(r[1], roots[1], 25);
-        CHECK_ULP_CLOSE(r[2], roots[2], (std::numeric_limits<Real>::digits > 100 ? 120 : 25));
+        CHECK_ULP_CLOSE(r[1], roots[1], 80);
+        CHECK_ULP_CLOSE(r[2], roots[2], (std::numeric_limits<Real>::digits > 100 ? 120 : 80));
         for (auto root : roots) {
             auto res = cubic_root_residual(a, b, c, d, root);
             CHECK_LE(abs(res[0]), res[1]);

diff --git a/test/linear_regression_test.cpp b/test/linear_regression_test.cpp
@@ -223,7 +223,7 @@ void test_scaling_relations()
     Real c1_lambda = std::get<1>(temp);
     Real Rsquared_lambda = std::get<2>(temp);
 
-    CHECK_ULP_CLOSE(lambda*c0, c0_lambda, 50);
+    CHECK_ULP_CLOSE(lambda*c0, c0_lambda, 70);
     CHECK_ULP_CLOSE(lambda*c1, c1_lambda, 30);
     CHECK_ULP_CLOSE(Rsquared, Rsquared_lambda, 3);
 
@@ -241,7 +241,7 @@ void test_scaling_relations()
     Real c1_ = std::get<1>(temp);
     Real Rsquared_ = std::get<2>(temp);
 
-    CHECK_ULP_CLOSE(c0, c0_, 70);
+    CHECK_ULP_CLOSE(c0, c0_, 100);
     CHECK_ULP_CLOSE(c1, c1_*lambda, 50);
     CHECK_ULP_CLOSE(Rsquared, Rsquared_, 50);
 

diff --git a/test/quartic_roots_test.cpp b/test/quartic_roots_test.cpp
@@ -117,8 +117,8 @@ void test_zero_coefficients()
 
         roots = quartic_roots(a, b, c, d, e);
         // I could check the condition number here, but this is fine right?
-        CHECK_ULP_CLOSE(r[0], roots[0], 160);
-        CHECK_ULP_CLOSE(r[1], roots[1], 260);
+        CHECK_ULP_CLOSE(r[0], roots[0], 340);
+        CHECK_ULP_CLOSE(r[1], roots[1], 440);
         CHECK_ULP_CLOSE(r[2], roots[2], 220);
         CHECK_ULP_CLOSE(r[3], roots[3], 160);
     }

diff --git a/test/test_bessel_i.hpp b/test/test_bessel_i.hpp
@@ -214,7 +214,7 @@ void test_bessel(T, const char* name)
        BOOST_CHECK_CLOSE_FRACTION(boost::math::cyl_bessel_i(T(0.5), T(11357)), SC_(7.173138695269929329584326974917488634629578339622112563648e4929), tolerance * mul);
     }
 #endif
-    BOOST_IF_CONSTEXPR (std::numeric_limits<T>::max_exponent > 1000)
+    BOOST_IF_CONSTEXPR (std::numeric_limits<T>::max_exponent10 > 304)
     {
        BOOST_IF_CONSTEXPR(std::is_floating_point<T>::value == false)
           tolerance *= 4; // multiprecision type.

diff --git a/test/test_bessel_k.hpp b/test/test_bessel_k.hpp
@@ -133,7 +133,7 @@ void test_bessel(T, const char* name)
         {{ SC_(-10.0), SC_(1.0), SC_(1.80713289901029454691597861302340015908245782948536080022119e8) }},
         {{ SC_(100.0), SC_(5.0), SC_(7.03986019306167654653386616796116726248616158936088056952477e115) }},
         {{ SC_(100.0), SC_(80.0), SC_(8.39287107246490782848985384895907681748152272748337807033319e-12) }},
-        {{ SC_(-1000.0), SC_(700.0), SC_(6.51561979144735818903553852606383312984409361984128221539405e-31) }},
+        {{ SC_(-129.0), SC_(200.0), SC_(3.61744436315860678558682169223740584132967454950379795115566e-71) }},
     }};
     static const std::array<std::array<typename table_type<T>::type, 3>, 11> kv_data = {{
         {{ SC_(0.5), SC_(0.875), SC_(0.558532231646608646115729767013630967055657943463362504577189) }},

diff --git a/test/test_bessel_k_prime.hpp b/test/test_bessel_k_prime.hpp
@@ -131,7 +131,7 @@ void test_bessel(T, const char* name)
         {{ SC_(-10.0), SC_(1.0), SC_(-1.8171379399979651461891429013401068319174853467388121e9) }},
         {{ SC_(100.0), SC_(5.0), SC_(-1.4097486373570936520327835736048715219413065916411893e117) }},
         {{ SC_(100.0), SC_(80.0), SC_(-1.34557011017664184003144916855685180771861680634827508e-11) }},
-        {{ SC_(-1000.0), SC_(700.0), SC_(-1.136342773238774160870536985092768591616106526374957e-30) }},
+        {{ SC_(-129.0), SC_(200.0), SC_(-4.3110345255133348027545113739271337415489367194240537230455182e-71) }},
     }};
     static const std::array<std::array<T, 3>, 11> kv_prime_data = {{
         {{ SC_(0.5), SC_(0.875), SC_(-0.8776935068732421581818610624499915196588910540138553643355820) }},

diff --git a/test/test_cauchy.cpp b/test/test_cauchy.cpp
@@ -124,6 +124,22 @@ void test_spots(RealType T)
          static_cast<RealType>(-10.0)),              // x
          static_cast<RealType>(0.031725517430553569514977118601302L),                // probability.
          tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::cdf(
+         cauchy_distribution<RealType>(),
+         static_cast<RealType>(-15000000.0)),
+         static_cast<RealType>(0.000000021220659078919346664504384865488560725L),
+         tolerance); // %
+   BOOST_CHECK_CLOSE(
+      // Test the CDF at -max_value()/4.
+      // For an input x of this magnitude, the reference value is 4/|x|/pi.
+      ::boost::math::cdf(
+         cauchy_distribution<RealType>(),
+         -boost::math::tools::max_value<RealType>()/4),
+         static_cast<RealType>(4)
+                      / boost::math::tools::max_value<RealType>()
+                      / boost::math::constants::pi<RealType>(),
+         tolerance); // %
 
    //
    // Complements:
@@ -188,6 +204,22 @@ void test_spots(RealType T)
          static_cast<RealType>(-10.0))),              // x
          static_cast<RealType>(0.9682744825694464304850228813987L),                // probability.
          tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::cdf(
+         complement(cauchy_distribution<RealType>(),
+         static_cast<RealType>(15000000.0))),
+         static_cast<RealType>(0.000000021220659078919346664504384865488560725L),
+         tolerance); // %
+   BOOST_CHECK_CLOSE(
+      // Test the complemented CDF at max_value()/4.
+      // For an input x of this magnitude, the reference value is 4/x/pi.
+      ::boost::math::cdf(
+         complement(cauchy_distribution<RealType>(),
+         boost::math::tools::max_value<RealType>()/4)),
+         static_cast<RealType>(4)
+                      / boost::math::tools::max_value<RealType>()
+                      / boost::math::constants::pi<RealType>(),
+         tolerance); // %
 
    //
    // Quantiles: