Update contrib/restricted/boost/math to 1.84.0

17 files changed, 262 insertions, 195 deletions

diff --git a/contrib/restricted/boost/math/include/boost/math/ccmath/abs.hpp b/contrib/restricted/boost/math/include/boost/math/ccmath/abs.hpp
index 84e01fba01..db59502ab6 100644
--- a/contrib/restricted/boost/math/include/boost/math/ccmath/abs.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/ccmath/abs.hpp
@@ -8,22 +8,16 @@
 #ifndef BOOST_MATH_CCMATH_ABS
 #define BOOST_MATH_CCMATH_ABS
 
-#include <cmath>
-#include <type_traits>
-#include <limits>
-#include <boost/math/tools/is_constant_evaluated.hpp>
+#include <boost/math/ccmath/detail/config.hpp>
+
+#ifdef BOOST_MATH_NO_CCMATH
+#error "The header <boost/math/abs.hpp> can only be used in C++17 and later."
+#endif
+
 #include <boost/math/tools/assert.hpp>
 #include <boost/math/ccmath/isnan.hpp>
 #include <boost/math/ccmath/isinf.hpp>
 
-#include <boost/math/tools/is_standalone.hpp>
-#ifndef BOOST_MATH_STANDALONE
-#include <boost/config.hpp>
-#ifdef BOOST_NO_CXX17_IF_CONSTEXPR
-#error "The header <boost/math/norms.hpp> can only be used in C++17 and later."
-#endif
-#endif
-
 namespace boost::math::ccmath {
 
 namespace detail {
diff --git a/contrib/restricted/boost/math/include/boost/math/ccmath/detail/config.hpp b/contrib/restricted/boost/math/include/boost/math/ccmath/detail/config.hpp
new file mode 100644
index 0000000000..61405e781a
--- /dev/null
+++ b/contrib/restricted/boost/math/include/boost/math/ccmath/detail/config.hpp
@@ -0,0 +1,52 @@
+//  (C) Copyright John Maddock 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)
+//
+//  Core configuration for ccmath functions, basically will they work or not?
+
+#ifndef BOOST_MATH_CCMATH_DETAIL_CONFIG
+#define BOOST_MATH_CCMATH_DETAIL_CONFIG
+
+#include <cmath>
+#include <type_traits>
+#include <limits>
+#include <boost/math/tools/is_constant_evaluated.hpp>
+#include <boost/math/tools/is_standalone.hpp>
+
+#ifndef BOOST_MATH_STANDALONE
+
+#include <boost/config.hpp>
+#ifdef BOOST_NO_CXX17_IF_CONSTEXPR
+#  define BOOST_MATH_NO_CCMATH
+#endif
+
+#else // BOOST_MATH_STANDALONE
+
+#if defined(_MSC_VER)
+
+#if defined(_MSVC_LANG) && (_MSVC_LANG < 201703)
+#  define BOOST_MATH_NO_CCMATH
+#endif
+
+#else // _MSC_VER
+
+#if (__cplusplus < 201703)
+#  define BOOST_MATH_NO_CCMATH
+#endif
+
+#endif
+
+#endif
+
+#ifndef _MSC_VER
+//
+// Don't check here for msvc as they didn't get std lib configuration macros at the same time as C++17 <type_traits>
+//
+#if (__cpp_lib_bool_constant < 201505L) && !defined(BOOST_MATH_NO_CCMATH)
+#  define BOOST_MATH_NO_CCMATH
+#endif
+#endif
+
+
+#endif
diff --git a/contrib/restricted/boost/math/include/boost/math/ccmath/isinf.hpp b/contrib/restricted/boost/math/include/boost/math/ccmath/isinf.hpp
index 7b942dcf72..dfbd1f00af 100644
--- a/contrib/restricted/boost/math/include/boost/math/ccmath/isinf.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/ccmath/isinf.hpp
@@ -6,18 +6,11 @@
 #ifndef BOOST_MATH_CCMATH_ISINF
 #define BOOST_MATH_CCMATH_ISINF
 
-#include <cmath>
-#include <limits>
-#include <type_traits>
-#include <boost/math/tools/is_constant_evaluated.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/ccmath/detail/config.hpp>
 
-#include <boost/math/tools/is_standalone.hpp>
-#ifndef BOOST_MATH_STANDALONE
-#include <boost/config.hpp>
-#ifdef BOOST_NO_CXX17_IF_CONSTEXPR
-#error "The header <boost/math/norms.hpp> can only be used in C++17 and later."
-#endif
+#ifdef BOOST_MATH_NO_CCMATH
+#error "The header <boost/math/isinf.hpp> can only be used in C++17 and later."
 #endif
 
 namespace boost::math::ccmath {
diff --git a/contrib/restricted/boost/math/include/boost/math/ccmath/isnan.hpp b/contrib/restricted/boost/math/include/boost/math/ccmath/isnan.hpp
index a5e121a97f..504fa90d2c 100644
--- a/contrib/restricted/boost/math/include/boost/math/ccmath/isnan.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/ccmath/isnan.hpp
@@ -6,17 +6,11 @@
 #ifndef BOOST_MATH_CCMATH_ISNAN
 #define BOOST_MATH_CCMATH_ISNAN
 
-#include <cmath>
-#include <type_traits>
-#include <boost/math/tools/is_constant_evaluated.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/ccmath/detail/config.hpp>
 
-#include <boost/math/tools/is_standalone.hpp>
-#ifndef BOOST_MATH_STANDALONE
-#include <boost/config.hpp>
-#ifdef BOOST_NO_CXX17_IF_CONSTEXPR
-#error "The header <boost/math/norms.hpp> can only be used in C++17 and later."
-#endif
+#ifdef BOOST_MATH_NO_CCMATH
+#error "The header <boost/math/isnan.hpp> can only be used in C++17 and later."
 #endif
 
 namespace boost::math::ccmath {
diff --git a/contrib/restricted/boost/math/include/boost/math/ccmath/ldexp.hpp b/contrib/restricted/boost/math/include/boost/math/ccmath/ldexp.hpp
index 31730b265f..3e2cd3610f 100644
--- a/contrib/restricted/boost/math/include/boost/math/ccmath/ldexp.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/ccmath/ldexp.hpp
@@ -7,11 +7,13 @@
 #ifndef BOOST_MATH_CCMATH_LDEXP_HPP
 #define BOOST_MATH_CCMATH_LDEXP_HPP
 
-#include <cmath>
-#include <limits>
-#include <type_traits>
+#include <boost/math/ccmath/detail/config.hpp>
+
+#ifdef BOOST_MATH_NO_CCMATH
+#error "The header <boost/math/ldexp.hpp> can only be used in C++17 and later."
+#endif
+
 #include <stdexcept>
-#include <boost/math/tools/is_constant_evaluated.hpp>
 #include <boost/math/ccmath/abs.hpp>
 #include <boost/math/ccmath/isinf.hpp>
 #include <boost/math/ccmath/isnan.hpp>
diff --git a/contrib/restricted/boost/math/include/boost/math/cstdfloat/cstdfloat_cmath.hpp b/contrib/restricted/boost/math/include/boost/math/cstdfloat/cstdfloat_cmath.hpp
index 049ac4091e..515f3052a2 100644
--- a/contrib/restricted/boost/math/include/boost/math/cstdfloat/cstdfloat_cmath.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/cstdfloat/cstdfloat_cmath.hpp
@@ -918,62 +918,6 @@ namespace std
    using boost::math::cstdfloat::detail::isunordered;
    //   end more functions
 
-   //
-   // Very basic iostream operator:
-   //
-   inline std::ostream& operator << (std::ostream& os, __float128 m_value)
-   {
-      std::streamsize digits = os.precision();
-      std::ios_base::fmtflags f = os.flags();
-      std::string s;
-
-      char buf[100];
-      std::unique_ptr<char[]> buf2;
-      std::string format = "%";
-      if (f & std::ios_base::showpos)
-         format += "+";
-      if (f & std::ios_base::showpoint)
-         format += "#";
-      format += ".*";
-      if (digits == 0)
-         digits = 36;
-      format += "Q";
-      if (f & std::ios_base::scientific)
-         format += "e";
-      else if (f & std::ios_base::fixed)
-         format += "f";
-      else
-         format += "g";
-
-      int v = quadmath_snprintf(buf, 100, format.c_str(), digits, m_value);
-
-      if ((v < 0) || (v >= 99))
-      {
-         int v_max = v;
-         buf2.reset(new char[v + 3]);
-         v = quadmath_snprintf(&buf2[0], v_max + 3, format.c_str(), digits, m_value);
-         if (v >= v_max + 3)
-         {
-            BOOST_MATH_THROW_EXCEPTION(std::runtime_error("Formatting of float128_type failed."));
-         }
-         s = &buf2[0];
-      }
-      else
-         s = buf;
-      std::streamsize ss = os.width();
-      if (ss > static_cast<std::streamsize>(s.size()))
-      {
-         char fill = os.fill();
-         if ((os.flags() & std::ios_base::left) == std::ios_base::left)
-            s.append(static_cast<std::string::size_type>(ss - s.size()), fill);
-         else
-            s.insert(static_cast<std::string::size_type>(0), static_cast<std::string::size_type>(ss - s.size()), fill);
-      }
-
-      return os << s;
-   }
-
-
 } // namespace std
 
 // We will now remove the preprocessor symbols representing quadruple-precision <cmath>
diff --git a/contrib/restricted/boost/math/include/boost/math/cstdfloat/cstdfloat_iostream.hpp b/contrib/restricted/boost/math/include/boost/math/cstdfloat/cstdfloat_iostream.hpp
index c01c236f71..7c2c396776 100644
--- a/contrib/restricted/boost/math/include/boost/math/cstdfloat/cstdfloat_iostream.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/cstdfloat/cstdfloat_iostream.hpp
@@ -32,8 +32,7 @@
   #include <boost/math/tools/nothrow.hpp>
   #include <boost/math/tools/throw_exception.hpp>
 
-//  #if (0)
-  #if defined(__GNUC__)
+#if defined(__GNUC__) && !defined(BOOST_MATH_TEST_IO_AS_INTEL_QUAD)
 
   // Forward declarations of quadruple-precision string functions.
   extern "C" int quadmath_snprintf(char *str, size_t size, const char *format, ...) BOOST_MATH_NOTHROW;
@@ -96,7 +95,7 @@
         // So we have to use dynamic memory allocation for the output
         // string buffer.
 
-        char* my_buffer2 = static_cast<char*>(0U);
+        char* my_buffer2 = nullptr;
 
 #ifndef BOOST_NO_EXCEPTIONS
         try
@@ -160,8 +159,7 @@
     }
   }
 
-//  #elif defined(__GNUC__)
-  #elif defined(__INTEL_COMPILER)
+#elif defined(__INTEL_COMPILER) || defined(BOOST_MATH_TEST_IO_AS_INTEL_QUAD)
 
   // The section for I/O stream support for the ICC compiler is particularly
   // long, because these functions must be painstakingly synthesized from
@@ -172,6 +170,7 @@
   // used in Boost.Multiprecision by John Maddock and Christopher Kormanyos.
   // This methodology has been slightly modified here for boost::float128_t.
 
+
   #include <cstring>
   #include <cctype>
   
@@ -266,7 +265,7 @@
     {
       // Pad out the end with zero's if we need to.
 
-      int chars = static_cast<int>(str.size());
+      std::ptrdiff_t chars = static_cast<std::ptrdiff_t>(str.size());
       chars = digits - chars;
 
       if(scientific)
@@ -442,7 +441,7 @@
       if(isneg) { x = -x; }
 
       float_type t;
-      float_type ten = 10;
+      constexpr float_type ten = 10;
 
       eval_log10(t, x);
       eval_floor(t, t);
@@ -507,6 +506,8 @@
         eval_subtract(t, digit);
         eval_multiply(t, ten);
       }
+      if (result.size() == 0)
+         result = "0";
 
       // Possibly round the result.
       if(digits >= 0)
@@ -522,11 +523,13 @@
           if((static_cast<int>(*result.rbegin() - '0') & 1) != 0)
           {
             round_string_up_at(result, static_cast<int>(result.size() - 1U), expon);
+            if (digits == 0) digits = 1;
           }
         }
         else if(cdigit >= 5)
         {
-          round_string_up_at(result, static_cast<int>(result.size() - 1), expon);
+          round_string_up_at(result, static_cast<int>(result.size() - 1u), expon);
+          if (digits == 0) digits = 1;
         }
       }
     }
