Merge pull request #523 from cppalliance/FMA

Author: mborland

include/boost/decimal/decimal128.hpp

@@ -1937,21 +1937,15 @@ constexpr auto operator*(decimal128 lhs, decimal128 rhs) noexcept -> decimal128
 
     auto lhs_sig {lhs.full_significand()};
     auto lhs_exp {lhs.biased_exponent()};
-
-    while (lhs_sig % 10 == 0 && lhs_sig != 0)
-    {
-        lhs_sig /= 10;
-        ++lhs_exp;
-    }
+    const auto lhs_zeros {detail::remove_trailing_zeros(lhs_sig)};
+    lhs_sig = lhs_zeros.trimmed_number;
+    lhs_exp += static_cast<std::int32_t>(lhs_zeros.number_of_removed_zeros);
 
     auto rhs_sig {rhs.full_significand()};
     auto rhs_exp {rhs.biased_exponent()};
-
-    while (rhs_sig % 10 == 0 && rhs_sig != 0)
-    {
-        rhs_sig /= 10;
-        ++rhs_exp;
-    }
+    const auto rhs_zeros {detail::remove_trailing_zeros(rhs_sig)};
+    rhs_sig = rhs_zeros.trimmed_number;
+    rhs_exp += static_cast<std::int32_t>(rhs_zeros.number_of_removed_zeros);
 
     const auto result {d128_mul_impl(lhs_sig, lhs_exp, lhs.isneg(),
                                      rhs_sig, rhs_exp, rhs.isneg())};
@@ -1970,20 +1964,16 @@ constexpr auto operator*(decimal128 lhs, Integer rhs) noexcept
 
     auto lhs_sig {lhs.full_significand()};
     auto lhs_exp {lhs.biased_exponent()};
-    while (lhs_sig % 10 == 0 && lhs_sig != 0)
-    {
-        lhs_sig /= 10;
-        ++lhs_exp;
-    }
+    const auto lhs_zeros {detail::remove_trailing_zeros(lhs_sig)};
+    lhs_sig = lhs_zeros.trimmed_number;
+    lhs_exp += static_cast<std::int32_t>(lhs_zeros.number_of_removed_zeros);
     auto lhs_components {detail::decimal128_components{lhs_sig, lhs_exp, lhs.isneg()}};
 
     auto rhs_sig {static_cast<detail::uint128>(detail::make_positive_unsigned(rhs))};
     std::int32_t rhs_exp {0};
-    while (rhs_sig % 10 == 0 && rhs_sig != 0)
-    {
-        rhs_sig /= 10;
-        ++rhs_exp;
-    }
+    const auto rhs_zeros {detail::remove_trailing_zeros(rhs_sig)};
+    rhs_sig = rhs_zeros.trimmed_number;
+    rhs_exp += static_cast<std::int32_t>(rhs_zeros.number_of_removed_zeros);
     auto unsigned_sig_rhs {detail::make_positive_unsigned(rhs_sig)};
     auto rhs_components {detail::decimal128_components{unsigned_sig_rhs, rhs_exp, (rhs < 0)}};
 
@@ -2433,25 +2423,18 @@ constexpr auto fmad128(decimal128 x, decimal128 y, decimal128 z) noexcept -> dec
 
     auto sig_lhs {x.full_significand()};
     auto exp_lhs {x.biased_exponent()};
-
-    while (sig_lhs % 10 == 0 && sig_lhs != 0)
-    {
-        sig_lhs /= 10;
-        ++exp_lhs;
-    }
+    detail::normalize<decimal128>(sig_lhs, exp_lhs);
 
     auto sig_rhs {y.full_significand()};
     auto exp_rhs {y.biased_exponent()};
-
-    while (sig_rhs % 10 == 0 && sig_rhs != 0)
-    {
-        sig_rhs /= 10;
-        ++exp_rhs;
-    }
+    detail::normalize<decimal128>(sig_rhs, exp_rhs);
 
     auto mul_result {d128_mul_impl(sig_lhs, exp_lhs, x.isneg(), sig_rhs, exp_rhs, y.isneg())};
     const decimal128 dec_result {mul_result.sig, mul_result.exp, mul_result.sign};
 
+    return dec_result + z;
+
+    /*
     const auto res_add {detail::check_non_finite(dec_result, z)};
     if (res_add != zero)
     {
@@ -2493,6 +2476,7 @@ constexpr auto fmad128(decimal128 x, decimal128 y, decimal128 z) noexcept -> dec
     }
 
     return {result.sig, result.exp, result.sign};
+    */
 }
 
 } //namespace decimal

