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