Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-04-19 08:19:37

0001 // Copyright 2020-2023 Daniel Lemire
0002 // Copyright 2023 Matt Borland
0003 // Distributed under the Boost Software License, Version 1.0.
0004 // https://www.boost.org/LICENSE_1_0.txt
0005 //
0006 // Derivative of: https://github.com/fastfloat/fast_float
0007 
0008 #ifndef BOOST_CHARCONV_DETAIL_FASTFLOAT_DECIMAL_TO_BINARY_HPP
0009 #define BOOST_CHARCONV_DETAIL_FASTFLOAT_DECIMAL_TO_BINARY_HPP
0010 
0011 #include <boost/charconv/detail/fast_float/float_common.hpp>
0012 #include <boost/charconv/detail/fast_float/fast_table.hpp>
0013 #include <cfloat>
0014 #include <cinttypes>
0015 #include <cmath>
0016 #include <cstdint>
0017 #include <cstdlib>
0018 #include <cstring>
0019 
0020 namespace boost { namespace charconv { namespace detail { namespace fast_float {
0021 
0022 // This will compute or rather approximate w * 5**q and return a pair of 64-bit words approximating
0023 // the result, with the "high" part corresponding to the most significant bits and the
0024 // low part corresponding to the least significant bits.
0025 //
0026 template <int bit_precision>
0027 BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
0028 value128 compute_product_approximation(int64_t q, uint64_t w) {
0029   const int index = 2 * int(q - powers::smallest_power_of_five);
0030   // For small values of q, e.g., q in [0,27], the answer is always exact because
0031   // The line value128 firstproduct = full_multiplication(w, power_of_five_128[index]);
0032   // gives the exact answer.
0033   value128 firstproduct = full_multiplication(w, powers::power_of_five_128[index]);
0034   static_assert((bit_precision >= 0) && (bit_precision <= 64), " precision should  be in (0,64]");
0035   constexpr uint64_t precision_mask = (bit_precision < 64) ?
0036                (uint64_t(0xFFFFFFFFFFFFFFFF) >> bit_precision)
0037                : uint64_t(0xFFFFFFFFFFFFFFFF);
0038   if((firstproduct.high & precision_mask) == precision_mask) { // could further guard with  (lower + w < lower)
0039     // regarding the second product, we only need secondproduct.high, but our expectation is that the compiler will optimize this extra work away if needed.
0040     value128 secondproduct = full_multiplication(w, powers::power_of_five_128[index + 1]);
0041     firstproduct.low += secondproduct.high;
0042     if(secondproduct.high > firstproduct.low) {
0043       firstproduct.high++;
0044     }
0045   }
0046   return firstproduct;
0047 }
0048 
0049 namespace detail {
0050 /**
0051  * For q in (0,350), we have that
0052  *  f = (((152170 + 65536) * q ) >> 16);
0053  * is equal to
0054  *   floor(p) + q
0055  * where
0056  *   p = log(5**q)/log(2) = q * log(5)/log(2)
0057  *
0058  * For negative values of q in (-400,0), we have that
0059  *  f = (((152170 + 65536) * q ) >> 16);
0060  * is equal to
0061  *   -ceil(p) + q
0062  * where
0063  *   p = log(5**-q)/log(2) = -q * log(5)/log(2)
0064  */
0065   constexpr BOOST_FORCEINLINE int32_t power(int32_t q)  noexcept  {
0066     return (((152170 + 65536) * q) >> 16) + 63;
0067   }
0068 } // namespace detail
0069 
0070 // create an adjusted mantissa, biased by the invalid power2
0071 // for significant digits already multiplied by 10 ** q.
0072 template <typename binary>
0073 BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR14
0074 adjusted_mantissa compute_error_scaled(int64_t q, uint64_t w, int lz) noexcept  {
0075   int hilz = int(w >> 63) ^ 1;
0076   adjusted_mantissa answer;
0077   answer.mantissa = w << hilz;
0078   int bias = binary::mantissa_explicit_bits() - binary::minimum_exponent();
0079   answer.power2 = int32_t(detail::power(int32_t(q)) + bias - hilz - lz - 62 + invalid_am_bias);
0080   return answer;
0081 }
0082 
0083 // w * 10 ** q, without rounding the representation up.
0084 // the power2 in the exponent will be adjusted by invalid_am_bias.
0085 template <typename binary>
0086 BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
0087 adjusted_mantissa compute_error(int64_t q, uint64_t w)  noexcept  {
0088   int lz = leading_zeroes(w);
0089   w <<= lz;
0090   value128 product = compute_product_approximation<binary::mantissa_explicit_bits() + 3>(q, w);
0091   return compute_error_scaled<binary>(q, product.high, lz);
0092 }
0093 
0094 // w * 10 ** q
0095 // The returned value should be a valid ieee64 number that simply need to be packed.
0096 // However, in some very rare cases, the computation will fail. In such cases, we
0097 // return an adjusted_mantissa with a negative power of 2: the caller should recompute
0098 // in such cases.
0099 template <typename binary>
0100 BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
0101 adjusted_mantissa compute_float(int64_t q, uint64_t w)  noexcept  {
0102   adjusted_mantissa answer;
0103   if ((w == 0) || (q < binary::smallest_power_of_ten())) {
0104     answer.power2 = 0;
0105     answer.mantissa = 0;
0106     // result should be zero
0107     return answer;
0108   }
0109   if (q > binary::largest_power_of_ten()) {
0110     // we want to get infinity:
0111     answer.power2 = binary::infinite_power();
0112     answer.mantissa = 0;
0113     return answer;
0114   }
0115   // At this point in time q is in [powers::smallest_power_of_five, powers::largest_power_of_five].
0116 
0117   // We want the most significant bit of i to be 1. Shift if needed.
0118   int lz = leading_zeroes(w);
0119   w <<= lz;
0120 
0121   // The required precision is binary::mantissa_explicit_bits() + 3 because
0122   // 1. We need the implicit bit
0123   // 2. We need an extra bit for rounding purposes
0124   // 3. We might lose a bit due to the "upperbit" routine (result too small, requiring a shift)
0125 
0126   value128 product = compute_product_approximation<binary::mantissa_explicit_bits() + 3>(q, w);
0127   // The computed 'product' is always sufficient.
0128   // Mathematical proof:
0129   // Noble Mushtak and Daniel Lemire, Fast Number Parsing Without Fallback (to appear)
0130   // See script/mushtak_lemire.py
0131 
0132   // The "compute_product_approximation" function can be slightly slower than a branchless approach:
0133   // value128 product = compute_product(q, w);
0134   // but in practice, we can win big with the compute_product_approximation if its additional branch
0135   // is easily predicted. Which is best is data specific.
0136   int upperbit = int(product.high >> 63);
0137 
0138   answer.mantissa = product.high >> (upperbit + 64 - binary::mantissa_explicit_bits() - 3);
0139 
0140   answer.power2 = int32_t(detail::power(int32_t(q)) + upperbit - lz - binary::minimum_exponent());
0141   if (answer.power2 <= 0) { // we have a subnormal?
0142     // Here have that answer.power2 <= 0 so -answer.power2 >= 0
0143     if(-answer.power2 + 1 >= 64) { // if we have more than 64 bits below the minimum exponent, you have a zero for sure.
0144       answer.power2 = 0;
0145       answer.mantissa = 0;
0146       // result should be zero
0147       return answer;
0148     }
0149     // next line is safe because -answer.power2 + 1 < 64
0150     answer.mantissa >>= -answer.power2 + 1;
0151     // Thankfully, we can't have both "round-to-even" and subnormals because
0152     // "round-to-even" only occurs for powers close to 0.
0153     answer.mantissa += (answer.mantissa & 1); // round up
0154     answer.mantissa >>= 1;
0155     // There is a weird scenario where we don't have a subnormal but just.
0156     // Suppose we start with 2.2250738585072013e-308, we end up
0157     // with 0x3fffffffffffff x 2^-1023-53 which is technically subnormal
0158     // whereas 0x40000000000000 x 2^-1023-53  is normal. Now, we need to round
0159     // up 0x3fffffffffffff x 2^-1023-53  and once we do, we are no longer
0160     // subnormal, but we can only know this after rounding.
0161     // So we only declare a subnormal if we are smaller than the threshold.
0162     answer.power2 = (answer.mantissa < (uint64_t(1) << binary::mantissa_explicit_bits())) ? 0 : 1;
0163     return answer;
0164   }
0165 
0166   // usually, we round *up*, but if we fall right in between and and we have an
0167   // even basis, we need to round down
0168   // We are only concerned with the cases where 5**q fits in single 64-bit word.
0169   if ((product.low <= 1) &&  (q >= binary::min_exponent_round_to_even()) && (q <= binary::max_exponent_round_to_even()) &&
0170       ((answer.mantissa & 3) == 1) ) { // we may fall between two floats!
0171     // To be in-between two floats we need that in doing
0172     //   answer.mantissa = product.high >> (upperbit + 64 - binary::mantissa_explicit_bits() - 3);
0173     // ... we dropped out only zeroes. But if this happened, then we can go back!!!
0174     if((answer.mantissa  << (upperbit + 64 - binary::mantissa_explicit_bits() - 3)) ==  product.high) {
0175       answer.mantissa &= ~uint64_t(1);          // flip it so that we do not round up
0176     }
0177   }
0178 
0179   answer.mantissa += (answer.mantissa & 1); // round up
0180   answer.mantissa >>= 1;
0181   if (answer.mantissa >= (uint64_t(2) << binary::mantissa_explicit_bits())) {
0182     answer.mantissa = (uint64_t(1) << binary::mantissa_explicit_bits());
0183     answer.power2++; // undo previous addition
0184   }
0185 
0186   answer.mantissa &= ~(uint64_t(1) << binary::mantissa_explicit_bits());
0187   if (answer.power2 >= binary::infinite_power()) { // infinity
0188     answer.power2 = binary::infinite_power();
0189     answer.mantissa = 0;
0190   }
0191   return answer;
0192 }
0193 
0194 }}}} // namespace fast_float
0195 
0196 #endif