Skip to content

Commit

Permalink
Streamline double-fp details and multiply
Browse files Browse the repository at this point in the history
  • Loading branch information
@ckormanyos
ckormanyos committed Dec 30, 2024
1 parent 37e28bf commit 7faf35e
Show file tree
Hide file tree
Showing 2 changed files with 120 additions and 42 deletions.
59 changes: 52 additions & 7 deletions include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,18 +14,63 @@
#include <boost/multiprecision/number.hpp>
#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath.hpp>

#include <cmath>
#include <limits>
#include <tuple>
#include <utility>

namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail {

inline constexpr float split (float) { return static_cast<float> (1U + static_cast<unsigned long long>(static_cast<unsigned long long>(1U) << static_cast<unsigned>((cpp_df_qf_detail::ccmath::numeric_limits<float >::digits + 1) / 2))); }
inline constexpr double split (double) { return static_cast<double> (1U + static_cast<unsigned long long>(static_cast<unsigned long long>(1U) << static_cast<unsigned>((cpp_df_qf_detail::ccmath::numeric_limits<double >::digits + 1) / 2))); }
inline constexpr long double split (long double) { return static_cast<long double>(1U + static_cast<unsigned long long>(static_cast<unsigned long long>(1U) << static_cast<unsigned>((cpp_df_qf_detail::ccmath::numeric_limits<long double >::digits + 1) / 2))); }
// Define a templated function with an EnableType
template <typename FloatType, typename EnableType = void>
struct split_maker { };

template <typename FloatType>
struct split_maker<FloatType,
typename ::std::enable_if<::std::is_floating_point<FloatType>::value>::type>
{
private:
using float_type = FloatType;

public:
static constexpr int
n_shl
{
static_cast<int>((ccmath::numeric_limits<float_type>::digits + 1) / 2)
};

static constexpr float_type
value
{
static_cast<float_type>
(
1ULL + static_cast<unsigned long long>(1ULL << static_cast<unsigned>(n_shl))
)
};
};

#if defined(BOOST_HAS_FLOAT128)
inline constexpr ::boost::float128_type split(::boost::float128_type) { return static_cast<::boost::float128_type>(1) + static_cast<::boost::float128_type>(static_cast<::boost::uint128_type>(1U) << static_cast<unsigned>((cpp_df_qf_detail::ccmath::numeric_limits<::boost::float128_type>::digits + 1) / 2)); }
template <typename FloatType>
struct split_maker<FloatType,
typename ::std::enable_if<::std::is_same<FloatType, ::boost::float128_type>::value>::type>
{
private:
using float_type = FloatType;

public:
static constexpr int
n_shl
{
static_cast<int>((ccmath::numeric_limits<float_type>::digits + 1) / 2)
};

static constexpr float_type
value
{
static_cast<float_type>
(
::boost::uint128_type { 1 }
+ ::boost::uint128_type { ::boost::uint128_type { 1 } << static_cast<unsigned>(n_shl) }
)
};
};
#endif

template <class FloatingPointType>
Expand Down
103 changes: 68 additions & 35 deletions include/boost/multiprecision/cpp_double_fp.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -452,6 +452,7 @@ class cpp_double_fp_backend
{
*this = cpp_double_fp_backend::my_value_nan();
}

return *this;
}

