Study mul256

doc/decimal/numbers.adoc

@@ -1,5 +1,5 @@
 ////
-Copyright 2023 Matt Borland
+Copyright 2023 - 2024 Matt Borland
 Distributed under the Boost Software License, Version 1.0.
 https://www.boost.org/LICENSE_1_0.txt
 ////
@@ -10,7 +10,7 @@ https://www.boost.org/LICENSE_1_0.txt
 
 == Overview
 
-Contains all constants provided by C+\+20's `<numbers>` specialized for the decimal types. These does not require C++20.
+Contains all constants provided by C+\+20's `<numbers>` specialized for the decimal types. These do not require C++20.
 
 - e_v - https://en.wikipedia.org/wiki/E_(mathematical_constant)[Euler's Number]
 - log2e_v - log2(e)
@@ -22,11 +22,14 @@ Contains all constants provided by C+\+20's `<numbers>` specialized for the deci
 - ln10_v - ln(10)
 - sqrt2_v - sqrt(2)
 - sqrt3_v - sqrt(3)
+- sqrt10_v - sqrt(10)
 - inv_sqrt3_v - 1 / sqrt(3)
+- cbrt2_v - cbrt(2)
+- cbrt10_v - cbrt(10)
 - egamma_v - https://en.wikipedia.org/wiki/Euler%27s_constant[Euler–Mascheroni constant]
 - phi_v - https://en.wikipedia.org/wiki/Golden_ratio[The Golden Ratio]
 
-There are also non-template variables that provide the constant as a decimal32 type.
+There are also non-template variables that provide the constant as a decimal64 type.
 
 == Reference
 
@@ -68,32 +71,44 @@ static constexpr Decimal sqrt2_v;
 template <typename Decimal>
 static constexpr Decimal sqrt3_v;
 
+template <typename Decimal>
+static constexpr Decimal sqrt10_v;
+
 template <typename Decimal>
 static constexpr Decimal inv_sqrt2_v;
 
 template <typename Decimal>
 static constexpr Decimal inv_sqrt3_v;
 
+template <typename Decimal>
+static constexpr Decimal cbrt2_v;
+
+template <typename Decimal>
+static constexpr Decimal cbrt10_v;
+
 template <typename Decimal>
 static constexpr Decimal egamma_v;
 
 template <typename Decimal>
 static constexpr Decimal phi_v;
 
-static constexpr auto e {e_v<decimal32>};
-static constexpr auto log2e {log2e_v<decimal32>};
-static constexpr auto log10e {log10e_v<decimal32>};
-static constexpr auto pi {pi_v<decimal32>};
-static constexpr auto inv_pi {inv_pi_v<decimal32>};
-static constexpr auto inv_sqrtpi {inv_sqrtpi_v<decimal32>};
-static constexpr auto ln2 {ln2_v<decimal32>};
-static constexpr auto ln10 {ln10_v<decimal32>};
-static constexpr auto sqrt2 {sqrt2_v<decimal32>};
-static constexpr auto sqrt3 {sqrt3_v<decimal32>};
-static constexpr auto inv_sqrt2 {inv_sqrt2_v<decimal32>};
-static constexpr auto inv_sqrt3 {inv_sqrt3_v<decimal32>};
-static constexpr auto egamma {egamma_v<decimal32>};
-static constexpr auto phi {phi_v<decimal32>};
+static constexpr auto e {e_v<decimal64>};
+static constexpr auto log2e {log2e_v<decimal64>};
+static constexpr auto log10e {log10e_v<decimal64>};
+static constexpr auto pi {pi_v<decimal64>};
+static constexpr auto inv_pi {inv_pi_v<decimal64>};
+static constexpr auto inv_sqrtpi {inv_sqrtpi_v<decimal64>};
+static constexpr auto ln2 {ln2_v<decimal64>};
+static constexpr auto ln10 {ln10_v<decimal64>};
+static constexpr auto sqrt2 {sqrt2_v<decimal64>};
+static constexpr auto sqrt3 {sqrt3_v<decimal64>};
+static constexpr auto sqrt10 {sqrt10_v<decimal64>};
+static constexpr auto inv_sqrt2 {inv_sqrt2_v<decimal64>};
+static constexpr auto inv_sqrt3 {inv_sqrt3_v<decimal64>};
+static constexpr auto cbrt2 {cbrt2_v<decimal64>};
+static constexpr auto cbrt10 {cbrt10_v<decimal64>};
+static constexpr auto egamma {egamma_v<decimal64>};
+static constexpr auto phi {phi_v<decimal64>};
 
 } //namespace decimal
 } //namespace boost

include/boost/decimal/detail/cmath/cbrt.hpp

@@ -1,5 +1,5 @@
-// Copyright 2023 Matt Borland
-// Copyright 2023 Christopher Kormanyos
+// Copyright 2023 - 2024 Matt Borland
+// Copyright 2023 - 2024 Christopher Kormanyos
 // Distributed under the Boost Software License, Version 1.0.
 // https://www.boost.org/LICENSE_1_0.txt
 
