Merge pull request #978 from boostorg/cpp23-float

C++23 Floats

.drone.star

@@ -16,13 +16,15 @@ windowsglobalimage="cppalliance/dronevs2019"
 def main(ctx):
 
   things_to_test = [ "special_fun", "distribution_tests", "mp", "misc", "interpolators", "quadrature", "autodiff", "long-running-tests", "float128_tests" ]
+  gcc13_things_to_test = [ "special_fun", "distribution_tests", "mp", "misc", "interpolators", "quadrature", "autodiff", "long-running-tests", "float128_tests", "new_floats" ]
   sanitizer_test = [ "special_fun", "distribution_tests", "misc", "interpolators", "quadrature", "float128_tests" ]
   gnu_5_stds = [ "gnu++14", "c++14" ]
   gnu_6_stds = [ "gnu++14", "c++14", "gnu++17", "c++17" ]
   clang_6_stds = [ "c++14", "c++17" ]
   gnu_9_stds = [ "gnu++14", "c++14", "gnu++17", "c++17", "gnu++2a", "c++2a" ]
   clang_10_stds = [ "c++14", "c++17", "c++2a" ]
   gnu_non_native = [ "gnu++17" ]
+  gcc13_stds = [ "c++23" ]
 
   result = []
 
@@ -59,6 +61,9 @@ def main(ctx):
       result.append(linux_cxx("Ubuntu g++ ARM64" + cxx + " " + suite, "g++", packages="g++", buildtype="boost", image="cppalliance/droneubuntu2204:multiarch", arch="arm64", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv))
     for cxx in gnu_non_native:
       result.append(osx_cxx("M1 Clang " + cxx + " " + suite, "clang++", buildscript="drone", buildtype="boost", xcode_version="14.1", environment={'TOOLSET': 'clang', 'CXXSTD': cxx, 'TEST_SUITE': suite, 'DEFINE': 'BOOST_MATH_NO_REAL_CONCEPT_TESTS,BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS,BOOST_MATH_MULTI_ARCH_CI_RUN', }, globalenv=globalenv))