Expand Down Expand Up @@ -553,13 +554,9 @@ class cpp_double_fp_backend
constexpr cpp_double_fp_backend& operator*=(const cpp_double_fp_backend& v)
{
// Evaluate the sign of the result.
const auto isneg_u = isneg();
const auto isneg_v = v.isneg();

const bool b_result_is_neg = (isneg_u != isneg_v);

const auto fpc_u = eval_fpclassify(*this);
const auto fpc_v = eval_fpclassify(v);
const int fpc_u { eval_fpclassify(*this) };
const int fpc_v { eval_fpclassify(v) };

// Handle special cases like zero, inf and NaN.
const bool isinf_u { (fpc_u == FP_INFINITE) };
Expand All @@ -576,10 +573,15 @@ class cpp_double_fp_backend

if (isinf_u || isinf_v)
{
const bool b_neg { (isneg() != v.isneg()) };

*this = cpp_double_fp_backend::my_value_inf();

if (b_result_is_neg)
if (b_neg)
{
negate();
}

return *this;
}

Expand All @@ -598,40 +600,59 @@ class cpp_double_fp_backend
// K. Briggs work. The algorithm has been significantly simplified
// while still atempting to retain proper rounding corrections.

float_type C { cpp_df_qf_detail::split(float_type()) * data.first };
float_type C { cpp_df_qf_detail::split_maker<float_type>::value * data.first };

float_type hu { C - float_type { C - data.first } };
float_type hu;

if (cpp_df_qf_detail::ccmath::isinf(C))
{
// In this case, data.first is evidently so large that multipllication
// with the split has overflowed to infinity. Take a last-chance try
// for multiplication by scaling with the split and subsequently
// taking the difference term.
// In this case, data.first is evidently so large that multipllication
// with the split has overflowed to infinity. Take a last-chance try
// for multiplication by scaling down (and then back up) with the split.

// TBD: Is this numericaly correct? Or is it simply the "best" we can
// do in this edge-case? Or is there an alternate multiplication scheme
// that is not susceptible to this particular sort of overflow?
C = data.first - cpp_df_qf_detail::ccmath::ldexp(data.first, -cpp_df_qf_detail::split_maker<float_type>::n_shl);

C = data.first - (data.first / cpp_df_qf_detail::split(float_type()));

hu = (data.first - C) * cpp_df_qf_detail::split(float_type());
hu = cpp_df_qf_detail::ccmath::ldexp(data.first - C, cpp_df_qf_detail::split_maker<float_type>::n_shl);
}
else
{
hu = C - float_type { C - data.first };
}

C = data.first * v.data.first;

if (cpp_df_qf_detail::ccmath::isinf(C))
{
const bool b_neg { (isneg() != v.isneg()) };

*this = cpp_double_fp_backend::my_value_inf();

if (b_result_is_neg)
if (b_neg)
{
negate();
}

return *this;
}

float_type c { cpp_df_qf_detail::split(float_type()) * v.data.first };
float_type c { cpp_df_qf_detail::split_maker<float_type>::value * v.data.first };

const float_type hv { c - float_type { c - v.data.first } };
float_type hv;

if (cpp_df_qf_detail::ccmath::isinf(c))
{
// In this case, v.data.first is evidently so large that multipllication
// with the split has overflowed to infinity. Take a last-chance try
// for multiplication by scaling down (and then back up) with the split.

c = v.data.first - cpp_df_qf_detail::ccmath::ldexp(v.data.first, -cpp_df_qf_detail::split_maker<float_type>::n_shl);

hv = cpp_df_qf_detail::ccmath::ldexp(v.data.first - c, cpp_df_qf_detail::split_maker<float_type>::n_shl);
}
else
{
hv = c - float_type { c - v.data.first };
}

const float_type tv { v.data.first - hv };

Expand Down Expand Up @@ -733,16 +754,22 @@ class cpp_double_fp_backend

const float_type C { data.first / v.data.first };

float_type c { cpp_df_qf_detail::split(float_type()) * C };
float_type c { cpp_df_qf_detail::split_maker<float_type>::value * C };

float_type u { cpp_df_qf_detail::split(float_type()) * v.data.first };
const float_type hc { c - float_type { c - C } };

float_type u { cpp_df_qf_detail::split_maker<float_type>::value * v.data.first };

const float_type hv { u - float_type { u - v.data.first } };

if (cpp_df_qf_detail::ccmath::isinf(u) || cpp_df_qf_detail::ccmath::isinf(c))
{
// Evidently we have some very large operands. Take a last-chance try
// for finite division. Use the ratio of square roots and subsequently
// square the ratio, handling the sign of the result separately.

// TBD: Is there a more sensible catch (and handling) for this case?

const bool u_neg { isneg() };
const bool v_neg { v.isneg() };
const bool b_neg { u_neg != v_neg };
Expand All @@ -766,14 +793,13 @@ class cpp_double_fp_backend
eval_multiply(*this, sqrt_ratio);

if (b_neg)
{
negate();
}

return *this;
}

const float_type hc { c - float_type { c - C } };

const float_type hv { u - float_type { u - v.data.first } };

{
const float_type U { C * v.data.first };

Expand Down Expand Up @@ -806,14 +832,18 @@ class cpp_double_fp_backend
constexpr cpp_double_fp_backend operator++(int)
{
cpp_double_fp_backend t(*this);
++*this;

++(*this);

return t;
}

constexpr cpp_double_fp_backend operator--(int)
{
cpp_double_fp_backend t(*this);
--*this;

--(*this);

return t;
}

Expand Down Expand Up @@ -910,7 +940,9 @@ class cpp_double_fp_backend
std::string str(std::streamsize number_of_digits, const std::ios::fmtflags format_flags) const
{
if (number_of_digits == 0)
{
number_of_digits = cpp_double_fp_backend::my_digits10;
}

// Use cpp_dec_float to write string (as is similarly done to read string).

Expand Down Expand Up @@ -1031,9 +1063,9 @@ class cpp_double_fp_backend
// TBD: Or should we use the static-initializer trick (as used in cpp_dec_float)
// to initialize this and possible other similar values?

cpp_double_fp_backend result { };
cpp_double_fp_backend result;

eval_ldexp(result, cpp_double_fp_backend(1), 3 - my_digits);
eval_ldexp(result, cpp_double_fp_backend { 1 }, 3 - my_digits);

return result;
}
Expand Down Expand Up @@ -1097,9 +1129,10 @@ bool cpp_double_fp_backend<FloatingPointType>::rd_string(const char* pstr)
}
else
{
cpp_dec_float_read_write_type dummy_frexp { };
cpp_dec_float_read_write_type dummy_frexp;

int e2_from_f_dec;

int e2_from_f_dec { };
eval_frexp(dummy_frexp, f_dec, &e2_from_f_dec);

const auto is_definitely_zero =
Expand Down Expand Up @@ -1480,7 +1513,7 @@ constexpr void eval_sqrt(cpp_double_fp_backend<FloatingPointType>& result, const

const local_float_type c = cpp_df_qf_detail::ccmath::sqrt(o.crep().first);

local_float_type p = cpp_df_qf_detail::split(local_float_type()) * c;
local_float_type p = cpp_df_qf_detail::split_maker<local_float_type>::value * c;
local_float_type hx = (c - p);
hx = hx + p;
local_float_type tx = c - hx;
Expand Down

0 comments on commit 7faf35e

Please sign in to comment.