include/boost/decimal/detail/cmath/impl/taylor_series_result.hpp

@@ -27,27 +27,6 @@ constexpr auto taylor_series_result(T x, const Array &coeffs) noexcept
     return result;
 }
 
-template <typename Array>
-constexpr auto taylor_series_result(boost::decimal::decimal128 x, const Array &coeffs) noexcept
-{
-    const std::size_t N = coeffs.size();
-
-    auto result = coeffs[N - 1];
-
-    auto my_own_fma =
-        [](boost::decimal::decimal128 x, boost::decimal::decimal128 y, boost::decimal::decimal128 z)
-        {
-            return (x * y) + z;
-        };
-
-    for (std::size_t i = N - 1; i-- > 0;)
-    {
-        result = my_own_fma(result, x, coeffs[i]);
-    }
-
-    return result;
-}
-
 } //namespace detail
 } //namespace decimal
 } //namespace boost

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

@@ -99,14 +99,8 @@ constexpr auto tgamma_impl(T x) noexcept
                     z = z - nx;
                 }
 
-                auto my_own_fma =
-                    [](T x, T y, T z)
-                    {
-                        return (x * y) + z;
-                    };
-
                 result = detail::tgamma_series_expansion(z);
-                result = one / (z * my_own_fma(result, z, one));
+                result = one / (z * fma(result, z, one));
 
                 if (x_is_gt_one)
                 {

include/boost/decimal/detail/remove_trailing_zeros.hpp

@@ -85,9 +85,15 @@ constexpr auto remove_trailing_zeros(std::uint64_t n) noexcept -> remove_trailin
     return {n, s};
 }
 
-// TODO(mborland): Make this better. Check lower word for equal to 0.
+// TODO(mborland): Make this better for the 2-word case
 constexpr auto remove_trailing_zeros(uint128 n) noexcept -> remove_trailing_zeros_return<uint128>
 {
+    if (n.high == UINT64_C(0))
+    {
+        const auto temp {remove_trailing_zeros(n.low)};
+        return {static_cast<uint128>(temp.trimmed_number), temp.number_of_removed_zeros};
+    }
+
     std::size_t s {};
 
     while (n % 10 == 0)

test/Jamfile

@@ -45,6 +45,7 @@ run-fail benchmarks.cpp ;
 compile-fail concepts_test.cpp ;
 run github_issue_426.cpp ;
 run github_issue_448.cpp ;
+run-fail github_issue_519.cpp ;
 run ibm_abs.cpp ;
 run ibm_add.cpp ;
 run link_1.cpp link_2.cpp link_3.cpp ;

test/github_issue_519.cpp

@@ -0,0 +1,47 @@
+// Copyright 2024 Matt Borland
+// Copyright 2024 Christopher Kormanyos
+// Distributed under the Boost Software License, Version 1.0.
+// https://www.boost.org/LICENSE_1_0.txt
+
+#include <boost/decimal.hpp>
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <sstream>
+
+auto main() -> int
+{
+    using local_decimal_type = boost::decimal::decimal128;
+
+    // N[Gamma[456/100], 34]
+    // 12.64819265438397922113369900828315
+    //
+    // For Decimal 128:
+    // 4.56
+    // 12.64819265438397922113369900828314
+
+    const local_decimal_type x = local_decimal_type { 456, -2 };
+
+    const auto tg = tgamma(x);
+
+    {
+        std::stringstream strm;
+
+        strm << std::setprecision(std::numeric_limits<local_decimal_type>::digits10) << x;
+
+        std::cout << "Expected: 4.56" << std::endl;
+        std::cout << "     Got: " << strm.str() << std::endl;
+    }
+
+    {
+        std::stringstream strm;
+
+        strm << std::setprecision(std::numeric_limits<local_decimal_type>::digits10) << tg;
+
+        std::cout << "Expected: 12.64819265438397922113369900828314" << std::endl;
+        std::cout << "     Got: " << strm.str() << std::endl;
+    }
+
+    return 1;
+}

test/test_tgamma.cpp

@@ -63,6 +63,14 @@ namespace local
       const auto delta = fabs(1 - (a / b));
 
       result_is_ok = (delta < tol);
+
+      if (!result_is_ok)
+      {
+          std::cerr << std::setprecision(std::numeric_limits<NumericType>::digits10) << "a: " << a
+                    << "\nb: " << b
+                    << "\ndelta: " << delta
+                    << "\ntol: " << tol << std::endl;
+      }
     }
 
     return result_is_ok;