From a38f983e434cc49dff80e3877c3f2b707e7c27a2 Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Sun, 2 Jun 2024 17:27:25 -0400 Subject: [PATCH 01/18] Decimal <--> Floating conversion. --- cpp/benchmarks/decimal/convert_floating.cpp | 17 - .../cudf/fixed_point/floating_conversion.hpp | 746 +++++++++++++++++- cpp/include/cudf/unary.hpp | 35 +- cpp/tests/fixed_point/fixed_point_tests.cpp | 128 ++- 4 files changed, 847 insertions(+), 79 deletions(-) diff --git a/cpp/benchmarks/decimal/convert_floating.cpp b/cpp/benchmarks/decimal/convert_floating.cpp index a367036c4940..ac09c3400cbd 100644 --- a/cpp/benchmarks/decimal/convert_floating.cpp +++ b/cpp/benchmarks/decimal/convert_floating.cpp @@ -32,8 +32,6 @@ void bench_cast_decimal(nvbench::state& state, nvbench::type_list || std::is_same_v; - static constexpr bool is_32bit = - std::is_same_v || std::is_same_v; static constexpr bool is_128bit = std::is_same_v || std::is_same_v; @@ -69,21 +67,6 @@ void bench_cast_decimal(nvbench::state& state, nvbench::type_list decimal conversion algorithm is limited - static constexpr bool is_64bit = !is_32bit && !is_128bit; - if (is_32bit && (exp_mode != 3)) { - state.skip("Decimal32 conversion only works up to scale factors of 10^9."); - return; - } - if (is_64bit && ((exp_mode < 2) || (exp_mode > 4))) { - state.skip("Decimal64 conversion only works up to scale factors of 10^18."); - return; - } - if (is_128bit && ((exp_mode == 0) || (exp_mode == 6))) { - state.skip("Decimal128 conversion only works up to scale factors of 10^38."); - return; - } - // Type IDs auto const input_id = cudf::type_to_id(); auto const output_id = cudf::type_to_id(); diff --git a/cpp/include/cudf/fixed_point/floating_conversion.hpp b/cpp/include/cudf/fixed_point/floating_conversion.hpp index 2c3a5c5629dd..62499fc078f4 100644 --- a/cpp/include/cudf/fixed_point/floating_conversion.hpp +++ b/cpp/include/cudf/fixed_point/floating_conversion.hpp @@ -81,8 +81,7 @@ struct floating_converter { // To store positive and negative exponents as unsigned values, the stored value for // the power-of-2 is exponent + bias. The bias is 127 for floats and 1023 for doubles. /// 127 / 1023 for float / double - static constexpr IntegralType exponent_bias = - cuda::std::numeric_limits::max_exponent - 1; + static constexpr int exponent_bias = cuda::std::numeric_limits::max_exponent - 1; /** * @brief Reinterpret the bits of a floating-point value as an integer @@ -113,15 +112,15 @@ struct floating_converter { } /** - * @brief Extracts the integral significand of a bit-casted floating-point number + * @brief Checks whether the bit-casted floating-point value is +/-0 * - * @param integer_rep The bit-casted floating value to extract the exponent from - * @return The integral significand, bit-shifted to a (large) whole number + * @param integer_rep The bit-casted floating value to check if is +/-0 + * @return True if is a zero, else false */ - CUDF_HOST_DEVICE inline static IntegralType get_base2_value(IntegralType integer_rep) + CUDF_HOST_DEVICE inline static bool is_zero(IntegralType integer_rep) { - // Extract the significand, setting the high bit for the understood 1/2 - return (integer_rep & mantissa_mask) | understood_bit_mask; + // It's a zero if every non-sign bit is zero + return ((integer_rep & ~sign_mask) == 0); } /** @@ -137,40 +136,52 @@ struct floating_converter { } /** - * @brief Extracts the exponent of a bit-casted floating-point number + * @brief Extracts the significand and exponent of a bit-casted floating-point number * - * @note This returns INT_MIN for +/-0, +/-inf, NaN's, and denormals - * For all of these cases, the decimal fixed_point number should be set to zero + * @note This returns (1 - exponent_bias) for denormals. Zeros/inf/NaN not handled. * * @param integer_rep The bit-casted floating value to extract the exponent from - * @return The stored base-2 exponent, or INT_MIN for special values + * @return The stored base-2 exponent, or (1 - exponent_bias) for denormals */ - CUDF_HOST_DEVICE inline static int get_exp2(IntegralType integer_rep) + CUDF_HOST_DEVICE inline static std::pair get_significand_and_exp2( + IntegralType integer_rep) { - // First extract the exponent bits and handle its special values. - // To minimize branching, all of these special cases will return INT_MIN. - // For all of these cases, the decimal fixed_point number should be set to zero. + // Extract the significand + auto significand = (integer_rep & mantissa_mask); + + // Extract the exponent bits. auto const exponent_bits = integer_rep & exponent_mask; + + // Notes on special values of exponent_bits: + // bits = exponent_mask is +/-inf or NaN, but those are handled prior to input. + // bits = 0 is either a denormal (handled below) or a zero (handled earlier by caller). + int floating_exp2; if (exponent_bits == 0) { - // Because of the understood set-bit not stored in the mantissa, it is not possible - // to store the value zero directly. Instead both +/-0 and denormals are represented with - // the exponent bits set to zero. - // Thus it's fastest to just floor (generally unwanted) denormals to zero. - return INT_MIN; - } else if (exponent_bits == exponent_mask) { - //+/-inf and NaN values are stored with all of the exponent bits set. - // As none of these are representable by integers, we'll return the same value for all cases. - return INT_MIN; + // Denormal values are 2^(1 - exponent_bias) * Sum_i(B_i * 2^-i) + // Where i is the i-th mantissa bit (counting from the LEFT, starting at 1), + // and B_i is the value of that bit (0 or 1) + // So e.g. for the minimum denormal, only the lowest bit is set: + // FLT_TRUE_MIN = 2^(1 - 127) * 2^-23 = 2^-149 + // DBL_TRUE_MIN = 2^(1 - 1023) * 2^-52 = 2^-1074 + floating_exp2 = 1 - exponent_bias; + } else { + // Extract the exponent value: shift the bits down and subtract the bias. + auto const shifted_exponent_bits = exponent_bits >> num_mantissa_bits; + floating_exp2 = static_cast(shifted_exponent_bits) - exponent_bias; + + // Set the high bit for the understood 1/2 + significand |= understood_bit_mask; } - // Extract the exponent value: shift the bits down and subtract the bias. - using SignedIntegralType = cuda::std::make_signed_t; - SignedIntegralType const shifted_exponent_bits = exponent_bits >> num_mantissa_bits; - return shifted_exponent_bits - static_cast(exponent_bias); + // To convert the mantissa to an integer, we effectively applied #-mantissa-bits + // powers of 2 to convert the fractional value to an integer, so subtract them off here + int const exp2 = floating_exp2 - num_mantissa_bits; + + return {significand, exp2}; } /** - * @brief Sets the sign bit of a positive floating-point number + * @brief Sets the sign bit of a floating-point number * * @param floating The floating-point value to set the sign of. Must be positive. * @param is_negative The sign bit to set for the floating-point number @@ -192,38 +203,58 @@ struct floating_converter { /** * @brief Adds to the base-2 exponent of a floating-point number * + * @note Where called, the input is guaranteed to be a positive whole number. + * * @param floating The floating value to add to the exponent of. Must be positive. * @param exp2 The power-of-2 to add to the floating-point number * @return The input floating-point value * 2^exp2 */ CUDF_HOST_DEVICE inline static FloatingType add_exp2(FloatingType floating, int exp2) { + // Note that the input floating-point number is positive (& whole), so we don't have to + // worry about the sign here; the sign will be set later in set_is_negative() + // Convert floating to integer auto integer_rep = bit_cast_to_integer(floating); // Extract the currently stored (biased) exponent + using SignedType = std::make_signed_t; auto exponent_bits = integer_rep & exponent_mask; - auto stored_exp2 = exponent_bits >> num_mantissa_bits; + auto stored_exp2 = static_cast(exponent_bits >> num_mantissa_bits); // Add the additional power-of-2 stored_exp2 += exp2; // Check for exponent over/under-flow. - // Note that the input floating-point number is always positive, so we don't have to - // worry about the sign here; the sign will be set later in set_is_negative() if (stored_exp2 <= 0) { - return 0.0; - } else if (stored_exp2 >= unshifted_exponent_mask) { + // Denormal (zero handled prior to input) + + // Early out if bit shift will zero it anyway. + // Note: We must handle this explicitly, as too-large a bit-shift is UB + auto const bit_shift = -stored_exp2 + 1; //+1 due to understood bit set below + if (bit_shift > num_mantissa_bits) { return 0.0; } + + // Clear the exponent bits (zero means 2^-126/2^-1022 w/ no understood bit) + integer_rep &= (~exponent_mask); + + // The input floating-point number has an "understood" bit that we need to set + // prior to bit-shifting. Set the understood bit. + integer_rep |= understood_bit_mask; + + // Convert to denormal: bit shift off the low bits + integer_rep >>= bit_shift; + } else if (stored_exp2 >= static_cast(unshifted_exponent_mask)) { + // Overflow: Set infinity return cuda::std::numeric_limits::infinity(); } else { - // Clear existing exponent bits and set new ones - exponent_bits = stored_exp2 << num_mantissa_bits; + // Normal number: Clear existing exponent bits and set new ones + exponent_bits = static_cast(stored_exp2) << num_mantissa_bits; integer_rep &= (~exponent_mask); integer_rep |= exponent_bits; - - // Convert back to float - return bit_cast_to_floating(integer_rep); } + + // Convert back to float + return bit_cast_to_floating(integer_rep); } }; @@ -609,6 +640,641 @@ CUDF_HOST_DEVICE inline constexpr T divide_power10(T value, int exp10) } } +/** + * @brief Perform a bit-shift left, guarding against undefined behavior + * + * @tparam IntegerType Type of input unsigned integer value + * @param value The integer whose bits are being shifted + * @param bit_shift The number of bits to shift left + * @return The bit-shifted integer, except max value if overflow would occur + */ +template )> +CUDF_HOST_DEVICE inline IntegerType guarded_left_shift(IntegerType value, int bit_shift) +{ + // Bit shifts larger than this are undefined behavior + static constexpr int max_safe_bit_shift = cuda::std::numeric_limits::digits - 1; + return (bit_shift <= max_safe_bit_shift) ? value << bit_shift + : cuda::std::numeric_limits::max(); +} + +/** + * @brief Perform a bit-shift right, guarding against undefined behavior + * + * @tparam IntegerType Type of input unsigned integer value + * @param value The integer whose bits are being shifted + * @param bit_shift The number of bits to shift right + * @return The bit-shifted integer, which is zero on underflow + */ +template )> +CUDF_HOST_DEVICE inline IntegerType guarded_right_shift(IntegerType value, int bit_shift) +{ + // Bit shifts larger than this are undefined behavior + static constexpr int max_safe_bit_shift = cuda::std::numeric_limits::digits - 1; + return (bit_shift <= max_safe_bit_shift) ? value >> bit_shift + : cuda::std::numeric_limits::max(); +} + +/** + * @brief Helper struct with common constants needed by the floating <--> decimal conversions + */ +template +struct shifting_constants { + /// Whether the type is double + static constexpr bool is_double = cuda::std::is_same_v; + + /// Integer type that can hold the value of the significand + using IntegerRep = std::conditional_t; + + /// Num bits needed to hold the significand + static constexpr auto num_significand_bits = cuda::std::numeric_limits::digits; + + /// Shift data back and forth in space of a type with 2x the starting bits, to give us enough room + using ShiftingRep = std::conditional_t; + + // The significand of a float / double is 24 / 53 bits + // However, to uniquely represent each double / float as different #'s in decimal + // you need 17 / 9 digits (from std::numeric_limits::max_digits10) + // To represent 10^17 / 10^9, you need 57 / 30 bits + // So we need to keep track of this # of bits during shifting to ensure no info is lost + /// # bits needed to represent the value + static constexpr int num_rep_bits = is_double ? 57 : 30; + + // We will be alternately shifting our data back and forth by powers of 2 and 10 to convert + // between floating and decimal (see shifting functions for details). + // For float -> decimal, we want to start with our significand bits at the top of the + // num_rep_bits range, so that we don't lose information we need on intermediary right-shifts. + // For normal numbers, this bit shift is a fixed distance, defined by the understood 2^0 bit. + // For denormals this bit is not set, and must be determined for each value. + /// Bit shift needed to line-up value to the top of the representation range + static constexpr int normal_lineup_shift = num_rep_bits - num_significand_bits; + + // To iteratively shift back and forth, our 2's (bit-) and 10's (divide-/multiply-) shifts must + // be of nearly the same magnitude, or else we'll over-/under-flow our shifting integer + + // 2^10 is approximately 10^3, so the largest shifts will have a 10/3 ratio + // The difference between 2^10 and 10^3 is 1024/1000: 2.4% + // So every time we shift by 10 bits and 3 decimal places, the 2s shift is an extra 2.4% + + // This 2.4% error compounds each time we do an iteration. + // The min (normal) float is 2^-126. + // Min denormal: 2^-126 * 2^-23 (mantissa bits): 2^-149 = ~1.4E-45 + // With our 10/3 shifting ratio, 149 (bit-shifts) * (3 / 10) = 44.7 (10s-shifts) + // 10^(-44.7) = 2E-45, which is off by ~1.4x from 1.4E-45 + + // Similarly, the min (normal) double is 2^-1022. + // Min denormal: 2^-1022 * 2^-52 (mantissa bits): 2^-1074 = 4.94E-324 + // With our 10/3 shifting ratio, 1074 (bit-shifts) * (3 / 10) = 322.2 (10s-shifts) + // 10^(-322.2) = 6.4E-323, which is off by ~13.2x from 4.94E-324 + + // To account for this compounding error, we can either complicate our loop code (slow), + // or use extra bits (in the direction we're shifting the 2s!) to compensate: + // 4 extra bits for doubles (2^4 = 16 > 13.2x error), 1 extra for floats (2 > 1.4x error) + /// # buffer bits to account for shifting error + static constexpr int num_2s_shift_buffer_bits = is_double ? 4 : 1; + + // How much room do we have for shifting? + // Float: 64-bit ShiftingRep - 31 (rep + buffer) = 33 bits. 2^33 = 8.6E9 + // Double: 128-bit ShiftingRep - 61 (rep + buffer) = 67 bits. 2^67 = 1.5E20 + // Thus for double / float we can shift up to 20 / 9 decimal places at once + + // But, we need to stick to our 10-bits / 3-decimals shift ratio to not over/under-flow. + // To simplify our loop code, we'll keep to this ratio by instead shifting a max of + // 18 / 9 decimal places, for double / float (60 / 30 bits) + /// Max at-once decimal place shift + static constexpr int max_digits_shift = is_double ? 18 : 9; + /// Max at-once bit shift + static constexpr int max_bits_shift = max_digits_shift * 10 / 3; + + // Pre-calculate 10^max_digits_shift. Note that 10^18 / 10^9 fits within IntegerRep + /// 10^max_digits_shift + static constexpr auto max_digits_shift_pow = + multiply_power10(IntegerRep(1), max_digits_shift); +}; + +/** + * @brief Increment integer rep of floating point if conversion causes truncation + * + * @note This fixes problems like 1.2 (value = 1.1999...) at scale -1 -> 11 + * + * @tparam T Type of integer holding the floating-point significand + * @param integral_mantissa The integer representation of the floating-point significand + * @param exp2 The power of 2 that needs to be applied to the significand + * @param exp10 The power of 10 that needs to be applied to the significand + * @return significand, incremented if the conversion to decimal causes truncation + */ +template >* = nullptr> +CUDF_HOST_DEVICE T increment_on_truncation(T const integral_mantissa, + int const exp2, + int const exp10) +{ + // The user-supplied scale may truncate information, so we need to talk about rounding. + // We have chosen not to round, so we want 1.23456f with scale -4 to be decimal 12345 + + // But if we don't round at all, 1.2 (double) with scale -1 is 11 instead of 12! + // Why? Because 1.2 (double) is actually stored as 1.1999999... which we truncate to 1.1 + // While correct (given our choice to truncate), this is surprising and undesirable. + // This problem happens because 1.2 is not perfectly representable in floating point, + // and the value 1.199999... happened to be closer to 1.2 than the next value (1.2000...1...) + + // If the scale truncates information (we didn't choose to keep exactly 1.1999...), how + // do we make sure we store 1.2? All we have to do is add 1 ulp! (unit in the last place) + // Then 1.1999... becomes 1.2000...1... which truncates to 1.2. + // And if it had been 1.2000...1..., adding 1 ulp still truncates to 1.2, the result is unchanged. + + // The only way that this produces the incorrect result is if, when we entered 1.19999..., + // we truly meant 1.19999... (exactly, out to the very last bit), but then decided to truncate + // anyway. By choosing to truncate, you're saying you don't actually care about that level of + // precision, so being off by < 1 ulp should be just fine, compared to screwing up 1.2 with scale + // -1 -> 11 + + // So when does the user-supplied scale truncate info? + // For powers > 0: When the 10s (scale) shift is larger than the corresponding bit-shift. + // For powers < 0: When the 10s shift is less than the corresponding bit-shift. + + // Corresponding bit-shift: + // 2^10 is approximately 10^3, but this is off by 1.024% + // 1.024^30 is 2.03704, so this is high by one bit for every 30*3 = 90 powers of 10 + // So 10^N = 2^(10*N/3 - N/90) = 2^(299*N/90) + int const corresponding_exp2 = 299 * exp10 / 90; + + // If exp10 > 0, truncate if divide by more 10s than we shift up by 2s + // If exp10 < 0, truncate if shift down by more OR THE SAME 2s than multiply by 10s + // Truncate on the same: because for our approximation 2^299 > 10^90 + // Note that this works for both +/- exponents + bool const conversion_truncates = + (exp2 < corresponding_exp2) || ((exp2 == corresponding_exp2) && (exp2 < 0)); + + // (Potentially) increment and return + return integral_mantissa + static_cast(conversion_truncates); +} + +/** + * @brief Perform lossless base-2 -> base-10 fixed-point conversion for exp10 > 0 + * + * @note Info is lost if the chosen scale factor truncates information. + * + * @tparam FloatingType The type of the original floating-point value we are converting from + * @param base2_value The base-2 fixed-point value we are converting from + * @param exp2 The number of powers of 2 to apply to convert from base-2 + * @param exp10 The number of powers of 10 to apply to reach the desired scale factor + * @return Magnitude of the converted-to decimal integer + */ + +template >* = nullptr> +CUDF_HOST_DEVICE inline typename shifting_constants::ShiftingRep +shift_to_decimal_posexp(typename shifting_constants::IntegerRep const base2_value, + int exp2, + int exp10) +{ + // To convert to decimal, we need to apply the input powers of 2 and 10 + // The result will be (integer) base2_value * (2^exp2) / (10^exp10) + // Output type is ShiftingRep + + // Here exp10 > 0 and exp2 > 0, so we need to shift left by 2s and divide by 10s. + // To do this losslessly, we will iterate back and forth between them, shifting + // up by 2s and down by 10s until all of the powers have been applied. + + // However the input base2_value type has virtually no spare room to shift our data + // without over- or under-flowing and losing precision. + // So we'll cast up to ShiftingRep: uint64 for float's, __uint128_t for double's + using Constants = shifting_constants; + using ShiftingRep = typename Constants::ShiftingRep; + auto shifting_rep = static_cast(base2_value); + + // We want to start by lining up our bits in num_rep_bits (see comments on normal_lineup_shift), + // but since we start by bit-shifting up anyway, combine the normal_lineup_shift & max_bits_shift. + // Note that since we're shifting 2s up, we need num_2s_shift_buffer_bits space on the high side, + // which we do (our max bit shift is low enough that we don't shift into the highest bits) + static constexpr int max_init_shift = Constants::normal_lineup_shift + Constants::max_bits_shift; + + // If our total bit shift is less than this, we don't need to iterate + if (exp2 <= max_init_shift) { + // Shift bits left, divide by 10s to apply the scale factor, and we're done. + return divide_power10(shifting_rep << exp2, exp10); + } + + // We need to iterate. Do the combined initial shift + shifting_rep <<= max_init_shift; + exp2 -= max_init_shift; + + // Iterate, dividing by 10s and shifting up by 2s until we're almost done + while (exp10 > Constants::max_digits_shift) { + // More decimal places to shift than we have room: Divide the max number of 10s + + // Note that the result of this division is guaranteed to fit within the low half of the bits. + // The highest set bit is num_rep_bits + num_2s_shift_buffer_bits + max_bits_shift + // For float this is 30 + 1 + 30 = 61, for double 57 + 4 + 60 = 121 + // 2^61 / 10^9 (~2^30) is ~2^31, and 2^121 / 10^18 (~2^60) is ~2^61 + // As a future optimization, we could use a faster division routine that takes this account. + shifting_rep /= Constants::max_digits_shift_pow; + exp10 -= Constants::max_digits_shift; + + // If our remaining bit shift is less than the max, we're finished iterating + if (exp2 <= Constants::max_bits_shift) { + // Shift bits left, divide by 10s to apply the scale factor, and we're done. + // Note: This divide result may not fit in the low half of the bit range + return divide_power10(shifting_rep << exp2, exp10); + } + + // Shift the max number of bits left again + shifting_rep <<= Constants::max_bits_shift; + exp2 -= Constants::max_bits_shift; + } + + // Last 10s-shift: Divdie all remaining decimal places, shift all remaining bits, then bail + // Note: This divide result may not fit in the low half of the bit range + // But the divisor is less than the max-shift, and thus fits within 64 / 32 bits + if constexpr (Constants::is_double) { + shifting_rep = divide_power10_64bit(shifting_rep, exp10); + } else { + shifting_rep = divide_power10_32bit(shifting_rep, exp10); + } + + // Final bit shift: Shift may be large, guard against UB + // NOTE: This can overflow! + return guarded_left_shift(shifting_rep, exp2); +} + +/** + * @brief Perform lossless base-2 -> base-10 fixed-point conversion for exp10 < 0 + * + * @note Info is lost if the chosen scale factor truncates information. + * + * @tparam FloatingType The type of the original floating-point value we are converting from + * @param base2_value The base-2 fixed-point value we are converting from + * @param exp2 The number of powers of 2 to apply to convert from base-2 + * @param exp10 The number of powers of 10 to apply to reach the desired scale factor + * @return Magnitude of the converted-to decimal integer + */ +template >* = nullptr> +CUDF_HOST_DEVICE inline typename shifting_constants::ShiftingRep +shift_to_decimal_negexp(typename shifting_constants::IntegerRep const base2_value, + int exp2, + int exp10) +{ + // This is similar to shift_to_decimal_posexp(), except exp10 < 0 & exp2 < 0 + // See comments in that function for details. + // Instead here we need to multiply by 10s and shift right by 2s + + // Convert to using positive values so we don't have keep negating each time we multiply + int exp2_mag = -exp2; + int exp10_mag = -exp10; + + // ShiftingRep: uint64 for float's, __uint128_t for double's + using Constants = shifting_constants; + using ShiftingRep = typename Constants::ShiftingRep; + ShiftingRep shifting_rep; + + // For performing final 10s-shift + auto final_shifts_low10s = [&]() { + // Last 10s-shift: multiply all remaining decimal places, shift all remaining bits, then bail + // The multiplier is less than the max-shift, and thus fits within 64 / 32 bits + if constexpr (Constants::is_double) { + shifting_rep = multiply_power10_64bit(shifting_rep, exp10_mag); + } else { + shifting_rep = multiply_power10_32bit(shifting_rep, exp10_mag); + } + + // Final bit shift: Shift may be large, guard against UB + return guarded_right_shift(shifting_rep, exp2_mag); + }; + + // If our total decimal shift is less than the max, we don't need to iterate + if (exp10_mag <= Constants::max_digits_shift) { + shifting_rep = base2_value; + return final_shifts_low10s(); + } + + // We want to start by lining up our bits to num_rep_bits, but since we'll be bit-shifting + // down, we need even more low bits as a buffer (see comments on these constants) + auto const lineup_shift = Constants::num_rep_bits - count_significant_bits(base2_value); + auto const num_init_bit_shift = lineup_shift + Constants::num_2s_shift_buffer_bits; + + // Constants::num_2s_shift_buffer_bits; Note: This shift is safe to do in the smaller IntegerRep + // as it is up to bit 61 / 31 + shifting_rep = base2_value << num_init_bit_shift; + exp2_mag += num_init_bit_shift; + + // Iterate, multiplying by 10s and shifting down by 2s until we're almost done + do { + // More decimal places to shift than we have room: Multiply the max number of 10s + shifting_rep *= Constants::max_digits_shift_pow; + exp10_mag -= Constants::max_digits_shift; + + // If our remaining bit shift is less than the max, we're finished iterating + if (exp2_mag <= Constants::max_bits_shift) { + // Last bit-shift: Shift all remaining bits, apply the remaining scale, then bail + shifting_rep >>= exp2_mag; + + // We need to convert to the output rep for the final scale-factor multiply, because if (e.g.) + // float -> dec128 and some large exp10_mag, it might overflow the 64bit shifting rep. + // It's not needed for exp10 > 0 because we're dividing by 10s there instead of multiplying. + using UnsignedRep = cuda::std::make_unsigned_t; + // NOTE: This can overflow! (Both multiply and cast) + return multiply_power10(static_cast(shifting_rep), exp10_mag); + } + + // More bits to shift than we have room: Shift the max number of 2s + shifting_rep >>= Constants::max_bits_shift; + exp2_mag -= Constants::max_bits_shift; + } while (exp10_mag > Constants::max_digits_shift); + + // Do our final shifts + return final_shifts_low10s(); +} + +/** + * @brief Perform lossless floating-point -> integer decimal conversion + * + * @note Info is lost if the chosen scale factor truncates information. + * + * @tparam Rep The type of integer we are converting to, to store the decimal value + * @tparam FloatingType The type of floating-point object we are converting from + * @param floating The floating point value to convert + * @param scale The desired base-10 scale factor: decimal value = returned value * 10^scale + * @return Integer representation of the floating-point value, given the desired scale + */ +template >* = nullptr> +CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& floating, + scale_type const& scale) +{ + // Extract components of the floating point number + using converter = floating_converter; + auto const integer_rep = converter::bit_cast_to_integer(floating); + if (converter::is_zero(integer_rep)) { return 0; } + + // Note that the base2_value here is an unsigned integer with sizeof(FloatingType) + auto const is_negative = converter::get_is_negative(integer_rep); + auto const [base2_value, floating_exp2] = converter::get_significand_and_exp2(integer_rep); + + auto const exp2 = floating_exp2; // Can't capture from a parameter pack + auto const exp10 = static_cast(scale); + + // Increment if truncating to yield expected value, see function for discussion + auto const incremented = increment_on_truncation(base2_value, exp2, exp10); + + // Apply the powers of 2 and 10 to convert to decimal. + // The result will be incremented * (2^exp2) / (10^exp10) + // + // Note that while this code is branchy, the decimal scale factor is part of the + // column type itself, so every thread will take the same branches on exp10. + // Also data within a column tends to be similar, so they will often take the + // same branches on exp2 as well. + // + // NOTE: All returns here can overflow (e.g. unsigned -> signed) + auto const magnitude = [&]() -> Rep { + using UnsignedRep = cuda::std::make_unsigned_t; + + if (exp10 == 0) { + // NOTE: Left Bit-shift can overflow! As can cast! (e.g. double -> decimal32) + // Bit shifts may be large, guard against UB + if (exp2 >= 0) { + return guarded_left_shift(static_cast(incremented), exp2); + } else { + return guarded_right_shift(incremented, -exp2); + } + } else if (exp10 > 0) { + if (exp2 <= 0) { + // Power-2/10 shifts both downward: order doesn't matter, apply and bail. + // Guard against shift being undefined behavior + auto const shifted = guarded_right_shift(incremented, -exp2); + return divide_power10(shifted, exp10); + } + return shift_to_decimal_posexp(incremented, exp2, exp10); + } else { // exp10 < 0 + if (exp2 >= 0) { + // Power-2/10 shifts both upward: order doesn't matter, apply and bail. + // NOTE: Either shift, multiply, or cast (e.g. double -> decimal32) can overflow! + auto const shifted = guarded_left_shift(static_cast(incremented), exp2); + return multiply_power10(shifted, -exp10); + } + return shift_to_decimal_negexp(incremented, exp2, exp10); + } + }(); + + // Reapply the sign and return + return is_negative ? -magnitude : magnitude; +} + +/** + * @brief Perform (nearly) lossless base-10 -> base-2 fixed-point conversion for exp10 > 0 + * + * @note Intended to only be called internally. + * + * @tparam DecimalRep The decimal integer type we are converting from + * @tparam FloatingType The type of floating point object we are converting to + * @param decimal_rep The decimal integer to convert + * @param exp10 The number of powers of 10 to apply to undo the scale factor. + * @return A pair of the base-2 value and the remaining powers of 2 to be applied. + */ +template >* = nullptr> +CUDF_HOST_DEVICE inline auto shift_to_binary_posexp(DecimalRep decimal_rep, int exp10) +{ + // This is the reverse of shift_to_decimal_posexp(), see that for more details. + + // ShiftingRep: uint64 for float's, __uint128_t for double's + using Constants = shifting_constants; + using ShiftingRep = typename Constants::ShiftingRep; + + // We would start by lining up our data to num_rep_bits, but since we'll be bit-shifting + // down, we need even more low bits as a buffer (see comments on these constants) + auto const num_significant_bits = count_significant_bits(decimal_rep); + int exp2 = num_significant_bits - (Constants::num_rep_bits + Constants::num_2s_shift_buffer_bits); + + // Perform the initial bit shift + ShiftingRep shifting_rep; + if constexpr (sizeof(ShiftingRep) < sizeof(DecimalRep)) { + // Shift within DecimalRep before dropping to the smaller ShiftingRep + decimal_rep = (exp2 >= 0) ? (decimal_rep >> exp2) : (decimal_rep << -exp2); + shifting_rep = static_cast(decimal_rep); + } else { + // Scale up to ShiftingRep before shifting + shifting_rep = static_cast(decimal_rep); + shifting_rep = (exp2 >= 0) ? (shifting_rep >> exp2) : (shifting_rep << -exp2); + } + + // Iterate, multiplying by 10s and shifting down by 2s until we're almost done + while (exp10 > Constants::max_digits_shift) { + // More decimal places to shift than we have room: Multiply the max number of 10s + shifting_rep *= Constants::max_digits_shift_pow; + exp10 -= Constants::max_digits_shift; + + // Then make more room by bit shifting down by the max # of 2s + shifting_rep >>= Constants::max_bits_shift; + exp2 += Constants::max_bits_shift; + } + + // Last 10s-shift: multiply all remaining decimal places + // The multiplier is less than the max-shift, and thus fits within 64 / 32 bits + if constexpr (Constants::is_double) { + shifting_rep = multiply_power10_64bit(shifting_rep, exp10); + } else { + shifting_rep = multiply_power10_32bit(shifting_rep, exp10); + } + + // Our shifting_rep is now the integer mantissa, return it and the powers of 2 + return std::pair{shifting_rep, exp2}; +} + +/** + * @brief Perform (nearly) lossless base-10 -> base-2 fixed-point conversion for exp10 < 0 + * + * @note Intended to only be called internally. + * @note A 1-ulp loss may occur, but only for magnitudes E-270 or smaller. + * + * @tparam DecimalRep The decimal integer type we are converting from + * @tparam FloatingType The type of floating point object we are converting to + * @param decimal_rep The decimal integer to convert + * @param exp10 The number of powers of 10 to apply to undo the scale factor. + * @return A pair of the base-2 value and the remaining powers of 2 to be applied. + */ +template >* = nullptr> +CUDF_HOST_DEVICE inline auto shift_to_binary_negexp(DecimalRep decimal_rep, int const exp10) +{ + // This is the reverse of shift_to_decimal_negexp(), see that for more details. + + // ShiftingRep: uint64 for float's, __uint128_t for double's + using Constants = shifting_constants; + using ShiftingRep = typename Constants::ShiftingRep; + + // We would start by lining up our data to num_rep_bits, but if we originated with a floating + // number, we had to keep track of extra bits on the low-side because we were bit-shifting down. + // + // Those bits were not rounded. If we didn't truncate those bits before, we don't want to now + // either to ensure that we end up at the same floating point value that we started with. + // + // We could try to round here instead, but we don't know if we came from a floating-point value + // or not, so any rounding may not be desired. + // + // Note that here we are bit-shifting up, so we also need num_2s_shift_buffer_bits + // room on the high side. We have barely enough room for this for floats, but we're one bit + // over for doubles. So for doubles we'll keep one less bit on the low-side. + // + // This MAY cause a discrepancy in the last bit of our double from the value that we started with. + // However we only need all 4 bits for extremely large exponents + // (one bit to start + one extra bit every 90 powers of 10, so < E-270). + // And it's only a partial bit, and the eventual cast to double rounds, so we + // are often (always?) fine anyway (e.g. DBL_MIN & DBL_TRUE_MIN work fine). + // + // See comments on these constants for more details. + auto const num_significant_bits = count_significant_bits(decimal_rep); + int exp2 = num_significant_bits - (Constants::num_rep_bits + Constants::num_2s_shift_buffer_bits); + if constexpr (Constants::is_double) { ++exp2; } + + // Max bit shift left to give us the most room for shifting 10s: Multiply by 2s + exp2 -= Constants::max_bits_shift; + + // Perform the initial bit shift + ShiftingRep shifting_rep; + if constexpr (sizeof(ShiftingRep) < sizeof(DecimalRep)) { + // Shift within DecimalRep before dropping to the smaller ShiftingRep + decimal_rep = (exp2 >= 0) ? (decimal_rep >> exp2) : (decimal_rep << -exp2); + shifting_rep = static_cast(decimal_rep); + } else { + // Scale up to ShiftingRep before shifting + shifting_rep = static_cast(decimal_rep); + shifting_rep = (exp2 >= 0) ? (shifting_rep >> exp2) : (shifting_rep << -exp2); + } + + // Convert to using positive values upfront, simpler than doing later. + int exp10_mag = -exp10; + + // Iterate, dividing by 10s and shifting up by 2s until we're almost done + while (exp10_mag > Constants::max_digits_shift) { + // More decimal places to shift than we have room: Divide the max number of 10s + // Note that the result of this division is guaranteed to fit within low 64/32 bits + // See discussion in shift_to_decimal_posexp() for more details + shifting_rep /= Constants::max_digits_shift_pow; + exp10_mag -= Constants::max_digits_shift; + + // Then make more room by bit shifting up by the max # of 2s + shifting_rep <<= Constants::max_bits_shift; + exp2 -= Constants::max_bits_shift; + } + + // Last 10s-shift: Divdie all remaining decimal places. + // This divide result may not fit in the low half of the bit range + // But the divisor is less than the max-shift, and thus fits within 64 / 32 bits + if constexpr (Constants::is_double) { + shifting_rep = divide_power10_64bit(shifting_rep, exp10_mag); + } else { + shifting_rep = divide_power10_32bit(shifting_rep, exp10_mag); + } + + // Our shifting_rep is now the integer mantissa, return it and the powers of 2 + return std::pair{shifting_rep, exp2}; +} + +/** + * @brief Perform (nearly) lossless integer decimal -> floating-point conversion + * + * @note A 1 ulp loss may occur, but only to doubles with magnitude <= 1E-270 + * + * @tparam FloatingType The type of floating-point object we are converting to + * @tparam Rep The decimal integer type we are converting from + * @param value The decimal integer to convert + * @param scale The base-10 scale factor for the input integer + * @return Floating-point representation of the scaled integral value + */ +template >* = nullptr> +CUDF_HOST_DEVICE inline FloatingType convert_integral_to_floating(Rep const& value, + scale_type const& scale) +{ + // Check the sign of the input + bool const is_negative = (value < 0); + + // Convert to unsigned for bit counting/shifting + using UnsignedType = cuda::std::make_unsigned_t; + auto const unsigned_value = [&]() -> UnsignedType { + // Use built-in abs functions where available + if constexpr (cuda::std::is_same_v) { + return cuda::std::llabs(value); + } else if constexpr (!cuda::std::is_same_v) { + return cuda::std::abs(value); + } + + // No abs function for 128bit types, so have to do it manually. + // Must guard against minimum value, as we can't just negate it: not representable. + if (value == cuda::std::numeric_limits<__int128_t>::min()) { + return static_cast(value); + } else { + return static_cast(is_negative ? -value : value); + } + }(); + + // Shift by powers of 2 and 10 to get our integer mantissa + auto const [mantissa, exp2] = [&]() { + auto const exp10 = static_cast(scale); + if (exp10 >= 0) { + return shift_to_binary_posexp(unsigned_value, exp10); + } else { // exp10 < 0 + return shift_to_binary_negexp(unsigned_value, exp10); + } + }(); + + // Zero has special exponent bits, just handle it here + if (mantissa == 0) { return FloatingType(0.0f); } + + // Cast our integer mantissa to floating point + auto const floating = static_cast(mantissa); // IEEE-754 rounds to even + + // Apply the sign and the remaining powers of 2 + using converter = floating_converter; + auto const magnitude = converter::add_exp2(floating, exp2); + return converter::set_is_negative(magnitude, is_negative); +} + } // namespace detail /** @} */ // end of group diff --git a/cpp/include/cudf/unary.hpp b/cpp/include/cudf/unary.hpp index 74c8bc67d3a8..8a515335351d 100644 --- a/cpp/include/cudf/unary.hpp +++ b/cpp/include/cudf/unary.hpp @@ -17,6 +17,7 @@ #pragma once #include +#include #include #include #include @@ -50,14 +51,19 @@ namespace cudf { */ template () && - cuda::std::is_floating_point_v>* = nullptr> + CUDF_ENABLE_IF(cuda::std::is_floating_point_v&& is_fixed_point())> CUDF_HOST_DEVICE Fixed convert_floating_to_fixed(Floating floating, numeric::scale_type scale) { - using Rep = typename Fixed::rep; - auto const shifted = numeric::detail::shift(floating, scale); - numeric::scaled_integer scaled{static_cast(shifted), scale}; - return Fixed(scaled); + using Rep = typename Fixed::rep; + auto const value = [&]() { + if constexpr (Fixed::rad == numeric::Radix::BASE_10) { + return numeric::detail::convert_floating_to_integral(floating, scale); + } else { + return static_cast(numeric::detail::shift(floating, scale)); + } + }(); + + return Fixed(numeric::scaled_integer{value, scale}); } /** @@ -75,14 +81,17 @@ CUDF_HOST_DEVICE Fixed convert_floating_to_fixed(Floating floating, numeric::sca */ template && - is_fixed_point()>* = nullptr> + CUDF_ENABLE_IF(cuda::std::is_floating_point_v&& is_fixed_point())> CUDF_HOST_DEVICE Floating convert_fixed_to_floating(Fixed fixed) { - using Rep = typename Fixed::rep; - auto const casted = static_cast(fixed.value()); - auto const scale = numeric::scale_type{-fixed.scale()}; - return numeric::detail::shift(casted, scale); + using Rep = typename Fixed::rep; + if constexpr (Fixed::rad == numeric::Radix::BASE_10) { + return numeric::detail::convert_integral_to_floating(fixed.value(), fixed.scale()); + } else { + auto const casted = static_cast(fixed.value()); + auto const scale = numeric::scale_type{-fixed.scale()}; + return numeric::detail::shift(casted, scale); + } } /** @@ -95,7 +104,7 @@ CUDF_HOST_DEVICE Floating convert_fixed_to_floating(Fixed fixed) */ template >* = nullptr> + CUDF_ENABLE_IF(cuda::std::is_floating_point_v)> CUDF_HOST_DEVICE Floating convert_to_floating(Input input) { if constexpr (is_fixed_point()) { diff --git a/cpp/tests/fixed_point/fixed_point_tests.cpp b/cpp/tests/fixed_point/fixed_point_tests.cpp index 73de1fbaa68b..86a218b7c8bf 100644 --- a/cpp/tests/fixed_point/fixed_point_tests.cpp +++ b/cpp/tests/fixed_point/fixed_point_tests.cpp @@ -38,7 +38,7 @@ struct FixedPointTest : public cudf::test::BaseFixture {}; template struct FixedPointTestAllReps : public cudf::test::BaseFixture {}; -using RepresentationTypes = ::testing::Types; +using RepresentationTypes = ::testing::Types; TYPED_TEST_SUITE(FixedPointTestAllReps, RepresentationTypes); @@ -53,6 +53,7 @@ TYPED_TEST(FixedPointTestAllReps, SimpleDecimalXXConstruction) auto num4 = cudf::convert_floating_to_fixed(1.234567, scale_type(-4)); auto num5 = cudf::convert_floating_to_fixed(1.234567, scale_type(-5)); auto num6 = cudf::convert_floating_to_fixed(1.234567, scale_type(-6)); + auto num7 = cudf::convert_floating_to_fixed(0.0, scale_type(-4)); EXPECT_EQ(1, cudf::convert_fixed_to_floating(num0)); EXPECT_EQ(1.2, cudf::convert_fixed_to_floating(num1)); @@ -61,6 +62,7 @@ TYPED_TEST(FixedPointTestAllReps, SimpleDecimalXXConstruction) EXPECT_EQ(1.2345, cudf::convert_fixed_to_floating(num4)); EXPECT_EQ(1.23456, cudf::convert_fixed_to_floating(num5)); EXPECT_EQ(1.234567, cudf::convert_fixed_to_floating(num6)); + EXPECT_EQ(0.0, cudf::convert_fixed_to_floating(num7)); } TYPED_TEST(FixedPointTestAllReps, SimpleNegativeDecimalXXConstruction) @@ -74,6 +76,7 @@ TYPED_TEST(FixedPointTestAllReps, SimpleNegativeDecimalXXConstruction) auto num4 = cudf::convert_floating_to_fixed(-1.234567, scale_type(-4)); auto num5 = cudf::convert_floating_to_fixed(-1.234567, scale_type(-5)); auto num6 = cudf::convert_floating_to_fixed(-1.234567, scale_type(-6)); + auto num7 = cudf::convert_floating_to_fixed(-0.0, scale_type(-4)); EXPECT_EQ(-1, cudf::convert_fixed_to_floating(num0)); EXPECT_EQ(-1.2, cudf::convert_fixed_to_floating(num1)); @@ -82,6 +85,7 @@ TYPED_TEST(FixedPointTestAllReps, SimpleNegativeDecimalXXConstruction) EXPECT_EQ(-1.2345, cudf::convert_fixed_to_floating(num4)); EXPECT_EQ(-1.23456, cudf::convert_fixed_to_floating(num5)); EXPECT_EQ(-1.234567, cudf::convert_fixed_to_floating(num6)); + EXPECT_EQ(-0.0, cudf::convert_fixed_to_floating(num7)); } TYPED_TEST(FixedPointTestAllReps, PaddedDecimalXXConstruction) @@ -99,14 +103,10 @@ TYPED_TEST(FixedPointTestAllReps, PaddedDecimalXXConstruction) EXPECT_EQ(1.1, cudf::convert_fixed_to_floating(a)); EXPECT_EQ(1.01, cudf::convert_fixed_to_floating(b)); - EXPECT_EQ(1, - cudf::convert_fixed_to_floating( - c)); // intentional (inherited problem from floating point) + EXPECT_EQ(1.001, cudf::convert_fixed_to_floating(c)); EXPECT_EQ(1.0001, cudf::convert_fixed_to_floating(d)); EXPECT_EQ(1.00001, cudf::convert_fixed_to_floating(e)); - EXPECT_EQ(1, - cudf::convert_fixed_to_floating( - f)); // intentional (inherited problem from floating point) + EXPECT_EQ(1.000001, cudf::convert_fixed_to_floating(f)); EXPECT_TRUE(1.000123 - cudf::convert_fixed_to_floating(x) < std::numeric_limits::epsilon()); @@ -153,6 +153,118 @@ TYPED_TEST(FixedPointTestAllReps, MoreSimpleBinaryFPConstruction) EXPECT_EQ(2.0625, cudf::convert_fixed_to_floating(num1)); } +TEST_F(FixedPointTest, PreciseFloatDecimal64Construction) +{ + // Need 9 decimal digits to uniquely represent all floats (numeric_limits::max_digits10()). + // Precise conversion: set the scale factor to 9 less than the order-of-magnitude. + // But with -9 scale factor decimal32 can overflow: use decimal64 instead. + + // Positive Exponent + { + auto num0 = cudf::convert_floating_to_fixed(3.141593E7f, scale_type(-2)); + auto num1 = cudf::convert_floating_to_fixed(3.141593E12f, scale_type(3)); + auto num2 = cudf::convert_floating_to_fixed(3.141593E17f, scale_type(8)); + auto num3 = cudf::convert_floating_to_fixed(3.141593E22f, scale_type(13)); + auto num4 = cudf::convert_floating_to_fixed(3.141593E27f, scale_type(18)); + auto num5 = cudf::convert_floating_to_fixed(3.141593E32f, scale_type(23)); + auto num6 = cudf::convert_floating_to_fixed(FLT_MAX, scale_type(29)); + + EXPECT_EQ(3.141593E7f, cudf::convert_fixed_to_floating(num0)); + EXPECT_EQ(3.141593E12f, cudf::convert_fixed_to_floating(num1)); + EXPECT_EQ(3.141593E17f, cudf::convert_fixed_to_floating(num2)); + EXPECT_EQ(3.141593E22f, cudf::convert_fixed_to_floating(num3)); + EXPECT_EQ(3.141593E27f, cudf::convert_fixed_to_floating(num4)); + EXPECT_EQ(3.141593E32f, cudf::convert_fixed_to_floating(num5)); + EXPECT_EQ(FLT_MAX, cudf::convert_fixed_to_floating(num6)); + } + + // Negative Exponent + { + auto num0 = cudf::convert_floating_to_fixed(3.141593E-7f, scale_type(-16)); + auto num1 = cudf::convert_floating_to_fixed(3.141593E-12f, scale_type(-21)); + auto num2 = cudf::convert_floating_to_fixed(3.141593E-17f, scale_type(-26)); + auto num3 = cudf::convert_floating_to_fixed(3.141593E-22f, scale_type(-31)); + auto num4 = cudf::convert_floating_to_fixed(3.141593E-27f, scale_type(-36)); + auto num5 = cudf::convert_floating_to_fixed(3.141593E-32f, scale_type(-41)); + auto num6 = cudf::convert_floating_to_fixed(FLT_MIN, scale_type(-47)); + + EXPECT_EQ(3.141593E-7f, cudf::convert_fixed_to_floating(num0)); + EXPECT_EQ(3.141593E-12f, cudf::convert_fixed_to_floating(num1)); + EXPECT_EQ(3.141593E-17f, cudf::convert_fixed_to_floating(num2)); + EXPECT_EQ(3.141593E-22f, cudf::convert_fixed_to_floating(num3)); + EXPECT_EQ(3.141593E-27f, cudf::convert_fixed_to_floating(num4)); + EXPECT_EQ(3.141593E-32f, cudf::convert_fixed_to_floating(num5)); + EXPECT_EQ(FLT_MIN, cudf::convert_fixed_to_floating(num6)); + + // Denormals + auto num7 = cudf::convert_floating_to_fixed(3.141593E-39f, scale_type(-48)); + auto num8 = cudf::convert_floating_to_fixed(3.141593E-41f, scale_type(-50)); + auto num9 = cudf::convert_floating_to_fixed(3.141593E-43f, scale_type(-52)); + auto num10 = cudf::convert_floating_to_fixed(FLT_TRUE_MIN, scale_type(-54)); + + EXPECT_EQ(3.141593E-39f, cudf::convert_fixed_to_floating(num7)); + EXPECT_EQ(3.141593E-41f, cudf::convert_fixed_to_floating(num8)); + EXPECT_EQ(3.141593E-43f, cudf::convert_fixed_to_floating(num9)); + EXPECT_EQ(FLT_TRUE_MIN, cudf::convert_fixed_to_floating(num10)); + } +} + +TEST_F(FixedPointTest, PreciseDoubleDecimal64Construction) +{ + // Need 17 decimal digits to uniquely represent all doubles (numeric_limits::max_digits10()). + // Precise conversion: set the scale factor to 17 less than the order-of-magnitude. + + using decimal64 = fixed_point; + + // Positive Exponent + { + auto num0 = cudf::convert_floating_to_fixed(3.141593E8, scale_type(-9)); + auto num1 = cudf::convert_floating_to_fixed(3.141593E58, scale_type(41)); + auto num2 = cudf::convert_floating_to_fixed(3.141593E108, scale_type(91)); + auto num3 = cudf::convert_floating_to_fixed(3.141593E158, scale_type(141)); + auto num4 = cudf::convert_floating_to_fixed(3.141593E208, scale_type(191)); + auto num5 = cudf::convert_floating_to_fixed(3.141593E258, scale_type(241)); + auto num6 = cudf::convert_floating_to_fixed(DBL_MAX, scale_type(291)); + + EXPECT_EQ(3.141593E8, cudf::convert_fixed_to_floating(num0)); + EXPECT_EQ(3.141593E58, cudf::convert_fixed_to_floating(num1)); + EXPECT_EQ(3.141593E108, cudf::convert_fixed_to_floating(num2)); + EXPECT_EQ(3.141593E158, cudf::convert_fixed_to_floating(num3)); + EXPECT_EQ(3.141593E208, cudf::convert_fixed_to_floating(num4)); + EXPECT_EQ(3.141593E258, cudf::convert_fixed_to_floating(num5)); + EXPECT_EQ(DBL_MAX, cudf::convert_fixed_to_floating(num6)); + } + + // Negative Exponent + { + auto num0 = cudf::convert_floating_to_fixed(3.141593E-8, scale_type(-25)); + auto num1 = cudf::convert_floating_to_fixed(3.141593E-58, scale_type(-75)); + auto num2 = cudf::convert_floating_to_fixed(3.141593E-108, scale_type(-125)); + auto num3 = cudf::convert_floating_to_fixed(3.141593E-158, scale_type(-175)); + auto num4 = cudf::convert_floating_to_fixed(3.141593E-208, scale_type(-225)); + auto num5 = cudf::convert_floating_to_fixed(3.141593E-258, scale_type(-275)); + auto num6 = cudf::convert_floating_to_fixed(DBL_MIN, scale_type(-325)); + + EXPECT_EQ(3.141593E-8, cudf::convert_fixed_to_floating(num0)); + EXPECT_EQ(3.141593E-58, cudf::convert_fixed_to_floating(num1)); + EXPECT_EQ(3.141593E-108, cudf::convert_fixed_to_floating(num2)); + EXPECT_EQ(3.141593E-158, cudf::convert_fixed_to_floating(num3)); + EXPECT_EQ(3.141593E-208, cudf::convert_fixed_to_floating(num4)); + EXPECT_EQ(3.141593E-258, cudf::convert_fixed_to_floating(num5)); + EXPECT_EQ(DBL_MIN, cudf::convert_fixed_to_floating(num6)); + + // Denormals + auto num7 = cudf::convert_floating_to_fixed(3.141593E-309, scale_type(-326)); + auto num8 = cudf::convert_floating_to_fixed(3.141593E-314, scale_type(-331)); + auto num9 = cudf::convert_floating_to_fixed(3.141593E-319, scale_type(-336)); + auto num10 = cudf::convert_floating_to_fixed(DBL_TRUE_MIN, scale_type(-341)); + EXPECT_EQ(3.141593E-309, cudf::convert_fixed_to_floating(num7)); + EXPECT_EQ(3.141593E-314, cudf::convert_fixed_to_floating(num8)); + EXPECT_EQ(3.141593E-319, cudf::convert_fixed_to_floating(num9)); + EXPECT_EQ(DBL_TRUE_MIN, cudf::convert_fixed_to_floating(num10)); + } +} + TYPED_TEST(FixedPointTestAllReps, SimpleDecimalXXMath) { using decimalXX = fixed_point; @@ -442,8 +554,6 @@ void float_vector_test(ValueType const initial_value, int32_t const scale, Binop binop) { - using decimal32 = fixed_point; - std::vector vec1(size); std::vector vec2(size); From 26ed8fef03148de2c6c041a3083b5d64f97d0849 Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Mon, 3 Jun 2024 09:43:50 -0400 Subject: [PATCH 02/18] Remove xfail from python conversion test --- python/cudf/cudf/tests/test_decimal.py | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/python/cudf/cudf/tests/test_decimal.py b/python/cudf/cudf/tests/test_decimal.py index 0745e5aba48b..f76e936b3162 100644 --- a/python/cudf/cudf/tests/test_decimal.py +++ b/python/cudf/cudf/tests/test_decimal.py @@ -1,4 +1,4 @@ -# Copyright (c) 2021-2023, NVIDIA CORPORATION. +# Copyright (c) 2021-2024, NVIDIA CORPORATION. import decimal from decimal import Decimal @@ -6,7 +6,6 @@ import numpy as np import pyarrow as pa import pytest -from packaging import version import cudf from cudf.core.column import Decimal32Column, Decimal64Column, NumericalColumn @@ -93,14 +92,6 @@ def test_from_arrow_max_precision_decimal32(): [Decimal64Dtype(7, 2), Decimal64Dtype(11, 4), Decimal64Dtype(18, 9)], ) def test_typecast_from_float_to_decimal(request, data, from_dtype, to_dtype): - request.applymarker( - pytest.mark.xfail( - condition=version.parse(pa.__version__) >= version.parse("13.0.0") - and from_dtype == np.dtype("float32") - and to_dtype.precision > 7, - reason="https://github.com/rapidsai/cudf/issues/14169", - ) - ) got = data.astype(from_dtype) pa_arr = got.to_arrow().cast( From b77472b30748e0092f2c5b66a8a2462fd4a42912 Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Mon, 24 Jun 2024 15:21:19 -0400 Subject: [PATCH 03/18] Corrected increment, use CUDF_ENABLE_IF, denormals working --- .../cudf/fixed_point/floating_conversion.hpp | 274 +++++++++--------- cpp/tests/fixed_point/fixed_point_tests.cpp | 17 +- python/cudf/cudf/tests/test_decimal.py | 10 + 3 files changed, 153 insertions(+), 148 deletions(-) diff --git a/cpp/include/cudf/fixed_point/floating_conversion.hpp b/cpp/include/cudf/fixed_point/floating_conversion.hpp index 62499fc078f4..ffcac44ac718 100644 --- a/cpp/include/cudf/fixed_point/floating_conversion.hpp +++ b/cpp/include/cudf/fixed_point/floating_conversion.hpp @@ -34,6 +34,49 @@ namespace numeric { namespace detail { +/** + * @brief Determine the number of significant bits in an integer + * + * @tparam T Type of input integer value. Must be either uint32_t, uint64_t, or __uint128_t + * @param value The integer whose bits are being counted + * @return The number of significant bits: the # of bits - # of leading zeroes + */ +template || std::is_same_v || + std::is_same_v)> +CUDF_HOST_DEVICE inline int count_significant_bits(T value) +{ +#ifdef __CUDA_ARCH__ + if constexpr (std::is_same_v) { + return 64 - __clzll(static_cast(value)); + } else if constexpr (std::is_same_v) { + return 32 - __clz(static_cast(value)); + } else if constexpr (std::is_same_v) { + // 128 bit type, must break up into high and low components + auto const high_bits = static_cast(value >> 64); + auto const low_bits = static_cast(value); + return 128 - (__clzll(high_bits) + static_cast(high_bits == 0) * __clzll(low_bits)); + } +#else + // Undefined behavior to call __builtin_clzll() with zero in gcc and clang + if (value == 0) { return 0; } + + if constexpr (std::is_same_v) { + return 64 - __builtin_clzll(value); + } else if constexpr (std::is_same_v) { + return 32 - __builtin_clz(value); + } else if constexpr (std::is_same_v) { + // 128 bit type, must break up into high and low components + auto const high_bits = static_cast(value >> 64); + if (high_bits == 0) { + return 64 - __builtin_clzll(static_cast(value)); + } else { + return 128 - __builtin_clzll(high_bits); + } + } +#endif +} + /** * @brief Helper struct for getting and setting the components of a floating-point value * @@ -62,21 +105,23 @@ struct floating_converter { // The low 23 / 52 bits (for float / double) are the mantissa. // The mantissa is normalized. There is an understood 1 bit to the left of the binary point. // The value of the mantissa is in the range [1, 2). + /// # significand bits (includes understood bit) + static constexpr int num_significand_bits = cuda::std::numeric_limits::digits; /// # mantissa bits (-1 for understood bit) - static constexpr int num_mantissa_bits = cuda::std::numeric_limits::digits - 1; + static constexpr int num_stored_mantissa_bits = num_significand_bits - 1; /// The mask for the understood bit - static constexpr IntegralType understood_bit_mask = (IntegralType(1) << num_mantissa_bits); + static constexpr IntegralType understood_bit_mask = (IntegralType(1) << num_stored_mantissa_bits); /// The mask to select the mantissa static constexpr IntegralType mantissa_mask = understood_bit_mask - 1; // And in between are the bits used to store the biased power-of-2 exponent. /// # exponents bits (-1 for sign bit) - static constexpr int num_exponent_bits = num_floating_bits - num_mantissa_bits - 1; + static constexpr int num_exponent_bits = num_floating_bits - num_stored_mantissa_bits - 1; /// The mask for the exponents, unshifted static constexpr IntegralType unshifted_exponent_mask = (IntegralType(1) << num_exponent_bits) - 1; /// The mask to select the exponents - static constexpr IntegralType exponent_mask = unshifted_exponent_mask << num_mantissa_bits; + static constexpr IntegralType exponent_mask = unshifted_exponent_mask << num_stored_mantissa_bits; // To store positive and negative exponents as unsigned values, the stored value for // the power-of-2 is exponent + bias. The bias is 127 for floats and 1023 for doubles. @@ -136,12 +181,13 @@ struct floating_converter { } /** - * @brief Extracts the significand and exponent of a bit-casted floating-point number + * @brief Extracts the significand and exponent of a bit-casted floating-point number, + * shifted for denormals. * - * @note This returns (1 - exponent_bias) for denormals. Zeros/inf/NaN not handled. + * @note Zeros/inf/NaN not handled. * * @param integer_rep The bit-casted floating value to extract the exponent from - * @return The stored base-2 exponent, or (1 - exponent_bias) for denormals + * @return The stored base-2 exponent and significand, shifted for denormals */ CUDF_HOST_DEVICE inline static std::pair get_significand_and_exp2( IntegralType integer_rep) @@ -164,9 +210,15 @@ struct floating_converter { // FLT_TRUE_MIN = 2^(1 - 127) * 2^-23 = 2^-149 // DBL_TRUE_MIN = 2^(1 - 1023) * 2^-52 = 2^-1074 floating_exp2 = 1 - exponent_bias; + + // Line-up denormal to same (understood) bit as normal numbers + // This is so bit-shifting starts at the same bit index + auto const lineup_shift = num_significand_bits - count_significant_bits(significand); + significand <<= lineup_shift; + floating_exp2 -= lineup_shift; } else { // Extract the exponent value: shift the bits down and subtract the bias. - auto const shifted_exponent_bits = exponent_bits >> num_mantissa_bits; + auto const shifted_exponent_bits = exponent_bits >> num_stored_mantissa_bits; floating_exp2 = static_cast(shifted_exponent_bits) - exponent_bias; // Set the high bit for the understood 1/2 @@ -175,7 +227,7 @@ struct floating_converter { // To convert the mantissa to an integer, we effectively applied #-mantissa-bits // powers of 2 to convert the fractional value to an integer, so subtract them off here - int const exp2 = floating_exp2 - num_mantissa_bits; + int const exp2 = floating_exp2 - num_stored_mantissa_bits; return {significand, exp2}; } @@ -203,7 +255,7 @@ struct floating_converter { /** * @brief Adds to the base-2 exponent of a floating-point number * - * @note Where called, the input is guaranteed to be a positive whole number. + * @note The caller must guarantee that the input is a positive (> 0) whole number. * * @param floating The floating value to add to the exponent of. Must be positive. * @param exp2 The power-of-2 to add to the floating-point number @@ -220,7 +272,7 @@ struct floating_converter { // Extract the currently stored (biased) exponent using SignedType = std::make_signed_t; auto exponent_bits = integer_rep & exponent_mask; - auto stored_exp2 = static_cast(exponent_bits >> num_mantissa_bits); + auto stored_exp2 = static_cast(exponent_bits >> num_stored_mantissa_bits); // Add the additional power-of-2 stored_exp2 += exp2; @@ -232,7 +284,7 @@ struct floating_converter { // Early out if bit shift will zero it anyway. // Note: We must handle this explicitly, as too-large a bit-shift is UB auto const bit_shift = -stored_exp2 + 1; //+1 due to understood bit set below - if (bit_shift > num_mantissa_bits) { return 0.0; } + if (bit_shift > num_stored_mantissa_bits) { return 0.0; } // Clear the exponent bits (zero means 2^-126/2^-1022 w/ no understood bit) integer_rep &= (~exponent_mask); @@ -248,7 +300,7 @@ struct floating_converter { return cuda::std::numeric_limits::infinity(); } else { // Normal number: Clear existing exponent bits and set new ones - exponent_bits = static_cast(stored_exp2) << num_mantissa_bits; + exponent_bits = static_cast(stored_exp2) << num_stored_mantissa_bits; integer_rep &= (~exponent_mask); integer_rep |= exponent_bits; } @@ -258,49 +310,6 @@ struct floating_converter { } }; -/** - * @brief Determine the number of significant bits in an integer - * - * @tparam T Type of input integer value. Must be either uint32_t, uint64_t, or __uint128_t - * @param value The integer whose bits are being counted - * @return The number of significant bits: the # of bits - # of leading zeroes - */ -template || std::is_same_v || - std::is_same_v)> -CUDF_HOST_DEVICE inline int count_significant_bits(T value) -{ -#ifdef __CUDA_ARCH__ - if constexpr (std::is_same_v) { - return 64 - __clzll(static_cast(value)); - } else if constexpr (std::is_same_v) { - return 32 - __clz(static_cast(value)); - } else if constexpr (std::is_same_v) { - // 128 bit type, must break up into high and low components - auto const high_bits = static_cast(value >> 64); - auto const low_bits = static_cast(value); - return 128 - (__clzll(high_bits) + static_cast(high_bits == 0) * __clzll(low_bits)); - } -#else - // Undefined behavior to call __builtin_clzll() with zero in gcc and clang - if (value == 0) { return 0; } - - if constexpr (std::is_same_v) { - return 64 - __builtin_clzll(value); - } else if constexpr (std::is_same_v) { - return 32 - __builtin_clz(value); - } else if constexpr (std::is_same_v) { - // 128 bit type, must break up into high and low components - auto const high_bits = static_cast(value >> 64); - if (high_bits == 0) { - return 64 - __builtin_clzll(static_cast(value)); - } else { - return 128 - __builtin_clzll(high_bits); - } - } -#endif -} - /** * @brief Recursively calculate a signed large power of 10 (>= 10^19) that can only be stored in an * 128bit integer @@ -329,7 +338,7 @@ constexpr __uint128_t large_power_of_10() * @param exp10 The power-of-10 of the denominator, from 0 to 9 inclusive. * @return Returns value / 10^exp10 */ -template >* = nullptr> +template )> CUDF_HOST_DEVICE inline T divide_power10_32bit(T value, int exp10) { // Computing division this way is much faster than the alternatives. @@ -379,7 +388,7 @@ CUDF_HOST_DEVICE inline T divide_power10_32bit(T value, int exp10) * @param exp10 The power-of-10 of the denominator, from 0 to 19 inclusive. * @return Returns value / 10^exp10 */ -template >* = nullptr> +template )> CUDF_HOST_DEVICE inline T divide_power10_64bit(T value, int exp10) { // See comments in divide_power10_32bit() for discussion. @@ -416,7 +425,7 @@ CUDF_HOST_DEVICE inline T divide_power10_64bit(T value, int exp10) * @param exp10 The power-of-10 of the denominator, from 0 to 38 inclusive. * @return Returns value / 10^exp10. */ -template >* = nullptr> +template )> CUDF_HOST_DEVICE inline constexpr T divide_power10_128bit(T value, int exp10) { // See comments in divide_power10_32bit() for an introduction. @@ -472,7 +481,7 @@ CUDF_HOST_DEVICE inline constexpr T divide_power10_128bit(T value, int exp10) * @param exp10 The power-of-10 of the multiplier, from 0 to 9 inclusive. * @return Returns value * 10^exp10 */ -template >* = nullptr> +template )> CUDF_HOST_DEVICE inline constexpr T multiply_power10_32bit(T value, int exp10) { // See comments in divide_power10_32bit() for discussion. @@ -499,7 +508,7 @@ CUDF_HOST_DEVICE inline constexpr T multiply_power10_32bit(T value, int exp10) * @param exp10 The power-of-10 of the multiplier, from 0 to 19 inclusive. * @return Returns value * 10^exp10 */ -template >* = nullptr> +template )> CUDF_HOST_DEVICE inline constexpr T multiply_power10_64bit(T value, int exp10) { // See comments in divide_power10_32bit() for discussion. @@ -536,7 +545,7 @@ CUDF_HOST_DEVICE inline constexpr T multiply_power10_64bit(T value, int exp10) * @param exp10 The power-of-10 of the multiplier, from 0 to 38 inclusive. * @return Returns value * 10^exp10. */ -template >* = nullptr> +template )> CUDF_HOST_DEVICE inline constexpr T multiply_power10_128bit(T value, int exp10) { // See comments in divide_power10_128bit() for discussion. @@ -596,9 +605,7 @@ CUDF_HOST_DEVICE inline constexpr T multiply_power10_128bit(T value, int exp10) * @param exp10 The power-of-10 of the multiplier. * @return Returns value * 10^exp10 */ -template )>* = nullptr> +template )> CUDF_HOST_DEVICE inline constexpr T multiply_power10(T value, int exp10) { // Use this function if you have no knowledge of what exp10 might be @@ -624,9 +631,7 @@ CUDF_HOST_DEVICE inline constexpr T multiply_power10(T value, int exp10) * @param exp10 The power-of-10 of the denominator. * @return Returns value / 10^exp10 */ -template )>* = nullptr> +template )> CUDF_HOST_DEVICE inline constexpr T divide_power10(T value, int exp10) { // Use this function if you have no knowledge of what exp10 might be @@ -646,13 +651,13 @@ CUDF_HOST_DEVICE inline constexpr T divide_power10(T value, int exp10) * @tparam IntegerType Type of input unsigned integer value * @param value The integer whose bits are being shifted * @param bit_shift The number of bits to shift left - * @return The bit-shifted integer, except max value if overflow would occur + * @return The bit-shifted integer, except max value if UB would occur */ template )> CUDF_HOST_DEVICE inline IntegerType guarded_left_shift(IntegerType value, int bit_shift) { // Bit shifts larger than this are undefined behavior - static constexpr int max_safe_bit_shift = cuda::std::numeric_limits::digits - 1; + constexpr int max_safe_bit_shift = cuda::std::numeric_limits::digits - 1; return (bit_shift <= max_safe_bit_shift) ? value << bit_shift : cuda::std::numeric_limits::max(); } @@ -669,9 +674,8 @@ template ::digits - 1; - return (bit_shift <= max_safe_bit_shift) ? value >> bit_shift - : cuda::std::numeric_limits::max(); + constexpr int max_safe_bit_shift = cuda::std::numeric_limits::digits - 1; + return (bit_shift <= max_safe_bit_shift) ? value >> bit_shift : 0; } /** @@ -695,18 +699,10 @@ struct shifting_constants { // However, to uniquely represent each double / float as different #'s in decimal // you need 17 / 9 digits (from std::numeric_limits::max_digits10) // To represent 10^17 / 10^9, you need 57 / 30 bits - // So we need to keep track of this # of bits during shifting to ensure no info is lost - /// # bits needed to represent the value - static constexpr int num_rep_bits = is_double ? 57 : 30; + // So we need to keep track of at least this # of bits during shifting to ensure no info is lost // We will be alternately shifting our data back and forth by powers of 2 and 10 to convert // between floating and decimal (see shifting functions for details). - // For float -> decimal, we want to start with our significand bits at the top of the - // num_rep_bits range, so that we don't lose information we need on intermediary right-shifts. - // For normal numbers, this bit shift is a fixed distance, defined by the understood 2^0 bit. - // For denormals this bit is not set, and must be determined for each value. - /// Bit shift needed to line-up value to the top of the representation range - static constexpr int normal_lineup_shift = num_rep_bits - num_significand_bits; // To iteratively shift back and forth, our 2's (bit-) and 10's (divide-/multiply-) shifts must // be of nearly the same magnitude, or else we'll over-/under-flow our shifting integer @@ -762,7 +758,7 @@ struct shifting_constants { * @param exp10 The power of 10 that needs to be applied to the significand * @return significand, incremented if the conversion to decimal causes truncation */ -template >* = nullptr> +template )> CUDF_HOST_DEVICE T increment_on_truncation(T const integral_mantissa, int const exp2, int const exp10) @@ -781,28 +777,21 @@ CUDF_HOST_DEVICE T increment_on_truncation(T const integral_mantissa, // Then 1.1999... becomes 1.2000...1... which truncates to 1.2. // And if it had been 1.2000...1..., adding 1 ulp still truncates to 1.2, the result is unchanged. - // The only way that this produces the incorrect result is if, when we entered 1.19999..., - // we truly meant 1.19999... (exactly, out to the very last bit), but then decided to truncate - // anyway. By choosing to truncate, you're saying you don't actually care about that level of - // precision, so being off by < 1 ulp should be just fine, compared to screwing up 1.2 with scale - // -1 -> 11 - - // So when does the user-supplied scale truncate info? - // For powers > 0: When the 10s (scale) shift is larger than the corresponding bit-shift. - // For powers < 0: When the 10s shift is less than the corresponding bit-shift. + // If we add 1 to the last bit, we are effectively adding 1/2 to the 2nd-to-last bit. + // This is like rounding to the 2nd-to-last bit (if we were to floor afterwards). + // However we don't want to add 1 to the last bit if we are keeping enough precision. + // So, only add 1 if exp10 shift is larger than corresponding exp2 shift to 2nd-to-last bit. // Corresponding bit-shift: // 2^10 is approximately 10^3, but this is off by 1.024% // 1.024^30 is 2.03704, so this is high by one bit for every 30*3 = 90 powers of 10 // So 10^N = 2^(10*N/3 - N/90) = 2^(299*N/90) - int const corresponding_exp2 = 299 * exp10 / 90; - - // If exp10 > 0, truncate if divide by more 10s than we shift up by 2s - // If exp10 < 0, truncate if shift down by more OR THE SAME 2s than multiply by 10s - // Truncate on the same: because for our approximation 2^299 > 10^90 - // Note that this works for both +/- exponents + // Do comparison without dividing, which loses information: + // Note: if shift is "equal," still truncates if exp2 < 0 (shifting UP by 2s, 2^10 > 10^3) + int const exp2_term = 90 * (exp2 + 1); //+1: effectively adding 1/2 to 2nd-to-last-bit + int const exp10_term = 299 * exp10; bool const conversion_truncates = - (exp2 < corresponding_exp2) || ((exp2 == corresponding_exp2) && (exp2 < 0)); + (exp10_term > exp2_term) || ((exp2_term == exp10_term) && (exp2 < 0)); // (Potentially) increment and return return integral_mantissa + static_cast(conversion_truncates); @@ -819,9 +808,7 @@ CUDF_HOST_DEVICE T increment_on_truncation(T const integral_mantissa, * @param exp10 The number of powers of 10 to apply to reach the desired scale factor * @return Magnitude of the converted-to decimal integer */ - -template >* = nullptr> +template )> CUDF_HOST_DEVICE inline typename shifting_constants::ShiftingRep shift_to_decimal_posexp(typename shifting_constants::IntegerRep const base2_value, int exp2, @@ -842,11 +829,13 @@ shift_to_decimal_posexp(typename shifting_constants::IntegerRep co using ShiftingRep = typename Constants::ShiftingRep; auto shifting_rep = static_cast(base2_value); - // We want to start by lining up our bits in num_rep_bits (see comments on normal_lineup_shift), - // but since we start by bit-shifting up anyway, combine the normal_lineup_shift & max_bits_shift. + // We want to start with our significand bits at the top of the shifting range, + // so that we don't lose information we need on intermediary right-shifts. // Note that since we're shifting 2s up, we need num_2s_shift_buffer_bits space on the high side, - // which we do (our max bit shift is low enough that we don't shift into the highest bits) - static constexpr int max_init_shift = Constants::normal_lineup_shift + Constants::max_bits_shift; + // For all numbers this bit shift is a fixed distance, due to the understood 2^0 bit. + static constexpr int shift_up_to = sizeof(ShiftingRep) * 8 - Constants::num_2s_shift_buffer_bits; + static constexpr int shift_from = Constants::num_significand_bits; + static constexpr int max_init_shift = shift_up_to - shift_from; // If our total bit shift is less than this, we don't need to iterate if (exp2 <= max_init_shift) { @@ -909,9 +898,9 @@ shift_to_decimal_posexp(typename shifting_constants::IntegerRep co */ template >* = nullptr> + CUDF_ENABLE_IF(cuda::std::is_floating_point_v)> CUDF_HOST_DEVICE inline typename shifting_constants::ShiftingRep -shift_to_decimal_negexp(typename shifting_constants::IntegerRep const base2_value, +shift_to_decimal_negexp(typename shifting_constants::IntegerRep base2_value, int exp2, int exp10) { @@ -919,14 +908,14 @@ shift_to_decimal_negexp(typename shifting_constants::IntegerRep co // See comments in that function for details. // Instead here we need to multiply by 10s and shift right by 2s - // Convert to using positive values so we don't have keep negating each time we multiply - int exp2_mag = -exp2; - int exp10_mag = -exp10; - // ShiftingRep: uint64 for float's, __uint128_t for double's using Constants = shifting_constants; using ShiftingRep = typename Constants::ShiftingRep; - ShiftingRep shifting_rep; + auto shifting_rep = static_cast(base2_value); + + // Convert to using positive values so we don't have keep negating + int exp10_mag = -exp10; + int exp2_mag = -exp2; // For performing final 10s-shift auto final_shifts_low10s = [&]() { @@ -943,19 +932,18 @@ shift_to_decimal_negexp(typename shifting_constants::IntegerRep co }; // If our total decimal shift is less than the max, we don't need to iterate - if (exp10_mag <= Constants::max_digits_shift) { - shifting_rep = base2_value; - return final_shifts_low10s(); - } + if (exp10_mag <= Constants::max_digits_shift) { return final_shifts_low10s(); } - // We want to start by lining up our bits to num_rep_bits, but since we'll be bit-shifting - // down, we need even more low bits as a buffer (see comments on these constants) - auto const lineup_shift = Constants::num_rep_bits - count_significant_bits(base2_value); - auto const num_init_bit_shift = lineup_shift + Constants::num_2s_shift_buffer_bits; + // We want to start by lining up our bits to the top of the shifting range, + // except our first operation is a multiply, so not quite that far + // We are bit-shifting down, so we need extra bits on the low-side, which this has. + static constexpr int shift_up_to = sizeof(ShiftingRep) * 8 - Constants::max_bits_shift; + static constexpr int shift_from = Constants::num_significand_bits; + static constexpr int num_init_bit_shift = shift_up_to - shift_from; // Constants::num_2s_shift_buffer_bits; Note: This shift is safe to do in the smaller IntegerRep // as it is up to bit 61 / 31 - shifting_rep = base2_value << num_init_bit_shift; + shifting_rep <<= num_init_bit_shift; exp2_mag += num_init_bit_shift; // Iterate, multiplying by 10s and shifting down by 2s until we're almost done @@ -999,7 +987,7 @@ shift_to_decimal_negexp(typename shifting_constants::IntegerRep co */ template >* = nullptr> + CUDF_ENABLE_IF(cuda::std::is_floating_point_v)> CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& floating, scale_type const& scale) { @@ -1026,10 +1014,9 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo // Also data within a column tends to be similar, so they will often take the // same branches on exp2 as well. // - // NOTE: All returns here can overflow (e.g. unsigned -> signed) - auto const magnitude = [&]() -> Rep { - using UnsignedRep = cuda::std::make_unsigned_t; - + // NOTE: some returns here can overflow (e.g. ShiftingRep -> UnsignedRep) + using UnsignedRep = cuda::std::make_unsigned_t; + auto const magnitude = [&]() -> UnsignedRep { if (exp10 == 0) { // NOTE: Left Bit-shift can overflow! As can cast! (e.g. double -> decimal32) // Bit shifts may be large, guard against UB @@ -1058,7 +1045,9 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo }(); // Reapply the sign and return - return is_negative ? -magnitude : magnitude; + // NOTE: Cast can overflow! + auto const signed_magnitude = static_cast(magnitude); + return is_negative ? -signed_magnitude : signed_magnitude; } /** @@ -1074,7 +1063,7 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo */ template >* = nullptr> + CUDF_ENABLE_IF(cuda::std::is_floating_point_v)> CUDF_HOST_DEVICE inline auto shift_to_binary_posexp(DecimalRep decimal_rep, int exp10) { // This is the reverse of shift_to_decimal_posexp(), see that for more details. @@ -1083,10 +1072,13 @@ CUDF_HOST_DEVICE inline auto shift_to_binary_posexp(DecimalRep decimal_rep, int using Constants = shifting_constants; using ShiftingRep = typename Constants::ShiftingRep; - // We would start by lining up our data to num_rep_bits, but since we'll be bit-shifting - // down, we need even more low bits as a buffer (see comments on these constants) - auto const num_significant_bits = count_significant_bits(decimal_rep); - int exp2 = num_significant_bits - (Constants::num_rep_bits + Constants::num_2s_shift_buffer_bits); + // We want to start by lining up our bits to the top of the shifting range, + // except our first operation is a multiply, so not quite that far + // We are bit-shifting down, so we need extra bits on the low-side, which this has. + static constexpr int shift_up_to = sizeof(ShiftingRep) * 8 - Constants::max_bits_shift; + int const shift_from = count_significant_bits(decimal_rep); + int const num_init_bit_shift = shift_up_to - shift_from; + int exp2 = -num_init_bit_shift; // Perform the initial bit shift ShiftingRep shifting_rep; @@ -1137,7 +1129,7 @@ CUDF_HOST_DEVICE inline auto shift_to_binary_posexp(DecimalRep decimal_rep, int */ template >* = nullptr> + CUDF_ENABLE_IF(cuda::std::is_floating_point_v)> CUDF_HOST_DEVICE inline auto shift_to_binary_negexp(DecimalRep decimal_rep, int const exp10) { // This is the reverse of shift_to_decimal_negexp(), see that for more details. @@ -1166,12 +1158,14 @@ CUDF_HOST_DEVICE inline auto shift_to_binary_negexp(DecimalRep decimal_rep, int // are often (always?) fine anyway (e.g. DBL_MIN & DBL_TRUE_MIN work fine). // // See comments on these constants for more details. - auto const num_significant_bits = count_significant_bits(decimal_rep); - int exp2 = num_significant_bits - (Constants::num_rep_bits + Constants::num_2s_shift_buffer_bits); - if constexpr (Constants::is_double) { ++exp2; } - // Max bit shift left to give us the most room for shifting 10s: Multiply by 2s - exp2 -= Constants::max_bits_shift; + // We want to start with our significand bits at the top of the shifting range, + // so that we don't lose information we need on intermediary right-shifts. + // Note that since we're shifting 2s up, we need num_2s_shift_buffer_bits space on the high side, + static constexpr int shift_up_to = sizeof(ShiftingRep) * 8 - Constants::num_2s_shift_buffer_bits; + int const shift_from = count_significant_bits(decimal_rep); + int const num_init_bit_shift = shift_up_to - shift_from; + int exp2 = -num_init_bit_shift; // Perform the initial bit shift ShiftingRep shifting_rep; @@ -1227,7 +1221,7 @@ CUDF_HOST_DEVICE inline auto shift_to_binary_negexp(DecimalRep decimal_rep, int */ template >* = nullptr> + CUDF_ENABLE_IF(cuda::std::is_floating_point_v)> CUDF_HOST_DEVICE inline FloatingType convert_integral_to_floating(Rep const& value, scale_type const& scale) { diff --git a/cpp/tests/fixed_point/fixed_point_tests.cpp b/cpp/tests/fixed_point/fixed_point_tests.cpp index 86a218b7c8bf..a24e42ac306b 100644 --- a/cpp/tests/fixed_point/fixed_point_tests.cpp +++ b/cpp/tests/fixed_point/fixed_point_tests.cpp @@ -167,7 +167,7 @@ TEST_F(FixedPointTest, PreciseFloatDecimal64Construction) auto num3 = cudf::convert_floating_to_fixed(3.141593E22f, scale_type(13)); auto num4 = cudf::convert_floating_to_fixed(3.141593E27f, scale_type(18)); auto num5 = cudf::convert_floating_to_fixed(3.141593E32f, scale_type(23)); - auto num6 = cudf::convert_floating_to_fixed(FLT_MAX, scale_type(29)); + auto num6 = cudf::convert_floating_to_fixed(3.141593E37f, scale_type(28)); EXPECT_EQ(3.141593E7f, cudf::convert_fixed_to_floating(num0)); EXPECT_EQ(3.141593E12f, cudf::convert_fixed_to_floating(num1)); @@ -175,7 +175,7 @@ TEST_F(FixedPointTest, PreciseFloatDecimal64Construction) EXPECT_EQ(3.141593E22f, cudf::convert_fixed_to_floating(num3)); EXPECT_EQ(3.141593E27f, cudf::convert_fixed_to_floating(num4)); EXPECT_EQ(3.141593E32f, cudf::convert_fixed_to_floating(num5)); - EXPECT_EQ(FLT_MAX, cudf::convert_fixed_to_floating(num6)); + EXPECT_EQ(3.141593E37f, cudf::convert_fixed_to_floating(num6)); } // Negative Exponent @@ -186,7 +186,7 @@ TEST_F(FixedPointTest, PreciseFloatDecimal64Construction) auto num3 = cudf::convert_floating_to_fixed(3.141593E-22f, scale_type(-31)); auto num4 = cudf::convert_floating_to_fixed(3.141593E-27f, scale_type(-36)); auto num5 = cudf::convert_floating_to_fixed(3.141593E-32f, scale_type(-41)); - auto num6 = cudf::convert_floating_to_fixed(FLT_MIN, scale_type(-47)); + auto num6 = cudf::convert_floating_to_fixed(3.141593E-37f, scale_type(-47)); EXPECT_EQ(3.141593E-7f, cudf::convert_fixed_to_floating(num0)); EXPECT_EQ(3.141593E-12f, cudf::convert_fixed_to_floating(num1)); @@ -194,7 +194,7 @@ TEST_F(FixedPointTest, PreciseFloatDecimal64Construction) EXPECT_EQ(3.141593E-22f, cudf::convert_fixed_to_floating(num3)); EXPECT_EQ(3.141593E-27f, cudf::convert_fixed_to_floating(num4)); EXPECT_EQ(3.141593E-32f, cudf::convert_fixed_to_floating(num5)); - EXPECT_EQ(FLT_MIN, cudf::convert_fixed_to_floating(num6)); + EXPECT_EQ(3.141593E-37f, cudf::convert_fixed_to_floating(num6)); // Denormals auto num7 = cudf::convert_floating_to_fixed(3.141593E-39f, scale_type(-48)); @@ -224,7 +224,7 @@ TEST_F(FixedPointTest, PreciseDoubleDecimal64Construction) auto num3 = cudf::convert_floating_to_fixed(3.141593E158, scale_type(141)); auto num4 = cudf::convert_floating_to_fixed(3.141593E208, scale_type(191)); auto num5 = cudf::convert_floating_to_fixed(3.141593E258, scale_type(241)); - auto num6 = cudf::convert_floating_to_fixed(DBL_MAX, scale_type(291)); + auto num6 = cudf::convert_floating_to_fixed(3.141593E307, scale_type(290)); EXPECT_EQ(3.141593E8, cudf::convert_fixed_to_floating(num0)); EXPECT_EQ(3.141593E58, cudf::convert_fixed_to_floating(num1)); @@ -232,7 +232,7 @@ TEST_F(FixedPointTest, PreciseDoubleDecimal64Construction) EXPECT_EQ(3.141593E158, cudf::convert_fixed_to_floating(num3)); EXPECT_EQ(3.141593E208, cudf::convert_fixed_to_floating(num4)); EXPECT_EQ(3.141593E258, cudf::convert_fixed_to_floating(num5)); - EXPECT_EQ(DBL_MAX, cudf::convert_fixed_to_floating(num6)); + EXPECT_EQ(3.141593E307, cudf::convert_fixed_to_floating(num6)); } // Negative Exponent @@ -243,7 +243,7 @@ TEST_F(FixedPointTest, PreciseDoubleDecimal64Construction) auto num3 = cudf::convert_floating_to_fixed(3.141593E-158, scale_type(-175)); auto num4 = cudf::convert_floating_to_fixed(3.141593E-208, scale_type(-225)); auto num5 = cudf::convert_floating_to_fixed(3.141593E-258, scale_type(-275)); - auto num6 = cudf::convert_floating_to_fixed(DBL_MIN, scale_type(-325)); + auto num6 = cudf::convert_floating_to_fixed(3.141593E-308, scale_type(-325)); EXPECT_EQ(3.141593E-8, cudf::convert_fixed_to_floating(num0)); EXPECT_EQ(3.141593E-58, cudf::convert_fixed_to_floating(num1)); @@ -251,13 +251,14 @@ TEST_F(FixedPointTest, PreciseDoubleDecimal64Construction) EXPECT_EQ(3.141593E-158, cudf::convert_fixed_to_floating(num3)); EXPECT_EQ(3.141593E-208, cudf::convert_fixed_to_floating(num4)); EXPECT_EQ(3.141593E-258, cudf::convert_fixed_to_floating(num5)); - EXPECT_EQ(DBL_MIN, cudf::convert_fixed_to_floating(num6)); + EXPECT_EQ(3.141593E-308, cudf::convert_fixed_to_floating(num6)); // Denormals auto num7 = cudf::convert_floating_to_fixed(3.141593E-309, scale_type(-326)); auto num8 = cudf::convert_floating_to_fixed(3.141593E-314, scale_type(-331)); auto num9 = cudf::convert_floating_to_fixed(3.141593E-319, scale_type(-336)); auto num10 = cudf::convert_floating_to_fixed(DBL_TRUE_MIN, scale_type(-341)); + EXPECT_EQ(3.141593E-309, cudf::convert_fixed_to_floating(num7)); EXPECT_EQ(3.141593E-314, cudf::convert_fixed_to_floating(num8)); EXPECT_EQ(3.141593E-319, cudf::convert_fixed_to_floating(num9)); diff --git a/python/cudf/cudf/tests/test_decimal.py b/python/cudf/cudf/tests/test_decimal.py index f76e936b3162..7d294e2fb70d 100644 --- a/python/cudf/cudf/tests/test_decimal.py +++ b/python/cudf/cudf/tests/test_decimal.py @@ -6,6 +6,7 @@ import numpy as np import pyarrow as pa import pytest +from packaging import version import cudf from cudf.core.column import Decimal32Column, Decimal64Column, NumericalColumn @@ -81,6 +82,7 @@ def test_from_arrow_max_precision_decimal32(): 94.31304, -112.2314, 0.3333333, + 10.03, np.nan, ] ), @@ -92,6 +94,14 @@ def test_from_arrow_max_precision_decimal32(): [Decimal64Dtype(7, 2), Decimal64Dtype(11, 4), Decimal64Dtype(18, 9)], ) def test_typecast_from_float_to_decimal(request, data, from_dtype, to_dtype): + request.applymarker( + pytest.mark.xfail( + condition=version.parse(pa.__version__) >= version.parse("13.0.0") + and from_dtype == np.dtype("float32") + and to_dtype.scale > 7, + reason="These are bits past the precision of float32. https://github.com/rapidsai/cudf/issues/14169", + ) + ) got = data.astype(from_dtype) pa_arr = got.to_arrow().cast( From 720b838493e97cef5b0753bdb3ff062805431b3e Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Mon, 24 Jun 2024 15:38:17 -0400 Subject: [PATCH 04/18] Update comments --- .../cudf/fixed_point/floating_conversion.hpp | 73 ++++--------------- 1 file changed, 15 insertions(+), 58 deletions(-) diff --git a/cpp/include/cudf/fixed_point/floating_conversion.hpp b/cpp/include/cudf/fixed_point/floating_conversion.hpp index ffcac44ac718..797b3aef0a75 100644 --- a/cpp/include/cudf/fixed_point/floating_conversion.hpp +++ b/cpp/include/cudf/fixed_point/floating_conversion.hpp @@ -798,9 +798,7 @@ CUDF_HOST_DEVICE T increment_on_truncation(T const integral_mantissa, } /** - * @brief Perform lossless base-2 -> base-10 fixed-point conversion for exp10 > 0 - * - * @note Info is lost if the chosen scale factor truncates information. + * @brief Perform base-2 -> base-10 fixed-point conversion for exp10 > 0 * * @tparam FloatingType The type of the original floating-point value we are converting from * @param base2_value The base-2 fixed-point value we are converting from @@ -819,8 +817,8 @@ shift_to_decimal_posexp(typename shifting_constants::IntegerRep co // Output type is ShiftingRep // Here exp10 > 0 and exp2 > 0, so we need to shift left by 2s and divide by 10s. - // To do this losslessly, we will iterate back and forth between them, shifting - // up by 2s and down by 10s until all of the powers have been applied. + // We'll iterate back and forth between them, shifting up by 2s + // and down by 10s until all of the powers have been applied. // However the input base2_value type has virtually no spare room to shift our data // without over- or under-flowing and losing precision. @@ -850,12 +848,6 @@ shift_to_decimal_posexp(typename shifting_constants::IntegerRep co // Iterate, dividing by 10s and shifting up by 2s until we're almost done while (exp10 > Constants::max_digits_shift) { // More decimal places to shift than we have room: Divide the max number of 10s - - // Note that the result of this division is guaranteed to fit within the low half of the bits. - // The highest set bit is num_rep_bits + num_2s_shift_buffer_bits + max_bits_shift - // For float this is 30 + 1 + 30 = 61, for double 57 + 4 + 60 = 121 - // 2^61 / 10^9 (~2^30) is ~2^31, and 2^121 / 10^18 (~2^60) is ~2^61 - // As a future optimization, we could use a faster division routine that takes this account. shifting_rep /= Constants::max_digits_shift_pow; exp10 -= Constants::max_digits_shift; @@ -886,9 +878,7 @@ shift_to_decimal_posexp(typename shifting_constants::IntegerRep co } /** - * @brief Perform lossless base-2 -> base-10 fixed-point conversion for exp10 < 0 - * - * @note Info is lost if the chosen scale factor truncates information. + * @brief Perform base-2 -> base-10 fixed-point conversion for exp10 < 0 * * @tparam FloatingType The type of the original floating-point value we are converting from * @param base2_value The base-2 fixed-point value we are converting from @@ -941,8 +931,7 @@ shift_to_decimal_negexp(typename shifting_constants::IntegerRep ba static constexpr int shift_from = Constants::num_significand_bits; static constexpr int num_init_bit_shift = shift_up_to - shift_from; - // Constants::num_2s_shift_buffer_bits; Note: This shift is safe to do in the smaller IntegerRep - // as it is up to bit 61 / 31 + // Perform initial shift shifting_rep <<= num_init_bit_shift; exp2_mag += num_init_bit_shift; @@ -975,9 +964,7 @@ shift_to_decimal_negexp(typename shifting_constants::IntegerRep ba } /** - * @brief Perform lossless floating-point -> integer decimal conversion - * - * @note Info is lost if the chosen scale factor truncates information. + * @brief Perform floating-point -> integer decimal conversion * * @tparam Rep The type of integer we are converting to, to store the decimal value * @tparam FloatingType The type of floating-point object we are converting from @@ -1051,15 +1038,13 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo } /** - * @brief Perform (nearly) lossless base-10 -> base-2 fixed-point conversion for exp10 > 0 - * - * @note Intended to only be called internally. + * @brief Perform base-10 -> base-2 fixed-point conversion for exp10 > 0 * * @tparam DecimalRep The decimal integer type we are converting from * @tparam FloatingType The type of floating point object we are converting to * @param decimal_rep The decimal integer to convert - * @param exp10 The number of powers of 10 to apply to undo the scale factor. - * @return A pair of the base-2 value and the remaining powers of 2 to be applied. + * @param exp10 The number of powers of 10 to apply to undo the scale factor + * @return A pair of the base-2 value and the remaining powers of 2 to be applied */ template base-2 fixed-point conversion for exp10 < 0 - * - * @note Intended to only be called internally. - * @note A 1-ulp loss may occur, but only for magnitudes E-270 or smaller. + * @brief Perform base-10 -> base-2 fixed-point conversion for exp10 < 0 * * @tparam DecimalRep The decimal integer type we are converting from * @tparam FloatingType The type of floating point object we are converting to * @param decimal_rep The decimal integer to convert - * @param exp10 The number of powers of 10 to apply to undo the scale factor. - * @return A pair of the base-2 value and the remaining powers of 2 to be applied. + * @param exp10 The number of powers of 10 to apply to undo the scale factor + * @return A pair of the base-2 value and the remaining powers of 2 to be applied */ template ; using ShiftingRep = typename Constants::ShiftingRep; - // We would start by lining up our data to num_rep_bits, but if we originated with a floating - // number, we had to keep track of extra bits on the low-side because we were bit-shifting down. - // - // Those bits were not rounded. If we didn't truncate those bits before, we don't want to now - // either to ensure that we end up at the same floating point value that we started with. - // - // We could try to round here instead, but we don't know if we came from a floating-point value - // or not, so any rounding may not be desired. - // - // Note that here we are bit-shifting up, so we also need num_2s_shift_buffer_bits - // room on the high side. We have barely enough room for this for floats, but we're one bit - // over for doubles. So for doubles we'll keep one less bit on the low-side. - // - // This MAY cause a discrepancy in the last bit of our double from the value that we started with. - // However we only need all 4 bits for extremely large exponents - // (one bit to start + one extra bit every 90 powers of 10, so < E-270). - // And it's only a partial bit, and the eventual cast to double rounds, so we - // are often (always?) fine anyway (e.g. DBL_MIN & DBL_TRUE_MIN work fine). - // - // See comments on these constants for more details. - // We want to start with our significand bits at the top of the shifting range, - // so that we don't lose information we need on intermediary right-shifts. - // Note that since we're shifting 2s up, we need num_2s_shift_buffer_bits space on the high side, + // so that we lose minimal information we need on intermediary right-shifts. + // Note that since we're shifting 2s up, we need num_2s_shift_buffer_bits space on the high side static constexpr int shift_up_to = sizeof(ShiftingRep) * 8 - Constants::num_2s_shift_buffer_bits; int const shift_from = count_significant_bits(decimal_rep); int const num_init_bit_shift = shift_up_to - shift_from; @@ -1185,8 +1146,6 @@ CUDF_HOST_DEVICE inline auto shift_to_binary_negexp(DecimalRep decimal_rep, int // Iterate, dividing by 10s and shifting up by 2s until we're almost done while (exp10_mag > Constants::max_digits_shift) { // More decimal places to shift than we have room: Divide the max number of 10s - // Note that the result of this division is guaranteed to fit within low 64/32 bits - // See discussion in shift_to_decimal_posexp() for more details shifting_rep /= Constants::max_digits_shift_pow; exp10_mag -= Constants::max_digits_shift; @@ -1209,9 +1168,7 @@ CUDF_HOST_DEVICE inline auto shift_to_binary_negexp(DecimalRep decimal_rep, int } /** - * @brief Perform (nearly) lossless integer decimal -> floating-point conversion - * - * @note A 1 ulp loss may occur, but only to doubles with magnitude <= 1E-270 + * @brief Perform integer decimal -> floating-point conversion * * @tparam FloatingType The type of floating-point object we are converting to * @tparam Rep The decimal integer type we are converting from From b9eb32d33d859beb44306184b50fe3ce5bda2eed Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Tue, 25 Jun 2024 13:27:21 -0400 Subject: [PATCH 05/18] Change ++ on trunc to add half-bit on trunc. --- .../cudf/fixed_point/floating_conversion.hpp | 79 ++++++++++++------- 1 file changed, 49 insertions(+), 30 deletions(-) diff --git a/cpp/include/cudf/fixed_point/floating_conversion.hpp b/cpp/include/cudf/fixed_point/floating_conversion.hpp index 797b3aef0a75..b072ad997683 100644 --- a/cpp/include/cudf/fixed_point/floating_conversion.hpp +++ b/cpp/include/cudf/fixed_point/floating_conversion.hpp @@ -657,7 +657,7 @@ template ::digits - 1; + static constexpr int max_safe_bit_shift = cuda::std::numeric_limits::digits - 1; return (bit_shift <= max_safe_bit_shift) ? value << bit_shift : cuda::std::numeric_limits::max(); } @@ -674,7 +674,7 @@ template ::digits - 1; + static constexpr int max_safe_bit_shift = cuda::std::numeric_limits::digits - 1; return (bit_shift <= max_safe_bit_shift) ? value >> bit_shift : 0; } @@ -748,20 +748,18 @@ struct shifting_constants { }; /** - * @brief Increment integer rep of floating point if conversion causes truncation + * @brief Add half a bit to integer rep of floating point if conversion causes truncation * * @note This fixes problems like 1.2 (value = 1.1999...) at scale -1 -> 11 * * @tparam T Type of integer holding the floating-point significand - * @param integral_mantissa The integer representation of the floating-point significand + * @param integer_rep The integer representation of the floating-point significand * @param exp2 The power of 2 that needs to be applied to the significand * @param exp10 The power of 10 that needs to be applied to the significand - * @return significand, incremented if the conversion to decimal causes truncation + * @return integer_rep, shifted 1 and ++'d if the conversion to decimal causes truncation */ template )> -CUDF_HOST_DEVICE T increment_on_truncation(T const integral_mantissa, - int const exp2, - int const exp10) +CUDF_HOST_DEVICE cuda::std::pair add_half_if_truncates(T integer_rep, int exp2, int exp10) { // The user-supplied scale may truncate information, so we need to talk about rounding. // We have chosen not to round, so we want 1.23456f with scale -4 to be decimal 12345 @@ -773,14 +771,24 @@ CUDF_HOST_DEVICE T increment_on_truncation(T const integral_mantissa, // and the value 1.199999... happened to be closer to 1.2 than the next value (1.2000...1...) // If the scale truncates information (we didn't choose to keep exactly 1.1999...), how - // do we make sure we store 1.2? All we have to do is add 1 ulp! (unit in the last place) + // do we make sure we store 1.2? We'll add half an ulp! (unit in the last place) // Then 1.1999... becomes 1.2000...1... which truncates to 1.2. - // And if it had been 1.2000...1..., adding 1 ulp still truncates to 1.2, the result is unchanged. + // And if it had been 1.2000...1..., adding half an ulp still truncates to 1.2 - // If we add 1 to the last bit, we are effectively adding 1/2 to the 2nd-to-last bit. - // This is like rounding to the 2nd-to-last bit (if we were to floor afterwards). - // However we don't want to add 1 to the last bit if we are keeping enough precision. - // So, only add 1 if exp10 shift is larger than corresponding exp2 shift to 2nd-to-last bit. + // Why 1/2 an ulp? Because that's all that is needed. The reason we have this problem in the + // first place is because the compiler rounded (e.g.) 1.2 to the nearest floating point number. + // The distance of this rounding is at most 1/2 ulp, otherwise we'd have rounded the other way. + + // How do we add 1/2 an ulp? Just shift the bits left (updating exp2) and add 1. + // We'll always shift up so every input to the conversion algorithm is aligned the same way. + + // If we add a full ulp we run into issues where we add too much and get the wrong result. + // This is because (e.g.) 2^23 = 8.4E6 which is not quite 7 digits of precision. + // So if we want 7 digits, that may "barely" truncate information; adding a 1 ulp is overkill. + + // So when does the user-supplied scale truncate info? + // For powers > 0: When the 10s (scale) shift is larger than the corresponding bit-shift. + // For powers < 0: When the 10s shift is less than the corresponding bit-shift. // Corresponding bit-shift: // 2^10 is approximately 10^3, but this is off by 1.024% @@ -788,13 +796,17 @@ CUDF_HOST_DEVICE T increment_on_truncation(T const integral_mantissa, // So 10^N = 2^(10*N/3 - N/90) = 2^(299*N/90) // Do comparison without dividing, which loses information: // Note: if shift is "equal," still truncates if exp2 < 0 (shifting UP by 2s, 2^10 > 10^3) - int const exp2_term = 90 * (exp2 + 1); //+1: effectively adding 1/2 to 2nd-to-last-bit + int const exp2_term = 90 * exp2; int const exp10_term = 299 * exp10; bool const conversion_truncates = (exp10_term > exp2_term) || ((exp2_term == exp10_term) && (exp2 < 0)); - // (Potentially) increment and return - return integral_mantissa + static_cast(conversion_truncates); + // Add half a bit on truncation (shift to make room and update exp2) + integer_rep <<= 1; + --exp2; + integer_rep += static_cast(conversion_truncates); + + return {integer_rep, exp2}; } /** @@ -831,8 +843,9 @@ shift_to_decimal_posexp(typename shifting_constants::IntegerRep co // so that we don't lose information we need on intermediary right-shifts. // Note that since we're shifting 2s up, we need num_2s_shift_buffer_bits space on the high side, // For all numbers this bit shift is a fixed distance, due to the understood 2^0 bit. + // Note that shift_from is +1 due to shift in add_half_if_truncates() static constexpr int shift_up_to = sizeof(ShiftingRep) * 8 - Constants::num_2s_shift_buffer_bits; - static constexpr int shift_from = Constants::num_significand_bits; + static constexpr int shift_from = Constants::num_significand_bits + 1; static constexpr int max_init_shift = shift_up_to - shift_from; // If our total bit shift is less than this, we don't need to iterate @@ -927,8 +940,9 @@ shift_to_decimal_negexp(typename shifting_constants::IntegerRep ba // We want to start by lining up our bits to the top of the shifting range, // except our first operation is a multiply, so not quite that far // We are bit-shifting down, so we need extra bits on the low-side, which this has. + // Note that shift_from is +1 due to shift in add_half_if_truncates() static constexpr int shift_up_to = sizeof(ShiftingRep) * 8 - Constants::max_bits_shift; - static constexpr int shift_from = Constants::num_significand_bits; + static constexpr int shift_from = Constants::num_significand_bits + 1; static constexpr int num_init_bit_shift = shift_up_to - shift_from; // Perform initial shift @@ -985,16 +999,21 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo // Note that the base2_value here is an unsigned integer with sizeof(FloatingType) auto const is_negative = converter::get_is_negative(integer_rep); - auto const [base2_value, floating_exp2] = converter::get_significand_and_exp2(integer_rep); + auto const [significand, floating_exp2] = converter::get_significand_and_exp2(integer_rep); - auto const exp2 = floating_exp2; // Can't capture from a parameter pack + auto exp2 = floating_exp2; // Can't capture from a parameter pack auto const exp10 = static_cast(scale); - // Increment if truncating to yield expected value, see function for discussion - auto const incremented = increment_on_truncation(base2_value, exp2, exp10); + // Add half a bit if truncating to yield expected value, see function for discussion. + auto const [base2_value_bound, exp2_bound] = + add_half_if_truncates(significand, floating_exp2, exp10); + + // Structured binding variables cannot be captured :/ + auto const base2_value = base2_value_bound; + auto const exp2 = exp2_bound; // Apply the powers of 2 and 10 to convert to decimal. - // The result will be incremented * (2^exp2) / (10^exp10) + // The result will be base2_value * (2^exp2) / (10^exp10) // // Note that while this code is branchy, the decimal scale factor is part of the // column type itself, so every thread will take the same branches on exp10. @@ -1008,26 +1027,26 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo // NOTE: Left Bit-shift can overflow! As can cast! (e.g. double -> decimal32) // Bit shifts may be large, guard against UB if (exp2 >= 0) { - return guarded_left_shift(static_cast(incremented), exp2); + return guarded_left_shift(static_cast(base2_value), exp2); } else { - return guarded_right_shift(incremented, -exp2); + return guarded_right_shift(base2_value, -exp2); } } else if (exp10 > 0) { if (exp2 <= 0) { // Power-2/10 shifts both downward: order doesn't matter, apply and bail. // Guard against shift being undefined behavior - auto const shifted = guarded_right_shift(incremented, -exp2); + auto const shifted = guarded_right_shift(base2_value, -exp2); return divide_power10(shifted, exp10); } - return shift_to_decimal_posexp(incremented, exp2, exp10); + return shift_to_decimal_posexp(base2_value, exp2, exp10); } else { // exp10 < 0 if (exp2 >= 0) { // Power-2/10 shifts both upward: order doesn't matter, apply and bail. // NOTE: Either shift, multiply, or cast (e.g. double -> decimal32) can overflow! - auto const shifted = guarded_left_shift(static_cast(incremented), exp2); + auto const shifted = guarded_left_shift(static_cast(base2_value), exp2); return multiply_power10(shifted, -exp10); } - return shift_to_decimal_negexp(incremented, exp2, exp10); + return shift_to_decimal_negexp(base2_value, exp2, exp10); } }(); From 027140f6608aa10bae879df5c7a0049dcca4aad9 Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Tue, 25 Jun 2024 13:29:10 -0400 Subject: [PATCH 06/18] Fix compile error --- cpp/include/cudf/fixed_point/floating_conversion.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/cpp/include/cudf/fixed_point/floating_conversion.hpp b/cpp/include/cudf/fixed_point/floating_conversion.hpp index b072ad997683..47eea2c1f315 100644 --- a/cpp/include/cudf/fixed_point/floating_conversion.hpp +++ b/cpp/include/cudf/fixed_point/floating_conversion.hpp @@ -1000,9 +1000,7 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo // Note that the base2_value here is an unsigned integer with sizeof(FloatingType) auto const is_negative = converter::get_is_negative(integer_rep); auto const [significand, floating_exp2] = converter::get_significand_and_exp2(integer_rep); - - auto exp2 = floating_exp2; // Can't capture from a parameter pack - auto const exp10 = static_cast(scale); + auto const exp10 = static_cast(scale); // Add half a bit if truncating to yield expected value, see function for discussion. auto const [base2_value_bound, exp2_bound] = From 6a47a0a4c900a611591b3c2e0a85175f5606a97f Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Thu, 27 Jun 2024 14:47:34 -0400 Subject: [PATCH 07/18] Fix java tests. Long.MAX_VALUE gets converted to a double of 2^64 on input so the test was bad. --- java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java b/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java index 1d6a3b3304a0..8f1499700166 100644 --- a/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java +++ b/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java @@ -3509,9 +3509,9 @@ void testCastFloatToDecimal() { @Test void testCastDoubleToDecimal() { testCastNumericToDecimalsAndBack(DType.FLOAT64, false, 0, - () -> ColumnVector.fromBoxedDoubles(1.0, 2.1, -3.23, null, 2.41281, (double) Long.MAX_VALUE), - () -> ColumnVector.fromBoxedDoubles(1.0, 2.0, -3.0, null, 2.0, (double) Long.MAX_VALUE), - new Long[]{1L, 2L, -3L, null, 2L, Long.MAX_VALUE} + () -> ColumnVector.fromBoxedDoubles(1.0, 2.1, -3.23, null, 2.41281, (double) Integer.MAX_VALUE), + () -> ColumnVector.fromBoxedDoubles(1.0, 2.0, -3.0, null, 2.0, (double) Integer.MAX_VALUE), + new Long[]{1L, 2L, -3L, null, 2L, Integer.MAX_VALUE} ); testCastNumericToDecimalsAndBack(DType.FLOAT64, false, -2, () -> ColumnVector.fromBoxedDoubles(1.0, 2.1, -3.23, null, 2.41281, -55.01999), From 28457b382072a03b8f72885832eff84754ce9dbd Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Thu, 27 Jun 2024 14:51:43 -0400 Subject: [PATCH 08/18] Update test comment --- python/cudf/cudf/tests/test_decimal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/cudf/cudf/tests/test_decimal.py b/python/cudf/cudf/tests/test_decimal.py index 7d294e2fb70d..b1bf0a042bc7 100644 --- a/python/cudf/cudf/tests/test_decimal.py +++ b/python/cudf/cudf/tests/test_decimal.py @@ -99,7 +99,7 @@ def test_typecast_from_float_to_decimal(request, data, from_dtype, to_dtype): condition=version.parse(pa.__version__) >= version.parse("13.0.0") and from_dtype == np.dtype("float32") and to_dtype.scale > 7, - reason="These are bits past the precision of float32. https://github.com/rapidsai/cudf/issues/14169", + reason="These fail on bits well past the precision of float32. https://github.com/rapidsai/cudf/issues/14169", ) ) got = data.astype(from_dtype) From ac9587db00ac115fdc339cab6a41484ab6941f06 Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Thu, 27 Jun 2024 16:22:39 -0400 Subject: [PATCH 09/18] Fix conversion from java int to long --- .../java/ai/rapids/cudf/ColumnVectorTest.java | 38 ++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java b/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java index 8f1499700166..1b5c6c757e24 100644 --- a/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java +++ b/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java @@ -20,6 +20,42 @@ import ai.rapids.cudf.ColumnView.FindOptions; import ai.rapids.cudf.HostColumnVector.*; +import src.main.java.ai.rapids.cudf.BaseDeviceMemoryBuffer; +import src.main.java.ai.rapids.cudf.BinaryOp; +import src.main.java.ai.rapids.cudf.BufferType; +import src.main.java.ai.rapids.cudf.CaptureGroups; +import src.main.java.ai.rapids.cudf.ColumnVector; +import src.main.java.ai.rapids.cudf.ColumnView; +import src.main.java.ai.rapids.cudf.ContiguousTable; +import src.main.java.ai.rapids.cudf.Cuda; +import src.main.java.ai.rapids.cudf.CudfException; +import src.main.java.ai.rapids.cudf.DType; +import src.main.java.ai.rapids.cudf.DefaultHostMemoryAllocator; +import src.main.java.ai.rapids.cudf.DeviceMemoryBuffer; +import src.main.java.ai.rapids.cudf.GetJsonObjectOptions; +import src.main.java.ai.rapids.cudf.HostColumnVector; +import src.main.java.ai.rapids.cudf.HostColumnVector.BasicType; +import src.main.java.ai.rapids.cudf.HostColumnVector.DataType; +import src.main.java.ai.rapids.cudf.HostColumnVector.ListType; +import src.main.java.ai.rapids.cudf.HostColumnVector.StructData; +import src.main.java.ai.rapids.cudf.HostColumnVector.StructType; +import src.main.java.ai.rapids.cudf.HostColumnVectorCore; +import src.main.java.ai.rapids.cudf.HostMemoryAllocator; +import src.main.java.ai.rapids.cudf.HostMemoryBuffer; +import src.main.java.ai.rapids.cudf.MemoryCleaner; +import src.main.java.ai.rapids.cudf.NullPolicy; +import src.main.java.ai.rapids.cudf.PadSide; +import src.main.java.ai.rapids.cudf.QuantileMethod; +import src.main.java.ai.rapids.cudf.RegexProgram; +import src.main.java.ai.rapids.cudf.ReplacePolicy; +import src.main.java.ai.rapids.cudf.RollingAggregation; +import src.main.java.ai.rapids.cudf.RoundMode; +import src.main.java.ai.rapids.cudf.Scalar; +import src.main.java.ai.rapids.cudf.ScanAggregation; +import src.main.java.ai.rapids.cudf.ScanType; +import src.main.java.ai.rapids.cudf.Table; +import src.main.java.ai.rapids.cudf.WindowOptions; + import com.google.common.collect.Lists; import org.junit.jupiter.api.Disabled; import org.junit.jupiter.api.Test; @@ -3511,7 +3547,7 @@ void testCastDoubleToDecimal() { testCastNumericToDecimalsAndBack(DType.FLOAT64, false, 0, () -> ColumnVector.fromBoxedDoubles(1.0, 2.1, -3.23, null, 2.41281, (double) Integer.MAX_VALUE), () -> ColumnVector.fromBoxedDoubles(1.0, 2.0, -3.0, null, 2.0, (double) Integer.MAX_VALUE), - new Long[]{1L, 2L, -3L, null, 2L, Integer.MAX_VALUE} + new Long[]{1L, 2L, -3L, null, 2L, Integer.MAX_VALUE.longValue()} ); testCastNumericToDecimalsAndBack(DType.FLOAT64, false, -2, () -> ColumnVector.fromBoxedDoubles(1.0, 2.1, -3.23, null, 2.41281, -55.01999), From ac042d31da80e2bc8ac1cc1a23de70d038d2b645 Mon Sep 17 00:00:00 2001 From: Paul Mattione <156858817+pmattione-nvidia@users.noreply.github.com> Date: Thu, 27 Jun 2024 16:35:48 -0400 Subject: [PATCH 10/18] Update java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java Co-authored-by: Nghia Truong <7416935+ttnghia@users.noreply.github.com> --- .../java/ai/rapids/cudf/ColumnVectorTest.java | 36 ------------------- 1 file changed, 36 deletions(-) diff --git a/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java b/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java index 1b5c6c757e24..ccdef2039f63 100644 --- a/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java +++ b/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java @@ -20,42 +20,6 @@ import ai.rapids.cudf.ColumnView.FindOptions; import ai.rapids.cudf.HostColumnVector.*; -import src.main.java.ai.rapids.cudf.BaseDeviceMemoryBuffer; -import src.main.java.ai.rapids.cudf.BinaryOp; -import src.main.java.ai.rapids.cudf.BufferType; -import src.main.java.ai.rapids.cudf.CaptureGroups; -import src.main.java.ai.rapids.cudf.ColumnVector; -import src.main.java.ai.rapids.cudf.ColumnView; -import src.main.java.ai.rapids.cudf.ContiguousTable; -import src.main.java.ai.rapids.cudf.Cuda; -import src.main.java.ai.rapids.cudf.CudfException; -import src.main.java.ai.rapids.cudf.DType; -import src.main.java.ai.rapids.cudf.DefaultHostMemoryAllocator; -import src.main.java.ai.rapids.cudf.DeviceMemoryBuffer; -import src.main.java.ai.rapids.cudf.GetJsonObjectOptions; -import src.main.java.ai.rapids.cudf.HostColumnVector; -import src.main.java.ai.rapids.cudf.HostColumnVector.BasicType; -import src.main.java.ai.rapids.cudf.HostColumnVector.DataType; -import src.main.java.ai.rapids.cudf.HostColumnVector.ListType; -import src.main.java.ai.rapids.cudf.HostColumnVector.StructData; -import src.main.java.ai.rapids.cudf.HostColumnVector.StructType; -import src.main.java.ai.rapids.cudf.HostColumnVectorCore; -import src.main.java.ai.rapids.cudf.HostMemoryAllocator; -import src.main.java.ai.rapids.cudf.HostMemoryBuffer; -import src.main.java.ai.rapids.cudf.MemoryCleaner; -import src.main.java.ai.rapids.cudf.NullPolicy; -import src.main.java.ai.rapids.cudf.PadSide; -import src.main.java.ai.rapids.cudf.QuantileMethod; -import src.main.java.ai.rapids.cudf.RegexProgram; -import src.main.java.ai.rapids.cudf.ReplacePolicy; -import src.main.java.ai.rapids.cudf.RollingAggregation; -import src.main.java.ai.rapids.cudf.RoundMode; -import src.main.java.ai.rapids.cudf.Scalar; -import src.main.java.ai.rapids.cudf.ScanAggregation; -import src.main.java.ai.rapids.cudf.ScanType; -import src.main.java.ai.rapids.cudf.Table; -import src.main.java.ai.rapids.cudf.WindowOptions; - import com.google.common.collect.Lists; import org.junit.jupiter.api.Disabled; import org.junit.jupiter.api.Test; From c7441eb60e6009f8c4dc9722a73cc6e04958179f Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Thu, 27 Jun 2024 18:58:08 -0400 Subject: [PATCH 11/18] Actually fix conversion from java int to long --- java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java b/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java index ccdef2039f63..7136b162c13a 100644 --- a/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java +++ b/java/src/test/java/ai/rapids/cudf/ColumnVectorTest.java @@ -3511,7 +3511,7 @@ void testCastDoubleToDecimal() { testCastNumericToDecimalsAndBack(DType.FLOAT64, false, 0, () -> ColumnVector.fromBoxedDoubles(1.0, 2.1, -3.23, null, 2.41281, (double) Integer.MAX_VALUE), () -> ColumnVector.fromBoxedDoubles(1.0, 2.0, -3.0, null, 2.0, (double) Integer.MAX_VALUE), - new Long[]{1L, 2L, -3L, null, 2L, Integer.MAX_VALUE.longValue()} + new Long[]{1L, 2L, -3L, null, 2L, (long) Integer.MAX_VALUE} ); testCastNumericToDecimalsAndBack(DType.FLOAT64, false, -2, () -> ColumnVector.fromBoxedDoubles(1.0, 2.1, -3.23, null, 2.41281, -55.01999), From 3e4265fe6858fc7e126719cfabcbacfada3d18b6 Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Mon, 1 Jul 2024 12:40:44 -0400 Subject: [PATCH 12/18] Last python tests pass, matching pandas --- .../cudf/fixed_point/floating_conversion.hpp | 24 +++++++++++++------ python/cudf/cudf/tests/test_decimal.py | 8 ------- 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/cpp/include/cudf/fixed_point/floating_conversion.hpp b/cpp/include/cudf/fixed_point/floating_conversion.hpp index 47eea2c1f315..07a176d3b2cf 100644 --- a/cpp/include/cudf/fixed_point/floating_conversion.hpp +++ b/cpp/include/cudf/fixed_point/floating_conversion.hpp @@ -19,6 +19,7 @@ #include #include +#include #include #include @@ -759,7 +760,9 @@ struct shifting_constants { * @return integer_rep, shifted 1 and ++'d if the conversion to decimal causes truncation */ template )> -CUDF_HOST_DEVICE cuda::std::pair add_half_if_truncates(T integer_rep, int exp2, int exp10) +CUDF_HOST_DEVICE cuda::std::tuple add_half_if_truncates(T integer_rep, + int exp2, + int exp10) { // The user-supplied scale may truncate information, so we need to talk about rounding. // We have chosen not to round, so we want 1.23456f with scale -4 to be decimal 12345 @@ -806,7 +809,7 @@ CUDF_HOST_DEVICE cuda::std::pair add_half_if_truncates(T integer_rep, in --exp2; integer_rep += static_cast(conversion_truncates); - return {integer_rep, exp2}; + return {integer_rep, exp2, conversion_truncates}; } /** @@ -897,6 +900,7 @@ shift_to_decimal_posexp(typename shifting_constants::IntegerRep co * @param base2_value The base-2 fixed-point value we are converting from * @param exp2 The number of powers of 2 to apply to convert from base-2 * @param exp10 The number of powers of 10 to apply to reach the desired scale factor + * @param truncating Whether the scale factor truncates floating-point information * @return Magnitude of the converted-to decimal integer */ template ::ShiftingRep shift_to_decimal_negexp(typename shifting_constants::IntegerRep base2_value, int exp2, - int exp10) + int exp10, + bool truncating) { // This is similar to shift_to_decimal_posexp(), except exp10 < 0 & exp2 < 0 // See comments in that function for details. @@ -930,8 +935,12 @@ shift_to_decimal_negexp(typename shifting_constants::IntegerRep ba shifting_rep = multiply_power10_32bit(shifting_rep, exp10_mag); } - // Final bit shift: Shift may be large, guard against UB - return guarded_right_shift(shifting_rep, exp2_mag); + // Final bit shifting: Shift may be large, guard against UB + // If we aren't truncating information from the original floating-point number, + // then always round up (to match pandas) by adding 1 to the last bit + shifting_rep = guarded_right_shift(shifting_rep, exp2_mag - 1); + shifting_rep += static_cast(!truncating); + return (shifting_rep >> 1); }; // If our total decimal shift is less than the max, we don't need to iterate @@ -1003,12 +1012,13 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo auto const exp10 = static_cast(scale); // Add half a bit if truncating to yield expected value, see function for discussion. - auto const [base2_value_bound, exp2_bound] = + auto const [base2_value_bound, exp2_bound, truncates_bound] = add_half_if_truncates(significand, floating_exp2, exp10); // Structured binding variables cannot be captured :/ auto const base2_value = base2_value_bound; auto const exp2 = exp2_bound; + auto const truncates = truncates_bound; // Apply the powers of 2 and 10 to convert to decimal. // The result will be base2_value * (2^exp2) / (10^exp10) @@ -1044,7 +1054,7 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo auto const shifted = guarded_left_shift(static_cast(base2_value), exp2); return multiply_power10(shifted, -exp10); } - return shift_to_decimal_negexp(base2_value, exp2, exp10); + return shift_to_decimal_negexp(base2_value, exp2, exp10, truncates); } }(); diff --git a/python/cudf/cudf/tests/test_decimal.py b/python/cudf/cudf/tests/test_decimal.py index a34b7a97e86e..5dc9a6c5ae52 100644 --- a/python/cudf/cudf/tests/test_decimal.py +++ b/python/cudf/cudf/tests/test_decimal.py @@ -94,14 +94,6 @@ def test_from_arrow_max_precision_decimal32(): [Decimal64Dtype(7, 2), Decimal64Dtype(11, 4), Decimal64Dtype(18, 9)], ) def test_typecast_from_float_to_decimal(request, data, from_dtype, to_dtype): - request.applymarker( - pytest.mark.xfail( - condition=version.parse(pa.__version__) >= version.parse("13.0.0") - and from_dtype == np.dtype("float32") - and to_dtype.scale > 7, - reason="These fail on bits well past the precision of float32. https://github.com/rapidsai/cudf/issues/14169", - ) - ) got = data.astype(from_dtype) pa_arr = got.to_arrow().cast( From f7a86638c8b63d2e7be7cc52d01aaba11d599c48 Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Mon, 1 Jul 2024 12:54:54 -0400 Subject: [PATCH 13/18] Remove unused python package --- python/cudf/cudf/tests/test_decimal.py | 1 - 1 file changed, 1 deletion(-) diff --git a/python/cudf/cudf/tests/test_decimal.py b/python/cudf/cudf/tests/test_decimal.py index 5dc9a6c5ae52..f3b2ac25fa73 100644 --- a/python/cudf/cudf/tests/test_decimal.py +++ b/python/cudf/cudf/tests/test_decimal.py @@ -6,7 +6,6 @@ import numpy as np import pyarrow as pa import pytest -from packaging import version import cudf from cudf.core.column import Decimal32Column, Decimal64Column, NumericalColumn From cf531d467d8eb6ab5a4f97eeca869be122a7f5e8 Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Mon, 1 Jul 2024 18:45:29 -0400 Subject: [PATCH 14/18] Address comments --- .../cudf/fixed_point/floating_conversion.hpp | 365 +++++++++--------- 1 file changed, 180 insertions(+), 185 deletions(-) diff --git a/cpp/include/cudf/fixed_point/floating_conversion.hpp b/cpp/include/cudf/fixed_point/floating_conversion.hpp index 07a176d3b2cf..6568461626ab 100644 --- a/cpp/include/cudf/fixed_point/floating_conversion.hpp +++ b/cpp/include/cudf/fixed_point/floating_conversion.hpp @@ -190,7 +190,7 @@ struct floating_converter { * @param integer_rep The bit-casted floating value to extract the exponent from * @return The stored base-2 exponent and significand, shifted for denormals */ - CUDF_HOST_DEVICE inline static std::pair get_significand_and_exp2( + CUDF_HOST_DEVICE inline static std::pair get_significand_and_pow2( IntegralType integer_rep) { // Extract the significand @@ -202,7 +202,7 @@ struct floating_converter { // Notes on special values of exponent_bits: // bits = exponent_mask is +/-inf or NaN, but those are handled prior to input. // bits = 0 is either a denormal (handled below) or a zero (handled earlier by caller). - int floating_exp2; + int floating_pow2; if (exponent_bits == 0) { // Denormal values are 2^(1 - exponent_bias) * Sum_i(B_i * 2^-i) // Where i is the i-th mantissa bit (counting from the LEFT, starting at 1), @@ -210,17 +210,17 @@ struct floating_converter { // So e.g. for the minimum denormal, only the lowest bit is set: // FLT_TRUE_MIN = 2^(1 - 127) * 2^-23 = 2^-149 // DBL_TRUE_MIN = 2^(1 - 1023) * 2^-52 = 2^-1074 - floating_exp2 = 1 - exponent_bias; + floating_pow2 = 1 - exponent_bias; // Line-up denormal to same (understood) bit as normal numbers // This is so bit-shifting starts at the same bit index auto const lineup_shift = num_significand_bits - count_significant_bits(significand); significand <<= lineup_shift; - floating_exp2 -= lineup_shift; + floating_pow2 -= lineup_shift; } else { // Extract the exponent value: shift the bits down and subtract the bias. auto const shifted_exponent_bits = exponent_bits >> num_stored_mantissa_bits; - floating_exp2 = static_cast(shifted_exponent_bits) - exponent_bias; + floating_pow2 = static_cast(shifted_exponent_bits) - exponent_bias; // Set the high bit for the understood 1/2 significand |= understood_bit_mask; @@ -228,9 +228,9 @@ struct floating_converter { // To convert the mantissa to an integer, we effectively applied #-mantissa-bits // powers of 2 to convert the fractional value to an integer, so subtract them off here - int const exp2 = floating_exp2 - num_stored_mantissa_bits; + int const pow2 = floating_pow2 - num_stored_mantissa_bits; - return {significand, exp2}; + return {significand, pow2}; } /** @@ -259,10 +259,10 @@ struct floating_converter { * @note The caller must guarantee that the input is a positive (> 0) whole number. * * @param floating The floating value to add to the exponent of. Must be positive. - * @param exp2 The power-of-2 to add to the floating-point number - * @return The input floating-point value * 2^exp2 + * @param pow2 The power-of-2 to add to the floating-point number + * @return The input floating-point value * 2^pow2 */ - CUDF_HOST_DEVICE inline static FloatingType add_exp2(FloatingType floating, int exp2) + CUDF_HOST_DEVICE inline static FloatingType add_pow2(FloatingType floating, int pow2) { // Note that the input floating-point number is positive (& whole), so we don't have to // worry about the sign here; the sign will be set later in set_is_negative() @@ -273,18 +273,18 @@ struct floating_converter { // Extract the currently stored (biased) exponent using SignedType = std::make_signed_t; auto exponent_bits = integer_rep & exponent_mask; - auto stored_exp2 = static_cast(exponent_bits >> num_stored_mantissa_bits); + auto stored_pow2 = static_cast(exponent_bits >> num_stored_mantissa_bits); // Add the additional power-of-2 - stored_exp2 += exp2; + stored_pow2 += pow2; // Check for exponent over/under-flow. - if (stored_exp2 <= 0) { + if (stored_pow2 <= 0) { // Denormal (zero handled prior to input) // Early out if bit shift will zero it anyway. // Note: We must handle this explicitly, as too-large a bit-shift is UB - auto const bit_shift = -stored_exp2 + 1; //+1 due to understood bit set below + auto const bit_shift = -stored_pow2 + 1; //+1 due to understood bit set below if (bit_shift > num_stored_mantissa_bits) { return 0.0; } // Clear the exponent bits (zero means 2^-126/2^-1022 w/ no understood bit) @@ -296,12 +296,12 @@ struct floating_converter { // Convert to denormal: bit shift off the low bits integer_rep >>= bit_shift; - } else if (stored_exp2 >= static_cast(unshifted_exponent_mask)) { + } else if (stored_pow2 >= static_cast(unshifted_exponent_mask)) { // Overflow: Set infinity return cuda::std::numeric_limits::infinity(); } else { // Normal number: Clear existing exponent bits and set new ones - exponent_bits = static_cast(stored_exp2) << num_stored_mantissa_bits; + exponent_bits = static_cast(stored_pow2) << num_stored_mantissa_bits; integer_rep &= (~exponent_mask); integer_rep |= exponent_bits; } @@ -317,18 +317,18 @@ struct floating_converter { * * @note Intended to be run at compile time. * - * @tparam Exp10 The power of 10 to calculate - * @return Returns 10^Exp10 + * @tparam Pow10 The power of 10 to calculate + * @return Returns 10^Pow10 */ -template +template constexpr __uint128_t large_power_of_10() { // Stop at 10^19 to speed up compilation; literals can be used for smaller powers of 10. - static_assert(Exp10 >= 19); - if constexpr (Exp10 == 19) + static_assert(Pow10 >= 19); + if constexpr (Pow10 == 19) return __uint128_t(10000000000000000000ULL); else - return large_power_of_10() * __uint128_t(10); + return large_power_of_10() * __uint128_t(10); } /** @@ -336,11 +336,11 @@ constexpr __uint128_t large_power_of_10() * * @tparam T Type of value to be divided-from. * @param value The number to be divided-from. - * @param exp10 The power-of-10 of the denominator, from 0 to 9 inclusive. - * @return Returns value / 10^exp10 + * @param pow10 The power-of-10 of the denominator, from 0 to 9 inclusive. + * @return Returns value / 10^pow10 */ template )> -CUDF_HOST_DEVICE inline T divide_power10_32bit(T value, int exp10) +CUDF_HOST_DEVICE inline T divide_power10_32bit(T value, int pow10) { // Computing division this way is much faster than the alternatives. // Division is not implemented in GPU hardware, and the compiler will often implement it as a @@ -350,7 +350,7 @@ CUDF_HOST_DEVICE inline T divide_power10_32bit(T value, int exp10) // Instead, if the compiler can see exactly what number it is dividing by, it can // produce much more optimal assembly, doing bit shifting, multiplies by a constant, etc. - // For the compiler to see the value though, array lookup (with exp10 as the index) + // For the compiler to see the value though, array lookup (with pow10 as the index) // is not sufficient: We have to use a switch statement. Although this introduces a branch, // it is still much faster than doing the divide any other way. // Perhaps an array can be used in C++23 with the assume attribute? @@ -366,7 +366,7 @@ CUDF_HOST_DEVICE inline T divide_power10_32bit(T value, int exp10) // introduces too much pressure on the kernels that use this code, slowing down their benchmarks. // It also dramatically slows down the compile time. - switch (exp10) { + switch (pow10) { case 0: return value; case 1: return value / 10U; case 2: return value / 100U; @@ -386,14 +386,14 @@ CUDF_HOST_DEVICE inline T divide_power10_32bit(T value, int exp10) * * @tparam T Type of value to be divided-from. * @param value The number to be divided-from. - * @param exp10 The power-of-10 of the denominator, from 0 to 19 inclusive. - * @return Returns value / 10^exp10 + * @param pow10 The power-of-10 of the denominator, from 0 to 19 inclusive. + * @return Returns value / 10^pow10 */ template )> -CUDF_HOST_DEVICE inline T divide_power10_64bit(T value, int exp10) +CUDF_HOST_DEVICE inline T divide_power10_64bit(T value, int pow10) { // See comments in divide_power10_32bit() for discussion. - switch (exp10) { + switch (pow10) { case 0: return value; case 1: return value / 10U; case 2: return value / 100U; @@ -423,14 +423,14 @@ CUDF_HOST_DEVICE inline T divide_power10_64bit(T value, int exp10) * * @tparam T Type of value to be divided-from. * @param value The number to be divided-from. - * @param exp10 The power-of-10 of the denominator, from 0 to 38 inclusive. - * @return Returns value / 10^exp10. + * @param pow10 The power-of-10 of the denominator, from 0 to 38 inclusive. + * @return Returns value / 10^pow10. */ template )> -CUDF_HOST_DEVICE inline constexpr T divide_power10_128bit(T value, int exp10) +CUDF_HOST_DEVICE inline constexpr T divide_power10_128bit(T value, int pow10) { // See comments in divide_power10_32bit() for an introduction. - switch (exp10) { + switch (pow10) { case 0: return value; case 1: return value / 10U; case 2: return value / 100U; @@ -479,14 +479,14 @@ CUDF_HOST_DEVICE inline constexpr T divide_power10_128bit(T value, int exp10) * * @tparam T Type of value to be multiplied. * @param value The number to be multiplied. - * @param exp10 The power-of-10 of the multiplier, from 0 to 9 inclusive. - * @return Returns value * 10^exp10 + * @param pow10 The power-of-10 of the multiplier, from 0 to 9 inclusive. + * @return Returns value * 10^pow10 */ template )> -CUDF_HOST_DEVICE inline constexpr T multiply_power10_32bit(T value, int exp10) +CUDF_HOST_DEVICE inline constexpr T multiply_power10_32bit(T value, int pow10) { // See comments in divide_power10_32bit() for discussion. - switch (exp10) { + switch (pow10) { case 0: return value; case 1: return value * 10U; case 2: return value * 100U; @@ -506,14 +506,14 @@ CUDF_HOST_DEVICE inline constexpr T multiply_power10_32bit(T value, int exp10) * * @tparam T Type of value to be multiplied. * @param value The number to be multiplied. - * @param exp10 The power-of-10 of the multiplier, from 0 to 19 inclusive. - * @return Returns value * 10^exp10 + * @param pow10 The power-of-10 of the multiplier, from 0 to 19 inclusive. + * @return Returns value * 10^pow10 */ template )> -CUDF_HOST_DEVICE inline constexpr T multiply_power10_64bit(T value, int exp10) +CUDF_HOST_DEVICE inline constexpr T multiply_power10_64bit(T value, int pow10) { // See comments in divide_power10_32bit() for discussion. - switch (exp10) { + switch (pow10) { case 0: return value; case 1: return value * 10U; case 2: return value * 100U; @@ -543,14 +543,14 @@ CUDF_HOST_DEVICE inline constexpr T multiply_power10_64bit(T value, int exp10) * * @tparam T Type of value to be multiplied. * @param value The number to be multiplied. - * @param exp10 The power-of-10 of the multiplier, from 0 to 38 inclusive. - * @return Returns value * 10^exp10. + * @param pow10 The power-of-10 of the multiplier, from 0 to 38 inclusive. + * @return Returns value * 10^pow10. */ template )> -CUDF_HOST_DEVICE inline constexpr T multiply_power10_128bit(T value, int exp10) +CUDF_HOST_DEVICE inline constexpr T multiply_power10_128bit(T value, int pow10) { // See comments in divide_power10_128bit() for discussion. - switch (exp10) { + switch (pow10) { case 0: return value; case 1: return value * 10U; case 2: return value * 100U; @@ -597,52 +597,52 @@ CUDF_HOST_DEVICE inline constexpr T multiply_power10_128bit(T value, int exp10) /** * @brief Multiply an integer by a power of 10. * - * @note Use this function if you have no a-priori knowledge of what exp10 might be. + * @note Use this function if you have no a-priori knowledge of what pow10 might be. * If you do, prefer calling the bit-size-specific versions * * @tparam Rep Representation type needed for integer exponentiation * @tparam T Integral type of value to be multiplied. * @param value The number to be multiplied. - * @param exp10 The power-of-10 of the multiplier. - * @return Returns value * 10^exp10 + * @param pow10 The power-of-10 of the multiplier. + * @return Returns value * 10^pow10 */ template )> -CUDF_HOST_DEVICE inline constexpr T multiply_power10(T value, int exp10) +CUDF_HOST_DEVICE inline constexpr T multiply_power10(T value, int pow10) { - // Use this function if you have no knowledge of what exp10 might be + // Use this function if you have no knowledge of what pow10 might be // If you do, prefer calling the bit-size-specific versions if constexpr (sizeof(Rep) <= 4) { - return multiply_power10_32bit(value, exp10); + return multiply_power10_32bit(value, pow10); } else if constexpr (sizeof(Rep) <= 8) { - return multiply_power10_64bit(value, exp10); + return multiply_power10_64bit(value, pow10); } else { - return multiply_power10_128bit(value, exp10); + return multiply_power10_128bit(value, pow10); } } /** * @brief Divide an integer by a power of 10. * - * @note Use this function if you have no a-priori knowledge of what exp10 might be. + * @note Use this function if you have no a-priori knowledge of what pow10 might be. * If you do, prefer calling the bit-size-specific versions * * @tparam Rep Representation type needed for integer exponentiation * @tparam T Integral type of value to be divided-from. * @param value The number to be divided-from. - * @param exp10 The power-of-10 of the denominator. - * @return Returns value / 10^exp10 + * @param pow10 The power-of-10 of the denominator. + * @return Returns value / 10^pow10 */ template )> -CUDF_HOST_DEVICE inline constexpr T divide_power10(T value, int exp10) +CUDF_HOST_DEVICE inline constexpr T divide_power10(T value, int pow10) { - // Use this function if you have no knowledge of what exp10 might be + // Use this function if you have no knowledge of what pow10 might be // If you do, prefer calling the bit-size-specific versions if constexpr (sizeof(Rep) <= 4) { - return divide_power10_32bit(value, exp10); + return divide_power10_32bit(value, pow10); } else if constexpr (sizeof(Rep) <= 8) { - return divide_power10_64bit(value, exp10); + return divide_power10_64bit(value, pow10); } else { - return divide_power10_128bit(value, exp10); + return divide_power10_128bit(value, pow10); } } @@ -658,7 +658,7 @@ template ::digits - 1; + constexpr int max_safe_bit_shift = cuda::std::numeric_limits::digits - 1; return (bit_shift <= max_safe_bit_shift) ? value << bit_shift : cuda::std::numeric_limits::max(); } @@ -675,7 +675,7 @@ template ::digits - 1; + constexpr int max_safe_bit_shift = cuda::std::numeric_limits::digits - 1; return (bit_shift <= max_safe_bit_shift) ? value >> bit_shift : 0; } @@ -755,14 +755,14 @@ struct shifting_constants { * * @tparam T Type of integer holding the floating-point significand * @param integer_rep The integer representation of the floating-point significand - * @param exp2 The power of 2 that needs to be applied to the significand - * @param exp10 The power of 10 that needs to be applied to the significand + * @param pow2 The power of 2 that needs to be applied to the significand + * @param pow10 The power of 10 that needs to be applied to the significand * @return integer_rep, shifted 1 and ++'d if the conversion to decimal causes truncation */ template )> CUDF_HOST_DEVICE cuda::std::tuple add_half_if_truncates(T integer_rep, - int exp2, - int exp10) + int pow2, + int pow10) { // The user-supplied scale may truncate information, so we need to talk about rounding. // We have chosen not to round, so we want 1.23456f with scale -4 to be decimal 12345 @@ -782,7 +782,7 @@ CUDF_HOST_DEVICE cuda::std::tuple add_half_if_truncates(T integer_ // first place is because the compiler rounded (e.g.) 1.2 to the nearest floating point number. // The distance of this rounding is at most 1/2 ulp, otherwise we'd have rounded the other way. - // How do we add 1/2 an ulp? Just shift the bits left (updating exp2) and add 1. + // How do we add 1/2 an ulp? Just shift the bits left (updating pow2) and add 1. // We'll always shift up so every input to the conversion algorithm is aligned the same way. // If we add a full ulp we run into issues where we add too much and get the wrong result. @@ -798,40 +798,40 @@ CUDF_HOST_DEVICE cuda::std::tuple add_half_if_truncates(T integer_ // 1.024^30 is 2.03704, so this is high by one bit for every 30*3 = 90 powers of 10 // So 10^N = 2^(10*N/3 - N/90) = 2^(299*N/90) // Do comparison without dividing, which loses information: - // Note: if shift is "equal," still truncates if exp2 < 0 (shifting UP by 2s, 2^10 > 10^3) - int const exp2_term = 90 * exp2; - int const exp10_term = 299 * exp10; + // Note: if shift is "equal," still truncates if pow2 < 0 (shifting UP by 2s, 2^10 > 10^3) + int const pow2_term = 90 * pow2; + int const pow10_term = 299 * pow10; bool const conversion_truncates = - (exp10_term > exp2_term) || ((exp2_term == exp10_term) && (exp2 < 0)); + (pow10_term > pow2_term) || ((pow2_term == pow10_term) && (pow2 < 0)); - // Add half a bit on truncation (shift to make room and update exp2) + // Add half a bit on truncation (shift to make room and update pow2) integer_rep <<= 1; - --exp2; + --pow2; integer_rep += static_cast(conversion_truncates); - return {integer_rep, exp2, conversion_truncates}; + return {integer_rep, pow2, conversion_truncates}; } /** - * @brief Perform base-2 -> base-10 fixed-point conversion for exp10 > 0 + * @brief Perform base-2 -> base-10 fixed-point conversion for pow10 > 0 * * @tparam FloatingType The type of the original floating-point value we are converting from * @param base2_value The base-2 fixed-point value we are converting from - * @param exp2 The number of powers of 2 to apply to convert from base-2 - * @param exp10 The number of powers of 10 to apply to reach the desired scale factor + * @param pow2 The number of powers of 2 to apply to convert from base-2 + * @param pow10 The number of powers of 10 to apply to reach the desired scale factor * @return Magnitude of the converted-to decimal integer */ template )> CUDF_HOST_DEVICE inline typename shifting_constants::ShiftingRep -shift_to_decimal_posexp(typename shifting_constants::IntegerRep const base2_value, - int exp2, - int exp10) +shift_to_decimal_pospow(typename shifting_constants::IntegerRep const base2_value, + int pow2, + int pow10) { // To convert to decimal, we need to apply the input powers of 2 and 10 - // The result will be (integer) base2_value * (2^exp2) / (10^exp10) + // The result will be (integer) base2_value * (2^pow2) / (10^pow10) // Output type is ShiftingRep - // Here exp10 > 0 and exp2 > 0, so we need to shift left by 2s and divide by 10s. + // Here pow10 > 0 and pow2 > 0, so we need to shift left by 2s and divide by 10s. // We'll iterate back and forth between them, shifting up by 2s // and down by 10s until all of the powers have been applied. @@ -852,54 +852,54 @@ shift_to_decimal_posexp(typename shifting_constants::IntegerRep co static constexpr int max_init_shift = shift_up_to - shift_from; // If our total bit shift is less than this, we don't need to iterate - if (exp2 <= max_init_shift) { + if (pow2 <= max_init_shift) { // Shift bits left, divide by 10s to apply the scale factor, and we're done. - return divide_power10(shifting_rep << exp2, exp10); + return divide_power10(shifting_rep << pow2, pow10); } // We need to iterate. Do the combined initial shift shifting_rep <<= max_init_shift; - exp2 -= max_init_shift; + pow2 -= max_init_shift; // Iterate, dividing by 10s and shifting up by 2s until we're almost done - while (exp10 > Constants::max_digits_shift) { + while (pow10 > Constants::max_digits_shift) { // More decimal places to shift than we have room: Divide the max number of 10s shifting_rep /= Constants::max_digits_shift_pow; - exp10 -= Constants::max_digits_shift; + pow10 -= Constants::max_digits_shift; // If our remaining bit shift is less than the max, we're finished iterating - if (exp2 <= Constants::max_bits_shift) { + if (pow2 <= Constants::max_bits_shift) { // Shift bits left, divide by 10s to apply the scale factor, and we're done. // Note: This divide result may not fit in the low half of the bit range - return divide_power10(shifting_rep << exp2, exp10); + return divide_power10(shifting_rep << pow2, pow10); } // Shift the max number of bits left again shifting_rep <<= Constants::max_bits_shift; - exp2 -= Constants::max_bits_shift; + pow2 -= Constants::max_bits_shift; } - // Last 10s-shift: Divdie all remaining decimal places, shift all remaining bits, then bail + // Last 10s-shift: Divide all remaining decimal places, shift all remaining bits, then bail // Note: This divide result may not fit in the low half of the bit range // But the divisor is less than the max-shift, and thus fits within 64 / 32 bits if constexpr (Constants::is_double) { - shifting_rep = divide_power10_64bit(shifting_rep, exp10); + shifting_rep = divide_power10_64bit(shifting_rep, pow10); } else { - shifting_rep = divide_power10_32bit(shifting_rep, exp10); + shifting_rep = divide_power10_32bit(shifting_rep, pow10); } // Final bit shift: Shift may be large, guard against UB // NOTE: This can overflow! - return guarded_left_shift(shifting_rep, exp2); + return guarded_left_shift(shifting_rep, pow2); } /** - * @brief Perform base-2 -> base-10 fixed-point conversion for exp10 < 0 + * @brief Perform base-2 -> base-10 fixed-point conversion for pow10 < 0 * * @tparam FloatingType The type of the original floating-point value we are converting from * @param base2_value The base-2 fixed-point value we are converting from - * @param exp2 The number of powers of 2 to apply to convert from base-2 - * @param exp10 The number of powers of 10 to apply to reach the desired scale factor + * @param pow2 The number of powers of 2 to apply to convert from base-2 + * @param pow10 The number of powers of 10 to apply to reach the desired scale factor * @param truncating Whether the scale factor truncates floating-point information * @return Magnitude of the converted-to decimal integer */ @@ -907,12 +907,12 @@ template )> CUDF_HOST_DEVICE inline typename shifting_constants::ShiftingRep -shift_to_decimal_negexp(typename shifting_constants::IntegerRep base2_value, - int exp2, - int exp10, +shift_to_decimal_negpow(typename shifting_constants::IntegerRep base2_value, + int pow2, + int pow10, bool truncating) { - // This is similar to shift_to_decimal_posexp(), except exp10 < 0 & exp2 < 0 + // This is similar to shift_to_decimal_pospow(), except pow10 < 0 & pow2 < 0 // See comments in that function for details. // Instead here we need to multiply by 10s and shift right by 2s @@ -922,29 +922,29 @@ shift_to_decimal_negexp(typename shifting_constants::IntegerRep ba auto shifting_rep = static_cast(base2_value); // Convert to using positive values so we don't have keep negating - int exp10_mag = -exp10; - int exp2_mag = -exp2; + int pow10_mag = -pow10; + int pow2_mag = -pow2; // For performing final 10s-shift auto final_shifts_low10s = [&]() { // Last 10s-shift: multiply all remaining decimal places, shift all remaining bits, then bail // The multiplier is less than the max-shift, and thus fits within 64 / 32 bits if constexpr (Constants::is_double) { - shifting_rep = multiply_power10_64bit(shifting_rep, exp10_mag); + shifting_rep = multiply_power10_64bit(shifting_rep, pow10_mag); } else { - shifting_rep = multiply_power10_32bit(shifting_rep, exp10_mag); + shifting_rep = multiply_power10_32bit(shifting_rep, pow10_mag); } // Final bit shifting: Shift may be large, guard against UB // If we aren't truncating information from the original floating-point number, // then always round up (to match pandas) by adding 1 to the last bit - shifting_rep = guarded_right_shift(shifting_rep, exp2_mag - 1); + shifting_rep = guarded_right_shift(shifting_rep, pow2_mag - 1); shifting_rep += static_cast(!truncating); return (shifting_rep >> 1); }; // If our total decimal shift is less than the max, we don't need to iterate - if (exp10_mag <= Constants::max_digits_shift) { return final_shifts_low10s(); } + if (pow10_mag <= Constants::max_digits_shift) { return final_shifts_low10s(); } // We want to start by lining up our bits to the top of the shifting range, // except our first operation is a multiply, so not quite that far @@ -956,31 +956,31 @@ shift_to_decimal_negexp(typename shifting_constants::IntegerRep ba // Perform initial shift shifting_rep <<= num_init_bit_shift; - exp2_mag += num_init_bit_shift; + pow2_mag += num_init_bit_shift; // Iterate, multiplying by 10s and shifting down by 2s until we're almost done do { // More decimal places to shift than we have room: Multiply the max number of 10s shifting_rep *= Constants::max_digits_shift_pow; - exp10_mag -= Constants::max_digits_shift; + pow10_mag -= Constants::max_digits_shift; // If our remaining bit shift is less than the max, we're finished iterating - if (exp2_mag <= Constants::max_bits_shift) { + if (pow2_mag <= Constants::max_bits_shift) { // Last bit-shift: Shift all remaining bits, apply the remaining scale, then bail - shifting_rep >>= exp2_mag; + shifting_rep >>= pow2_mag; // We need to convert to the output rep for the final scale-factor multiply, because if (e.g.) - // float -> dec128 and some large exp10_mag, it might overflow the 64bit shifting rep. - // It's not needed for exp10 > 0 because we're dividing by 10s there instead of multiplying. + // float -> dec128 and some large pow10_mag, it might overflow the 64bit shifting rep. + // It's not needed for pow10 > 0 because we're dividing by 10s there instead of multiplying. using UnsignedRep = cuda::std::make_unsigned_t; // NOTE: This can overflow! (Both multiply and cast) - return multiply_power10(static_cast(shifting_rep), exp10_mag); + return multiply_power10(static_cast(shifting_rep), pow10_mag); } // More bits to shift than we have room: Shift the max number of 2s shifting_rep >>= Constants::max_bits_shift; - exp2_mag -= Constants::max_bits_shift; - } while (exp10_mag > Constants::max_digits_shift); + pow2_mag -= Constants::max_bits_shift; + } while (pow10_mag > Constants::max_digits_shift); // Do our final shifts return final_shifts_low10s(); @@ -1008,53 +1008,53 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo // Note that the base2_value here is an unsigned integer with sizeof(FloatingType) auto const is_negative = converter::get_is_negative(integer_rep); - auto const [significand, floating_exp2] = converter::get_significand_and_exp2(integer_rep); - auto const exp10 = static_cast(scale); + auto const [significand, floating_pow2] = converter::get_significand_and_pow2(integer_rep); + auto const pow10 = static_cast(scale); // Add half a bit if truncating to yield expected value, see function for discussion. - auto const [base2_value_bound, exp2_bound, truncates_bound] = - add_half_if_truncates(significand, floating_exp2, exp10); + auto const [base2_value_bound, pow2_bound, truncates_bound] = + add_half_if_truncates(significand, floating_pow2, pow10); // Structured binding variables cannot be captured :/ auto const base2_value = base2_value_bound; - auto const exp2 = exp2_bound; + auto const pow2 = pow2_bound; auto const truncates = truncates_bound; // Apply the powers of 2 and 10 to convert to decimal. - // The result will be base2_value * (2^exp2) / (10^exp10) + // The result will be base2_value * (2^pow2) / (10^pow10) // // Note that while this code is branchy, the decimal scale factor is part of the - // column type itself, so every thread will take the same branches on exp10. + // column type itself, so every thread will take the same branches on pow10. // Also data within a column tends to be similar, so they will often take the - // same branches on exp2 as well. + // same branches on pow2 as well. // // NOTE: some returns here can overflow (e.g. ShiftingRep -> UnsignedRep) using UnsignedRep = cuda::std::make_unsigned_t; auto const magnitude = [&]() -> UnsignedRep { - if (exp10 == 0) { + if (pow10 == 0) { // NOTE: Left Bit-shift can overflow! As can cast! (e.g. double -> decimal32) // Bit shifts may be large, guard against UB - if (exp2 >= 0) { - return guarded_left_shift(static_cast(base2_value), exp2); + if (pow2 >= 0) { + return guarded_left_shift(static_cast(base2_value), pow2); } else { - return guarded_right_shift(base2_value, -exp2); + return guarded_right_shift(base2_value, -pow2); } - } else if (exp10 > 0) { - if (exp2 <= 0) { + } else if (pow10 > 0) { + if (pow2 <= 0) { // Power-2/10 shifts both downward: order doesn't matter, apply and bail. // Guard against shift being undefined behavior - auto const shifted = guarded_right_shift(base2_value, -exp2); - return divide_power10(shifted, exp10); + auto const shifted = guarded_right_shift(base2_value, -pow2); + return divide_power10(shifted, pow10); } - return shift_to_decimal_posexp(base2_value, exp2, exp10); - } else { // exp10 < 0 - if (exp2 >= 0) { + return shift_to_decimal_pospow(base2_value, pow2, pow10); + } else { // pow10 < 0 + if (pow2 >= 0) { // Power-2/10 shifts both upward: order doesn't matter, apply and bail. // NOTE: Either shift, multiply, or cast (e.g. double -> decimal32) can overflow! - auto const shifted = guarded_left_shift(static_cast(base2_value), exp2); - return multiply_power10(shifted, -exp10); + auto const shifted = guarded_left_shift(static_cast(base2_value), pow2); + return multiply_power10(shifted, -pow10); } - return shift_to_decimal_negexp(base2_value, exp2, exp10, truncates); + return shift_to_decimal_negpow(base2_value, pow2, pow10, truncates); } }(); @@ -1065,20 +1065,20 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo } /** - * @brief Perform base-10 -> base-2 fixed-point conversion for exp10 > 0 + * @brief Perform base-10 -> base-2 fixed-point conversion for pow10 > 0 * * @tparam DecimalRep The decimal integer type we are converting from * @tparam FloatingType The type of floating point object we are converting to * @param decimal_rep The decimal integer to convert - * @param exp10 The number of powers of 10 to apply to undo the scale factor + * @param pow10 The number of powers of 10 to apply to undo the scale factor * @return A pair of the base-2 value and the remaining powers of 2 to be applied */ template )> -CUDF_HOST_DEVICE inline auto shift_to_binary_posexp(DecimalRep decimal_rep, int exp10) +CUDF_HOST_DEVICE inline auto shift_to_binary_pospow(DecimalRep decimal_rep, int pow10) { - // This is the reverse of shift_to_decimal_posexp(), see that for more details. + // This is the reverse of shift_to_decimal_pospow(), see that for more details. // ShiftingRep: uint64 for float's, __uint128_t for double's using Constants = shifting_constants; @@ -1090,58 +1090,58 @@ CUDF_HOST_DEVICE inline auto shift_to_binary_posexp(DecimalRep decimal_rep, int static constexpr int shift_up_to = sizeof(ShiftingRep) * 8 - Constants::max_bits_shift; int const shift_from = count_significant_bits(decimal_rep); int const num_init_bit_shift = shift_up_to - shift_from; - int exp2 = -num_init_bit_shift; + int pow2 = -num_init_bit_shift; // Perform the initial bit shift ShiftingRep shifting_rep; if constexpr (sizeof(ShiftingRep) < sizeof(DecimalRep)) { // Shift within DecimalRep before dropping to the smaller ShiftingRep - decimal_rep = (exp2 >= 0) ? (decimal_rep >> exp2) : (decimal_rep << -exp2); + decimal_rep = (pow2 >= 0) ? (decimal_rep >> pow2) : (decimal_rep << -pow2); shifting_rep = static_cast(decimal_rep); } else { // Scale up to ShiftingRep before shifting shifting_rep = static_cast(decimal_rep); - shifting_rep = (exp2 >= 0) ? (shifting_rep >> exp2) : (shifting_rep << -exp2); + shifting_rep = (pow2 >= 0) ? (shifting_rep >> pow2) : (shifting_rep << -pow2); } // Iterate, multiplying by 10s and shifting down by 2s until we're almost done - while (exp10 > Constants::max_digits_shift) { + while (pow10 > Constants::max_digits_shift) { // More decimal places to shift than we have room: Multiply the max number of 10s shifting_rep *= Constants::max_digits_shift_pow; - exp10 -= Constants::max_digits_shift; + pow10 -= Constants::max_digits_shift; // Then make more room by bit shifting down by the max # of 2s shifting_rep >>= Constants::max_bits_shift; - exp2 += Constants::max_bits_shift; + pow2 += Constants::max_bits_shift; } // Last 10s-shift: multiply all remaining decimal places // The multiplier is less than the max-shift, and thus fits within 64 / 32 bits if constexpr (Constants::is_double) { - shifting_rep = multiply_power10_64bit(shifting_rep, exp10); + shifting_rep = multiply_power10_64bit(shifting_rep, pow10); } else { - shifting_rep = multiply_power10_32bit(shifting_rep, exp10); + shifting_rep = multiply_power10_32bit(shifting_rep, pow10); } // Our shifting_rep is now the integer mantissa, return it and the powers of 2 - return std::pair{shifting_rep, exp2}; + return std::pair{shifting_rep, pow2}; } /** - * @brief Perform base-10 -> base-2 fixed-point conversion for exp10 < 0 + * @brief Perform base-10 -> base-2 fixed-point conversion for pow10 < 0 * * @tparam DecimalRep The decimal integer type we are converting from * @tparam FloatingType The type of floating point object we are converting to * @param decimal_rep The decimal integer to convert - * @param exp10 The number of powers of 10 to apply to undo the scale factor + * @param pow10 The number of powers of 10 to apply to undo the scale factor * @return A pair of the base-2 value and the remaining powers of 2 to be applied */ template )> -CUDF_HOST_DEVICE inline auto shift_to_binary_negexp(DecimalRep decimal_rep, int const exp10) +CUDF_HOST_DEVICE inline auto shift_to_binary_negpow(DecimalRep decimal_rep, int const pow10) { - // This is the reverse of shift_to_decimal_negexp(), see that for more details. + // This is the reverse of shift_to_decimal_negpow(), see that for more details. // ShiftingRep: uint64 for float's, __uint128_t for double's using Constants = shifting_constants; @@ -1153,45 +1153,45 @@ CUDF_HOST_DEVICE inline auto shift_to_binary_negexp(DecimalRep decimal_rep, int static constexpr int shift_up_to = sizeof(ShiftingRep) * 8 - Constants::num_2s_shift_buffer_bits; int const shift_from = count_significant_bits(decimal_rep); int const num_init_bit_shift = shift_up_to - shift_from; - int exp2 = -num_init_bit_shift; + int pow2 = -num_init_bit_shift; // Perform the initial bit shift ShiftingRep shifting_rep; if constexpr (sizeof(ShiftingRep) < sizeof(DecimalRep)) { // Shift within DecimalRep before dropping to the smaller ShiftingRep - decimal_rep = (exp2 >= 0) ? (decimal_rep >> exp2) : (decimal_rep << -exp2); + decimal_rep = (pow2 >= 0) ? (decimal_rep >> pow2) : (decimal_rep << -pow2); shifting_rep = static_cast(decimal_rep); } else { // Scale up to ShiftingRep before shifting shifting_rep = static_cast(decimal_rep); - shifting_rep = (exp2 >= 0) ? (shifting_rep >> exp2) : (shifting_rep << -exp2); + shifting_rep = (pow2 >= 0) ? (shifting_rep >> pow2) : (shifting_rep << -pow2); } // Convert to using positive values upfront, simpler than doing later. - int exp10_mag = -exp10; + int pow10_mag = -pow10; // Iterate, dividing by 10s and shifting up by 2s until we're almost done - while (exp10_mag > Constants::max_digits_shift) { + while (pow10_mag > Constants::max_digits_shift) { // More decimal places to shift than we have room: Divide the max number of 10s shifting_rep /= Constants::max_digits_shift_pow; - exp10_mag -= Constants::max_digits_shift; + pow10_mag -= Constants::max_digits_shift; // Then make more room by bit shifting up by the max # of 2s shifting_rep <<= Constants::max_bits_shift; - exp2 -= Constants::max_bits_shift; + pow2 -= Constants::max_bits_shift; } // Last 10s-shift: Divdie all remaining decimal places. // This divide result may not fit in the low half of the bit range // But the divisor is less than the max-shift, and thus fits within 64 / 32 bits if constexpr (Constants::is_double) { - shifting_rep = divide_power10_64bit(shifting_rep, exp10_mag); + shifting_rep = divide_power10_64bit(shifting_rep, pow10_mag); } else { - shifting_rep = divide_power10_32bit(shifting_rep, exp10_mag); + shifting_rep = divide_power10_32bit(shifting_rep, pow10_mag); } // Our shifting_rep is now the integer mantissa, return it and the powers of 2 - return std::pair{shifting_rep, exp2}; + return std::pair{shifting_rep, pow2}; } /** @@ -1215,29 +1215,24 @@ CUDF_HOST_DEVICE inline FloatingType convert_integral_to_floating(Rep const& val // Convert to unsigned for bit counting/shifting using UnsignedType = cuda::std::make_unsigned_t; auto const unsigned_value = [&]() -> UnsignedType { - // Use built-in abs functions where available - if constexpr (cuda::std::is_same_v) { - return cuda::std::llabs(value); - } else if constexpr (!cuda::std::is_same_v) { - return cuda::std::abs(value); - } + // Must guard against minimum value, as we can't just negate it: not representable. + if (value == cuda::std::numeric_limits::min()) { return static_cast(value); } // No abs function for 128bit types, so have to do it manually. - // Must guard against minimum value, as we can't just negate it: not representable. - if (value == cuda::std::numeric_limits<__int128_t>::min()) { - return static_cast(value); - } else { + if constexpr (cuda::std::is_same_v) { return static_cast(is_negative ? -value : value); + } else { + return cuda::std::abs(value); } }(); // Shift by powers of 2 and 10 to get our integer mantissa - auto const [mantissa, exp2] = [&]() { - auto const exp10 = static_cast(scale); - if (exp10 >= 0) { - return shift_to_binary_posexp(unsigned_value, exp10); - } else { // exp10 < 0 - return shift_to_binary_negexp(unsigned_value, exp10); + auto const [mantissa, pow2] = [&]() { + auto const pow10 = static_cast(scale); + if (pow10 >= 0) { + return shift_to_binary_pospow(unsigned_value, pow10); + } else { // pow10 < 0 + return shift_to_binary_negpow(unsigned_value, pow10); } }(); @@ -1249,7 +1244,7 @@ CUDF_HOST_DEVICE inline FloatingType convert_integral_to_floating(Rep const& val // Apply the sign and the remaining powers of 2 using converter = floating_converter; - auto const magnitude = converter::add_exp2(floating, exp2); + auto const magnitude = converter::add_pow2(floating, pow2); return converter::set_is_negative(magnitude, is_negative); } From 55859d98fdf5f66b1e022f4d26021ba5414bb909 Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Mon, 8 Jul 2024 17:24:17 -0400 Subject: [PATCH 15/18] Remove last-bit-shift rounding needed to match pyarrow; not desired --- .../cudf/fixed_point/floating_conversion.hpp | 22 +++++-------------- python/cudf/cudf/tests/test_decimal.py | 10 ++++++++- 2 files changed, 15 insertions(+), 17 deletions(-) diff --git a/cpp/include/cudf/fixed_point/floating_conversion.hpp b/cpp/include/cudf/fixed_point/floating_conversion.hpp index 6568461626ab..44b8913ab761 100644 --- a/cpp/include/cudf/fixed_point/floating_conversion.hpp +++ b/cpp/include/cudf/fixed_point/floating_conversion.hpp @@ -19,7 +19,6 @@ #include #include -#include #include #include @@ -760,9 +759,7 @@ struct shifting_constants { * @return integer_rep, shifted 1 and ++'d if the conversion to decimal causes truncation */ template )> -CUDF_HOST_DEVICE cuda::std::tuple add_half_if_truncates(T integer_rep, - int pow2, - int pow10) +CUDF_HOST_DEVICE cuda::std::pair add_half_if_truncates(T integer_rep, int pow2, int pow10) { // The user-supplied scale may truncate information, so we need to talk about rounding. // We have chosen not to round, so we want 1.23456f with scale -4 to be decimal 12345 @@ -809,7 +806,7 @@ CUDF_HOST_DEVICE cuda::std::tuple add_half_if_truncates(T integer_ --pow2; integer_rep += static_cast(conversion_truncates); - return {integer_rep, pow2, conversion_truncates}; + return {integer_rep, pow2}; } /** @@ -900,7 +897,6 @@ shift_to_decimal_pospow(typename shifting_constants::IntegerRep co * @param base2_value The base-2 fixed-point value we are converting from * @param pow2 The number of powers of 2 to apply to convert from base-2 * @param pow10 The number of powers of 10 to apply to reach the desired scale factor - * @param truncating Whether the scale factor truncates floating-point information * @return Magnitude of the converted-to decimal integer */ template ::ShiftingRep shift_to_decimal_negpow(typename shifting_constants::IntegerRep base2_value, int pow2, - int pow10, - bool truncating) + int pow10) { // This is similar to shift_to_decimal_pospow(), except pow10 < 0 & pow2 < 0 // See comments in that function for details. @@ -936,11 +931,7 @@ shift_to_decimal_negpow(typename shifting_constants::IntegerRep ba } // Final bit shifting: Shift may be large, guard against UB - // If we aren't truncating information from the original floating-point number, - // then always round up (to match pandas) by adding 1 to the last bit - shifting_rep = guarded_right_shift(shifting_rep, pow2_mag - 1); - shifting_rep += static_cast(!truncating); - return (shifting_rep >> 1); + return guarded_right_shift(shifting_rep, pow2_mag); }; // If our total decimal shift is less than the max, we don't need to iterate @@ -1012,13 +1003,12 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo auto const pow10 = static_cast(scale); // Add half a bit if truncating to yield expected value, see function for discussion. - auto const [base2_value_bound, pow2_bound, truncates_bound] = + auto const [base2_value_bound, pow2_bound] = add_half_if_truncates(significand, floating_pow2, pow10); // Structured binding variables cannot be captured :/ auto const base2_value = base2_value_bound; auto const pow2 = pow2_bound; - auto const truncates = truncates_bound; // Apply the powers of 2 and 10 to convert to decimal. // The result will be base2_value * (2^pow2) / (10^pow10) @@ -1054,7 +1044,7 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo auto const shifted = guarded_left_shift(static_cast(base2_value), pow2); return multiply_power10(shifted, -pow10); } - return shift_to_decimal_negpow(base2_value, pow2, pow10, truncates); + return shift_to_decimal_negpow(base2_value, pow2, pow10); } }(); diff --git a/python/cudf/cudf/tests/test_decimal.py b/python/cudf/cudf/tests/test_decimal.py index f3b2ac25fa73..65f739bc74a2 100644 --- a/python/cudf/cudf/tests/test_decimal.py +++ b/python/cudf/cudf/tests/test_decimal.py @@ -6,6 +6,7 @@ import numpy as np import pyarrow as pa import pytest +from packaging import version import cudf from cudf.core.column import Decimal32Column, Decimal64Column, NumericalColumn @@ -81,7 +82,6 @@ def test_from_arrow_max_precision_decimal32(): 94.31304, -112.2314, 0.3333333, - 10.03, np.nan, ] ), @@ -93,6 +93,14 @@ def test_from_arrow_max_precision_decimal32(): [Decimal64Dtype(7, 2), Decimal64Dtype(11, 4), Decimal64Dtype(18, 9)], ) def test_typecast_from_float_to_decimal(request, data, from_dtype, to_dtype): + request.applymarker( + pytest.mark.xfail( + condition=version.parse(pa.__version__) >= version.parse("13.0.0") + and from_dtype == np.dtype("float32") + and to_dtype.precision > 12, + reason="https://github.com/rapidsai/cudf/issues/14169", + ) + ) got = data.astype(from_dtype) pa_arr = got.to_arrow().cast( From 5954b6b499c7999f980ceac83f2e923678cbe69a Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Tue, 9 Jul 2024 11:08:32 -0400 Subject: [PATCH 16/18] Fix overflow bug --- .../cudf/fixed_point/floating_conversion.hpp | 45 ++++++++++--------- 1 file changed, 25 insertions(+), 20 deletions(-) diff --git a/cpp/include/cudf/fixed_point/floating_conversion.hpp b/cpp/include/cudf/fixed_point/floating_conversion.hpp index 44b8913ab761..03cc6d2cecc3 100644 --- a/cpp/include/cudf/fixed_point/floating_conversion.hpp +++ b/cpp/include/cudf/fixed_point/floating_conversion.hpp @@ -812,17 +812,18 @@ CUDF_HOST_DEVICE cuda::std::pair add_half_if_truncates(T integer_rep, in /** * @brief Perform base-2 -> base-10 fixed-point conversion for pow10 > 0 * + * @tparam Rep The type of the storage for the decimal value * @tparam FloatingType The type of the original floating-point value we are converting from * @param base2_value The base-2 fixed-point value we are converting from * @param pow2 The number of powers of 2 to apply to convert from base-2 * @param pow10 The number of powers of 10 to apply to reach the desired scale factor * @return Magnitude of the converted-to decimal integer */ -template )> -CUDF_HOST_DEVICE inline typename shifting_constants::ShiftingRep -shift_to_decimal_pospow(typename shifting_constants::IntegerRep const base2_value, - int pow2, - int pow10) +template )> +CUDF_HOST_DEVICE inline cuda::std::make_unsigned_t shift_to_decimal_pospow( + typename shifting_constants::IntegerRep const base2_value, int pow2, int pow10) { // To convert to decimal, we need to apply the input powers of 2 and 10 // The result will be (integer) base2_value * (2^pow2) / (10^pow10) @@ -849,9 +850,12 @@ shift_to_decimal_pospow(typename shifting_constants::IntegerRep co static constexpr int max_init_shift = shift_up_to - shift_from; // If our total bit shift is less than this, we don't need to iterate + using UnsignedRep = cuda::std::make_unsigned_t; if (pow2 <= max_init_shift) { // Shift bits left, divide by 10s to apply the scale factor, and we're done. - return divide_power10(shifting_rep << pow2, pow10); + shifting_rep = divide_power10(shifting_rep << pow2, pow10); + // NOTE: Cast can overflow! + return static_cast(shifting_rep); } // We need to iterate. Do the combined initial shift @@ -867,8 +871,10 @@ shift_to_decimal_pospow(typename shifting_constants::IntegerRep co // If our remaining bit shift is less than the max, we're finished iterating if (pow2 <= Constants::max_bits_shift) { // Shift bits left, divide by 10s to apply the scale factor, and we're done. - // Note: This divide result may not fit in the low half of the bit range - return divide_power10(shifting_rep << pow2, pow10); + shifting_rep = divide_power10(shifting_rep << pow2, pow10); + + // NOTE: Cast can overflow! + return static_cast(shifting_rep); } // Shift the max number of bits left again @@ -886,13 +892,14 @@ shift_to_decimal_pospow(typename shifting_constants::IntegerRep co } // Final bit shift: Shift may be large, guard against UB - // NOTE: This can overflow! - return guarded_left_shift(shifting_rep, pow2); + // NOTE: This can overflow (both cast and shift)! + return guarded_left_shift(static_cast(shifting_rep), pow2); } /** * @brief Perform base-2 -> base-10 fixed-point conversion for pow10 < 0 * + * @tparam Rep The type of the storage for the decimal value * @tparam FloatingType The type of the original floating-point value we are converting from * @param base2_value The base-2 fixed-point value we are converting from * @param pow2 The number of powers of 2 to apply to convert from base-2 @@ -902,10 +909,8 @@ shift_to_decimal_pospow(typename shifting_constants::IntegerRep co template )> -CUDF_HOST_DEVICE inline typename shifting_constants::ShiftingRep -shift_to_decimal_negpow(typename shifting_constants::IntegerRep base2_value, - int pow2, - int pow10) +CUDF_HOST_DEVICE inline cuda::std::make_unsigned_t shift_to_decimal_negpow( + typename shifting_constants::IntegerRep base2_value, int pow2, int pow10) { // This is similar to shift_to_decimal_pospow(), except pow10 < 0 & pow2 < 0 // See comments in that function for details. @@ -921,6 +926,7 @@ shift_to_decimal_negpow(typename shifting_constants::IntegerRep ba int pow2_mag = -pow2; // For performing final 10s-shift + using UnsignedRep = cuda::std::make_unsigned_t; auto final_shifts_low10s = [&]() { // Last 10s-shift: multiply all remaining decimal places, shift all remaining bits, then bail // The multiplier is less than the max-shift, and thus fits within 64 / 32 bits @@ -931,7 +937,7 @@ shift_to_decimal_negpow(typename shifting_constants::IntegerRep ba } // Final bit shifting: Shift may be large, guard against UB - return guarded_right_shift(shifting_rep, pow2_mag); + return static_cast(guarded_right_shift(shifting_rep, pow2_mag)); }; // If our total decimal shift is less than the max, we don't need to iterate @@ -963,7 +969,6 @@ shift_to_decimal_negpow(typename shifting_constants::IntegerRep ba // We need to convert to the output rep for the final scale-factor multiply, because if (e.g.) // float -> dec128 and some large pow10_mag, it might overflow the 64bit shifting rep. // It's not needed for pow10 > 0 because we're dividing by 10s there instead of multiplying. - using UnsignedRep = cuda::std::make_unsigned_t; // NOTE: This can overflow! (Both multiply and cast) return multiply_power10(static_cast(shifting_rep), pow10_mag); } @@ -997,7 +1002,7 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo auto const integer_rep = converter::bit_cast_to_integer(floating); if (converter::is_zero(integer_rep)) { return 0; } - // Note that the base2_value here is an unsigned integer with sizeof(FloatingType) + // Note that the significand here is an unsigned integer with sizeof(FloatingType) auto const is_negative = converter::get_is_negative(integer_rep); auto const [significand, floating_pow2] = converter::get_significand_and_pow2(integer_rep); auto const pow10 = static_cast(scale); @@ -1027,16 +1032,16 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo if (pow2 >= 0) { return guarded_left_shift(static_cast(base2_value), pow2); } else { - return guarded_right_shift(base2_value, -pow2); + return static_cast(guarded_right_shift(base2_value, -pow2)); } } else if (pow10 > 0) { if (pow2 <= 0) { // Power-2/10 shifts both downward: order doesn't matter, apply and bail. // Guard against shift being undefined behavior auto const shifted = guarded_right_shift(base2_value, -pow2); - return divide_power10(shifted, pow10); + return static_cast(divide_power10(shifted, pow10)); } - return shift_to_decimal_pospow(base2_value, pow2, pow10); + return shift_to_decimal_pospow(base2_value, pow2, pow10); } else { // pow10 < 0 if (pow2 >= 0) { // Power-2/10 shifts both upward: order doesn't matter, apply and bail. From 738f1a5372e39efd1f3f54477147e87d061e34a3 Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Tue, 9 Jul 2024 11:33:55 -0400 Subject: [PATCH 17/18] Split shifting into separate function. --- .../cudf/fixed_point/floating_conversion.hpp | 102 ++++++++++-------- 1 file changed, 58 insertions(+), 44 deletions(-) diff --git a/cpp/include/cudf/fixed_point/floating_conversion.hpp b/cpp/include/cudf/fixed_point/floating_conversion.hpp index 03cc6d2cecc3..36822248add4 100644 --- a/cpp/include/cudf/fixed_point/floating_conversion.hpp +++ b/cpp/include/cudf/fixed_point/floating_conversion.hpp @@ -107,7 +107,7 @@ struct floating_converter { // The value of the mantissa is in the range [1, 2). /// # significand bits (includes understood bit) static constexpr int num_significand_bits = cuda::std::numeric_limits::digits; - /// # mantissa bits (-1 for understood bit) + /// # stored mantissa bits (-1 for understood bit) static constexpr int num_stored_mantissa_bits = num_significand_bits - 1; /// The mask for the understood bit static constexpr IntegralType understood_bit_mask = (IntegralType(1) << num_stored_mantissa_bits); @@ -982,6 +982,59 @@ CUDF_HOST_DEVICE inline cuda::std::make_unsigned_t shift_to_decimal_negpow( return final_shifts_low10s(); } +/** + * @brief Perform base-2 -> base-10 fixed-point conversion + * + * @tparam Rep The type of integer we are converting to, to store the decimal value + * @tparam FloatingType The type of floating-point object we are converting from + * @param base2_value The base-2 fixed-point value we are converting from + * @param pow2 The number of powers of 2 to apply to convert from base-2 + * @param pow10 The number of powers of 10 to apply to reach the desired scale factor + * @return Integer representation of the floating-point value, given the desired scale + */ +template )> +CUDF_HOST_DEVICE inline cuda::std::make_unsigned_t convert_floating_to_integral_shifting( + typename floating_converter::IntegralType base2_value, int pow10, int pow2) +{ + // Apply the powers of 2 and 10 to convert to decimal. + // The result will be base2_value * (2^pow2) / (10^pow10) + + // Note that while this code is branchy, the decimal scale factor is part of the + // column type itself, so every thread will take the same branches on pow10. + // Also data within a column tends to be similar, so they will often take the + // same branches on pow2 as well. + + // NOTE: some returns here can overflow (e.g. ShiftingRep -> UnsignedRep) + using UnsignedRep = cuda::std::make_unsigned_t; + if (pow10 == 0) { + // NOTE: Left Bit-shift can overflow! As can cast! (e.g. double -> decimal32) + // Bit shifts may be large, guard against UB + if (pow2 >= 0) { + return guarded_left_shift(static_cast(base2_value), pow2); + } else { + return static_cast(guarded_right_shift(base2_value, -pow2)); + } + } else if (pow10 > 0) { + if (pow2 <= 0) { + // Power-2/10 shifts both downward: order doesn't matter, apply and bail. + // Guard against shift being undefined behavior + auto const shifted = guarded_right_shift(base2_value, -pow2); + return static_cast(divide_power10(shifted, pow10)); + } + return shift_to_decimal_pospow(base2_value, pow2, pow10); + } else { // pow10 < 0 + if (pow2 >= 0) { + // Power-2/10 shifts both upward: order doesn't matter, apply and bail. + // NOTE: Either shift, multiply, or cast (e.g. double -> decimal32) can overflow! + auto const shifted = guarded_left_shift(static_cast(base2_value), pow2); + return multiply_power10(shifted, -pow10); + } + return shift_to_decimal_negpow(base2_value, pow2, pow10); + } +} + /** * @brief Perform floating-point -> integer decimal conversion * @@ -1005,53 +1058,14 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo // Note that the significand here is an unsigned integer with sizeof(FloatingType) auto const is_negative = converter::get_is_negative(integer_rep); auto const [significand, floating_pow2] = converter::get_significand_and_pow2(integer_rep); - auto const pow10 = static_cast(scale); // Add half a bit if truncating to yield expected value, see function for discussion. - auto const [base2_value_bound, pow2_bound] = - add_half_if_truncates(significand, floating_pow2, pow10); - - // Structured binding variables cannot be captured :/ - auto const base2_value = base2_value_bound; - auto const pow2 = pow2_bound; + auto const pow10 = static_cast(scale); + auto const [base2_value, pow2] = add_half_if_truncates(significand, floating_pow2, pow10); // Apply the powers of 2 and 10 to convert to decimal. - // The result will be base2_value * (2^pow2) / (10^pow10) - // - // Note that while this code is branchy, the decimal scale factor is part of the - // column type itself, so every thread will take the same branches on pow10. - // Also data within a column tends to be similar, so they will often take the - // same branches on pow2 as well. - // - // NOTE: some returns here can overflow (e.g. ShiftingRep -> UnsignedRep) - using UnsignedRep = cuda::std::make_unsigned_t; - auto const magnitude = [&]() -> UnsignedRep { - if (pow10 == 0) { - // NOTE: Left Bit-shift can overflow! As can cast! (e.g. double -> decimal32) - // Bit shifts may be large, guard against UB - if (pow2 >= 0) { - return guarded_left_shift(static_cast(base2_value), pow2); - } else { - return static_cast(guarded_right_shift(base2_value, -pow2)); - } - } else if (pow10 > 0) { - if (pow2 <= 0) { - // Power-2/10 shifts both downward: order doesn't matter, apply and bail. - // Guard against shift being undefined behavior - auto const shifted = guarded_right_shift(base2_value, -pow2); - return static_cast(divide_power10(shifted, pow10)); - } - return shift_to_decimal_pospow(base2_value, pow2, pow10); - } else { // pow10 < 0 - if (pow2 >= 0) { - // Power-2/10 shifts both upward: order doesn't matter, apply and bail. - // NOTE: Either shift, multiply, or cast (e.g. double -> decimal32) can overflow! - auto const shifted = guarded_left_shift(static_cast(base2_value), pow2); - return multiply_power10(shifted, -pow10); - } - return shift_to_decimal_negpow(base2_value, pow2, pow10); - } - }(); + auto const magnitude = + convert_floating_to_integral_shifting(base2_value, pow10, pow2); // Reapply the sign and return // NOTE: Cast can overflow! From 060667305f7c8082c147c82b13faa39661e88d7d Mon Sep 17 00:00:00 2001 From: Paul Mattione Date: Wed, 10 Jul 2024 18:30:16 -0400 Subject: [PATCH 18/18] Don't add half bit if converting a whole number! --- .../cudf/fixed_point/floating_conversion.hpp | 24 ++++++++++++++----- 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/cpp/include/cudf/fixed_point/floating_conversion.hpp b/cpp/include/cudf/fixed_point/floating_conversion.hpp index 36822248add4..c64ae8877d49 100644 --- a/cpp/include/cudf/fixed_point/floating_conversion.hpp +++ b/cpp/include/cudf/fixed_point/floating_conversion.hpp @@ -18,6 +18,7 @@ #include +#include #include #include @@ -752,14 +753,19 @@ struct shifting_constants { * * @note This fixes problems like 1.2 (value = 1.1999...) at scale -1 -> 11 * - * @tparam T Type of integer holding the floating-point significand + * @tparam FloatingType Type of integer holding the floating-point significand + * @param floating The floating-point number to convert * @param integer_rep The integer representation of the floating-point significand * @param pow2 The power of 2 that needs to be applied to the significand * @param pow10 The power of 10 that needs to be applied to the significand * @return integer_rep, shifted 1 and ++'d if the conversion to decimal causes truncation */ -template )> -CUDF_HOST_DEVICE cuda::std::pair add_half_if_truncates(T integer_rep, int pow2, int pow10) +template )> +CUDF_HOST_DEVICE cuda::std::pair::IntegralType, int> +add_half_if_truncates(FloatingType floating, + typename floating_converter::IntegralType integer_rep, + int pow2, + int pow10) { // The user-supplied scale may truncate information, so we need to talk about rounding. // We have chosen not to round, so we want 1.23456f with scale -4 to be decimal 12345 @@ -801,10 +807,15 @@ CUDF_HOST_DEVICE cuda::std::pair add_half_if_truncates(T integer_rep, in bool const conversion_truncates = (pow10_term > pow2_term) || ((pow2_term == pow10_term) && (pow2 < 0)); + // However, don't add a half-bit if the input is a whole number! + // This is only for errors introduced by rounding decimal fractions! + bool const is_whole_number = (cuda::std::floor(floating) == floating); + bool const add_half_bit = conversion_truncates && !is_whole_number; + // Add half a bit on truncation (shift to make room and update pow2) integer_rep <<= 1; --pow2; - integer_rep += static_cast(conversion_truncates); + integer_rep += static_cast(add_half_bit); return {integer_rep, pow2}; } @@ -1060,8 +1071,9 @@ CUDF_HOST_DEVICE inline Rep convert_floating_to_integral(FloatingType const& flo auto const [significand, floating_pow2] = converter::get_significand_and_pow2(integer_rep); // Add half a bit if truncating to yield expected value, see function for discussion. - auto const pow10 = static_cast(scale); - auto const [base2_value, pow2] = add_half_if_truncates(significand, floating_pow2, pow10); + auto const pow10 = static_cast(scale); + auto const [base2_value, pow2] = + add_half_if_truncates(floating, significand, floating_pow2, pow10); // Apply the powers of 2 and 10 to convert to decimal. auto const magnitude =