Add cpp_double_fp and exercise arithmetic_tests

.github/workflows/multiprecision.yml

@@ -149,7 +149,7 @@ jobs:
       matrix:
         compiler: [ g++-10 ]
         standard: [ c++17, c++2a ]
-        suite: [ specfun_mpfr, specfun_gmp, specfun_cpp_dec_float, specfun_cpp_bin_float, specfun_float128 ]
+        suite: [ specfun_mpfr, specfun_gmp, specfun_cpp_dec_float, specfun_cpp_bin_float, specfun_float128, specfun_cpp_double_double ]
     steps:
       - uses: actions/checkout@v3
         with:
@@ -398,7 +398,7 @@ jobs:
       fail-fast: false
       matrix:
         toolset: [ gcc ]
-        standard: [ 11, 14, 17 ]
+        standard: [ 14, 17 ]
         suite: [ github_ci_block_1, github_ci_block_2 ]
     steps:
       - uses: actions/checkout@v3

example/cpp_double_double_del_v_jv.cpp

@@ -0,0 +1,76 @@
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright Christopher Kormanyos 2023.
+//  Distributed under 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 <iomanip>
+#include <iostream>
+
+#include <boost/math/special_functions/bessel.hpp>
+#include <boost/multiprecision/cpp_double_fp.hpp>
+
+namespace local
+{
+  template<typename RealValueType,
+           typename RealFunctionType>
+  RealValueType derivative(const RealValueType& x,
+                           const RealValueType& dx,
+                           RealFunctionType real_function)
+  {
+    using real_value_type = RealValueType;
+
+    // Compute the first derivative of the input function
+    // using a three-point central difference rule of O(dx^6).
+
+    const auto dx2 = static_cast<real_value_type>(dx  + dx);
+    const auto dx3 = static_cast<real_value_type>(dx2 + dx);
+
+    const auto m1 = static_cast<real_value_type>((  real_function(x + dx)
+                                                  - real_function(x - dx))  / 2U);
+    const auto m2 = static_cast<real_value_type>((  real_function(x + dx2)
+                                                  - real_function(x - dx2)) / 4U);
+    const auto m3 = static_cast<real_value_type>((  real_function(x + dx3)
+                                                  - real_function(x - dx3)) / 6U);
+
+    const auto fifteen_m1 = static_cast<real_value_type>(m1 * 15U);
+    const auto six_m2     = static_cast<real_value_type>(m2 *  6U);
+    const auto ten_dx     = static_cast<real_value_type>(dx * 10U);
+
+    return ((fifteen_m1 - six_m2) + m3) / ten_dx;
+  }
+
+  using wide_double = boost::multiprecision::cpp_double_double;
+
+  wide_double cyl_bessel_j_of_v(const wide_double& v)
+  {
+    const auto x = static_cast<wide_double>(static_cast<wide_double>(34U) / 10U);
+
+    return boost::math::cyl_bessel_j(v, x);
+  }
+} // namespace local
+
+int main()
+{
+  using local::wide_double;
+
+  // D[BesselJ[v, 34/10], v]
+  // v := 12/10
+  // N[Out[1], 30]
+  // 0.439649800900385297241807133820
+
+  const auto  v = static_cast<wide_double>(static_cast<wide_double>(12U) /     10U); // 1.2
+  const auto dv = static_cast<wide_double>(static_cast<wide_double> (1U) / 100000U); // 1E-5
+
+  const auto del_v_jv = local::derivative(v, dv, local::cyl_bessel_j_of_v);
+
+  const auto flg = std::cout.flags();
+
+  std::cout << std::fixed
+            << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<double>::max_digits10))
+            << static_cast<double>(del_v_jv)
+            << std::endl;
+
+  std::cout.flags(flg);
+}

include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail.hpp

