Line data    Source code

       1             : #ifndef FASTFLOAT_DECIMAL_TO_BINARY_H
       2             : #define FASTFLOAT_DECIMAL_TO_BINARY_H
       3             : 
       4             : #include "float_common.h"
       5             : #include "fast_table.h"
       6             : #include <cfloat>
       7             : #include <cinttypes>
       8             : #include <cmath>
       9             : #include <cstdint>
      10             : #include <cstdlib>
      11             : #include <cstring>
      12             : 
      13             : namespace fast_float {
      14             : 
      15             : // This will compute or rather approximate w * 5**q and return a pair of 64-bit
      16             : // words approximating the result, with the "high" part corresponding to the
      17             : // most significant bits and the low part corresponding to the least significant
      18             : // bits.
      19             : //
      20             : template <int bit_precision>
      21             : fastfloat_really_inline FASTFLOAT_CONSTEXPR20 value128
      22             : compute_product_approximation(int64_t q, uint64_t w) {
      23      330041 :   const int index = 2 * int(q - powers::smallest_power_of_five);
      24             :   // For small values of q, e.g., q in [0,27], the answer is always exact
      25             :   // because The line value128 firstproduct = full_multiplication(w,
      26             :   // power_of_five_128[index]); gives the exact answer.
      27             :   value128 firstproduct =
      28      660082 :       full_multiplication(w, powers::power_of_five_128[index]);
      29             :   static_assert((bit_precision >= 0) && (bit_precision <= 64),
      30             :                 " precision should  be in (0,64]");
      31      330041 :   constexpr uint64_t precision_mask =
      32             :       (bit_precision < 64) ? (uint64_t(0xFFFFFFFFFFFFFFFF) >> bit_precision)
      33             :                            : uint64_t(0xFFFFFFFFFFFFFFFF);
      34      330041 :   if ((firstproduct.high & precision_mask) ==
      35             :       precision_mask) { // could further guard with  (lower + w < lower)
      36             :     // regarding the second product, we only need secondproduct.high, but our
      37             :     // expectation is that the compiler will optimize this extra work away if
      38             :     // needed.
      39             :     value128 secondproduct =
      40       34981 :         full_multiplication(w, powers::power_of_five_128[index + 1]);
      41       34981 :     firstproduct.low += secondproduct.high;
      42       34981 :     if (secondproduct.high > firstproduct.low) {
      43       17874 :       firstproduct.high++;
      44             :     }
      45             :   }
      46      330041 :   return firstproduct;
      47             : }
      48             : 
      49             : namespace detail {
      50             : /**
      51             :  * For q in (0,350), we have that
      52             :  *  f = (((152170 + 65536) * q ) >> 16);
      53             :  * is equal to
      54             :  *   floor(p) + q
      55             :  * where
      56             :  *   p = log(5**q)/log(2) = q * log(5)/log(2)
      57             :  *
      58             :  * For negative values of q in (-400,0), we have that
      59             :  *  f = (((152170 + 65536) * q ) >> 16);
      60             :  * is equal to
      61             :  *   -ceil(p) + q
      62             :  * where
      63             :  *   p = log(5**-q)/log(2) = -q * log(5)/log(2)
      64             :  */
      65             : constexpr fastfloat_really_inline int32_t power(int32_t q) noexcept {
      66      330041 :   return (((152170 + 65536) * q) >> 16) + 63;
      67             : }
      68             : } // namespace detail
      69             : 
      70             : // create an adjusted mantissa, biased by the invalid power2
      71             : // for significant digits already multiplied by 10 ** q.
      72             : template <typename binary>
      73             : fastfloat_really_inline FASTFLOAT_CONSTEXPR14 adjusted_mantissa
      74             : compute_error_scaled(int64_t q, uint64_t w, int lz) noexcept {
      75           0 :   int hilz = int(w >> 63) ^ 1;
      76           0 :   adjusted_mantissa answer;
      77           0 :   answer.mantissa = w << hilz;
      78           0 :   int bias = binary::mantissa_explicit_bits() - binary::minimum_exponent();
      79           0 :   answer.power2 = int32_t(detail::power(int32_t(q)) + bias - hilz - lz - 62 +
      80             :                           invalid_am_bias);
      81           0 :   return answer;
      82             : }
      83             : 
      84             : // w * 10 ** q, without rounding the representation up.
      85             : // the power2 in the exponent will be adjusted by invalid_am_bias.
      86             : template <typename binary>
      87             : fastfloat_really_inline FASTFLOAT_CONSTEXPR20 adjusted_mantissa
      88             : compute_error(int64_t q, uint64_t w) noexcept {
      89           0 :   int lz = leading_zeroes(w);
      90           0 :   w <<= lz;
      91             :   value128 product =
      92             :       compute_product_approximation<binary::mantissa_explicit_bits() + 3>(q, w);
      93           0 :   return compute_error_scaled<binary>(q, product.high, lz);
      94             : }
      95             : 
      96             : // w * 10 ** q
      97             : // The returned value should be a valid ieee64 number that simply need to be
      98             : // packed. However, in some very rare cases, the computation will fail. In such
      99             : // cases, we return an adjusted_mantissa with a negative power of 2: the caller
     100             : // should recompute in such cases.
     101             : template <typename binary>
     102             : fastfloat_really_inline FASTFLOAT_CONSTEXPR20 adjusted_mantissa
     103             : compute_float(int64_t q, uint64_t w) noexcept {
     104      330044 :   adjusted_mantissa answer;
     105      330044 :   if ((w == 0) || (q < binary::smallest_power_of_ten())) {
     106           0 :     answer.power2 = 0;
     107           0 :     answer.mantissa = 0;
     108             :     // result should be zero
     109           0 :     return answer;
     110             :   }
     111      330044 :   if (q > binary::largest_power_of_ten()) {
     112             :     // we want to get infinity:
     113           2 :     answer.power2 = binary::infinite_power();
     114           2 :     answer.mantissa = 0;
     115           2 :     return answer;
     116             :   }
     117             :   // At this point in time q is in [powers::smallest_power_of_five,
     118             :   // powers::largest_power_of_five].
     119             : 
     120             :   // We want the most significant bit of i to be 1. Shift if needed.
     121      330041 :   int lz = leading_zeroes(w);
     122      330041 :   w <<= lz;
     123             : 
     124             :   // The required precision is binary::mantissa_explicit_bits() + 3 because
     125             :   // 1. We need the implicit bit
     126             :   // 2. We need an extra bit for rounding purposes
     127             :   // 3. We might lose a bit due to the "upperbit" routine (result too small,
     128             :   // requiring a shift)
     129             : 
     130             :   value128 product =
     131             :       compute_product_approximation<binary::mantissa_explicit_bits() + 3>(q, w);
     132             :   // The computed 'product' is always sufficient.
     133             :   // Mathematical proof:
     134             :   // Noble Mushtak and Daniel Lemire, Fast Number Parsing Without Fallback (to
     135             :   // appear) See script/mushtak_lemire.py
     136             : 
     137             :   // The "compute_product_approximation" function can be slightly slower than a
     138             :   // branchless approach: value128 product = compute_product(q, w); but in
     139             :   // practice, we can win big with the compute_product_approximation if its
     140             :   // additional branch is easily predicted. Which is best is data specific.
     141      330041 :   int upperbit = int(product.high >> 63);
     142      330041 :   int shift = upperbit + 64 - binary::mantissa_explicit_bits() - 3;
     143             : 
     144      330041 :   answer.mantissa = product.high >> shift;
     145             : 
     146      330041 :   answer.power2 = int32_t(detail::power(int32_t(q)) + upperbit - lz -
     147      330041 :                           binary::minimum_exponent());
     148      330041 :   if (answer.power2 <= 0) { // we have a subnormal?
     149             :     // Here have that answer.power2 <= 0 so -answer.power2 >= 0
     150          46 :     if (-answer.power2 + 1 >=
     151             :         64) { // if we have more than 64 bits below the minimum exponent, you
     152             :               // have a zero for sure.
     153           0 :       answer.power2 = 0;
     154           0 :       answer.mantissa = 0;
     155             :       // result should be zero
     156           0 :       return answer;
     157             :     }
     158             :     // next line is safe because -answer.power2 + 1 < 64
     159          46 :     answer.mantissa >>= -answer.power2 + 1;
     160             :     // Thankfully, we can't have both "round-to-even" and subnormals because
     161             :     // "round-to-even" only occurs for powers close to 0.
     162          46 :     answer.mantissa += (answer.mantissa & 1); // round up
     163          46 :     answer.mantissa >>= 1;
     164             :     // There is a weird scenario where we don't have a subnormal but just.
     165             :     // Suppose we start with 2.2250738585072013e-308, we end up
     166             :     // with 0x3fffffffffffff x 2^-1023-53 which is technically subnormal
     167             :     // whereas 0x40000000000000 x 2^-1023-53  is normal. Now, we need to round
     168             :     // up 0x3fffffffffffff x 2^-1023-53  and once we do, we are no longer
     169             :     // subnormal, but we can only know this after rounding.
     170             :     // So we only declare a subnormal if we are smaller than the threshold.
     171          46 :     answer.power2 =
     172          46 :         (answer.mantissa < (uint64_t(1) << binary::mantissa_explicit_bits()))
     173          46 :             ? 0
     174             :             : 1;
     175          46 :     return answer;
     176             :   }
     177             : 
     178             :   // usually, we round *up*, but if we fall right in between and and we have an
     179             :   // even basis, we need to round down
     180             :   // We are only concerned with the cases where 5**q fits in single 64-bit word.
     181        8692 :   if ((product.low <= 1) && (q >= binary::min_exponent_round_to_even()) &&
     182      338702 :       (q <= binary::max_exponent_round_to_even()) &&
     183          15 :       ((answer.mantissa & 3) == 1)) { // we may fall between two floats!
     184             :     // To be in-between two floats we need that in doing
     185             :     //   answer.mantissa = product.high >> (upperbit + 64 -
     186             :     //   binary::mantissa_explicit_bits() - 3);
     187             :     // ... we dropped out only zeroes. But if this happened, then we can go
     188             :     // back!!!
     189           0 :     if ((answer.mantissa << shift) == product.high) {
     190           0 :       answer.mantissa &= ~uint64_t(1); // flip it so that we do not round up
     191             :     }
     192             :   }
     193             : 
     194      329995 :   answer.mantissa += (answer.mantissa & 1); // round up
     195      329995 :   answer.mantissa >>= 1;
     196      329995 :   if (answer.mantissa >= (uint64_t(2) << binary::mantissa_explicit_bits())) {
     197           3 :     answer.mantissa = (uint64_t(1) << binary::mantissa_explicit_bits());
     198           3 :     answer.power2++; // undo previous addition
     199             :   }
     200             : 
     201      329995 :   answer.mantissa &= ~(uint64_t(1) << binary::mantissa_explicit_bits());
     202      329995 :   if (answer.power2 >= binary::infinite_power()) { // infinity
     203           1 :     answer.power2 = binary::infinite_power();
     204           1 :     answer.mantissa = 0;
     205             :   }
     206      329995 :   return answer;
     207             : }
     208             : 
     209             : } // namespace fast_float
     210             : 
     211             : #endif