+  for suite in gcc13_things_to_test:
+    for cxx in gcc13_stds:
+      result.append(linux_cxx("Ubuntu g++-13 " + cxx + " " + suite, "g++-13", packages="g++-13", buildtype="boost", image="cppalliance/droneubuntu2304:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-13', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv))
 
   return result
 

include/boost/math/concepts/real_concept.hpp

@@ -45,6 +45,10 @@
 #  include <cstdio>
 #endif
 
+#if __has_include(<stdfloat>)
+#  include <stdfloat>
+#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<real_concept_base_type>(c)){}
+#endif
+#ifdef __STDCPP_FLOAT64_T__
+   real_concept(std::float64_t c) : m_value(static_cast<real_concept_base_type>(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; }

include/boost/math/differentiation/autodiff.hpp

@@ -1491,7 +1491,7 @@ fvar<RealType, Order> sqrt(fvar<RealType, Order> const& cr) {
   BOOST_IF_CONSTEXPR (order == 0)
     return fvar<RealType, Order>(*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;

include/boost/math/differentiation/finite_difference.hpp

@@ -153,7 +153,7 @@ namespace detail {
 
       const Real eps = (numeric_limits<Real>::epsilon)();
       // Error bound ~eps^4/5
-      Real h = pow(11.25*eps, static_cast<Real>(1) / static_cast<Real>(5));
+      Real h = pow(Real(11.25)*eps, static_cast<Real>(1) / static_cast<Real>(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<Real>(1) / static_cast<Real>(9));
+      Real h = pow(Real(551.25)*eps, static_cast<Real>(1) / static_cast<Real>(9));
       h = detail::make_xph_representable(x, h);
 
       Real yh = f(x + h);

include/boost/math/differentiation/lanczos_smoothing.hpp

@@ -10,6 +10,7 @@
 #include <limits> // to nan initialize
 #include <vector>
 #include <string>
+#include <cstdint>
 #include <stdexcept>
 #include <type_traits>
 #include <boost/math/tools/assert.hpp>

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<RealType>(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<RealType>(ibeta_invb(alpha, x, probability, Policy()));
       } //  RealType find_beta(alpha, x, probability)
 
     private:
@@ -396,7 +397,7 @@ namespace boost
       {
         if (a == 1)
         {
-          return 1 / beta(a, b);
+          return static_cast<RealType>(1 / beta(a, b));
         }
         else if (a < 1)
         {
@@ -411,7 +412,7 @@ namespace boost
       {
         if (b == 1)
         {
-          return 1 / beta(a, b);
+          return static_cast<RealType>(1 / beta(a, b));
         }
         else if (b < 1)
         {
@@ -423,7 +424,7 @@ namespace boost
         }
       }
 
-      return ibeta_derivative(a, b, x, Policy());
+      return static_cast<RealType>(ibeta_derivative(a, b, x, Policy()));
     } // pdf
 
     template <class RealType, class Policy>
@@ -454,7 +455,7 @@ namespace boost
       {
         return 1;
       }
-      return ibeta(a, b, x, Policy());
+      return static_cast<RealType>(ibeta(a, b, x, Policy()));
     } // beta cdf
 
     template <class RealType, class Policy>
@@ -481,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<RealType>(ibetac(a, b, x, Policy()));
     } // beta cdf
 
     template <class RealType, class Policy>
@@ -519,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<RealType*>(nullptr), Policy());
+      return static_cast<RealType>(ibeta_inv(a, b, p, static_cast<RealType*>(nullptr), Policy()));
     } // quantile
 
     template <class RealType, class Policy>
@@ -555,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<RealType*>(nullptr), Policy());
+      return static_cast<RealType>(ibetac_inv(a, b, q, static_cast<RealType*>(nullptr), Policy()));
     } // Quantile Complement
 
   } // namespace math

include/boost/math/distributions/cauchy.hpp

@@ -81,7 +81,7 @@ RealType cdf_imp(const cauchy_distribution<RealType, Policy>& dist, const RealTy
    RealType mx = -fabs((x - location) / scale); // scale is > 0
    if(mx > -tools::epsilon<RealType>() / 8)
    {  // special case first: x extremely close to location.
-      return 0.5;
+      return static_cast<RealType>(0.5f);
    }
    result = -atan(1 / mx) / constants::pi<RealType>();
    return (((x > location) != complement) ? 1 - result : result);

include/boost/math/distributions/detail/hypergeometric_pdf.hpp

@@ -262,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<T>(pow(x, ex));
 #endif
 }
 template <class T>

include/boost/math/distributions/inverse_gaussian.hpp

@@ -382,7 +382,7 @@ inline RealType quantile(const inverse_gaussian_distribution<RealType, Policy>&
   RealType guess = detail::guess_ig(p, dist.mean(), dist.scale());
   using boost::math::tools::max_value;
 
-  RealType min = 0.; // Minimum possible value is bottom of range of distribution.
+  RealType min = static_cast<RealType>(0); // Minimum possible value is bottom of range of distribution.
   RealType max = max_value<RealType>();// Maximum possible value is top of range. 
   // int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.
   // digits used to control how accurate to try to make the result.
@@ -454,7 +454,7 @@ inline RealType quantile(const complemented2_type<inverse_gaussian_distribution<
    // Complement.
    using boost::math::tools::max_value;
 
-  RealType min = 0.; // Minimum possible value is bottom of range of distribution.
+  RealType min = static_cast<RealType>(0); // Minimum possible value is bottom of range of distribution.
   RealType max = max_value<RealType>();// Maximum possible value is top of range. 
   // int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.
   // digits used to control how accurate to try to make the result.

include/boost/math/distributions/skew_normal.hpp

@@ -486,7 +486,7 @@ namespace boost{ namespace math{
 
     // 21 elements
     static const RealType shapes[] = {
-      0.0,
+      static_cast<RealType>(0.0),
       static_cast<RealType>(1.000000000000000e-004),
       static_cast<RealType>(2.069138081114790e-004),
       static_cast<RealType>(4.281332398719396e-004),
@@ -511,7 +511,7 @@ namespace boost{ namespace math{
 
     // 21 elements
     static const RealType guess[] = {
-      0.0,
+      static_cast<RealType>(0.0),
       static_cast<RealType>(5.000050000525391e-005),
       static_cast<RealType>(1.500015000148736e-004),
       static_cast<RealType>(3.500035000350010e-004),

include/boost/math/interpolators/detail/cardinal_cubic_b_spline_detail.hpp

@@ -225,7 +225,7 @@ cardinal_cubic_b_spline_imp<Real>::cardinal_cubic_b_spline_imp(BidiIterator f, B
     // mapsto
     // 1 0 -1 | r0
     // 0 1 1/2| (r1 - r0)/4
-    super_diagonal[1] = 0.5;
+    super_diagonal[1] = static_cast<Real>(0.5);
     rhs[1] = (rhs[1] - rhs[0])/4;
 
     // Now do a tridiagonal row reduction the standard way, until just before the last row:

include/boost/math/policies/error_handling.hpp

@@ -82,7 +82,7 @@ namespace detail
 {
 
 template <class T>
-std::string prec_format(const T& val)
+inline std::string prec_format(const T& val)
 {
    typedef typename boost::math::policies::precision<T, boost::math::policies::policy<> >::type prec_type;
    std::stringstream ss;
@@ -95,6 +95,18 @@ std::string prec_format(const T& val)
    return ss.str();
 }
 
+#ifdef BOOST_MATH_USE_CHARCONV_FOR_CONVERSION
+
+template <>
+inline std::string prec_format<std::float128_t>(const std::float128_t& val)
+{
+   char buffer[128] {};
+   const auto r = std::to_chars(buffer, buffer + sizeof(buffer), val);
+   return std::string(buffer, r.ptr);
+}
+
+#endif
+
 inline void replace_all_in_string(std::string& result, const char* what, const char* with)
 {
    std::string::size_type pos = 0;

include/boost/math/quadrature/naive_monte_carlo.hpp

@@ -21,6 +21,11 @@
 #include <map>
 #include <type_traits>
 #include <boost/math/policies/error_handling.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
+
+#ifdef BOOST_NAIVE_MONTE_CARLO_DEBUG_FAILURES
+#  include <iostream>
+#endif
 
 namespace boost { namespace math { namespace quadrature {
 
@@ -45,6 +50,8 @@ class naive_monte_carlo
     {
         using std::numeric_limits;
         using std::sqrt;
+        using boost::math::isinf;
+
         uint64_t n = bounds.size();
         m_lbs.resize(n);
         m_dxs.resize(n);
@@ -58,9 +65,9 @@ class naive_monte_carlo
                 boost::math::policies::raise_domain_error(function, "The upper bound is <= the lower bound.\n", bounds[i].second, Policy());
                 return;
             }
-            if (bounds[i].first == -numeric_limits<Real>::infinity())
+            if (isinf(bounds[i].first))
             {
-                if (bounds[i].second == numeric_limits<Real>::infinity())
+                if (isinf(bounds[i].second))
                 {
                     m_limit_types[i] = detail::limit_classification::DOUBLE_INFINITE;
                 }
@@ -72,7 +79,7 @@ class naive_monte_carlo
                     m_dxs[i] = numeric_limits<Real>::quiet_NaN();
                 }
             }
-            else if (bounds[i].second == numeric_limits<Real>::infinity())
+            else if (isinf(bounds[i].second))
             {
                 m_limit_types[i] = detail::limit_classification::UPPER_BOUND_INFINITE;
                 if (singular)
@@ -285,17 +292,17 @@ class naive_monte_carlo
             m_done = false;
 
 #ifdef BOOST_NAIVE_MONTE_CARLO_DEBUG_FAILURES
-            std::cout << "Failed to achieve required tolerance first time through..\n";
-            std::cout << "  variance =    " << m_variance << std::endl;
-            std::cout << "  average =     " << m_avg << std::endl;
-            std::cout << "  total calls = " << m_total_calls << std::endl;
+            std::cerr << "Failed to achieve required tolerance first time through..\n";
+            std::cerr << "  variance =    " << m_variance << std::endl;
+            std::cerr << "  average =     " << m_avg << std::endl;
+            std::cerr << "  total calls = " << m_total_calls << std::endl;
 
             for (std::size_t i = 0; i < m_num_threads; ++i)
-               std::cout << "  thread_calls[" << i << "] = " << m_thread_calls[i] << std::endl;
+               std::cerr << "  thread_calls[" << i << "] = " << m_thread_calls[i] << std::endl;
             for (std::size_t i = 0; i < m_num_threads; ++i)
-               std::cout << "  thread_averages[" << i << "] = " << m_thread_averages[i] << std::endl;
+               std::cerr << "  thread_averages[" << i << "] = " << m_thread_averages[i] << std::endl;
             for (std::size_t i = 0; i < m_num_threads; ++i)
-               std::cout << "  thread_Ss[" << i << "] = " << m_thread_Ss[i] << std::endl;
+               std::cerr << "  thread_Ss[" << i << "] = " << m_thread_Ss[i] << std::endl;
 #endif
          }
 

include/boost/math/special_functions/bessel_prime.hpp

@@ -38,9 +38,9 @@ inline T cyl_bessel_j_prime_imp(T v, T x, const Policy& pol)
    if (x == 0)
    {
       if (v == 1)
-         return 0.5;
+         return static_cast<T>(0.5);
       else if (v == -1)
-         return -0.5;
+         return static_cast<T>(-0.5);
       else if (floor(v) == v || v > 1)
          return 0;
       else return boost::math::policies::raise_domain_error<T>(
@@ -126,7 +126,7 @@ inline T cyl_bessel_i_prime_imp(T v, T x, const Policy& pol)
    if (x == 0)
    {
       if (v == 1 || v == -1)
-         return 0.5;
+         return static_cast<T>(0.5);
       else if (floor(v) == v || v > 1)
          return 0;
       else return boost::math::policies::raise_domain_error<T>(