@@ -0,0 +1,239 @@
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2021 Fahad Syed.
+//  Copyright 2021 - 2023 Christopher Kormanyos.
+//  Copyright 2021 Janek Kozicki.
+//  Distributed under 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)
+//
+
+#ifndef BOOST_MP_CPP_DF_QF_DETAIL_2023_01_02_HPP
+#define BOOST_MP_CPP_DF_QF_DETAIL_2023_01_02_HPP
+
+#include <cmath>
+#include <limits>
+#include <tuple>
+#include <utility>
+
+#include <boost/config.hpp>
+#include <boost/multiprecision/number.hpp>
+#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath.hpp>
+
+#ifdef BOOST_HAS_FLOAT128
+#include <quadmath.h>
+#endif
+
+namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail {
+
+inline float                  floor_of_constituent(float                  x) { return ::floorf(x); }
+inline double                 floor_of_constituent(double                 x) { return ::floor (x); }
+inline long double            floor_of_constituent(long double            x) { return ::floorl(x); }
+#if defined(BOOST_HAS_FLOAT128)
+inline ::boost::float128_type floor_of_constituent(::boost::float128_type x) { return ::floorq(x); }
+#endif
+
+inline float                  log_of_constituent(float                  x) { return ::logf(x); }
+inline double                 log_of_constituent(double                 x) { return ::log (x); }
+inline long double            log_of_constituent(long double            x) { return ::logl(x); }
+#if defined(BOOST_HAS_FLOAT128)
+inline ::boost::float128_type log_of_constituent(::boost::float128_type x) { return ::logq(x); }
+#endif
+
+template <class FloatingPointType>
+struct is_floating_point_or_float128
+{
+   static constexpr auto value =    std::is_same<FloatingPointType, float>::value
+                                 || std::is_same<FloatingPointType, double>::value
+                                 || std::is_same<FloatingPointType, long double>::value
+#if defined(BOOST_HAS_FLOAT128)
+                                 || std::is_same<FloatingPointType, ::boost::float128_type>::value
+#endif
+                                 ;
+};
+
+template <typename FloatingPointType>
+struct exact_arithmetic
+{
+   // The exact_arithmetic<> struct implements extended precision
+   // techniques that are used in cpp_double_fp_backend and cpp_quad_float.
+
+   static_assert(is_floating_point_or_float128<FloatingPointType>::value,
+                 "Error: exact_arithmetic<> invoked with unknown floating-point type");
+
+   using float_type  = FloatingPointType;
+   using float_pair  = std::pair<float_type, float_type>;
+   using float_tuple = std::tuple<float_type, float_type, float_type, float_type>;
+
+   static
+   #if (defined(_MSC_VER) && (_MSC_VER <= 1900))
+   BOOST_MP_CXX14_CONSTEXPR
+   #else
+   constexpr
+   #endif
+   float_pair split(const float_type& a)
+   {
+      // Split a floating point number in two (high and low) parts approximating the
+      // upper-half and lower-half bits of the float
+
+      static_assert(is_floating_point_or_float128<FloatingPointType>::value,
+                    "Error: exact_arithmetic<>::split invoked with unknown floating-point type");
+
+      // TODO Replace bit shifts with constexpr funcs or ldexp for better compaitibility
+      constexpr int MantissaBits = cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits;
+      constexpr int SplitBits    = MantissaBits / 2 + 1;
+
+      // Check if the integer is wide enough to hold the Splitter.
+      static_assert(std::numeric_limits<std::uintmax_t>::digits > SplitBits,
+                    "Inadequate integer width for binary shifting needed in split(), try using ldexp instead");
+
+      // If the above line gives an compilation error, replace the
+      // line below it with the commented line
+
+      constexpr float_type Splitter       = FloatingPointType((static_cast<std::uintmax_t>(UINT8_C(1)) << SplitBits) + 1);
+      const     float_type SplitThreshold = (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max)() / (Splitter * 2);
+
+      float_type hi { };
+      float_type lo { };
+
+      // Handle if multiplication with the splitter would cause overflow
+      if (a > SplitThreshold || a < -SplitThreshold)
+      {
+         constexpr float_type Normalizer = float_type(1ULL << (SplitBits + 1));
+
+         const float_type a_ = a / Normalizer;
+
+         const float_type temp = Splitter * a_;
+
+         hi   = temp - (temp - a_);
+         lo   = a_ - hi;
+
+         hi *= Normalizer;
+         lo *= Normalizer;
+      }
+      else
+      {
+         const float_type temp = Splitter * a;
+
+         hi   = temp - (temp - a);
+         lo   = a - hi;
+      }
+
+      return std::make_pair(hi, lo);
+   }
+
+   static
+   #if (defined(_MSC_VER) && (_MSC_VER <= 1900))
+   BOOST_MP_CXX14_CONSTEXPR
+   #else
+   constexpr
+   #endif
+   float_pair fast_sum(float_type a, float_type b)
+   {
+      // Exact addition of two floating point numbers, given |a| > |b|
+      const float_type a_plus_b = a + b;
+
+      const float_pair result(a_plus_b, b - (a_plus_b - a));
+
+      return result;
+   }
+
+   static
+   #if (defined(_MSC_VER) && (_MSC_VER <= 1900))
+   BOOST_MP_CXX14_CONSTEXPR
+   #else
+   constexpr
+   #endif
+   void sum(float_pair& result, float_type a, float_type b)
+   {
+      // Exact addition of two floating point numbers
+      const float_type a_plus_b = a + b;
+      const float_type v        = a_plus_b - a;
+
+      const float_pair tmp(a_plus_b, (a - (a_plus_b - v)) + (b - v));
+
+      result.first  = tmp.first;
+      result.second = tmp.second;
+   }
+
+   #if 0
+   static
+   #if (defined(_MSC_VER) && (_MSC_VER <= 1900))
+   BOOST_MP_CXX14_CONSTEXPR
+   #else
+   constexpr
+   #endif
+   void sum(float_pair& p, float_type& e)
+   {
+      using std::tie;
+
+      float_pair t;
+      float_type t_;
+
+      t                = sum(p.first, p.second);
+      tie(p.first, t_) = sum(e, t.first);
+      tie(p.second, e) = sum(t.second, t_);
+   }
+   #endif
+
+   static
+   #if (defined(_MSC_VER) && (_MSC_VER <= 1900))
+   BOOST_MP_CXX14_CONSTEXPR
+   #else
+   constexpr
+   #endif
+   float_pair difference(const float_type& a, const float_type& b)
+   {
+      // Exact subtraction of two floating point numbers
+      const float_type a_minus_b = a - b;
+      const float_type v         = a_minus_b - a;
+
+      return float_pair(a_minus_b, (a - (a_minus_b - v)) - (b + v));
+   }
+
+   static
+   #if (defined(_MSC_VER) && (_MSC_VER <= 1900))
+   BOOST_MP_CXX14_CONSTEXPR
+   #else
+   constexpr
+   #endif
+   float_pair product(const float_type& a, const float_type& b)
+   {
+      // Exact product of two floating point numbers
+      const float_pair a_split = split(a);
+      const float_pair b_split = split(b);
+
+      const volatile float_type pf = a * b;
+
+      return
+        std::make_pair
+        (
+          const_cast<const float_type&>(pf),
+            (
+                ((a_split.first  * b_split.first) - const_cast<const float_type&>(pf))
+              +  (a_split.first  * b_split.second)
+              +  (a_split.second * b_split.first)
+            )
+          +
+            (a_split.second * b_split.second)
+        );
+   }
+
+   static
+   #if (defined(_MSC_VER) && (_MSC_VER <= 1900))
+   BOOST_MP_CXX14_CONSTEXPR
+   #else
+   constexpr
+   #endif
+   void normalize(float_pair& result, float_type a, float_type b)
+   {
+      // Converts a pair of floats to standard form.
+      const float_pair tmp = fast_sum(a, b);
+
+      result.first  = tmp.first;
+      result.second = tmp.second;
+   }
+};
+
+} } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail
+
+#endif // BOOST_MP_CPP_DF_QF_DETAIL_2023_01_02_HPP

