@@ -519,322 +519,6 @@ constexpr std::array<double, 2 * kPrecomputedPowersOfTen + 1> kDoublePowersOfTen
519519 1e56 , 1e57 , 1e58 , 1e59 , 1e60 , 1e61 , 1e62 , 1e63 , 1e64 , 1e65 , 1e66 , 1e67 ,
520520 1e68 , 1e69 , 1e70 , 1e71 , 1e72 , 1e73 , 1e74 , 1e75 , 1e76 };
521521
522- // ceil(log2(10 ^ k)) for k in [0...76]
523- constexpr std::array<int32_t , 76 + 1 > kCeilLog2PowersOfTen = {
524- 0 , 4 , 7 , 10 , 14 , 17 , 20 , 24 , 27 , 30 , 34 , 37 , 40 , 44 , 47 , 50 ,
525- 54 , 57 , 60 , 64 , 67 , 70 , 74 , 77 , 80 , 84 , 87 , 90 , 94 , 97 , 100 , 103 ,
526- 107 , 110 , 113 , 117 , 120 , 123 , 127 , 130 , 133 , 137 , 140 , 143 , 147 , 150 , 153 , 157 ,
527- 160 , 163 , 167 , 170 , 173 , 177 , 180 , 183 , 187 , 190 , 193 , 196 , 200 , 203 , 206 , 210 ,
528- 213 , 216 , 220 , 223 , 226 , 230 , 233 , 236 , 240 , 243 , 246 , 250 , 253 };
529-
530- template <typename Real>
531- struct RealTraits {};
532-
533- template <>
534- struct RealTraits <float > {
535- static constexpr const float * powers_of_ten () { return kFloatPowersOfTen .data (); }
536-
537- static constexpr float two_to_64 (float x) { return x * 1 .8446744e+19f ; }
538-
539- static constexpr int32_t kMantissaBits = 24 ;
540- // ceil(log10(2 ^ kMantissaBits))
541- static constexpr int32_t kMantissaDigits = 8 ;
542- // Integers between zero and kMaxPreciseInteger can be precisely represented
543- static constexpr uint64_t kMaxPreciseInteger = (1ULL << kMantissaBits ) - 1 ;
544- };
545-
546- template <>
547- struct RealTraits <double > {
548- static constexpr const double * powers_of_ten () { return kDoublePowersOfTen .data (); }
549-
550- static constexpr double two_to_64 (double x) { return x * 1.8446744073709552e+19 ; }
551-
552- static constexpr int32_t kMantissaBits = 53 ;
553- // ceil(log10(2 ^ kMantissaBits))
554- static constexpr int32_t kMantissaDigits = 16 ;
555- // Integers between zero and kMaxPreciseInteger can be precisely represented
556- static constexpr uint64_t kMaxPreciseInteger = (1ULL << kMantissaBits ) - 1 ;
557- };
558-
559- struct DecimalRealConversion {
560- // Return 10**exp, with a fast lookup, assuming `exp` is within bounds
561- template <typename Real>
562- static Real PowerOfTen (int32_t exp) {
563- constexpr int32_t N = kPrecomputedPowersOfTen ;
564- ICEBERG_DCHECK (exp >= -N && exp <= N, " " );
565- return RealTraits<Real>::powers_of_ten ()[N + exp];
566- }
567-
568- // Return 10**exp, with a fast lookup if possible
569- template <typename Real>
570- static Real LargePowerOfTen (int32_t exp) {
571- constexpr int32_t N = kPrecomputedPowersOfTen ;
572- if (exp >= -N && exp <= N) {
573- return RealTraits<Real>::powers_of_ten ()[N + exp];
574- } else {
575- return std::pow (static_cast <Real>(10.0 ), static_cast <Real>(exp));
576- }
577- }
578-
579- static constexpr int32_t kMaxPrecision = Decimal::kMaxPrecision ;
580- static constexpr int32_t kMaxScale = Decimal::kMaxScale ;
581-
582- static constexpr auto & DecimalPowerOfTen (int32_t exp) {
583- ICEBERG_DCHECK (exp >= 0 && exp <= kMaxPrecision , " " );
584- return kDecimal128PowersOfTen [exp];
585- }
586-
587- // Right shift positive `x` by positive `bits`, rounded half to even
588- static Decimal RoundedRightShift (const Decimal& x, int32_t bits) {
589- if (bits == 0 ) {
590- return x;
591- }
592- int64_t result_hi = x.high ();
593- uint64_t result_lo = x.low ();
594- uint64_t shifted = 0 ;
595- while (bits >= 64 ) {
596- // Retain the information that set bits were shifted right.
597- // This is important to detect an exact half.
598- shifted = result_lo | (shifted > 0 );
599- result_lo = result_hi;
600- result_hi >>= 63 ; // for sign
601- bits -= 64 ;
602- }
603- if (bits > 0 ) {
604- shifted = (result_lo << (64 - bits)) | (shifted > 0 );
605- result_lo >>= bits;
606- result_lo |= static_cast <uint64_t >(result_hi) << (64 - bits);
607- result_hi >>= bits;
608- }
609- // We almost have our result, but now do the rounding.
610- constexpr uint64_t kHalf = 0x8000000000000000ULL ;
611- if (shifted > kHalf ) {
612- // Strictly more than half => round up
613- result_lo += 1 ;
614- result_hi += (result_lo == 0 );
615- } else if (shifted == kHalf ) {
616- // Exactly half => round to even
617- if ((result_lo & 1 ) != 0 ) {
618- result_lo += 1 ;
619- result_hi += (result_lo == 0 );
620- }
621- } else {
622- // Strictly less than half => round down
623- }
624- return Decimal{result_hi, result_lo};
625- }
626-
627- template <typename Real>
628- static Result<Decimal> FromPositiveRealApprox (Real real, int32_t precision,
629- int32_t scale) {
630- // Approximate algorithm that operates in the FP domain (thus subject
631- // to precision loss).
632- const auto x = std::nearbyint (real * PowerOfTen<Real>(scale));
633- const auto max_abs = PowerOfTen<Real>(precision);
634- if (x <= -max_abs || x >= max_abs) {
635- return Invalid (" Cannot convert {} to Decimal(precision = {}, scale = {}): overflow" ,
636- real, precision, scale);
637- }
638-
639- // Extract high and low bits
640- const auto high = std::floor (std::ldexp (x, -64 ));
641- const auto low = x - std::ldexp (high, 64 );
642-
643- ICEBERG_DCHECK (high >= 0 , " " );
644- ICEBERG_DCHECK (high < 9.223372036854776e+18 , " " ); // 2**63
645- ICEBERG_DCHECK (low >= 0 , " " );
646- ICEBERG_DCHECK (low < 1.8446744073709552e+19 , " " ); // 2**64
647- return Decimal (static_cast <int64_t >(high), static_cast <uint64_t >(low));
648- }
649-
650- template <typename Real>
651- static Result<Decimal> FromPositiveReal (Real real, int32_t precision, int32_t scale) {
652- constexpr int32_t kMantissaBits = RealTraits<Real>::kMantissaBits ;
653- constexpr int32_t kMantissaDigits = RealTraits<Real>::kMantissaDigits ;
654-
655- // Problem statement: construct the Decimal with the value
656- // closest to `real * 10^scale`.
657- if (scale < 0 ) {
658- // Negative scales are not handled below, fall back to approx algorithm
659- return FromPositiveRealApprox (real, precision, scale);
660- }
661-
662- // 1. Check that `real` is within acceptable bounds.
663- const Real limit = PowerOfTen<Real>(precision - scale);
664- if (real > limit) {
665- // Checking the limit early helps ensure the computations below do not
666- // overflow.
667- // NOTE: `limit` is allowed here as rounding can make it smaller than
668- // the theoretical limit (for example, 1.0e23 < 10^23).
669- return Invalid (" Cannot convert {} to Decimal(precision = {}, scale = {}): overflow" ,
670- real, precision, scale);
671- }
672-
673- // 2. Losslessly convert `real` to `mant * 2**k`
674- int32_t binary_exp = 0 ;
675- const Real real_mant = std::frexp (real, &binary_exp);
676- // `real_mant` is within 0.5 and 1 and has M bits of precision.
677- // Multiply it by 2^M to get an exact integer.
678- const auto mant = static_cast <uint64_t >(std::ldexp (real_mant, kMantissaBits ));
679- const int32_t k = binary_exp - kMantissaBits ;
680- // (note that `real = mant * 2^k`)
681-
682- // 3. Start with `mant`.
683- // We want to end up with `real * 10^scale` i.e. `mant * 2^k * 10^scale`.
684- Decimal x (mant);
685-
686- if (k < 0 ) {
687- // k < 0 (i.e. binary_exp < kMantissaBits), is probably the common case
688- // when converting to decimal. It implies right-shifting by -k bits,
689- // while multiplying by 10^scale. We also must avoid overflow (losing
690- // bits on the left) and precision loss (losing bits on the right).
691- int32_t right_shift_by = -k;
692- int32_t mul_by_ten_to = scale;
693-
694- // At this point, `x` has kMantissaDigits significant digits but it can
695- // fit kMaxPrecision (excluding sign). We can therefore multiply by up
696- // to 10^(kMaxPrecision - kMantissaDigits).
697- constexpr int32_t kSafeMulByTenTo = kMaxPrecision - kMantissaDigits ;
698-
699- if (mul_by_ten_to <= kSafeMulByTenTo ) {
700- // Scale is small enough, so we can do it all at once.
701- x *= DecimalPowerOfTen (mul_by_ten_to);
702- x = RoundedRightShift (x, right_shift_by);
703- } else {
704- // Scale is too large, we cannot multiply at once without overflow.
705- // We use an iterative algorithm which alternately shifts left by
706- // multiplying by a power of ten, and shifts right by a number of bits.
707-
708- // First multiply `x` by as large a power of ten as possible
709- // without overflowing.
710- x *= DecimalPowerOfTen (kSafeMulByTenTo );
711- mul_by_ten_to -= kSafeMulByTenTo ;
712-
713- // `x` now has full precision. However, we know we'll only
714- // keep `precision` digits at the end. Extraneous bits/digits
715- // on the right can be safely shifted away, before multiplying
716- // again.
717- // NOTE: if `precision` is the full precision then the algorithm will
718- // lose the last digit. If `precision` is almost the full precision,
719- // there can be an off-by-one error due to rounding.
720- const int32_t mul_step = std::max (1 , kMaxPrecision - precision);
721-
722- // The running exponent, useful to compute by how much we must
723- // shift right to make place on the left before the next multiply.
724- int32_t total_exp = 0 ;
725- int32_t total_shift = 0 ;
726- while (mul_by_ten_to > 0 && right_shift_by > 0 ) {
727- const int32_t exp = std::min (mul_by_ten_to, mul_step);
728- total_exp += exp;
729- // The supplementary right shift required so that
730- // `x * 10^total_exp / 2^total_shift` fits in the decimal.
731- ICEBERG_DCHECK (static_cast <size_t >(total_exp) < sizeof (kCeilLog2PowersOfTen ),
732- " " );
733- const int32_t bits =
734- std::min (right_shift_by, kCeilLog2PowersOfTen [total_exp] - total_shift);
735- total_shift += bits;
736- // Right shift to make place on the left, then multiply
737- x = RoundedRightShift (x, bits);
738- right_shift_by -= bits;
739- // Should not overflow thanks to the precautions taken
740- x *= DecimalPowerOfTen (exp);
741- mul_by_ten_to -= exp;
742- }
743- if (mul_by_ten_to > 0 ) {
744- x *= DecimalPowerOfTen (mul_by_ten_to);
745- }
746- if (right_shift_by > 0 ) {
747- x = RoundedRightShift (x, right_shift_by);
748- }
749- }
750- } else {
751- // k >= 0 implies left-shifting by k bits and multiplying by 10^scale.
752- // The order of these operations therefore doesn't matter. We know
753- // we won't overflow because of the limit check above, and we also
754- // won't lose any significant bits on the right.
755- x *= DecimalPowerOfTen (scale);
756- x <<= k;
757- }
758-
759- // Rounding might have pushed `x` just above the max precision, check again
760- if (!x.FitsInPrecision (precision)) {
761- return Invalid (" Cannot convert {} to Decimal(precision = {}, scale = {}): overflow" ,
762- real, precision, scale);
763- }
764- return x;
765- }
766-
767- template <typename Real>
768- static Real ToRealPositiveNoSplit (const Decimal& decimal, int32_t scale) {
769- Real x = RealTraits<Real>::two_to_64 (static_cast <Real>(decimal.high ()));
770- x += static_cast <Real>(decimal.low ());
771- x *= LargePowerOfTen<Real>(-scale);
772- return x;
773- }
774-
775- // / An approximate conversion from Decimal to Real that guarantees:
776- // / 1. If the decimal is an integer, the conversion is exact.
777- // / 2. If the number of fractional digits is <= RealTraits<Real>::kMantissaDigits (e.g.
778- // / 8 for float and 16 for double), the conversion is within 1 ULP of the exact
779- // / value.
780- // / 3. Otherwise, the conversion is within 2^(-RealTraits<Real>::kMantissaDigits+1)
781- // / (e.g. 2^-23 for float and 2^-52 for double) of the exact value.
782- // / Here "exact value" means the closest representable value by Real.
783- template <typename Real>
784- static Real ToRealPositive (const Decimal& decimal, int32_t scale) {
785- if (scale <= 0 ||
786- (decimal.high () == 0 && decimal.low () <= RealTraits<Real>::kMaxPreciseInteger )) {
787- // No need to split the decimal if it is already an integer (scale <= 0) or if it
788- // can be precisely represented by Real
789- return ToRealPositiveNoSplit<Real>(decimal, scale);
790- }
791-
792- // Split decimal into whole and fractional parts to avoid precision loss
793- auto s = decimal.GetWholeAndFraction (scale);
794- ICEBERG_DCHECK (s, " Decimal::GetWholeAndFraction failed" );
795-
796- Real whole = ToRealPositiveNoSplit<Real>(s->first , 0 );
797- Real fraction = ToRealPositiveNoSplit<Real>(s->second , scale);
798-
799- return whole + fraction;
800- }
801-
802- template <typename Real>
803- static Result<Decimal> FromReal (Real value, int32_t precision, int32_t scale) {
804- ICEBERG_DCHECK (precision >= 1 && precision <= kMaxPrecision , " " );
805- ICEBERG_DCHECK (scale >= -kMaxScale && scale <= kMaxScale , " " );
806-
807- if (!std::isfinite (value)) {
808- return InvalidArgument (" Cannot convert {} to Decimal" , value);
809- }
810-
811- if (value == 0 ) {
812- return Decimal{};
813- }
814-
815- if (value < 0 ) {
816- ICEBERG_ASSIGN_OR_RAISE (auto decimal, FromPositiveReal (-value, precision, scale));
817- return decimal.Negate ();
818- } else {
819- return FromPositiveReal (value, precision, scale);
820- }
821- }
822-
823- template <typename Real>
824- static Real ToReal (const Decimal& decimal, int32_t scale) {
825- ICEBERG_DCHECK (scale >= -kMaxScale && scale <= kMaxScale , " " );
826-
827- if (decimal.IsNegative ()) {
828- // Convert the absolute value to avoid precision loss
829- auto abs = decimal;
830- abs.Negate ();
831- return -ToRealPositive<Real>(abs, scale);
832- } else {
833- return ToRealPositive<Real>(decimal, scale);
834- }
835- }
836- };
837-
838522// Helper function used by Decimal::FromBigEndian
839523static inline uint64_t UInt64FromBigEndian (const uint8_t * bytes, int32_t length) {
840524 // We don't bounds check the length here because this is called by
@@ -869,21 +553,6 @@ static bool RescaleWouldCauseDataLoss(const Decimal& value, int32_t delta_scale,
869553
870554} // namespace
871555
872- Result<Decimal> Decimal::FromReal (float x, int32_t precision, int32_t scale) {
873- return DecimalRealConversion::FromReal (x, precision, scale);
874- }
875-
876- Result<Decimal> Decimal::FromReal (double x, int32_t precision, int32_t scale) {
877- return DecimalRealConversion::FromReal (x, precision, scale);
878- }
879-
880- Result<std::pair<Decimal, Decimal>> Decimal::GetWholeAndFraction (int32_t scale) const {
881- ICEBERG_DCHECK (scale >= 0 && scale <= kMaxScale , " " );
882-
883- Decimal multiplier (kDecimal128PowersOfTen [scale]);
884- return Divide (multiplier);
885- }
886-
887556Result<Decimal> Decimal::FromBigEndian (const uint8_t * bytes, int32_t length) {
888557 static constexpr int32_t kMinDecimalBytes = 1 ;
889558 static constexpr int32_t kMaxDecimalBytes = 16 ;
@@ -966,14 +635,6 @@ bool Decimal::FitsInPrecision(int32_t precision) const {
966635 return Decimal::Abs (*this ) < kDecimal128PowersOfTen [precision];
967636}
968637
969- float Decimal::ToFloat (int32_t scale) const {
970- return DecimalRealConversion::ToReal<float >(*this , scale);
971- }
972-
973- double Decimal::ToDouble (int32_t scale) const {
974- return DecimalRealConversion::ToReal<double >(*this , scale);
975- }
976-
977638std::array<uint8_t , Decimal::kByteWidth > Decimal::ToBytes () const {
978639 std::array<uint8_t , kByteWidth > out{{0 }};
979640 std::memcpy (out.data (), &data_, kByteWidth );
0 commit comments