From 954af38c4a354c8bb5e093acaefb157185b0b2de Mon Sep 17 00:00:00 2001 From: "Tarn W. Burton" Date: Thu, 22 Aug 2024 14:08:55 -0400 Subject: [PATCH] Add float_convert::to_quadruple/from_quadruple --- include/clasp/core/bignum.h | 27 +++++- include/clasp/core/float_util.h | 165 ++++++++++++++++++++++++++++++++ include/clasp/core/foundation.h | 36 ------- include/clasp/core/object.h | 1 + src/core/bignum.cc | 50 +++++----- src/core/foundation.cc | 4 - src/core/num_co.cc | 91 ++++++------------ src/core/numbers.cc | 52 ++++------ 8 files changed, 265 insertions(+), 161 deletions(-) create mode 100644 include/clasp/core/float_util.h diff --git a/include/clasp/core/bignum.h b/include/clasp/core/bignum.h index 1a563c5055..5017d2efa2 100644 --- a/include/clasp/core/bignum.h +++ b/include/clasp/core/bignum.h @@ -86,6 +86,12 @@ class Bignum_O : public Integer_O { static Bignum_sp create(int64_t v) { return create_from_limbs((v < 0) ? -1 : 1, std::abs(v), true); } #endif static Bignum_sp create(uint64_t v) { return create_from_limbs(1, v, true); } + static Bignum_sp create(__uint128_t v) { + Bignum_sp b = create_from_limbs(2); + b->_limbs[0] = static_cast(v); + b->_limbs[1] = static_cast(v >> 64); + return b; + } static Bignum_sp create(double d) { // KLUDGE: there is no mpn function for converting from a double. // However, this conses, which we shouldn't need to do. @@ -133,7 +139,26 @@ class Bignum_O : public Integer_O { virtual bool evenp_() const override { return !((this->limbs())[0] & 1); } virtual bool oddp_() const override { return (this->limbs())[0] & 1; } - template integral to_integral() const { + template integral to_integral() const { + constexpr auto limb_width = sizeof(mp_limb_t) * 8; + constexpr auto integral_width = std::bit_width(std::numeric_limits::max()); + constexpr auto max_limb_count = 1 + ((sizeof(integral) - 1) / sizeof(mp_limb_t)); + + const mp_limb_t* limbs = this->limbs(); + mp_size_t len = this->length(); + + if ((len > -1) && ((limb_width * (len - 1) + std::bit_width(limbs[len - 1])) <= integral_width)) { + integral value = 0; + for (mp_size_t i = 0; i < length(); i++) + value |= (integral)limbs[i] << (i * limb_width); + return value; + } + + TYPE_ERROR(this->asSmartPtr(), Cons_O::createList(cl::_sym_Integer_O, Integer_O::create(std::numeric_limits::min()), + Integer_O::create(std::numeric_limits::max()))); + } + + template integral to_integral() const { integral mn = std::numeric_limits::min(); integral mx = std::numeric_limits::max(); // First, if integral can only hold fixnums, conversion will always fail. diff --git a/include/clasp/core/float_util.h b/include/clasp/core/float_util.h new file mode 100644 index 0000000000..da6e4c794b --- /dev/null +++ b/include/clasp/core/float_util.h @@ -0,0 +1,165 @@ +#pragma once + +#include + +namespace core { + +template struct float_convert { + static constexpr uint16_t significand_width = std::numeric_limits::digits; + static constexpr uint16_t exponent_width = std::bit_width((unsigned int)std::numeric_limits::max_exponent); + static constexpr uint16_t sign_width = 1; + static constexpr bool has_hidden_bit = ((sign_width + exponent_width + significand_width) % 8) != 0; + static constexpr uint16_t storage_width = sign_width + exponent_width + significand_width + ((has_hidden_bit) ? -1 : 0); + static constexpr int32_t exponent_bias = std::numeric_limits::max_exponent + significand_width - 2; + static constexpr int32_t max_exponent = std::numeric_limits::max_exponent - significand_width; + static constexpr int32_t min_exponent = 2 - exponent_bias - significand_width; + static constexpr int32_t min_normalized_exponent = 1 - exponent_bias; + using uint_t = + std::conditional_t>>>; + static constexpr uint_t hidden_bit = has_hidden_bit ? uint_t{1} << (significand_width - 1) : uint_t{0}; + static constexpr uint16_t exponent_shift = storage_width - sign_width - exponent_width; + static constexpr uint16_t sign_shift = storage_width - sign_width; + static constexpr uint_t significand_mask = (uint_t{1} << (significand_width + ((has_hidden_bit) ? -1 : 0))) - uint_t{1}; + static constexpr uint_t payload_mask = (uint_t{1} << (significand_width + ((has_hidden_bit) ? -2 : -1))) - uint_t{1}; + static constexpr int32_t non_finite_exponent = static_cast(((uint_t{1} << exponent_width) - uint_t{1})); + static constexpr uint_t exponent_mask = ((uint_t{1} << exponent_width) - uint_t{1}) << exponent_shift; + static constexpr uint_t sign_mask = ((uint_t{1} << sign_width) - uint_t{1}) << sign_shift; + static constexpr uint_t nan_type_mask = uint_t{1} << (significand_width + ((has_hidden_bit) ? -2 : -1)); + + enum class category { finite, quiet_nan, signaling_nan, infinity }; + + typedef union { + Float f; + uint_t b; + } convert_t; + + static inline uint_t to_bits(Float f) { + convert_t convert = {.f = f}; + return convert.b; + } + + static inline Float from_bits(uint_t b) { + convert_t convert = {.b = b}; + return convert.f; + } + + struct quadruple { + category category; + uint_t significand; + int32_t exponent; + int16_t sign; + }; + + static quadruple to_quadruple(Float f) { + quadruple q; + uint_t b = to_bits(f); + + q.significand = b & significand_mask; + q.exponent = static_cast((b & exponent_mask) >> exponent_shift); + q.sign = (b & sign_mask) ? int32_t{-1} : int32_t{1}; + + if (q.exponent != non_finite_exponent) { + q.category = category::finite; + if (q.exponent != 0) { + // Normal non-zero + q.significand |= hidden_bit; + q.exponent -= exponent_bias; + } else if (q.significand != 0) { + // Subnormal + int32_t shift = significand_width - std::bit_width(q.significand); + q.significand = q.significand << shift; + q.exponent = 1 - exponent_bias - shift; + } + } else if (q.significand == 0) { + q.category = category::infinity; + } else if (q.significand & nan_type_mask) { + q.category = category::quiet_nan; + q.significand &= ~nan_type_mask; + } else { + q.category = category::signaling_nan; + } + + return q; + } + + static Float from_quadruple(quadruple q) { + uint_t b = 0; + + if (q.sign < 0) + b |= sign_mask; + + switch (q.category) { + case category::infinity: + b |= exponent_mask; + break; + case category::quiet_nan: + b |= exponent_mask | nan_type_mask | (q.significand & payload_mask); + break; + case category::signaling_nan: + b |= exponent_mask | ((q.significand == 0) ? uint_t{1} : (q.significand & payload_mask)); + break; + default: + if (q.significand != 0) { + uint_t significand = q.significand; + int32_t exponent = q.exponent; + int32_t shift = std::bit_width(significand) - significand_width; + + // If we don't have enough bits then right shift. + if (shift < 0) { + significand = significand << -shift; + exponent += shift; + shift = 0; + } + + // Check for subnormals and set the shift needed to normalize. + if ((exponent + shift) < min_normalized_exponent) { + shift = min_normalized_exponent - exponent; + } + + // If we shift away all of the bits that is an underflow. + if (shift > std::bit_width(significand)) { + feraiseexcept(FE_UNDERFLOW); + // Return +/- zero if traps masked. + return from_bits(b); + } + + // Round if we have extra bits. + if (shift > 0) { + if (significand & (1 << shift) - 1)) + feraiseexcept(FE_INEXACT); + significand = ((significand & (1 << (shift - 1))) ? (significand + (1 << shift)) : significand) >> shift; + exponent += shift; + } + + // Check one more time to ensure rounding hasn't increased the width. + shift = std::max(std::bit_width(significand) - significand_width, 0); + significand = significand >> shift; + exponent += shift; + + // Check for overflow. + if (exponent > max_exponent) { + feraiseexcept(FE_OVERFLOW); + // Return +/- infinity if traps masked. + return from_bits(b | exponent_mask); + } + + if (std::bit_width(significand) < significand_width) { + // Subnormals + b |= significand & significand_mask; + } else { + // Normal + b |= (significand & significand_mask) | + ((static_cast(exponent + exponent_bias) << exponent_shift) & exponent_mask); + } + } + break; + } + + return from_bits(b); + }; +}; + +} // namespace core diff --git a/include/clasp/core/foundation.h b/include/clasp/core/foundation.h index dd48065d62..98eebb03bd 100644 --- a/include/clasp/core/foundation.h +++ b/include/clasp/core/foundation.h @@ -51,46 +51,10 @@ THE SOFTWARE. namespace core { [[noreturn]] void lisp_throwLispError(const std::string& str); -[[noreturn]] void lisp_floating_point_invalid_operation(); [[noreturn]] void lisp_error_simple(const char* functionName, const char* fileName, int lineNumber, const string& fmt); [[noreturn]] void lisp_error_simple(const char* functionName, const char* fileName, int lineNumber, const std::string& str); void lisp_debugLogWrite(const char* fileName, const char* funcName, uint lineNumber, uint column, const std::string& message, uint debugFlags = DEBUG_CPP_FUNCTION); - -template struct float_convert { - static constexpr uint16_t significand_width = std::numeric_limits::digits; - static constexpr uint16_t exponent_width = std::bit_width((unsigned int)std::numeric_limits::max_exponent); - static constexpr uint16_t sign_width = 1; - static constexpr bool has_hidden_bit = ((sign_width + exponent_width + significand_width) % 8) != 0; - static constexpr uint16_t storage_width = sign_width + exponent_width + significand_width + ((has_hidden_bit) ? -1 : 0); - static constexpr int32_t exponent_bias = std::numeric_limits::max_exponent + significand_width - 2; - using uint_t = - std::conditional_t>>>; - static constexpr uint16_t exponent_shift = storage_width - sign_width - exponent_width; - static constexpr uint16_t sign_shift = storage_width - sign_width; - static constexpr uint_t significand_mask = (uint_t{1} << (significand_width + ((has_hidden_bit) ? -1 : 0))) - uint_t{1}; - static constexpr uint_t exponent_mask = ((uint_t{1} << exponent_width) - uint_t{1}) << exponent_shift; - static constexpr uint_t sign_mask = ((uint_t{1} << sign_width) - uint_t{1}) << sign_shift; - - typedef union { - Float f; - uint_t b; - } convert_t; - - static inline uint_t to_bits(Float f) { - convert_t convert = {.f = f}; - return convert.b; - } - - static inline Float from_bits(uint_t b) { - convert_t convert = {.b = b}; - return convert.f; - } -}; - }; // namespace core template struct fmt::formatter : fmt::formatter> { diff --git a/include/clasp/core/object.h b/include/clasp/core/object.h index dcb4abd9b6..cfb32890db 100644 --- a/include/clasp/core/object.h +++ b/include/clasp/core/object.h @@ -36,6 +36,7 @@ THE SOFTWARE. #include #include +#include #include #include #include diff --git a/src/core/bignum.cc b/src/core/bignum.cc index 393e40ab13..b0cc494d64 100644 --- a/src/core/bignum.cc +++ b/src/core/bignum.cc @@ -799,35 +799,39 @@ Number_sp Bignum_O::oneMinus_() const { return next_fadd(this->limbs(), this->le Number_sp Bignum_O::onePlus_() const { return next_fadd(this->limbs(), this->length(), 1); } -double next_to_double(mp_size_t len, const mp_limb_t* limbs) { - // MPN does not seem to export a float conversion function, - // so we roll our own. FIXME: Could get bad with float rounding - // or otherwise as I'm making this up as I go. +template Float limbs_to_float(mp_size_t len, const mp_limb_t* limbs) { + constexpr size_t limb_width = sizeof(mp_limb_t) * 8; mp_size_t size = std::abs(len); + struct float_convert::quadruple q = { + .category = float_convert::category::finite, .significand = 0, .exponent = (size - 1) * limb_width, + .sign = (len < 0) ? -1 : 1 + }; + size_t shift = float_convert::significand_width + 1; + size_t width = std::bit_width(limbs[size - 1]); + + if (width >= shift) { + q.significand = limbs[size - 1] >> (width - shift); + q.exponent += width - shift; + } else { + q.significand = limbs[size - 1]; + shift -= width; + + for (mp_size_t i = size - 2; i > -1 && shift > 0; i--) { + size_t limb_shift = std::min(shift, limb_width); + q.significand = (q.significand << limb_shift) | (limbs[i] >> (limb_width - limb_shift)); + q.exponent -= limb_shift; + shift -= limb_shift; + } + } - // There are no length zero bignums. - // If the bignum is only one long, we just convert directly. - if (len == 1) - return static_cast(limbs[0]); - else if (len == -1) - return -static_cast(limbs[0]); - - // Otherwise, we just use the most significant two limbs. - // Assuming the float format has at most 128 bits of significand, - // the lower limbs ought to be irrelevant. - // (In fact one 64-bit number would probably be enough, but - // it's possible the high bits of this number are almost all zero.) - double ultimate = static_cast(limbs[size - 1]); - double penultimate = static_cast(limbs[size - 2]); - double soon = penultimate + std::ldexp(ultimate, 64); - return std::ldexp(((len < 0) ? -soon : soon), 64 * (size - 2)); + return float_convert::from_quadruple(q); } -float Bignum_O::as_float_() const { return static_cast(mpz_get_d(this->mpz().get_mpz_t())); } +float Bignum_O::as_float_() const { return limbs_to_float(this->length(), this->limbs()); } -double Bignum_O::as_double_() const { return mpz_get_d(this->mpz().get_mpz_t()); } +double Bignum_O::as_double_() const { return limbs_to_float(this->length(), this->limbs()); } -LongFloat Bignum_O::as_long_float_() const { return static_cast(next_to_double(this->length(), this->limbs())); } +LongFloat Bignum_O::as_long_float_() const { return limbs_to_float(this->length(), this->limbs()); } DOCGROUP(clasp); CL_DEFUN int core__next_compare(Bignum_sp left, Bignum_sp right) { diff --git a/src/core/foundation.cc b/src/core/foundation.cc index 6c6aa7b90b..acfe1b6d81 100644 --- a/src/core/foundation.cc +++ b/src/core/foundation.cc @@ -214,10 +214,6 @@ CL_DEFUN void core__dump_class_ids() { reg::dump_class_ids(); } SYMBOL_EXPORT_SC_(ClPkg,floating_point_invalid_operation); -void lisp_floating_point_invalid_operation() { - lisp_error(cl::_sym_floating_point_invalid_operation,nil()); -} - void lisp_errorExpectedList(core::T_O* v) { T_sp tv((gctools::Tagged)v); TYPE_ERROR(tv, cl::_sym_list); diff --git a/src/core/num_co.cc b/src/core/num_co.cc index 89536215a7..68ba57afbf 100644 --- a/src/core/num_co.cc +++ b/src/core/num_co.cc @@ -942,82 +942,45 @@ CL_DEFUN Integer_sp cl__float_precision(Float_sp x) { return clasp_make_fixnum(precision); } +template inline Real_mv integer_decode_float(Float f) { + struct float_convert::quadruple q; + + // IEEE-754-2019 7.2.j specifies that NaN or infinity should result in FE_INVALID for integer producing operations. + switch (std::fpclassify(f)) { + case FP_INFINITE: + feraiseexcept(FE_INVALID); + q = float_convert::to_quadruple(std::signbit(f) ? std::numeric_limits::min() : std::numeric_limits::max()); + break; + case FP_NAN: + feraiseexcept(FE_INVALID); + q = float_convert::to_quadruple(std::signbit(f) ? Float{-0.0} : Float{0.0}); + break; + default: + q = float_convert::to_quadruple(f); + break; + } + + return Values(Integer_O::create(q.significand), clasp_make_fixnum(q.exponent), clasp_make_fixnum(q.sign)); +} + CL_LAMBDA(x); CL_DECLARE(); CL_UNWIND_COOP(true); CL_DOCSTRING(R"dx(integer_decode_float)dx"); DOCGROUP(clasp); CL_DEFUN Real_mv cl__integer_decode_float(Float_sp x) { - int e = 0, s = 1; - Real_sp rx(nil()); switch (clasp_t_of(x)) { #ifdef CLASP_LONG_FLOAT - case number_LongFloat: { - LongFloat d = clasp_long_float(x); - if (std::isfinite(d)) { - if (std::signbit(d)) { - s = -1; - d = -d; - } - if (d == 0.0) { - e = 0; - rx = clasp_make_fixnum(0); - } else { - d = frexpl(d, &e); - rx = _clasp_long_double_to_integer(ldexpl(d, LDBL_MANT_DIG)); - e -= LDBL_MANT_DIG; - } - } else { - SIMPLE_ERROR("Can't decode NaN or infinity {}", _rep_(x)); - } - break; - } + case number_LongFloat: + return integer_decode_float(gc::As_unsafe(x)->get()); #endif - case number_DoubleFloat: { - double d = gc::As_unsafe(x)->get(); - if (std::isfinite(d)) { - if (std::signbit(d)) { - s = -1; - d = -d; - } - if (d == 0.0) { - e = 0; - rx = clasp_make_fixnum(0); - } else { - d = frexp(d, &e); - rx = _clasp_double_to_integer(ldexp(d, DBL_MANT_DIG)); - e -= DBL_MANT_DIG; - } - } else { - SIMPLE_ERROR("Can't decode NaN or infinity {}", _rep_(x)); - } - break; - } - case number_SingleFloat: { - float d = x.unsafe_single_float(); - if (std::isfinite(d)) { - if (std::signbit(d)) { - s = -1; - d = -d; - } - if (d == 0.0) { - e = 0; - rx = clasp_make_fixnum(0); - } else { - d = frexpf(d, &e); - rx = _clasp_double_to_integer(ldexp(d, FLT_MANT_DIG)); - e -= FLT_MANT_DIG; - } - } else { - SIMPLE_ERROR("Can't decode NaN or infinity {}", _rep_(x)); - } - break; - } + case number_DoubleFloat: + return integer_decode_float(gc::As_unsafe(x)->get()); + case number_SingleFloat: + return integer_decode_float(x.unsafe_single_float()); default: ERROR_WRONG_TYPE_NTH_ARG(cl::_sym_integer_decode_float, 1, x, cl::_sym_float); } - ASSERT(rx.notnilp()); - return Values(rx, clasp_make_fixnum(e), clasp_make_fixnum(s)); } CL_LAMBDA(r &optional (i 0)); diff --git a/src/core/numbers.cc b/src/core/numbers.cc index fe1c2470f7..c9d6b2c195 100644 --- a/src/core/numbers.cc +++ b/src/core/numbers.cc @@ -1713,42 +1713,28 @@ static Integer_sp mantissa_and_exponent_from_ratio(Integer_sp num, Integer_sp de return quotient; } -float Ratio_O::as_float_() const { - gc::Fixnum exponent; - Integer_sp mantissa = - mantissa_and_exponent_from_ratio(this->_numerator, this->_denominator, std::numeric_limits::digits, &exponent); - return std::ldexp(mantissa.fixnump() ? static_cast(mantissa.unsafe_fixnum()) : mantissa->as_float_(), exponent); -} - -double Ratio_O::as_double_() const { - gc::Fixnum exponent; - Integer_sp mantissa = - mantissa_and_exponent_from_ratio(this->_numerator, this->_denominator, std::numeric_limits::digits, &exponent); - return std::ldexp(mantissa.fixnump() ? static_cast(mantissa.unsafe_fixnum()) : mantissa->as_double_(), exponent); - - /*if ((this->_numerator).fixnump() && (this->_denominator).fixnump()) { - double d = clasp_to_double(this->_numerator); - d /= clasp_to_double(this->_denominator); - return d; - } else { - gc::Fixnum exponent; - Integer_sp mantissa = mantissa_and_exponent_from_ratio(this->_numerator, this->_denominator, DBL_MANT_DIG, &exponent); - double output; - if (mantissa.fixnump()) - output = mantissa.unsafe_fixnum(); - else - output = mantissa->as_double_(); - return std::ldexp(output, exponent); - }*/ -} +template inline Float ratio_to_float(Integer_sp num, Integer_sp den) { + struct float_convert::quadruple q = { + .sign = 1 + }; -LongFloat Ratio_O::as_long_float_() const { - gc::Fixnum exponent; - Integer_sp mantissa = - mantissa_and_exponent_from_ratio(this->_numerator, this->_denominator, std::numeric_limits::digits, &exponent); - return std::ldexp(mantissa.fixnump() ? static_cast(mantissa.unsafe_fixnum()) : mantissa->as_long_float_(), exponent); + if (clasp_minusp(num)) { + q.sign = -1; + num = gc::As_unsafe(clasp_negate(num)); + } + + q.exponent = clasp_integer_length(num) - clasp_integer_length(den) - float_convert::significand_width - 1; + q.significand = clasp_to_integral::uint_t>(clasp_integer_divide(clasp_ash(num, -q.exponent), den)); + + return float_convert::from_quadruple(q); } +float Ratio_O::as_float_() const { return ratio_to_float(this->_numerator, this->_denominator); } + +double Ratio_O::as_double_() const { return ratio_to_float(this->_numerator, this->_denominator); } + +LongFloat Ratio_O::as_long_float_() const { return ratio_to_float(this->_numerator, this->_denominator); } + string Ratio_O::__repr__() const { stringstream ss; ss << _rep_(this->_numerator) << "/" << _rep_(this->_denominator);