include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath.hpp

@@ -0,0 +1,19 @@
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright Christopher Kormanyos 2023.
+//  Distributed under 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)
+//
+
+#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_2023_01_07_HPP
+#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_2023_01_07_HPP
+
+#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_fabs.hpp>
+#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_frexp.hpp>
+#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_isinf.hpp>
+#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_isnan.hpp>
+#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_ldexp.hpp>
+#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_limits.hpp>
+#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_sqrt.hpp>
+
+#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_2023_01_07_HPP

include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_fabs.hpp

@@ -0,0 +1,26 @@
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright Christopher Kormanyos 2023.
+//  Distributed under 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)
+//
+
+#ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FABS_2023_01_07_HPP
+#define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FABS_2023_01_07_HPP
+
+#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_isnan.hpp>
+#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_limits.hpp>
+
+namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath {
+
+template <class FloatingPointType>
+constexpr FloatingPointType fabs(FloatingPointType x)
+{
+   return (cpp_df_qf_detail::ccmath::isnan(x))      ? cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::quiet_NaN() :
+          (x == static_cast<FloatingPointType>(-0)) ? static_cast<FloatingPointType>(0) :
+          (x >= 0)                                  ? x : -x;
+};
+
+} } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath
+
+#endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_FABS_2023_01_07_HPP