@@ -569,9 +572,9 @@
   {
     value = 0;
 
-    if((p == static_cast<const char*>(0U)) || (*p == static_cast<char>(0)))
+    if((p == nullptr) || (*p == '\0'))
     {
-      return;
+      return false;
     }
 
     bool is_neg       = false;
@@ -584,11 +587,11 @@
 
     constexpr int max_digits = std::numeric_limits<float_type>::max_digits10 + 1;
 
-    if(*p == static_cast<char>('+'))
+    if(*p == '+')
     {
       ++p;
     }
-    else if(*p == static_cast<char>('-'))
+    else if(*p == '-')
     {
       is_neg = true;
       ++p;
@@ -632,7 +635,7 @@
       ++digits_seen;
     }
 
-    if(*p == static_cast<char>('.'))
+    if(*p == '.')
     {
       // Grab everything after the point, stop when we've seen
       // enough digits, even if there are actually more available.
@@ -659,15 +662,15 @@
     }
 
     // Parse the exponent.
-    if((*p == static_cast<char>('e')) || (*p == static_cast<char>('E')))
+    if((*p == 'e') || (*p == 'E'))
     {
       ++p;
 
-      if(*p == static_cast<char>('+'))
+      if(*p == '+')
       {
         ++p;
       }
-      else if(*p == static_cast<char>('-'))
+      else if(*p == '-')
       {
         is_neg_expon = true;
         ++p;
@@ -718,7 +721,7 @@
       value = -value;
     }
 
-    return (*p == static_cast<char>(0));
+    return (*p == '\0');
   }
   } } } } // boost::math::cstdfloat::detail
 
