mfbt/double-conversion/bignum-dtoa.cc

changeset 0
6474c204b198
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/mfbt/double-conversion/bignum-dtoa.cc	Wed Dec 31 06:09:35 2014 +0100
     1.3 @@ -0,0 +1,640 @@
     1.4 +// Copyright 2010 the V8 project authors. All rights reserved.
     1.5 +// Redistribution and use in source and binary forms, with or without
     1.6 +// modification, are permitted provided that the following conditions are
     1.7 +// met:
     1.8 +//
     1.9 +//     * Redistributions of source code must retain the above copyright
    1.10 +//       notice, this list of conditions and the following disclaimer.
    1.11 +//     * Redistributions in binary form must reproduce the above
    1.12 +//       copyright notice, this list of conditions and the following
    1.13 +//       disclaimer in the documentation and/or other materials provided
    1.14 +//       with the distribution.
    1.15 +//     * Neither the name of Google Inc. nor the names of its
    1.16 +//       contributors may be used to endorse or promote products derived
    1.17 +//       from this software without specific prior written permission.
    1.18 +//
    1.19 +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
    1.20 +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
    1.21 +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
    1.22 +// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
    1.23 +// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
    1.24 +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
    1.25 +// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
    1.26 +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
    1.27 +// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
    1.28 +// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
    1.29 +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
    1.30 +
    1.31 +#include <math.h>
    1.32 +
    1.33 +#include "bignum-dtoa.h"
    1.34 +
    1.35 +#include "bignum.h"
    1.36 +#include "ieee.h"
    1.37 +
    1.38 +namespace double_conversion {
    1.39 +
    1.40 +static int NormalizedExponent(uint64_t significand, int exponent) {
    1.41 +  ASSERT(significand != 0);
    1.42 +  while ((significand & Double::kHiddenBit) == 0) {
    1.43 +    significand = significand << 1;
    1.44 +    exponent = exponent - 1;
    1.45 +  }
    1.46 +  return exponent;
    1.47 +}
    1.48 +
    1.49 +
    1.50 +// Forward declarations:
    1.51 +// Returns an estimation of k such that 10^(k-1) <= v < 10^k.
    1.52 +static int EstimatePower(int exponent);
    1.53 +// Computes v / 10^estimated_power exactly, as a ratio of two bignums, numerator
    1.54 +// and denominator.
    1.55 +static void InitialScaledStartValues(uint64_t significand,
    1.56 +                                     int exponent,
    1.57 +                                     bool lower_boundary_is_closer,
    1.58 +                                     int estimated_power,
    1.59 +                                     bool need_boundary_deltas,
    1.60 +                                     Bignum* numerator,
    1.61 +                                     Bignum* denominator,
    1.62 +                                     Bignum* delta_minus,
    1.63 +                                     Bignum* delta_plus);
    1.64 +// Multiplies numerator/denominator so that its values lies in the range 1-10.
    1.65 +// Returns decimal_point s.t.
    1.66 +//  v = numerator'/denominator' * 10^(decimal_point-1)
    1.67 +//     where numerator' and denominator' are the values of numerator and
    1.68 +//     denominator after the call to this function.
    1.69 +static void FixupMultiply10(int estimated_power, bool is_even,
    1.70 +                            int* decimal_point,
    1.71 +                            Bignum* numerator, Bignum* denominator,
    1.72 +                            Bignum* delta_minus, Bignum* delta_plus);
    1.73 +// Generates digits from the left to the right and stops when the generated
    1.74 +// digits yield the shortest decimal representation of v.
    1.75 +static void GenerateShortestDigits(Bignum* numerator, Bignum* denominator,
    1.76 +                                   Bignum* delta_minus, Bignum* delta_plus,
    1.77 +                                   bool is_even,
    1.78 +                                   Vector<char> buffer, int* length);
    1.79 +// Generates 'requested_digits' after the decimal point.
    1.80 +static void BignumToFixed(int requested_digits, int* decimal_point,
    1.81 +                          Bignum* numerator, Bignum* denominator,
    1.82 +                          Vector<char>(buffer), int* length);
    1.83 +// Generates 'count' digits of numerator/denominator.
    1.84 +// Once 'count' digits have been produced rounds the result depending on the
    1.85 +// remainder (remainders of exactly .5 round upwards). Might update the
    1.86 +// decimal_point when rounding up (for example for 0.9999).
    1.87 +static void GenerateCountedDigits(int count, int* decimal_point,
    1.88 +                                  Bignum* numerator, Bignum* denominator,
    1.89 +                                  Vector<char>(buffer), int* length);
    1.90 +
    1.91 +
    1.92 +void BignumDtoa(double v, BignumDtoaMode mode, int requested_digits,
    1.93 +                Vector<char> buffer, int* length, int* decimal_point) {
    1.94 +  ASSERT(v > 0);
    1.95 +  ASSERT(!Double(v).IsSpecial());
    1.96 +  uint64_t significand;
    1.97 +  int exponent;
    1.98 +  bool lower_boundary_is_closer;
    1.99 +  if (mode == BIGNUM_DTOA_SHORTEST_SINGLE) {
   1.100 +    float f = static_cast<float>(v);
   1.101 +    ASSERT(f == v);
   1.102 +    significand = Single(f).Significand();
   1.103 +    exponent = Single(f).Exponent();
   1.104 +    lower_boundary_is_closer = Single(f).LowerBoundaryIsCloser();
   1.105 +  } else {
   1.106 +    significand = Double(v).Significand();
   1.107 +    exponent = Double(v).Exponent();
   1.108 +    lower_boundary_is_closer = Double(v).LowerBoundaryIsCloser();
   1.109 +  }
   1.110 +  bool need_boundary_deltas =
   1.111 +      (mode == BIGNUM_DTOA_SHORTEST || mode == BIGNUM_DTOA_SHORTEST_SINGLE);
   1.112 +
   1.113 +  bool is_even = (significand & 1) == 0;
   1.114 +  int normalized_exponent = NormalizedExponent(significand, exponent);
   1.115 +  // estimated_power might be too low by 1.
   1.116 +  int estimated_power = EstimatePower(normalized_exponent);
   1.117 +
   1.118 +  // Shortcut for Fixed.
   1.119 +  // The requested digits correspond to the digits after the point. If the
   1.120 +  // number is much too small, then there is no need in trying to get any
   1.121 +  // digits.
   1.122 +  if (mode == BIGNUM_DTOA_FIXED && -estimated_power - 1 > requested_digits) {
   1.123 +    buffer[0] = '\0';
   1.124 +    *length = 0;
   1.125 +    // Set decimal-point to -requested_digits. This is what Gay does.
   1.126 +    // Note that it should not have any effect anyways since the string is
   1.127 +    // empty.
   1.128 +    *decimal_point = -requested_digits;
   1.129 +    return;
   1.130 +  }
   1.131 +
   1.132 +  Bignum numerator;
   1.133 +  Bignum denominator;
   1.134 +  Bignum delta_minus;
   1.135 +  Bignum delta_plus;
   1.136 +  // Make sure the bignum can grow large enough. The smallest double equals
   1.137 +  // 4e-324. In this case the denominator needs fewer than 324*4 binary digits.
   1.138 +  // The maximum double is 1.7976931348623157e308 which needs fewer than
   1.139 +  // 308*4 binary digits.
   1.140 +  ASSERT(Bignum::kMaxSignificantBits >= 324*4);
   1.141 +  InitialScaledStartValues(significand, exponent, lower_boundary_is_closer,
   1.142 +                           estimated_power, need_boundary_deltas,
   1.143 +                           &numerator, &denominator,
   1.144 +                           &delta_minus, &delta_plus);
   1.145 +  // We now have v = (numerator / denominator) * 10^estimated_power.
   1.146 +  FixupMultiply10(estimated_power, is_even, decimal_point,
   1.147 +                  &numerator, &denominator,
   1.148 +                  &delta_minus, &delta_plus);
   1.149 +  // We now have v = (numerator / denominator) * 10^(decimal_point-1), and
   1.150 +  //  1 <= (numerator + delta_plus) / denominator < 10
   1.151 +  switch (mode) {
   1.152 +    case BIGNUM_DTOA_SHORTEST:
   1.153 +    case BIGNUM_DTOA_SHORTEST_SINGLE:
   1.154 +      GenerateShortestDigits(&numerator, &denominator,
   1.155 +                             &delta_minus, &delta_plus,
   1.156 +                             is_even, buffer, length);
   1.157 +      break;
   1.158 +    case BIGNUM_DTOA_FIXED:
   1.159 +      BignumToFixed(requested_digits, decimal_point,
   1.160 +                    &numerator, &denominator,
   1.161 +                    buffer, length);
   1.162 +      break;
   1.163 +    case BIGNUM_DTOA_PRECISION:
   1.164 +      GenerateCountedDigits(requested_digits, decimal_point,
   1.165 +                            &numerator, &denominator,
   1.166 +                            buffer, length);
   1.167 +      break;
   1.168 +    default:
   1.169 +      UNREACHABLE();
   1.170 +  }
   1.171 +  buffer[*length] = '\0';
   1.172 +}
   1.173 +
   1.174 +
   1.175 +// The procedure starts generating digits from the left to the right and stops
   1.176 +// when the generated digits yield the shortest decimal representation of v. A
   1.177 +// decimal representation of v is a number lying closer to v than to any other
   1.178 +// double, so it converts to v when read.
   1.179 +//
   1.180 +// This is true if d, the decimal representation, is between m- and m+, the
   1.181 +// upper and lower boundaries. d must be strictly between them if !is_even.
   1.182 +//           m- := (numerator - delta_minus) / denominator
   1.183 +//           m+ := (numerator + delta_plus) / denominator
   1.184 +//
   1.185 +// Precondition: 0 <= (numerator+delta_plus) / denominator < 10.
   1.186 +//   If 1 <= (numerator+delta_plus) / denominator < 10 then no leading 0 digit
   1.187 +//   will be produced. This should be the standard precondition.
   1.188 +static void GenerateShortestDigits(Bignum* numerator, Bignum* denominator,
   1.189 +                                   Bignum* delta_minus, Bignum* delta_plus,
   1.190 +                                   bool is_even,
   1.191 +                                   Vector<char> buffer, int* length) {
   1.192 +  // Small optimization: if delta_minus and delta_plus are the same just reuse
   1.193 +  // one of the two bignums.
   1.194 +  if (Bignum::Equal(*delta_minus, *delta_plus)) {
   1.195 +    delta_plus = delta_minus;
   1.196 +  }
   1.197 +  *length = 0;
   1.198 +  while (true) {
   1.199 +    uint16_t digit;
   1.200 +    digit = numerator->DivideModuloIntBignum(*denominator);
   1.201 +    ASSERT(digit <= 9);  // digit is a uint16_t and therefore always positive.
   1.202 +    // digit = numerator / denominator (integer division).
   1.203 +    // numerator = numerator % denominator.
   1.204 +    buffer[(*length)++] = digit + '0';
   1.205 +
   1.206 +    // Can we stop already?
   1.207 +    // If the remainder of the division is less than the distance to the lower
   1.208 +    // boundary we can stop. In this case we simply round down (discarding the
   1.209 +    // remainder).
   1.210 +    // Similarly we test if we can round up (using the upper boundary).
   1.211 +    bool in_delta_room_minus;
   1.212 +    bool in_delta_room_plus;
   1.213 +    if (is_even) {
   1.214 +      in_delta_room_minus = Bignum::LessEqual(*numerator, *delta_minus);
   1.215 +    } else {
   1.216 +      in_delta_room_minus = Bignum::Less(*numerator, *delta_minus);
   1.217 +    }
   1.218 +    if (is_even) {
   1.219 +      in_delta_room_plus =
   1.220 +          Bignum::PlusCompare(*numerator, *delta_plus, *denominator) >= 0;
   1.221 +    } else {
   1.222 +      in_delta_room_plus =
   1.223 +          Bignum::PlusCompare(*numerator, *delta_plus, *denominator) > 0;
   1.224 +    }
   1.225 +    if (!in_delta_room_minus && !in_delta_room_plus) {
   1.226 +      // Prepare for next iteration.
   1.227 +      numerator->Times10();
   1.228 +      delta_minus->Times10();
   1.229 +      // We optimized delta_plus to be equal to delta_minus (if they share the
   1.230 +      // same value). So don't multiply delta_plus if they point to the same
   1.231 +      // object.
   1.232 +      if (delta_minus != delta_plus) {
   1.233 +        delta_plus->Times10();
   1.234 +      }
   1.235 +    } else if (in_delta_room_minus && in_delta_room_plus) {
   1.236 +      // Let's see if 2*numerator < denominator.
   1.237 +      // If yes, then the next digit would be < 5 and we can round down.
   1.238 +      int compare = Bignum::PlusCompare(*numerator, *numerator, *denominator);
   1.239 +      if (compare < 0) {
   1.240 +        // Remaining digits are less than .5. -> Round down (== do nothing).
   1.241 +      } else if (compare > 0) {
   1.242 +        // Remaining digits are more than .5 of denominator. -> Round up.
   1.243 +        // Note that the last digit could not be a '9' as otherwise the whole
   1.244 +        // loop would have stopped earlier.
   1.245 +        // We still have an assert here in case the preconditions were not
   1.246 +        // satisfied.
   1.247 +        ASSERT(buffer[(*length) - 1] != '9');
   1.248 +        buffer[(*length) - 1]++;
   1.249 +      } else {
   1.250 +        // Halfway case.
   1.251 +        // TODO(floitsch): need a way to solve half-way cases.
   1.252 +        //   For now let's round towards even (since this is what Gay seems to
   1.253 +        //   do).
   1.254 +
   1.255 +        if ((buffer[(*length) - 1] - '0') % 2 == 0) {
   1.256 +          // Round down => Do nothing.
   1.257 +        } else {
   1.258 +          ASSERT(buffer[(*length) - 1] != '9');
   1.259 +          buffer[(*length) - 1]++;
   1.260 +        }
   1.261 +      }
   1.262 +      return;
   1.263 +    } else if (in_delta_room_minus) {
   1.264 +      // Round down (== do nothing).
   1.265 +      return;
   1.266 +    } else {  // in_delta_room_plus
   1.267 +      // Round up.
   1.268 +      // Note again that the last digit could not be '9' since this would have
   1.269 +      // stopped the loop earlier.
   1.270 +      // We still have an ASSERT here, in case the preconditions were not
   1.271 +      // satisfied.
   1.272 +      ASSERT(buffer[(*length) -1] != '9');
   1.273 +      buffer[(*length) - 1]++;
   1.274 +      return;
   1.275 +    }
   1.276 +  }
   1.277 +}
   1.278 +
   1.279 +
   1.280 +// Let v = numerator / denominator < 10.
   1.281 +// Then we generate 'count' digits of d = x.xxxxx... (without the decimal point)
   1.282 +// from left to right. Once 'count' digits have been produced we decide wether
   1.283 +// to round up or down. Remainders of exactly .5 round upwards. Numbers such
   1.284 +// as 9.999999 propagate a carry all the way, and change the
   1.285 +// exponent (decimal_point), when rounding upwards.
   1.286 +static void GenerateCountedDigits(int count, int* decimal_point,
   1.287 +                                  Bignum* numerator, Bignum* denominator,
   1.288 +                                  Vector<char>(buffer), int* length) {
   1.289 +  ASSERT(count >= 0);
   1.290 +  for (int i = 0; i < count - 1; ++i) {
   1.291 +    uint16_t digit;
   1.292 +    digit = numerator->DivideModuloIntBignum(*denominator);
   1.293 +    ASSERT(digit <= 9);  // digit is a uint16_t and therefore always positive.
   1.294 +    // digit = numerator / denominator (integer division).
   1.295 +    // numerator = numerator % denominator.
   1.296 +    buffer[i] = digit + '0';
   1.297 +    // Prepare for next iteration.
   1.298 +    numerator->Times10();
   1.299 +  }
   1.300 +  // Generate the last digit.
   1.301 +  uint16_t digit;
   1.302 +  digit = numerator->DivideModuloIntBignum(*denominator);
   1.303 +  if (Bignum::PlusCompare(*numerator, *numerator, *denominator) >= 0) {
   1.304 +    digit++;
   1.305 +  }
   1.306 +  buffer[count - 1] = digit + '0';
   1.307 +  // Correct bad digits (in case we had a sequence of '9's). Propagate the
   1.308 +  // carry until we hat a non-'9' or til we reach the first digit.
   1.309 +  for (int i = count - 1; i > 0; --i) {
   1.310 +    if (buffer[i] != '0' + 10) break;
   1.311 +    buffer[i] = '0';
   1.312 +    buffer[i - 1]++;
   1.313 +  }
   1.314 +  if (buffer[0] == '0' + 10) {
   1.315 +    // Propagate a carry past the top place.
   1.316 +    buffer[0] = '1';
   1.317 +    (*decimal_point)++;
   1.318 +  }
   1.319 +  *length = count;
   1.320 +}
   1.321 +
   1.322 +
   1.323 +// Generates 'requested_digits' after the decimal point. It might omit
   1.324 +// trailing '0's. If the input number is too small then no digits at all are
   1.325 +// generated (ex.: 2 fixed digits for 0.00001).
   1.326 +//
   1.327 +// Input verifies:  1 <= (numerator + delta) / denominator < 10.
   1.328 +static void BignumToFixed(int requested_digits, int* decimal_point,
   1.329 +                          Bignum* numerator, Bignum* denominator,
   1.330 +                          Vector<char>(buffer), int* length) {
   1.331 +  // Note that we have to look at more than just the requested_digits, since
   1.332 +  // a number could be rounded up. Example: v=0.5 with requested_digits=0.
   1.333 +  // Even though the power of v equals 0 we can't just stop here.
   1.334 +  if (-(*decimal_point) > requested_digits) {
   1.335 +    // The number is definitively too small.
   1.336 +    // Ex: 0.001 with requested_digits == 1.
   1.337 +    // Set decimal-point to -requested_digits. This is what Gay does.
   1.338 +    // Note that it should not have any effect anyways since the string is
   1.339 +    // empty.
   1.340 +    *decimal_point = -requested_digits;
   1.341 +    *length = 0;
   1.342 +    return;
   1.343 +  } else if (-(*decimal_point) == requested_digits) {
   1.344 +    // We only need to verify if the number rounds down or up.
   1.345 +    // Ex: 0.04 and 0.06 with requested_digits == 1.
   1.346 +    ASSERT(*decimal_point == -requested_digits);
   1.347 +    // Initially the fraction lies in range (1, 10]. Multiply the denominator
   1.348 +    // by 10 so that we can compare more easily.
   1.349 +    denominator->Times10();
   1.350 +    if (Bignum::PlusCompare(*numerator, *numerator, *denominator) >= 0) {
   1.351 +      // If the fraction is >= 0.5 then we have to include the rounded
   1.352 +      // digit.
   1.353 +      buffer[0] = '1';
   1.354 +      *length = 1;
   1.355 +      (*decimal_point)++;
   1.356 +    } else {
   1.357 +      // Note that we caught most of similar cases earlier.
   1.358 +      *length = 0;
   1.359 +    }
   1.360 +    return;
   1.361 +  } else {
   1.362 +    // The requested digits correspond to the digits after the point.
   1.363 +    // The variable 'needed_digits' includes the digits before the point.
   1.364 +    int needed_digits = (*decimal_point) + requested_digits;
   1.365 +    GenerateCountedDigits(needed_digits, decimal_point,
   1.366 +                          numerator, denominator,
   1.367 +                          buffer, length);
   1.368 +  }
   1.369 +}
   1.370 +
   1.371 +
   1.372 +// Returns an estimation of k such that 10^(k-1) <= v < 10^k where
   1.373 +// v = f * 2^exponent and 2^52 <= f < 2^53.
   1.374 +// v is hence a normalized double with the given exponent. The output is an
   1.375 +// approximation for the exponent of the decimal approimation .digits * 10^k.
   1.376 +//
   1.377 +// The result might undershoot by 1 in which case 10^k <= v < 10^k+1.
   1.378 +// Note: this property holds for v's upper boundary m+ too.
   1.379 +//    10^k <= m+ < 10^k+1.
   1.380 +//   (see explanation below).
   1.381 +//
   1.382 +// Examples:
   1.383 +//  EstimatePower(0)   => 16
   1.384 +//  EstimatePower(-52) => 0
   1.385 +//
   1.386 +// Note: e >= 0 => EstimatedPower(e) > 0. No similar claim can be made for e<0.
   1.387 +static int EstimatePower(int exponent) {
   1.388 +  // This function estimates log10 of v where v = f*2^e (with e == exponent).
   1.389 +  // Note that 10^floor(log10(v)) <= v, but v <= 10^ceil(log10(v)).
   1.390 +  // Note that f is bounded by its container size. Let p = 53 (the double's
   1.391 +  // significand size). Then 2^(p-1) <= f < 2^p.
   1.392 +  //
   1.393 +  // Given that log10(v) == log2(v)/log2(10) and e+(len(f)-1) is quite close
   1.394 +  // to log2(v) the function is simplified to (e+(len(f)-1)/log2(10)).
   1.395 +  // The computed number undershoots by less than 0.631 (when we compute log3
   1.396 +  // and not log10).
   1.397 +  //
   1.398 +  // Optimization: since we only need an approximated result this computation
   1.399 +  // can be performed on 64 bit integers. On x86/x64 architecture the speedup is
   1.400 +  // not really measurable, though.
   1.401 +  //
   1.402 +  // Since we want to avoid overshooting we decrement by 1e10 so that
   1.403 +  // floating-point imprecisions don't affect us.
   1.404 +  //
   1.405 +  // Explanation for v's boundary m+: the computation takes advantage of
   1.406 +  // the fact that 2^(p-1) <= f < 2^p. Boundaries still satisfy this requirement
   1.407 +  // (even for denormals where the delta can be much more important).
   1.408 +
   1.409 +  const double k1Log10 = 0.30102999566398114;  // 1/lg(10)
   1.410 +
   1.411 +  // For doubles len(f) == 53 (don't forget the hidden bit).
   1.412 +  const int kSignificandSize = Double::kSignificandSize;
   1.413 +  double estimate = ceil((exponent + kSignificandSize - 1) * k1Log10 - 1e-10);
   1.414 +  return static_cast<int>(estimate);
   1.415 +}
   1.416 +
   1.417 +
   1.418 +// See comments for InitialScaledStartValues.
   1.419 +static void InitialScaledStartValuesPositiveExponent(
   1.420 +    uint64_t significand, int exponent,
   1.421 +    int estimated_power, bool need_boundary_deltas,
   1.422 +    Bignum* numerator, Bignum* denominator,
   1.423 +    Bignum* delta_minus, Bignum* delta_plus) {
   1.424 +  // A positive exponent implies a positive power.
   1.425 +  ASSERT(estimated_power >= 0);
   1.426 +  // Since the estimated_power is positive we simply multiply the denominator
   1.427 +  // by 10^estimated_power.
   1.428 +
   1.429 +  // numerator = v.
   1.430 +  numerator->AssignUInt64(significand);
   1.431 +  numerator->ShiftLeft(exponent);
   1.432 +  // denominator = 10^estimated_power.
   1.433 +  denominator->AssignPowerUInt16(10, estimated_power);
   1.434 +
   1.435 +  if (need_boundary_deltas) {
   1.436 +    // Introduce a common denominator so that the deltas to the boundaries are
   1.437 +    // integers.
   1.438 +    denominator->ShiftLeft(1);
   1.439 +    numerator->ShiftLeft(1);
   1.440 +    // Let v = f * 2^e, then m+ - v = 1/2 * 2^e; With the common
   1.441 +    // denominator (of 2) delta_plus equals 2^e.
   1.442 +    delta_plus->AssignUInt16(1);
   1.443 +    delta_plus->ShiftLeft(exponent);
   1.444 +    // Same for delta_minus. The adjustments if f == 2^p-1 are done later.
   1.445 +    delta_minus->AssignUInt16(1);
   1.446 +    delta_minus->ShiftLeft(exponent);
   1.447 +  }
   1.448 +}
   1.449 +
   1.450 +
   1.451 +// See comments for InitialScaledStartValues
   1.452 +static void InitialScaledStartValuesNegativeExponentPositivePower(
   1.453 +    uint64_t significand, int exponent,
   1.454 +    int estimated_power, bool need_boundary_deltas,
   1.455 +    Bignum* numerator, Bignum* denominator,
   1.456 +    Bignum* delta_minus, Bignum* delta_plus) {
   1.457 +  // v = f * 2^e with e < 0, and with estimated_power >= 0.
   1.458 +  // This means that e is close to 0 (have a look at how estimated_power is
   1.459 +  // computed).
   1.460 +
   1.461 +  // numerator = significand
   1.462 +  //  since v = significand * 2^exponent this is equivalent to
   1.463 +  //  numerator = v * / 2^-exponent
   1.464 +  numerator->AssignUInt64(significand);
   1.465 +  // denominator = 10^estimated_power * 2^-exponent (with exponent < 0)
   1.466 +  denominator->AssignPowerUInt16(10, estimated_power);
   1.467 +  denominator->ShiftLeft(-exponent);
   1.468 +
   1.469 +  if (need_boundary_deltas) {
   1.470 +    // Introduce a common denominator so that the deltas to the boundaries are
   1.471 +    // integers.
   1.472 +    denominator->ShiftLeft(1);
   1.473 +    numerator->ShiftLeft(1);
   1.474 +    // Let v = f * 2^e, then m+ - v = 1/2 * 2^e; With the common
   1.475 +    // denominator (of 2) delta_plus equals 2^e.
   1.476 +    // Given that the denominator already includes v's exponent the distance
   1.477 +    // to the boundaries is simply 1.
   1.478 +    delta_plus->AssignUInt16(1);
   1.479 +    // Same for delta_minus. The adjustments if f == 2^p-1 are done later.
   1.480 +    delta_minus->AssignUInt16(1);
   1.481 +  }
   1.482 +}
   1.483 +
   1.484 +
   1.485 +// See comments for InitialScaledStartValues
   1.486 +static void InitialScaledStartValuesNegativeExponentNegativePower(
   1.487 +    uint64_t significand, int exponent,
   1.488 +    int estimated_power, bool need_boundary_deltas,
   1.489 +    Bignum* numerator, Bignum* denominator,
   1.490 +    Bignum* delta_minus, Bignum* delta_plus) {
   1.491 +  // Instead of multiplying the denominator with 10^estimated_power we
   1.492 +  // multiply all values (numerator and deltas) by 10^-estimated_power.
   1.493 +
   1.494 +  // Use numerator as temporary container for power_ten.
   1.495 +  Bignum* power_ten = numerator;
   1.496 +  power_ten->AssignPowerUInt16(10, -estimated_power);
   1.497 +
   1.498 +  if (need_boundary_deltas) {
   1.499 +    // Since power_ten == numerator we must make a copy of 10^estimated_power
   1.500 +    // before we complete the computation of the numerator.
   1.501 +    // delta_plus = delta_minus = 10^estimated_power
   1.502 +    delta_plus->AssignBignum(*power_ten);
   1.503 +    delta_minus->AssignBignum(*power_ten);
   1.504 +  }
   1.505 +
   1.506 +  // numerator = significand * 2 * 10^-estimated_power
   1.507 +  //  since v = significand * 2^exponent this is equivalent to
   1.508 +  // numerator = v * 10^-estimated_power * 2 * 2^-exponent.
   1.509 +  // Remember: numerator has been abused as power_ten. So no need to assign it
   1.510 +  //  to itself.
   1.511 +  ASSERT(numerator == power_ten);
   1.512 +  numerator->MultiplyByUInt64(significand);
   1.513 +
   1.514 +  // denominator = 2 * 2^-exponent with exponent < 0.
   1.515 +  denominator->AssignUInt16(1);
   1.516 +  denominator->ShiftLeft(-exponent);
   1.517 +
   1.518 +  if (need_boundary_deltas) {
   1.519 +    // Introduce a common denominator so that the deltas to the boundaries are
   1.520 +    // integers.
   1.521 +    numerator->ShiftLeft(1);
   1.522 +    denominator->ShiftLeft(1);
   1.523 +    // With this shift the boundaries have their correct value, since
   1.524 +    // delta_plus = 10^-estimated_power, and
   1.525 +    // delta_minus = 10^-estimated_power.
   1.526 +    // These assignments have been done earlier.
   1.527 +    // The adjustments if f == 2^p-1 (lower boundary is closer) are done later.
   1.528 +  }
   1.529 +}
   1.530 +
   1.531 +
   1.532 +// Let v = significand * 2^exponent.
   1.533 +// Computes v / 10^estimated_power exactly, as a ratio of two bignums, numerator
   1.534 +// and denominator. The functions GenerateShortestDigits and
   1.535 +// GenerateCountedDigits will then convert this ratio to its decimal
   1.536 +// representation d, with the required accuracy.
   1.537 +// Then d * 10^estimated_power is the representation of v.
   1.538 +// (Note: the fraction and the estimated_power might get adjusted before
   1.539 +// generating the decimal representation.)
   1.540 +//
   1.541 +// The initial start values consist of:
   1.542 +//  - a scaled numerator: s.t. numerator/denominator == v / 10^estimated_power.
   1.543 +//  - a scaled (common) denominator.
   1.544 +//  optionally (used by GenerateShortestDigits to decide if it has the shortest
   1.545 +//  decimal converting back to v):
   1.546 +//  - v - m-: the distance to the lower boundary.
   1.547 +//  - m+ - v: the distance to the upper boundary.
   1.548 +//
   1.549 +// v, m+, m-, and therefore v - m- and m+ - v all share the same denominator.
   1.550 +//
   1.551 +// Let ep == estimated_power, then the returned values will satisfy:
   1.552 +//  v / 10^ep = numerator / denominator.
   1.553 +//  v's boundarys m- and m+:
   1.554 +//    m- / 10^ep == v / 10^ep - delta_minus / denominator
   1.555 +//    m+ / 10^ep == v / 10^ep + delta_plus / denominator
   1.556 +//  Or in other words:
   1.557 +//    m- == v - delta_minus * 10^ep / denominator;
   1.558 +//    m+ == v + delta_plus * 10^ep / denominator;
   1.559 +//
   1.560 +// Since 10^(k-1) <= v < 10^k    (with k == estimated_power)
   1.561 +//  or       10^k <= v < 10^(k+1)
   1.562 +//  we then have 0.1 <= numerator/denominator < 1
   1.563 +//           or    1 <= numerator/denominator < 10
   1.564 +//
   1.565 +// It is then easy to kickstart the digit-generation routine.
   1.566 +//
   1.567 +// The boundary-deltas are only filled if the mode equals BIGNUM_DTOA_SHORTEST
   1.568 +// or BIGNUM_DTOA_SHORTEST_SINGLE.
   1.569 +
   1.570 +static void InitialScaledStartValues(uint64_t significand,
   1.571 +                                     int exponent,
   1.572 +                                     bool lower_boundary_is_closer,
   1.573 +                                     int estimated_power,
   1.574 +                                     bool need_boundary_deltas,
   1.575 +                                     Bignum* numerator,
   1.576 +                                     Bignum* denominator,
   1.577 +                                     Bignum* delta_minus,
   1.578 +                                     Bignum* delta_plus) {
   1.579 +  if (exponent >= 0) {
   1.580 +    InitialScaledStartValuesPositiveExponent(
   1.581 +        significand, exponent, estimated_power, need_boundary_deltas,
   1.582 +        numerator, denominator, delta_minus, delta_plus);
   1.583 +  } else if (estimated_power >= 0) {
   1.584 +    InitialScaledStartValuesNegativeExponentPositivePower(
   1.585 +        significand, exponent, estimated_power, need_boundary_deltas,
   1.586 +        numerator, denominator, delta_minus, delta_plus);
   1.587 +  } else {
   1.588 +    InitialScaledStartValuesNegativeExponentNegativePower(
   1.589 +        significand, exponent, estimated_power, need_boundary_deltas,
   1.590 +        numerator, denominator, delta_minus, delta_plus);
   1.591 +  }
   1.592 +
   1.593 +  if (need_boundary_deltas && lower_boundary_is_closer) {
   1.594 +    // The lower boundary is closer at half the distance of "normal" numbers.
   1.595 +    // Increase the common denominator and adapt all but the delta_minus.
   1.596 +    denominator->ShiftLeft(1);  // *2
   1.597 +    numerator->ShiftLeft(1);    // *2
   1.598 +    delta_plus->ShiftLeft(1);   // *2
   1.599 +  }
   1.600 +}
   1.601 +
   1.602 +
   1.603 +// This routine multiplies numerator/denominator so that its values lies in the
   1.604 +// range 1-10. That is after a call to this function we have:
   1.605 +//    1 <= (numerator + delta_plus) /denominator < 10.
   1.606 +// Let numerator the input before modification and numerator' the argument
   1.607 +// after modification, then the output-parameter decimal_point is such that
   1.608 +//  numerator / denominator * 10^estimated_power ==
   1.609 +//    numerator' / denominator' * 10^(decimal_point - 1)
   1.610 +// In some cases estimated_power was too low, and this is already the case. We
   1.611 +// then simply adjust the power so that 10^(k-1) <= v < 10^k (with k ==
   1.612 +// estimated_power) but do not touch the numerator or denominator.
   1.613 +// Otherwise the routine multiplies the numerator and the deltas by 10.
   1.614 +static void FixupMultiply10(int estimated_power, bool is_even,
   1.615 +                            int* decimal_point,
   1.616 +                            Bignum* numerator, Bignum* denominator,
   1.617 +                            Bignum* delta_minus, Bignum* delta_plus) {
   1.618 +  bool in_range;
   1.619 +  if (is_even) {
   1.620 +    // For IEEE doubles half-way cases (in decimal system numbers ending with 5)
   1.621 +    // are rounded to the closest floating-point number with even significand.
   1.622 +    in_range = Bignum::PlusCompare(*numerator, *delta_plus, *denominator) >= 0;
   1.623 +  } else {
   1.624 +    in_range = Bignum::PlusCompare(*numerator, *delta_plus, *denominator) > 0;
   1.625 +  }
   1.626 +  if (in_range) {
   1.627 +    // Since numerator + delta_plus >= denominator we already have
   1.628 +    // 1 <= numerator/denominator < 10. Simply update the estimated_power.
   1.629 +    *decimal_point = estimated_power + 1;
   1.630 +  } else {
   1.631 +    *decimal_point = estimated_power;
   1.632 +    numerator->Times10();
   1.633 +    if (Bignum::Equal(*delta_minus, *delta_plus)) {
   1.634 +      delta_minus->Times10();
   1.635 +      delta_plus->AssignBignum(*delta_minus);
   1.636 +    } else {
   1.637 +      delta_minus->Times10();
   1.638 +      delta_plus->Times10();
   1.639 +    }
   1.640 +  }
   1.641 +}
   1.642 +
   1.643 +}  // namespace double_conversion

mercurial