@@ -24,65 +24,149 @@ namespace decimal {
 namespace detail {
 
 template <typename T>
-constexpr auto cbrt_impl(T val) noexcept
+constexpr auto cbrt_impl(T x) noexcept
     BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T)
 {
-    constexpr T zero {0, 0};
-    constexpr T one {1, 0};
+    const auto fpc = fpclassify(x);
 
     T result { };
 
-    if (isnan(val) || abs(val) == zero)
+    if ((fpc == FP_NAN) || (fpc == FP_ZERO))
     {
-        result = val;
-    }
-    else if (isinf(val))
-    {
-        if (signbit(val))
-        {
-            result = std::numeric_limits<T>::quiet_NaN();
-        }
-        else
-        {
-            result = val;
-        }
+        result = x;
     }
-    else if (val < zero)
+    else if (signbit(x))
     {
-        result = std::numeric_limits<T>::quiet_NaN();
+        result = -cbrt(-x);
     }
-    else if (val == one)
+    else if (fpc == FP_INFINITE)
     {
-        result = one;
+        result = std::numeric_limits<T>::infinity();
     }
     else
     {
-        constexpr T epsilon = std::numeric_limits<T>::epsilon() * 100;
-        T error = one / epsilon;
+        int exp10val { };
+
+        const auto gn { frexp10(x, &exp10val) };
 
-        T x {};
-        if (val > one)
+        const auto
+            zeros_removal
+            {
+                remove_trailing_zeros(gn)
+            };
+
+        const bool is_pure { static_cast<unsigned>(zeros_removal.trimmed_number) == 1U };
+
+        if(is_pure)
         {
-            // Scale down if val is large by dividing the exp by 3
-            int exp {};
-            auto sig = frexp10(val, &exp);
-            x = T{sig, exp / 3};
+            // Here, a pure power-of-10 argument gets a straightforward result.
+            // For argument 10^n where n is a multiple of 3, the result is exact.
+
+            const int p10 { exp10val + static_cast<int>(zeros_removal.number_of_removed_zeros) };
+
+            if (p10 == 0)
+            {
+                result = T { 1 };
+            }
+            else
+            {
+                const int p10_mod3 = (p10 % 3);
+                const int p10_div3 = (p10 / 3);
+
+                result = T { 1, p10_div3 };
+
+                switch (p10_mod3)
+                {
+                    case 2:
+                        result *= numbers::cbrt10_v<T>;
+                        // fallthrough
+
+                    case 1:
+                        result *= numbers::cbrt10_v<T>;
+                        break;
+
+                    case -2:
+                        result /= numbers::cbrt10_v<T>;
+                        // fallthrough
+
+                    case -1:
+                        result /= numbers::cbrt10_v<T>;
+                        break;
+                }
+            }
         }
         else
         {
-            // Trivial heuristic
-            x = val * 2;
-        }
-
-        while (error > epsilon)
-        {
-            const T new_x {(2 * x + val / (x * x)) / 3};
-
-            error = fabs(new_x - x);
-            x = new_x;
+            // Scale the argument to the interval 1/10 <= x < 1.
+            T gx { gn, -std::numeric_limits<T>::digits10 };
+
+            exp10val += std::numeric_limits<T>::digits10;
+
+            // For this work we perform an order-2 Pade approximation of the cube-root
+            // at argument x = 1/2. This results in slightly more than 2 decimal digits
+            // of accuracy over the interval 1/10 <= x < 1.
+
+            // PadeApproximant[x^(1/3), {x, 1/2, {2, 2}}]
+            // FullSimplify[%]
+
+            // HornerForm[Numerator[Out[2]]]
+            // Results in:
+            //   5 + x (70 + 56 x)
+
+            // HornerForm[Denominator[Out[2]]]
+            // Results in:
+            //   2^(1/3) (14 + x (70 + 20 x))
+
+            constexpr T five     {  5 };
+            constexpr T fourteen { 14 };
+            constexpr T seventy  {  7, 1 };
+
+            result =
+                  (five + gx * (seventy + gx * 56))
+                / (numbers::cbrt2_v<T> * (fourteen + gx * (seventy + gx * 20)));
+
+            // Perform 2, 3 or 4 Newton-Raphson iterations depending on precision.
+            // Note from above, we start with slightly more than 2 decimal digits
+            // of accuracy.
+
+            constexpr int iter_loops
+            {
+                  std::numeric_limits<T>::digits10 < 10 ? 2
+                : std::numeric_limits<T>::digits10 < 20 ? 3 : 4
+            };
+
+            for (int idx = 0; idx < iter_loops; ++idx)
+            {
+                result = ((result + result) + gx / (result * result)) / 3;
+            }
+
+            if (exp10val != 0)
+            {
+                const int exp10val_mod3 = (exp10val % 3);
+                const int exp10val_div3 = (exp10val / 3);
+
+                result *= T { 1, exp10val_div3 };
+
+                switch (exp10val_mod3)
+                {
+                    case 2:
+                        result *= numbers::cbrt10_v<T>;
+                        // fallthrough
+
+                    case 1:
+                        result *= numbers::cbrt10_v<T>;
+                        break;
+
+                    case -2:
+                        result /= numbers::cbrt10_v<T>;
+                        // fallthrough
+
+                    case -1:
+                        result /= numbers::cbrt10_v<T>;
+                        break;
+                }
+            }
         }
-
-        result = x;
     }
 
     return result;