diff --git a/contrib/restricted/boost/math/include/boost/math/special_functions/beta.hpp b/contrib/restricted/boost/math/include/boost/math/special_functions/beta.hpp
index ec824daf72..f15e157324 100644
--- a/contrib/restricted/boost/math/include/boost/math/special_functions/beta.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/special_functions/beta.hpp
@@ -230,7 +230,10 @@ T ibeta_power_terms(T a,
    T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
    T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
    T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
-   result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
+   if ((a < tools::min_value<T>()) || (b < tools::min_value<T>()))
+      result = 0;  // denominator overflows in this case
+   else
+      result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
    result *= prefix;
    // combine with the leftover terms from the Lanczos approximation:
    result *= sqrt(bgh / boost::math::constants::e<T>());
@@ -389,14 +392,15 @@ T ibeta_power_terms(T a,
          }
          else
          {
-            T p1 = pow(b1, a / b);
+            // This protects against spurious overflow in a/b:
+            T p1 = (b1 < 1) && (b < 1) && (tools::max_value<T>() * b < a) ? static_cast<T>(0) : static_cast<T>(pow(b1, a / b));
             T l3 = (p1 != 0) && (b2 != 0) ? (log(p1) + log(b2)) * b : tools::max_value<T>();  // arbitrary large value if the logs would fail!
             if((l3 < tools::log_max_value<T>())
                && (l3 > tools::log_min_value<T>()))
             {
                result *= pow(p1 * b2, b);
             }
-            else
+            else if(result != 0)  // we can elude the calculation below if we're already going to be zero
             {
                l2 += l1 + log(result);
                if(l2 >= tools::log_max_value<T>())
@@ -473,20 +477,108 @@ T ibeta_power_terms(T a,
    if ((shift_a == 0) && (shift_b == 0))
    {
       T power1, power2;
+      bool need_logs = false;
       if (a < b)
       {
-         power1 = pow((x * y * c * c) / (a * b), a);
-         power2 = pow((y * c) / b, b - a);
+         BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
+         {
+            power1 = pow((x * y * c * c) / (a * b), a);
+            power2 = pow((y * c) / b, b - a);
+         }
+         else
+         {
+            // We calculate these logs purely so we can check for overflow in the power functions
+            T l1 = log((x * y * c * c) / (a * b));
+            T l2 = log((y * c) / b);
+            if ((l1 * a > tools::log_min_value<T>()) && (l1 * a < tools::log_max_value<T>()) && (l2 * (b - a) < tools::log_max_value<T>()) && (l2 * (b - a) > tools::log_min_value<T>()))
+            {
+               power1 = pow((x * y * c * c) / (a * b), a);
+               power2 = pow((y * c) / b, b - a);
+            }
+            else
+            {
+               need_logs = true;
+            }
+         }
       }
       else
       {
-         power1 = pow((x * y * c * c) / (a * b), b);
-         power2 = pow((x * c) / a, a - b);
+         BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
+         {
+            power1 = pow((x * y * c * c) / (a * b), b);
+            power2 = pow((x * c) / a, a - b);
+         }
+         else
+         {
+            // We calculate these logs purely so we can check for overflow in the power functions
+            T l1 = log((x * y * c * c) / (a * b)) * b;
+            T l2 = log((x * c) / a) * (a - b);
+            if ((l1 * a > tools::log_min_value<T>()) && (l1 * a < tools::log_max_value<T>()) && (l2 * (b - a) < tools::log_max_value<T>()) && (l2 * (b - a) > tools::log_min_value<T>()))
+            {
+               power1 = pow((x * y * c * c) / (a * b), b);
+               power2 = pow((x * c) / a, a - b);
+            }
+            else
+               need_logs = true;
+         }
+      }
+      BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
+      {
+         if (!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2))
+         {
+            need_logs = true;
+         }
       }
-      if (!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2))
+      if (need_logs)
       {
-         // We have to use logs :(
-         return prefix * exp(a * log(x * c / a) + b * log(y * c / b)) * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol));
+         //
+         // We want:
+         //
+         // (xc / a)^a (yc / b)^b
+         //
+         // But we know that one or other term will over / underflow and combining the logs will be next to useless as that will cause significant cancellation.
+         // If we assume b > a and express z ^ b as(z ^ b / a) ^ a with z = (yc / b) then we can move one power term inside the other :
+         //
+         // ((xc / a) * (yc / b)^(b / a))^a
+         //
+         // However, we're not quite there yet, as the term being exponentiated is quite likely to be close to unity, so let:
+         //
+         // xc / a = 1 + (xb - ya) / a
+         //
+         // analogously let :
+         //
+         // 1 + p = (yc / b) ^ (b / a) = 1 + expm1((b / a) * log1p((ya - xb) / b))
+         //
+         // so putting the two together we have :
+         //
+         // exp(a * log1p((xb - ya) / a + p + p(xb - ya) / a))
+         //
+         // Analogously, when a > b we can just swap all the terms around.
+         // 
+         // Finally, there are a few cases (x or y is unity) when the above logic can't be used
+         // or where there is no logarithmic cancellation and accuracy is better just using
+         // the regular formula:
+         //
+         T xc_a = x * c / a;
+         T yc_b = y * c / b;
+         if ((x == 1) || (y == 1) || (fabs(xc_a - 1) > 0.25) || (fabs(yc_b - 1) > 0.25))
+         {
+            // The above logic fails, the result is almost certainly zero:
+            power1 = exp(log(xc_a) * a + log(yc_b) * b);
+            power2 = 1;
+         }
+         else if (b > a)
+         {
+            T p = boost::math::expm1((b / a) * boost::math::log1p((y * a - x * b) / b));
+            power1 = exp(a * boost::math::log1p((x * b - y * a) / a + p * (x * c / a)));
+            power2 = 1;
+         }
+         else
+         {
+            T p = boost::math::expm1((a / b) * boost::math::log1p((x * b - y * a) / a));
+            power1 = exp(b * boost::math::log1p((y * a - x * b) / b + p * (y * c / b)));
+            power2 = 1;
+         }
       }
       return prefix * power1 * power2 * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol));
    }
@@ -568,7 +660,10 @@ T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_deriva
       T agh = static_cast<T>(a + Lanczos::g() - 0.5f);
       T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
       T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
-      result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
+      if ((a < tools::min_value<T>()) || (b < tools::min_value<T>()))
+         result = 0;  // denorms cause overflow in the Lanzos series, result will be zero anyway
+      else
+         result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));
 
       if (!(boost::math::isfinite)(result))
          result = 0;
@@ -601,10 +696,13 @@ T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_deriva
          //
          // Oh dear, we need logs, and this *will* cancel:
          //
-         result = log(result) + l1 + l2 + (log(agh) - 1) / 2;
-         if(p_derivative)
-            *p_derivative = exp(result + b * log(y));
-         result = exp(result);
+         if (result != 0)  // elude calculation when result will be zero.
+         {
+            result = log(result) + l1 + l2 + (log(agh) - 1) / 2;
+            if (p_derivative)
+               *p_derivative = exp(result + b * log(y));
+            result = exp(result);
+         }
       }
    }
    else
@@ -1021,18 +1119,15 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de
    BOOST_MATH_ASSERT((p_derivative == 0) || normalised);
 
    if(!(boost::math::isfinite)(a))
-      return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
+      return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol);
    if(!(boost::math::isfinite)(b))
-      return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
-   if(!(boost::math::isfinite)(x))
+      return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol);
+   if (!(0 <= x && x <= 1))
       return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol);
 
    if(p_derivative)
       *p_derivative = -1; // value not set.
 
-   if((x < 0) || (x > 1))
-      return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol);
-
    if(normalised)
    {
       if(a < 0)
@@ -1422,18 +1517,16 @@ T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
    // start with the usual error checks:
    //
    if (!(boost::math::isfinite)(a))
-      return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
+      return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol);
    if (!(boost::math::isfinite)(b))
-      return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
-   if (!(boost::math::isfinite)(x))
+      return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol);
+   if (!(0 <= x && x <= 1))
       return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol);
 
    if(a <= 0)
       return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol);
    if(b <= 0)
       return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol);
-   if((x < 0) || (x > 1))
-      return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol);
    //
    // Now the corner cases:
    //
diff --git a/contrib/restricted/boost/math/include/boost/math/special_functions/detail/ibeta_inverse.hpp b/contrib/restricted/boost/math/include/boost/math/special_functions/detail/ibeta_inverse.hpp
index be4bd95399..98cc1d7485 100644
--- a/contrib/restricted/boost/math/include/boost/math/special_functions/detail/ibeta_inverse.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/special_functions/detail/ibeta_inverse.hpp
@@ -25,23 +25,16 @@ namespace boost{ namespace math{ namespace detail{
 template <class T>
 struct temme_root_finder
 {
-   temme_root_finder(const T t_, const T a_) : t(t_), a(a_) {}
+   temme_root_finder(const T t_, const T a_) : t(t_), a(a_) {
+      const T x_extrema = 1 / (1 + a);
+      BOOST_MATH_ASSERT(0 < x_extrema && x_extrema < 1);
+   }
 
    boost::math::tuple<T, T> operator()(T x)
    {
       BOOST_MATH_STD_USING // ADL of std names
 
       T y = 1 - x;
-      if(y == 0)
-      {
-         T big = tools::max_value<T>() / 4;
-         return boost::math::make_tuple(static_cast<T>(-big), static_cast<T>(-big));
-      }
-      if(x == 0)
-      {
-         T big = tools::max_value<T>() / 4;
-         return boost::math::make_tuple(static_cast<T>(-big), big);
-      }
       T f = log(x) + a * log(y) + t;
       T f1 = (1 / x) - (a / (y));
       return boost::math::make_tuple(f, f1);
@@ -410,6 +403,10 @@ T temme_method_3_ibeta_inverse(T a, T b, T p, T q, const Policy& pol)
    T lower = eta < mu ? cross : 0;
    T upper = eta < mu ? 1 : cross;
    T x = (lower + upper) / 2;
+
+   // Early exit for cases with numerical precision issues.
+   if (cross == 0 || cross == 1) { return cross; }
+   
    x = tools::newton_raphson_iterate(
       temme_root_finder<T>(u, mu), x, lower, upper, policies::digits<T, Policy>() / 2);
 #ifdef BOOST_INSTRUMENT
@@ -828,7 +825,21 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
          std::swap(p, q);
          invert = !invert;
       }
-      if(pow(p, 1/a) < 0.5)
+      if (a < tools::min_value<T>())
+      {
+         // Avoid spurious overflows for denorms:
+         if (p < 1)
+         {
+            x = 1;
+            y = 0;
+         }
+         else
+         {
+            x = 0;
+            y = 1;
+         }
+      }
+      else if(pow(p, 1/a) < 0.5)
       {
 #ifndef BOOST_NO_EXCEPTIONS
          try 
diff --git a/contrib/restricted/boost/math/include/boost/math/special_functions/next.hpp b/contrib/restricted/boost/math/include/boost/math/special_functions/next.hpp
index 07c55e8d4a..1b56932423 100644
--- a/contrib/restricted/boost/math/include/boost/math/special_functions/next.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/special_functions/next.hpp
@@ -82,8 +82,8 @@ inline T normalize_value(const T& val, const std::true_type&)
 }
 
 template <class T>
-inline T get_smallest_value(std::true_type const&)
-{
+inline T get_smallest_value(std::true_type const&) {
+   static_assert(std::numeric_limits<T>::is_specialized, "Type T must be specialized.");
    //
    // numeric_limits lies about denorms being present - particularly
    // when this can be turned on or off at runtime, as is the case
@@ -106,11 +106,12 @@ inline T get_smallest_value(std::false_type const&)
 template <class T>
 inline T get_smallest_value()
 {
-#if defined(BOOST_MSVC) && (BOOST_MSVC <= 1310)
-   return get_smallest_value<T>(std::integral_constant<bool, std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == 1)>());
-#else
-   return get_smallest_value<T>(std::integral_constant<bool, std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::has_denorm == std::denorm_present)>());
-#endif
+   return get_smallest_value<T>(std::integral_constant<bool, std::numeric_limits<T>::is_specialized>());
+}
+
+template <class T>
+inline bool has_denorm_now() {
+   return get_smallest_value<T>() < tools::min_value<T>();
 }
 
 //
diff --git a/contrib/restricted/boost/math/include/boost/math/special_functions/round.hpp b/contrib/restricted/boost/math/include/boost/math/special_functions/round.hpp
index b32c8e4347..e74acba85b 100644
--- a/contrib/restricted/boost/math/include/boost/math/special_functions/round.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/special_functions/round.hpp
@@ -12,6 +12,7 @@
 #endif
 
 #include <boost/math/tools/config.hpp>
+#include <boost/math/ccmath/detail/config.hpp>
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/math/special_functions/math_fwd.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
@@ -19,11 +20,9 @@
 #include <limits>
 #include <cmath>
 
-#ifndef BOOST_NO_CXX17_IF_CONSTEXPR
+#if !defined(BOOST_MATH_NO_CCMATH) && !defined(BOOST_MATH_NO_CONSTEXPR_DETECTION)
 #include <boost/math/ccmath/ldexp.hpp>
-#  if !defined(BOOST_MATH_NO_CONSTEXPR_DETECTION)
 #    define BOOST_MATH_HAS_CONSTEXPR_LDEXP
-#  endif
 #endif
 
 namespace boost{ namespace math{
diff --git a/contrib/restricted/boost/math/include/boost/math/special_functions/trunc.hpp b/contrib/restricted/boost/math/include/boost/math/special_functions/trunc.hpp
index 7ac33a92d3..a084de560b 100644
--- a/contrib/restricted/boost/math/include/boost/math/special_functions/trunc.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/special_functions/trunc.hpp
@@ -14,14 +14,14 @@
 #include <type_traits>
 #include <boost/math/special_functions/math_fwd.hpp>
 #include <boost/math/tools/config.hpp>
+#include <boost/math/ccmath/detail/config.hpp>
 #include <boost/math/policies/error_handling.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/tools/is_constant_evaluated.hpp>
 
-#ifndef BOOST_NO_CXX17_IF_CONSTEXPR
+#if !defined(BOOST_MATH_NO_CCMATH) && !defined(BOOST_MATH_NO_CONSTEXPR_DETECTION)
 #include <boost/math/ccmath/ldexp.hpp>
-#  if !defined(BOOST_MATH_NO_CONSTEXPR_DETECTION)
 #    define BOOST_MATH_HAS_CONSTEXPR_LDEXP
-#  endif
 #endif
 
 namespace boost{ namespace math{ namespace detail{
diff --git a/contrib/restricted/boost/math/include/boost/math/tools/big_constant.hpp b/contrib/restricted/boost/math/include/boost/math/tools/big_constant.hpp
index bec93a1165..b26d8cf86d 100644
--- a/contrib/restricted/boost/math/include/boost/math/tools/big_constant.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/tools/big_constant.hpp
@@ -60,7 +60,7 @@ inline T make_big_value(largest_float, const char* s, std::false_type const&, st
 }
 #else
 template <typename T>
-inline T make_big_value(largest_float, const char* s, std::false_type const&, std::false_type const&)
+inline T make_big_value(largest_float, const char*, std::false_type const&, std::false_type const&)
 {
    static_assert(sizeof(T) == 0, "Type is unsupported in standalone mode. Please disable and try again.");
 }
diff --git a/contrib/restricted/boost/math/include/boost/math/tools/config.hpp b/contrib/restricted/boost/math/include/boost/math/tools/config.hpp
index 4b4a8a66ba..584cbd5485 100644
--- a/contrib/restricted/boost/math/include/boost/math/tools/config.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/tools/config.hpp
@@ -146,34 +146,17 @@
 #  define BOOST_MATH_EXEC_COMPATIBLE
 #endif
 
-// Attributes from C++14 and newer
-#ifdef __has_cpp_attribute
-
-// C++17
-#if (__cplusplus >= 201703L || _MSVC_LANG >= 201703L)
-#  if __has_cpp_attribute(maybe_unused)
-#    define BOOST_MATH_MAYBE_UNUSED [[maybe_unused]]
-#  endif
-#endif
-
-#endif // Standalone config
-
-// If attributes are not defined make sure we don't have compiler errors
-#ifndef BOOST_MATH_MAYBE_UNUSED
-#  define BOOST_MATH_MAYBE_UNUSED 
-#endif
-
 // C++23
 #if __cplusplus > 202002L || _MSVC_LANG > 202002L
 #  if __GNUC__ >= 13
      // libstdc++3 only defines to/from_chars for std::float128_t when one of these defines are set
      // otherwise we're right out of luck...
 #    if defined(_GLIBCXX_LDOUBLE_IS_IEEE_BINARY128) || defined(_GLIBCXX_HAVE_FLOAT128_MATH)
-#    include <cstring> // std::strlen is used with from_chars
-#    include <charconv>
-#    error #include <stdfloat>
-#    define BOOST_MATH_USE_CHARCONV_FOR_CONVERSION
-#endif
+#      include <cstring> // std::strlen is used with from_chars
+#      include <charconv>
+#      error #include <stdfloat>
+#      define BOOST_MATH_USE_CHARCONV_FOR_CONVERSION
+#    endif
 #  endif
 #endif
 
diff --git a/contrib/restricted/boost/math/include/boost/math/tools/promotion.hpp b/contrib/restricted/boost/math/include/boost/math/tools/promotion.hpp
index 80b3d89a10..263f3a6205 100644
--- a/contrib/restricted/boost/math/include/boost/math/tools/promotion.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/tools/promotion.hpp
@@ -26,8 +26,12 @@
 #include <boost/math/tools/config.hpp>
 #include <type_traits>
 
-#if __has_include(<stdfloat>)
-#  error #include <stdfloat>
+#if defined __has_include
+#  if __cplusplus > 202002L || _MSVC_LANG > 202002L 
+#    if __has_include (<stdfloat>)
+#    error #include <stdfloat>
+#    endif
+#  endif
 #endif
 
 namespace boost
diff --git a/contrib/restricted/boost/math/include/boost/math/tools/roots.hpp b/contrib/restricted/boost/math/include/boost/math/tools/roots.hpp
index 9fd7bee488..af46e2fcdd 100644
--- a/contrib/restricted/boost/math/include/boost/math/tools/roots.hpp
+++ b/contrib/restricted/boost/math/include/boost/math/tools/roots.hpp
@@ -275,14 +275,7 @@ T newton_raphson_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t&
       if (fabs(delta * 2) > fabs(delta2))
       {
          // Last two steps haven't converged.
-         T shift = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
-         if ((result != 0) && (fabs(shift) > fabs(result)))
-         {
-            delta = sign(delta) * fabs(result) * 1.1f; // Protect against huge jumps!
-            //delta = sign(delta) * result; // Protect against huge jumps! Failed for negative result. https://github.com/boostorg/math/issues/216
-         }
-         else
-            delta = shift;
+         delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
          // reset delta1/2 so we don't take this branch next time round:
          delta1 = 3 * delta;
          delta2 = 3 * delta;
@@ -596,7 +589,8 @@ namespace detail {
    #ifdef BOOST_MATH_INSTRUMENT
          std::cout << "Second order root iteration, delta = " << delta << ", residual = " << f0 << "\n";
    #endif
-         T convergence = fabs(delta / delta2);
+         // We need to avoid delta/delta2 overflowing here:
+         T convergence = (fabs(delta2) > 1) || (fabs(tools::max_value<T>() * delta2) > fabs(delta)) ? fabs(delta / delta2) : tools::max_value<T>();
          if ((convergence > 0.8) && (convergence < 2))
          {
             // last two steps haven't converged.
diff --git a/contrib/restricted/boost/math/ya.make b/contrib/restricted/boost/math/ya.make
index a692f03ae5..660c3ca286 100644
--- a/contrib/restricted/boost/math/ya.make
+++ b/contrib/restricted/boost/math/ya.make
@@ -6,9 +6,9 @@ LICENSE(BSL-1.0)
 
 LICENSE_TEXTS(.yandex_meta/licenses.list.txt)
 
-VERSION(1.83.0)
+VERSION(1.84.0)
 
-ORIGINAL_SOURCE(https://github.com/boostorg/math/archive/boost-1.83.0.tar.gz)
+ORIGINAL_SOURCE(https://github.com/boostorg/math/archive/boost-1.84.0.tar.gz)
 
 PEERDIR(
     contrib/restricted/boost/assert